diff --git a/include/tensorwrapper/generate/add_noise.hpp b/include/tensorwrapper/generate/add_noise.hpp new file mode 100644 index 00000000..ee10a448 --- /dev/null +++ b/include/tensorwrapper/generate/add_noise.hpp @@ -0,0 +1,52 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once +#include +#include +#include + +namespace tensorwrapper::generate { + +/** @brief Adds clamped normal noise to each element of @p matrix. + * + * Draws `delta ~ Normal(0, t)` and clamps to `[-t, t]` before adding to each + * element. When @p t is zero the input is copied unchanged. + * + * @param[in] matrix The tensor to perturb. + * @param[in] t Non-negative noise scale (standard deviation and clamp bound). + * @param[in,out] gen Random number generator used for the normal draws. + * + * @return A new tensor with the same shape as @p matrix. + * + * @throw std::invalid_argument if @p t is negative. + */ +Tensor add_noise(const Tensor& matrix, double t, std::mt19937& gen); + +/** @brief Overload of add_noise that creates its own RNG from @p seed. + * + * @param[in] matrix The tensor to perturb. + * @param[in] t Non-negative noise scale (standard deviation and clamp bound). + * @param[in] seed Seed for the internal random number generator. A value of + * zero selects a non-deterministic seed. + * + * @return A new tensor with the same shape as @p matrix. + * + * @throw std::invalid_argument if @p t is negative. + */ +Tensor add_noise(const Tensor& matrix, double t, std::uint64_t seed = 42); + +} // namespace tensorwrapper::generate diff --git a/include/tensorwrapper/generate/generate.hpp b/include/tensorwrapper/generate/generate.hpp new file mode 100644 index 00000000..7b672c65 --- /dev/null +++ b/include/tensorwrapper/generate/generate.hpp @@ -0,0 +1,31 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once +#include +#include +#include +#include +#include +#include + +/** @brief Namespace for utilities that generate synthetic tensors and matrices. + * + * The functions in this namespace create reproducible test matrices, + * eigenvalue + * spectra, and related quantities for unit testing and benchmarking. + */ +namespace tensorwrapper::generate {} diff --git a/include/tensorwrapper/generate/generate_eigen_system.hpp b/include/tensorwrapper/generate/generate_eigen_system.hpp new file mode 100644 index 00000000..cbaeb653 --- /dev/null +++ b/include/tensorwrapper/generate/generate_eigen_system.hpp @@ -0,0 +1,56 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once +#include +#include +#include + +namespace tensorwrapper::generate { + +/** @brief A symmetric matrix together with its eigen-decomposition. + * + * The member @p matrix satisfies @f$M = Q D Q^T@f$ where @p eigenvectors is + * @f$Q@f$ and @p eigenvalues holds the diagonal entries of @f$D@f$. + */ +struct EigenSystem { + /// Dimension of the square matrix. + std::size_t n = 0; + /// Rank-1 tensor of length @p n containing the eigenvalues. + Tensor eigenvalues; + /// Rank-2 tensor of shape `(n, n)` whose columns are the eigenvectors. + Tensor eigenvectors; + /// Rank-2 tensor of shape `(n, n)` representing the symmetric matrix. + Tensor matrix; +}; + +/** @brief Generates a reproducible symmetric matrix and its + * eigen-decomposition. + * + * Constructs @f$M = Q D Q^T@f$ where @f$D@f$ is built from eigenvalues + * generated according to @p spec and @f$Q@f$ is a random orthogonal matrix. + * + * @param[in] spec Parameters controlling the matrix dimension, spectrum, and + * random seed. + * + * @return An @ref EigenSystem populated with @p spec.n, the eigenvalues, + * eigenvectors, and the assembled matrix. + * + * @throw std::invalid_argument if @p spec.n is outside `[1, kMaxMatrixDim]`. + */ +EigenSystem generate_eigen_system(const SymmetricMatrixSpec& spec); + +} // namespace tensorwrapper::generate diff --git a/include/tensorwrapper/generate/generate_eigenvalues.hpp b/include/tensorwrapper/generate/generate_eigenvalues.hpp new file mode 100644 index 00000000..2ae3d46c --- /dev/null +++ b/include/tensorwrapper/generate/generate_eigenvalues.hpp @@ -0,0 +1,43 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once +#include +#include +#include + +namespace tensorwrapper::generate { + +/** @brief Generates a sorted rank-1 tensor of eigenvalues from @p spec. + * + * The spectrum spans + * @f$[\texttt{spec.min\_eigenvalue}, + * \texttt{spec.min\_eigenvalue} \times \texttt{spec.condition\_number}]@f$ + * using the spacing strategy given by @p spec.spacing. + * + * @param[in] spec Parameters controlling the eigenvalue distribution. + * @param[in,out] gen Random number generator used when @p spec.spacing + * requires random draws. + * + * @return A rank-1 tensor of length @p spec.n containing the eigenvalues in + * ascending order. + * + * @throw std::invalid_argument if @p spec.n is outside `[1, kMaxMatrixDim]` or + * if @p spec.spacing is invalid. + */ +Tensor generate_eigenvalues(const SymmetricMatrixSpec& spec, std::mt19937& gen); + +} // namespace tensorwrapper::generate diff --git a/include/tensorwrapper/generate/generate_utils.hpp b/include/tensorwrapper/generate/generate_utils.hpp new file mode 100644 index 00000000..15a2bfb2 --- /dev/null +++ b/include/tensorwrapper/generate/generate_utils.hpp @@ -0,0 +1,103 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once +#include +#include +#include + +namespace tensorwrapper::generate { + +/// Maximum supported dimension for generated square matrices. +constexpr std::size_t kMaxMatrixDim = 10; + +/** @brief Specifies how eigenvalues are distributed between the endpoints. + * + * Each spacing mode fills the interval + * @f$[\lambda_{\min}, \lambda_{\max}]@f$ with @f$n@f$ eigenvalues, where + * @f$\lambda_{\max} = \lambda_{\min} \times@f$ the condition number. + */ +enum class EigenvalueSpacing { + /// Uniform spacing with @f$\Delta\lambda = (\lambda_{\max} - + /// \lambda_{\min}) / (n - 1)@f$. + Linear, + /// Uniform spacing in log space with + /// @f$\Delta\log\lambda = \log(\lambda_{\max} / \lambda_{\min}) / (n - + /// 1)@f$. + Logarithmic, + /// Eigenvalues are grouped into clusters of width @p cluster_width that are + /// separated by @f$(\lambda_{\max} - \lambda_{\min}) / (n_{\text{clusters}} + /// - 1)@f$. + Clustered, + /// Same cluster centers as @ref Clustered, but all eigenvalues in a cluster + /// are identical. + Degenerate +}; + +/** @brief Parameters controlling the generation of a symmetric test matrix. + * + * The resulting matrix has eigenvalues in + * @f$[\texttt{min\_eigenvalue}, + * \texttt{min\_eigenvalue} \times \texttt{condition\_number}]@f$ with the + * spacing determined by @p spacing. + */ +struct SymmetricMatrixSpec { + /// Dimension of the square matrix. + std::size_t n = 2; + /// Ratio of the largest to smallest eigenvalue. + double condition_number = 10.0; + /// Smallest eigenvalue in the spectrum. + double min_eigenvalue = 1.0; + /// Strategy used to distribute eigenvalues between the endpoints. + EigenvalueSpacing spacing = EigenvalueSpacing::Linear; + /// Number of eigenvalue clusters when @p spacing is @ref Clustered or + /// @ref Degenerate. + std::size_t n_clusters = 1; + /// Half-width of each cluster when @p spacing is @ref Clustered. + double cluster_width = 1e-6; + /// Seed for random number generation. A value of zero selects a + /// non-deterministic seed. + std::uint64_t seed = 42; +}; + +/** @brief Creates a Mersenne Twister RNG from @p seed. + * + * @param[in] seed The seed value. When zero, a seed is drawn from + * `std::random_device`. + * + * @return A `std::mt19937` generator initialized with @p seed. + */ +inline std::mt19937 make_rng(std::uint64_t seed) { + if(seed == 0) { + std::random_device rd; + return std::mt19937(rd()); + } + return std::mt19937(static_cast(seed)); +} + +/** @brief Validates that @p n is an allowed matrix dimension. + * + * @param[in] n The dimension to validate. + * + * @throw std::invalid_argument if @p n is not in `[1, kMaxMatrixDim]`. + */ +inline void require_valid_n(std::size_t n) { + if(n < 1 || n > kMaxMatrixDim) { + throw std::invalid_argument("n must be in [1, kMaxMatrixDim]"); + } +} + +} // namespace tensorwrapper::generate diff --git a/include/tensorwrapper/generate/identity_matrix.hpp b/include/tensorwrapper/generate/identity_matrix.hpp new file mode 100644 index 00000000..9e55f45e --- /dev/null +++ b/include/tensorwrapper/generate/identity_matrix.hpp @@ -0,0 +1,41 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once +#include +#include +#include +#include +#include + +namespace tensorwrapper::generate { + +/** @brief Creates an @p n x @p n identity matrix. + * + * @param[in] n Matrix dimension. Must be in `[1, kMaxMatrixDim]`. + * + * @return A rank-2 tensor with ones on the diagonal and zeros elsewhere. + * + * @throw std::invalid_argument if @p n is outside the allowed range. + */ +inline Tensor identity_matrix(std::size_t n) { + require_valid_n(n); + std::vector values(n, 1.0); + auto values_tensor = utilities::make_tensor({n}, values); + return utilities::diagonal_matrix(values_tensor); +} + +} // namespace tensorwrapper::generate diff --git a/include/tensorwrapper/generate/random_orthogonal_matrix.hpp b/include/tensorwrapper/generate/random_orthogonal_matrix.hpp new file mode 100644 index 00000000..0be2a3f9 --- /dev/null +++ b/include/tensorwrapper/generate/random_orthogonal_matrix.hpp @@ -0,0 +1,37 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once +#include +#include + +namespace tensorwrapper::generate { + +/** @brief Creates a random @p n x @p n orthogonal matrix. + * + * Draws entries from a standard normal distribution and applies Householder QR + * factorization to obtain an orthogonal matrix @f$Q@f$. + * + * @param[in] n Matrix dimension. Must be in `[1, kMaxMatrixDim]`. + * @param[in,out] gen Random number generator used for the normal draws. + * + * @return A rank-2 tensor whose columns form an orthonormal basis. + * + * @throw std::invalid_argument if @p n is outside the allowed range. + */ +Tensor random_orthogonal_matrix(std::size_t n, std::mt19937& gen); + +} // namespace tensorwrapper::generate diff --git a/include/tensorwrapper/tensorwrapper.hpp b/include/tensorwrapper/tensorwrapper.hpp index 861b0131..b6505f26 100644 --- a/include/tensorwrapper/tensorwrapper.hpp +++ b/include/tensorwrapper/tensorwrapper.hpp @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include diff --git a/src/tensorwrapper/generate/add_noise.cpp b/src/tensorwrapper/generate/add_noise.cpp new file mode 100644 index 00000000..f6cf6b08 --- /dev/null +++ b/src/tensorwrapper/generate/add_noise.cpp @@ -0,0 +1,58 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace tensorwrapper::generate { +namespace { + +std::vector shape_extents(const auto& shape) { + std::vector rv(shape.rank()); + for(std::size_t d = 0; d < shape.rank(); ++d) rv[d] = shape.extent(d); + return rv; +} + +} // namespace + +Tensor add_noise(const Tensor& matrix, double t, std::mt19937& gen) { + if(t < 0.0) { throw std::invalid_argument("t must be non-negative"); } + + auto& in_buf = buffer::make_contiguous(matrix.buffer()); + auto in_data = buffer::get_raw_data(in_buf); + std::vector data(in_data.begin(), in_data.end()); + + if(t > 0.0) { + std::normal_distribution dist(0.0, t); + for(auto& x : data) { x += std::clamp(dist(gen), -t, t); } + } + + return utilities::make_tensor(shape_extents(in_buf.shape()), data.begin(), + data.end()); +} + +Tensor add_noise(const Tensor& matrix, double t, std::uint64_t seed) { + auto gen = make_rng(seed); + return add_noise(matrix, t, gen); +} + +} // namespace tensorwrapper::generate diff --git a/src/tensorwrapper/generate/generate_eigen_system.cpp b/src/tensorwrapper/generate/generate_eigen_system.cpp new file mode 100644 index 00000000..c3b67d7f --- /dev/null +++ b/src/tensorwrapper/generate/generate_eigen_system.cpp @@ -0,0 +1,45 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include + +namespace tensorwrapper::generate { + +EigenSystem generate_eigen_system(const SymmetricMatrixSpec& spec) { + require_valid_n(spec.n); + auto gen = make_rng(spec.seed); + const auto n = spec.n; + + EigenSystem rv; + rv.n = n; + rv.eigenvalues = generate_eigenvalues(spec, gen); + rv.eigenvectors = random_orthogonal_matrix(n, gen); + + const auto D = utilities::diagonal_matrix(rv.eigenvalues); + + Tensor qd; + qd("i,k") = rv.eigenvectors("i,l") * D("l,k"); + + Tensor matrix; + matrix("i,j") = qd("i,k") * rv.eigenvectors("j,k"); + rv.matrix = std::move(matrix); + return rv; +} + +} // namespace tensorwrapper::generate diff --git a/src/tensorwrapper/generate/generate_eigenvalues.cpp b/src/tensorwrapper/generate/generate_eigenvalues.cpp new file mode 100644 index 00000000..1a95eca1 --- /dev/null +++ b/src/tensorwrapper/generate/generate_eigenvalues.cpp @@ -0,0 +1,172 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include +#include + +namespace tensorwrapper::generate { +namespace { + +/** @brief Generates @p n eigenvalues with uniform spacing in @f$\lambda@f$. + * + * @param[in] n Number of eigenvalues. + * @param[in] lambda_min Smallest eigenvalue. + * @param[in] lambda_max Largest eigenvalue. + * + * @return A vector of length @p n with values from @p lambda_min to + * @p lambda_max. + */ +auto linear_spacing(std::size_t n, double lambda_min, double lambda_max) { + const auto dx = (lambda_max - lambda_min) / (n - 1); + std::vector values(n); + for(std::size_t i = 0; i < n; ++i) { + values[i] = lambda_min + static_cast(i) * dx; + } + return values; +} + +/** @brief Generates @p n eigenvalues with uniform spacing in @f$\log\lambda@f$. + * + * @param[in] n Number of eigenvalues. + * @param[in] lambda_min Smallest eigenvalue. Must be positive. + * @param[in] lambda_max Largest eigenvalue. Must be positive. + * + * @return A vector of length @p n with values from @p lambda_min to + * @p lambda_max on a log scale. + */ +auto logarithmic_spacing(std::size_t n, double lambda_min, double lambda_max) { + const double log_min = std::log(lambda_min); + const double log_max = std::log(lambda_max); + const double dlog = (log_max - log_min) / (n - 1); + std::vector values(n); + for(std::size_t i = 0; i < n; ++i) { + values[i] = std::exp(log_min + static_cast(i) * dlog); + } + return values; +} + +/** @brief Generates @p n eigenvalues grouped into @p n_clusters jittered + * clusters. + * + * @param[in] n Number of eigenvalues. + * @param[in] n_clusters Number of cluster centers. + * @param[in] lambda_min Smallest cluster center. + * @param[in] lambda_max Largest cluster center. + * @param[in] cluster_width Half-width of the uniform jitter around each + * center. + * @param[in,out] gen Random number generator used for the jitter draws. + * + * @return A vector of length @p n with eigenvalues assigned cyclically to + * clusters. + */ +auto clustered_spacing(std::size_t n, std::size_t n_clusters, double lambda_min, + double lambda_max, double cluster_width, + std::mt19937& gen) { + std::uniform_real_distribution jitter(-cluster_width, + cluster_width); + std::vector values(n); + std::vector cluster_centers(n_clusters); + if(n_clusters == 1) { + cluster_centers[0] = lambda_min; + } else { + const double dx = (lambda_max - lambda_min) / (n_clusters - 1); + for(std::size_t c = 0; c < n_clusters; ++c) { + cluster_centers[c] = lambda_min + static_cast(c) * dx; + } + } + + for(std::size_t i = 0; i < n; ++i) { + const auto cluster_id = i % n_clusters; + values[i] = cluster_centers[cluster_id] + jitter(gen); + } + return values; +} + +/** @brief Generates @p n eigenvalues with exact degeneracies within each + * cluster. + * + * @param[in] n Number of eigenvalues. + * @param[in] n_clusters Number of distinct eigenvalue plateaus. + * @param[in] lambda_min Value of the smallest plateau. + * @param[in] lambda_max Value of the largest plateau. + * + * @return A vector of length @p n with eigenvalues assigned cyclically to + * @p n_clusters plateaus. + */ +auto degenerate_spacing(std::size_t n, std::size_t n_clusters, + double lambda_min, double lambda_max) { + std::vector values(n); + if(n_clusters <= 1) { + std::fill(values.begin(), values.end(), lambda_min); + } else { + const auto n_plateaus = std::min(n_clusters, n); + const auto nm1 = std::max(1, n_plateaus - 1); + const double dx = (lambda_max - lambda_min) / static_cast(nm1); + for(std::size_t i = 0; i < n; ++i) { + const auto plateau = i % n_plateaus; + values[i] = lambda_min + static_cast(plateau) * dx; + } + } + return values; +} +} // namespace + +Tensor generate_eigenvalues(const SymmetricMatrixSpec& spec, + std::mt19937& gen) { + require_valid_n(spec.n); + const auto n = spec.n; + + const double lambda_min = spec.min_eigenvalue; + const double lambda_max = spec.min_eigenvalue * spec.condition_number; + + std::vector values; + + if(n == 1) return utilities::make_tensor({n}, values); + + switch(spec.spacing) { + case EigenvalueSpacing::Linear: { + values = linear_spacing(n, lambda_min, lambda_max); + break; + } + case EigenvalueSpacing::Logarithmic: { + values = logarithmic_spacing(n, lambda_min, lambda_max); + break; + } + case EigenvalueSpacing::Clustered: { + const auto n_clusters = std::max(1, spec.n_clusters); + const auto n_clusters_clamped = std::min(n_clusters, n); + values = clustered_spacing(n, n_clusters_clamped, lambda_min, + lambda_max, spec.cluster_width, gen); + break; + } + case EigenvalueSpacing::Degenerate: { + values = + degenerate_spacing(n, spec.n_clusters, lambda_min, lambda_max); + break; + } + default: { + throw std::invalid_argument("Invalid eigenvalue spacing"); + } + } + + std::sort(values.begin(), values.end()); + return utilities::make_tensor({n}, values); +} + +} // namespace tensorwrapper::generate diff --git a/src/tensorwrapper/generate/random_orthogonal_matrix.cpp b/src/tensorwrapper/generate/random_orthogonal_matrix.cpp new file mode 100644 index 00000000..3822c62f --- /dev/null +++ b/src/tensorwrapper/generate/random_orthogonal_matrix.cpp @@ -0,0 +1,49 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include +#include + +namespace tensorwrapper::generate { + +Tensor random_orthogonal_matrix(std::size_t n, std::mt19937& gen) { + require_valid_n(n); + + Eigen::MatrixXd M(n, n); + std::normal_distribution dist(0.0, 1.0); + for(std::size_t i = 0; i < n; ++i) { + for(std::size_t j = 0; j < n; ++j) { + M(static_cast(i), static_cast(j)) = + dist(gen); + } + } + Eigen::HouseholderQR qr(M); + const Eigen::MatrixXd Q = qr.householderQ(); + + std::vector data(n * n); + for(std::size_t i = 0; i < n; ++i) { + for(std::size_t j = 0; j < n; ++j) { + data[i * n + j] = + Q(static_cast(i), static_cast(j)); + } + } + return utilities::make_tensor({n, n}, data); +} + +} // namespace tensorwrapper::generate diff --git a/tests/cxx/unit_tests/tensorwrapper/generate/add_noise.cpp b/tests/cxx/unit_tests/tensorwrapper/generate/add_noise.cpp new file mode 100644 index 00000000..9aaca48f --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/generate/add_noise.cpp @@ -0,0 +1,87 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include +#include +#include + +using namespace tensorwrapper; +using namespace tensorwrapper::buffer; +using namespace tensorwrapper::generate; +using namespace tensorwrapper::operations; +using namespace tensorwrapper::utilities; + +namespace { +double elem_as_double(const Contiguous::const_reference& elem) { + using wtf::fp::float_cast; + try { + return float_cast(elem); + } catch(const std::runtime_error&) { + return static_cast(float_cast(elem)); + } +} +} // namespace + +TEST_CASE("add_noise") { + const auto matrix = make_tensor({2, 2}, std::vector{1, 2, 3, 4}); + + SECTION("t equals zero returns copy") { + auto gen = make_rng(1); + auto out = add_noise(matrix, 0.0, gen); + REQUIRE(approximately_equal(out, matrix)); + } + + SECTION("within t of input") { + const double t = 0.05; + auto gen = make_rng(17); + auto out = add_noise(matrix, t, gen); + + auto in_buf = make_contiguous(matrix.buffer()); + auto out_buf = make_contiguous(out.buffer()); + for(std::size_t i = 0; i < 2; ++i) { + for(std::size_t j = 0; j < 2; ++j) { + const auto x = elem_as_double(in_buf.get_elem({i, j})); + const auto y = elem_as_double(out_buf.get_elem({i, j})); + REQUIRE(std::abs(y - x) <= t); + } + } + } + + SECTION("deterministic for fixed seed") { + const double t = 0.01; + auto out1 = add_noise(matrix, t, 99); + auto out2 = add_noise(matrix, t, 99); + REQUIRE(approximately_equal(out1, out2)); + } + + SECTION("shape preserved") { + auto gen = make_rng(3); + auto out = add_noise(matrix, 0.01, gen); + auto in_shape = make_contiguous(matrix.buffer()).shape(); + auto out_shape = make_contiguous(out.buffer()).shape(); + REQUIRE(out_shape.rank() == in_shape.rank()); + REQUIRE(out_shape.extent(0) == in_shape.extent(0)); + REQUIRE(out_shape.extent(1) == in_shape.extent(1)); + } + + SECTION("invalid t throws") { + auto gen = make_rng(1); + REQUIRE_THROWS_AS(add_noise(matrix, -1.0, gen), std::invalid_argument); + } +} diff --git a/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigen_system.cpp b/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigen_system.cpp new file mode 100644 index 00000000..3b72828b --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigen_system.cpp @@ -0,0 +1,127 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include +#include +#include + +using namespace tensorwrapper; +using namespace tensorwrapper::buffer; +using namespace tensorwrapper::generate; +using namespace tensorwrapper::operations; +using namespace tensorwrapper::utilities; + +TEST_CASE("generate_eigen_system") { + SECTION("shapes") { + SymmetricMatrixSpec spec; + spec.n = 4; + spec.condition_number = 1e3; + spec.spacing = EigenvalueSpacing::Linear; + spec.seed = 11; + auto system = generate_eigen_system(spec); + + REQUIRE(system.n == 4); + REQUIRE( + make_contiguous(system.eigenvalues.buffer()).shape().extent(0) == 4); + REQUIRE( + make_contiguous(system.eigenvectors.buffer()).shape().extent(0) == 4); + REQUIRE( + make_contiguous(system.eigenvectors.buffer()).shape().extent(1) == 4); + REQUIRE(make_contiguous(system.matrix.buffer()).shape().extent(0) == 4); + REQUIRE(make_contiguous(system.matrix.buffer()).shape().extent(1) == 4); + } + + SECTION("symmetry") { + SymmetricMatrixSpec spec; + spec.n = 3; + spec.seed = 5; + auto system = generate_eigen_system(spec); + + Tensor sym; + sym("i,j") = system.matrix("i,j"); + Tensor tran; + tran("i,j") = system.matrix("j,i"); + REQUIRE(approximately_equal(sym, tran, 1e-12)); + } + + SECTION("orthonormal eigenvectors") { + SymmetricMatrixSpec spec; + spec.n = 3; + spec.seed = 7; + auto system = generate_eigen_system(spec); + + Tensor product; + product("i,k") = + system.eigenvectors("i,j") * system.eigenvectors("k,j"); + + auto ones = make_tensor({3}, std::vector{1, 1, 1}); + auto ident = diagonal_matrix(ones); + REQUIRE(approximately_equal(product, ident, 1e-12)); + } + + SECTION("exact eigenpair") { + SymmetricMatrixSpec spec; + spec.n = 4; + spec.seed = 13; + auto system = generate_eigen_system(spec); + + Tensor av; + av("i,k") = system.matrix("i,j") * system.eigenvectors("j,k"); + + const auto D = diagonal_matrix(system.eigenvalues); + Tensor scaled; + scaled("i,k") = system.eigenvectors("i,l") * D("l,k"); + REQUIRE(approximately_equal(av, scaled, 1e-12)); + } + + SECTION("deterministic for fixed seed") { + SymmetricMatrixSpec spec; + spec.n = 4; + spec.seed = 17; + auto system1 = generate_eigen_system(spec); + auto system2 = generate_eigen_system(spec); + REQUIRE(approximately_equal(system1.eigenvalues, system2.eigenvalues)); + REQUIRE( + approximately_equal(system1.eigenvectors, system2.eigenvectors)); + REQUIRE(approximately_equal(system1.matrix, system2.matrix)); + } + + SECTION("degenerate") { + SymmetricMatrixSpec spec; + spec.n = 2; + spec.condition_number = 1.0; + spec.min_eigenvalue = 1.0; + spec.spacing = EigenvalueSpacing::Degenerate; + spec.n_clusters = 1; + spec.seed = 7; + auto system = generate_eigen_system(spec); + + auto evals = make_tensor({2}, std::vector{1.0, 1.0}); + REQUIRE(approximately_equal(system.eigenvalues, evals)); + + auto ident = make_tensor({2, 2}, std::vector{1, 0, 0, 1}); + REQUIRE(approximately_equal(system.matrix, ident, 1e-12)); + } + + SECTION("invalid n throws") { + SymmetricMatrixSpec spec; + spec.n = 0; + REQUIRE_THROWS_AS(generate_eigen_system(spec), std::invalid_argument); + } +} diff --git a/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigenvalues.cpp b/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigenvalues.cpp new file mode 100644 index 00000000..98cc4e0f --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigenvalues.cpp @@ -0,0 +1,129 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace tensorwrapper; +using namespace tensorwrapper::buffer; +using namespace tensorwrapper::generate; +using namespace tensorwrapper::operations; +using namespace tensorwrapper::utilities; + +namespace { +double elem_as_double(const Contiguous::const_reference& elem) { + using wtf::fp::float_cast; + try { + return float_cast(elem); + } catch(const std::runtime_error&) { + return static_cast(float_cast(elem)); + } +} + +std::vector tensor_to_vector(const Tensor& t) { + auto buf = make_contiguous(t.buffer()); + std::vector rv(buf.shape().extent(0)); + for(std::size_t i = 0; i < rv.size(); ++i) { + rv[i] = elem_as_double(buf.get_elem({i})); + } + return rv; +} +} // namespace + +TEST_CASE("generate_eigenvalues") { + SECTION("linear spacing") { + SymmetricMatrixSpec spec; + spec.n = 4; + spec.min_eigenvalue = 1.0; + spec.condition_number = 10.0; + spec.spacing = EigenvalueSpacing::Linear; + auto gen = make_rng(1); + auto result = generate_eigenvalues(spec, gen); + auto corr = make_tensor({4}, std::vector{1, 4, 7, 10}); + REQUIRE(approximately_equal(result, corr)); + } + + SECTION("logarithmic spacing") { + SymmetricMatrixSpec spec; + spec.n = 3; + spec.min_eigenvalue = 1.0; + spec.condition_number = 100.0; + spec.spacing = EigenvalueSpacing::Logarithmic; + auto gen = make_rng(1); + auto values = tensor_to_vector(generate_eigenvalues(spec, gen)); + REQUIRE(values.size() == 3); + REQUIRE(values[0] == Catch::Approx(1.0)); + REQUIRE(values[1] == Catch::Approx(10.0)); + REQUIRE(values[2] == Catch::Approx(100.0)); + REQUIRE(std::is_sorted(values.begin(), values.end())); + } + + SECTION("degenerate spacing") { + SymmetricMatrixSpec spec; + spec.n = 3; + spec.min_eigenvalue = 2.5; + spec.condition_number = 10.0; + spec.spacing = EigenvalueSpacing::Degenerate; + spec.n_clusters = 1; + auto gen = make_rng(1); + auto values = tensor_to_vector(generate_eigenvalues(spec, gen)); + REQUIRE(values.size() == 3); + for(const auto v : values) { REQUIRE(v == Catch::Approx(2.5)); } + } + + SECTION("clustered spacing") { + SymmetricMatrixSpec spec; + spec.n = 6; + spec.min_eigenvalue = 1.0; + spec.condition_number = 100.0; + spec.spacing = EigenvalueSpacing::Clustered; + spec.n_clusters = 3; + spec.cluster_width = 1e-6; + spec.seed = 23; + auto gen = make_rng(spec.seed); + auto values = tensor_to_vector(generate_eigenvalues(spec, gen)); + REQUIRE(values.size() == 6); + REQUIRE(std::is_sorted(values.begin(), values.end())); + + const double lambda_max = spec.min_eigenvalue * spec.condition_number; + const double dx = (lambda_max - spec.min_eigenvalue) / 2.0; + std::vector centers(3); + for(std::size_t c = 0; c < 3; ++c) { + centers[c] = spec.min_eigenvalue + static_cast(c) * dx; + } + for(const auto v : values) { + const auto near_center = + std::any_of(centers.begin(), centers.end(), [&](double center) { + return std::abs(v - center) <= spec.cluster_width; + }); + REQUIRE(near_center); + } + } + + SECTION("invalid n throws") { + SymmetricMatrixSpec spec; + spec.n = 0; + auto gen = make_rng(1); + REQUIRE_THROWS_AS(generate_eigenvalues(spec, gen), + std::invalid_argument); + } +} diff --git a/tests/cxx/unit_tests/tensorwrapper/generate/generate_utils.cpp b/tests/cxx/unit_tests/tensorwrapper/generate/generate_utils.cpp new file mode 100644 index 00000000..04b6b98f --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/generate/generate_utils.cpp @@ -0,0 +1,35 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include + +using namespace tensorwrapper::generate; + +TEST_CASE("generate_utils") { + SECTION("make_rng is deterministic for fixed seed") { + auto gen1 = make_rng(42); + auto gen2 = make_rng(42); + REQUIRE(gen1() == gen2()); + } + + SECTION("require_valid_n throws for invalid n") { + REQUIRE_THROWS_AS(require_valid_n(0), std::invalid_argument); + REQUIRE_THROWS_AS(require_valid_n(11), std::invalid_argument); + REQUIRE_NOTHROW(require_valid_n(1)); + REQUIRE_NOTHROW(require_valid_n(10)); + } +} diff --git a/tests/cxx/unit_tests/tensorwrapper/generate/identity_matrix.cpp b/tests/cxx/unit_tests/tensorwrapper/generate/identity_matrix.cpp new file mode 100644 index 00000000..c092a8de --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/generate/identity_matrix.cpp @@ -0,0 +1,44 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include + +using namespace tensorwrapper; +using namespace tensorwrapper::generate; +using namespace tensorwrapper::operations; +using namespace tensorwrapper::utilities; + +TEST_CASE("identity_matrix") { + SECTION("2 by 2") { + auto result = identity_matrix(2); + auto corr = make_tensor({2, 2}, std::vector{1, 0, 0, 1}); + REQUIRE(approximately_equal(result, corr)); + } + + SECTION("3 by 3") { + auto result = identity_matrix(3); + auto corr = + make_tensor({3, 3}, std::vector{1, 0, 0, 0, 1, 0, 0, 0, 1}); + REQUIRE(approximately_equal(result, corr)); + } + + SECTION("invalid n throws") { + REQUIRE_THROWS_AS(identity_matrix(0), std::invalid_argument); + } +} diff --git a/tests/cxx/unit_tests/tensorwrapper/generate/random_orthogonal_matrix.cpp b/tests/cxx/unit_tests/tensorwrapper/generate/random_orthogonal_matrix.cpp new file mode 100644 index 00000000..6761c082 --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/generate/random_orthogonal_matrix.cpp @@ -0,0 +1,58 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include +#include +#include + +using namespace tensorwrapper; +using namespace tensorwrapper::buffer; +using namespace tensorwrapper::generate; +using namespace tensorwrapper::operations; +using namespace tensorwrapper::utilities; + +TEST_CASE("random_orthogonal_matrix") { + SECTION("shape") { + auto gen = make_rng(7); + auto Q = random_orthogonal_matrix(3, gen); + auto buf = make_contiguous(Q.buffer()); + REQUIRE(buf.shape().rank() == 2); + REQUIRE(buf.shape().extent(0) == 3); + REQUIRE(buf.shape().extent(1) == 3); + } + + SECTION("orthogonality") { + auto gen = make_rng(11); + auto Q = random_orthogonal_matrix(3, gen); + Tensor product; + product("i,k") = Q("i,j") * Q("k,j"); + + auto ones = make_tensor({3}, std::vector{1, 1, 1}); + auto ident = diagonal_matrix(ones); + REQUIRE(approximately_equal(product, ident, 1e-12)); + } + + SECTION("deterministic for fixed seed") { + auto gen1 = make_rng(99); + auto gen2 = make_rng(99); + auto Q1 = random_orthogonal_matrix(4, gen1); + auto Q2 = random_orthogonal_matrix(4, gen2); + REQUIRE(approximately_equal(Q1, Q2)); + } +}