From eb0376a809417adc634373bd992de9a774a9a84c Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Sat, 30 May 2026 10:56:51 -0500 Subject: [PATCH 1/2] some updates to choose floating point type --- .../generate/generate_eigenvalues.hpp | 29 +++++ .../generate/identity_matrix.hpp | 21 +++- .../generate/random_orthogonal_matrix.hpp | 27 +++++ .../generate/generate_eigenvalues.cpp | 96 +++++++++++++---- .../generate/random_orthogonal_matrix.cpp | 36 +++++-- .../generate/generate_eigenvalues.cpp | 102 +++++++++--------- .../generate/identity_matrix.cpp | 19 ++-- .../generate/random_orthogonal_matrix.cpp | 24 +++-- 8 files changed, 254 insertions(+), 100 deletions(-) diff --git a/include/tensorwrapper/generate/generate_eigenvalues.hpp b/include/tensorwrapper/generate/generate_eigenvalues.hpp index 2ae3d46c..86377498 100644 --- a/include/tensorwrapper/generate/generate_eigenvalues.hpp +++ b/include/tensorwrapper/generate/generate_eigenvalues.hpp @@ -16,8 +16,10 @@ #pragma once #include +#include #include #include +#include namespace tensorwrapper::generate { @@ -28,6 +30,11 @@ namespace tensorwrapper::generate { * \texttt{spec.min\_eigenvalue} \times \texttt{spec.condition\_number}]@f$ * using the spacing strategy given by @p spec.spacing. * + * Explicit instantiations are provided only for types in + * @ref types::floating_point_types. + * + * @tparam T Element type of the returned tensor. + * * @param[in] spec Parameters controlling the eigenvalue distribution. * @param[in,out] gen Random number generator used when @p spec.spacing * requires random draws. @@ -38,6 +45,28 @@ namespace tensorwrapper::generate { * @throw std::invalid_argument if @p spec.n is outside `[1, kMaxMatrixDim]` or * if @p spec.spacing is invalid. */ +template +Tensor generate_eigenvalues(const SymmetricMatrixSpec& spec, std::mt19937& gen); + +/** @brief Generates eigenvalues with element type `double`. + * + * Equivalent to `generate_eigenvalues(spec, gen)`. + * + * @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. + */ Tensor generate_eigenvalues(const SymmetricMatrixSpec& spec, std::mt19937& gen); +#define DECLARE_GENERATE_EIGENVALUES(TYPE) \ + extern template Tensor generate_eigenvalues( \ + const SymmetricMatrixSpec& spec, std::mt19937& gen); + +TW_APPLY_FLOATING_POINT_TYPES(DECLARE_GENERATE_EIGENVALUES); + +#undef DECLARE_GENERATE_EIGENVALUES + } // namespace tensorwrapper::generate diff --git a/include/tensorwrapper/generate/identity_matrix.hpp b/include/tensorwrapper/generate/identity_matrix.hpp index 9e55f45e..6996cf0c 100644 --- a/include/tensorwrapper/generate/identity_matrix.hpp +++ b/include/tensorwrapper/generate/identity_matrix.hpp @@ -16,6 +16,7 @@ #pragma once #include +#include #include #include #include @@ -31,11 +32,23 @@ namespace tensorwrapper::generate { * * @throw std::invalid_argument if @p n is outside the allowed range. */ -inline Tensor identity_matrix(std::size_t n) { +template +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); + std::vector values(n, 1.0); + return utilities::diagonal_matrix(utilities::make_tensor({n}, values)); +} + +/** @brief Creates an identity matrix with element type `double`. + * + * Equivalent to `identity_matrix(n)`. + * + * @param[in] n Matrix dimension. Must be in `[1, kMaxMatrixDim]`. + * + * @return A rank-2 tensor with ones on the diagonal and zeros elsewhere. + */ +inline Tensor identity_matrix(std::size_t n) { + return identity_matrix(n); } } // namespace tensorwrapper::generate diff --git a/include/tensorwrapper/generate/random_orthogonal_matrix.hpp b/include/tensorwrapper/generate/random_orthogonal_matrix.hpp index 0be2a3f9..ecdb7a4e 100644 --- a/include/tensorwrapper/generate/random_orthogonal_matrix.hpp +++ b/include/tensorwrapper/generate/random_orthogonal_matrix.hpp @@ -16,7 +16,9 @@ #pragma once #include +#include #include +#include namespace tensorwrapper::generate { @@ -25,6 +27,11 @@ namespace tensorwrapper::generate { * Draws entries from a standard normal distribution and applies Householder QR * factorization to obtain an orthogonal matrix @f$Q@f$. * + * Explicit instantiations are provided only for types in + * @ref types::floating_point_types. + * + * @tparam T Element type of the returned tensor. + * * @param[in] n Matrix dimension. Must be in `[1, kMaxMatrixDim]`. * @param[in,out] gen Random number generator used for the normal draws. * @@ -32,6 +39,26 @@ namespace tensorwrapper::generate { * * @throw std::invalid_argument if @p n is outside the allowed range. */ +template +Tensor random_orthogonal_matrix(std::size_t n, std::mt19937& gen); + +/** @brief Creates a random orthogonal matrix with element type `double`. + * + * Equivalent to `random_orthogonal_matrix(n, gen)`. + * + * @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. + */ Tensor random_orthogonal_matrix(std::size_t n, std::mt19937& gen); +#define DECLARE_RANDOM_ORTHOGONAL_MATRIX(TYPE) \ + extern template Tensor random_orthogonal_matrix(std::size_t n, \ + std::mt19937& gen); + +TW_APPLY_FLOATING_POINT_TYPES(DECLARE_RANDOM_ORTHOGONAL_MATRIX); + +#undef DECLARE_RANDOM_ORTHOGONAL_MATRIX + } // namespace tensorwrapper::generate diff --git a/src/tensorwrapper/generate/generate_eigenvalues.cpp b/src/tensorwrapper/generate/generate_eigenvalues.cpp index 1a95eca1..5b4f722d 100644 --- a/src/tensorwrapper/generate/generate_eigenvalues.cpp +++ b/src/tensorwrapper/generate/generate_eigenvalues.cpp @@ -16,7 +16,9 @@ #include #include +#include #include +#include #include #include @@ -32,11 +34,12 @@ namespace { * @return A vector of length @p n with values from @p lambda_min to * @p lambda_max. */ +template 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); + std::vector values(n); for(std::size_t i = 0; i < n; ++i) { - values[i] = lambda_min + static_cast(i) * dx; + values[i] = static_cast(lambda_min + static_cast(i) * dx); } return values; } @@ -50,13 +53,14 @@ auto linear_spacing(std::size_t n, double lambda_min, double lambda_max) { * @return A vector of length @p n with values from @p lambda_min to * @p lambda_max on a log scale. */ +template 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); + const T log_min = types::log(static_cast(lambda_min)); + const T log_max = types::log(static_cast(lambda_max)); + const auto dlog = (log_max - log_min) / static_cast(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); + values[i] = types::exp(log_min + static_cast(i) * dlog); } return values; } @@ -75,12 +79,13 @@ auto logarithmic_spacing(std::size_t n, double lambda_min, double lambda_max) { * @return A vector of length @p n with eigenvalues assigned cyclically to * clusters. */ +template 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 values(n); std::vector cluster_centers(n_clusters); if(n_clusters == 1) { cluster_centers[0] = lambda_min; @@ -93,7 +98,7 @@ auto clustered_spacing(std::size_t n, std::size_t n_clusters, double lambda_min, for(std::size_t i = 0; i < n; ++i) { const auto cluster_id = i % n_clusters; - values[i] = cluster_centers[cluster_id] + jitter(gen); + values[i] = static_cast(cluster_centers[cluster_id] + jitter(gen)); } return values; } @@ -109,24 +114,59 @@ auto clustered_spacing(std::size_t n, std::size_t n_clusters, double lambda_min, * @return A vector of length @p n with eigenvalues assigned cyclically to * @p n_clusters plateaus. */ +template auto degenerate_spacing(std::size_t n, std::size_t n_clusters, double lambda_min, double lambda_max) { - std::vector values(n); + std::vector values(n); if(n_clusters <= 1) { - std::fill(values.begin(), values.end(), lambda_min); + std::fill(values.begin(), values.end(), static_cast(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; + values[i] = + static_cast(lambda_min + static_cast(plateau) * dx); } } return values; } + +template +double eigenvalue_sort_key(const T& value) { + if constexpr(types::is_interval_v) { + return static_cast(value.median()); + } else if constexpr(types::is_uncertain_v) { + return static_cast(value.mean()); + } else { + return static_cast(value); + } +} + +template +void sort_eigenvalues(std::vector& values) { + if constexpr(types::is_interval_v || types::is_uncertain_v) { + std::vector indices(values.size()); + std::iota(indices.begin(), indices.end(), 0); + std::sort(indices.begin(), indices.end(), + [&](std::size_t a, std::size_t b) { + return eigenvalue_sort_key(values[a]) < + eigenvalue_sort_key(values[b]); + }); + std::vector sorted(values.size()); + for(std::size_t i = 0; i < values.size(); ++i) { + sorted[i] = values[indices[i]]; + } + values = std::move(sorted); + } else { + std::sort(values.begin(), values.end()); + } +} + } // namespace +template Tensor generate_eigenvalues(const SymmetricMatrixSpec& spec, std::mt19937& gen) { require_valid_n(spec.n); @@ -135,29 +175,32 @@ Tensor generate_eigenvalues(const SymmetricMatrixSpec& spec, 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}, std::vector{static_cast(lambda_min)}); + } - if(n == 1) return utilities::make_tensor({n}, values); + std::vector values; switch(spec.spacing) { case EigenvalueSpacing::Linear: { - values = linear_spacing(n, lambda_min, lambda_max); + values = linear_spacing(n, lambda_min, lambda_max); break; } case EigenvalueSpacing::Logarithmic: { - values = logarithmic_spacing(n, lambda_min, lambda_max); + 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); + 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); + degenerate_spacing(n, spec.n_clusters, lambda_min, lambda_max); break; } default: { @@ -165,8 +208,21 @@ Tensor generate_eigenvalues(const SymmetricMatrixSpec& spec, } } - std::sort(values.begin(), values.end()); + sort_eigenvalues(values); return utilities::make_tensor({n}, values); } +Tensor generate_eigenvalues(const SymmetricMatrixSpec& spec, + std::mt19937& gen) { + return generate_eigenvalues(spec, gen); +} + +#define DEFINE_GENERATE_EIGENVALUES(TYPE) \ + template Tensor generate_eigenvalues( \ + const SymmetricMatrixSpec& spec, std::mt19937& gen); + +TW_APPLY_FLOATING_POINT_TYPES(DEFINE_GENERATE_EIGENVALUES); + +#undef DEFINE_GENERATE_EIGENVALUES + } // namespace tensorwrapper::generate diff --git a/src/tensorwrapper/generate/random_orthogonal_matrix.cpp b/src/tensorwrapper/generate/random_orthogonal_matrix.cpp index 3822c62f..491d3ea1 100644 --- a/src/tensorwrapper/generate/random_orthogonal_matrix.cpp +++ b/src/tensorwrapper/generate/random_orthogonal_matrix.cpp @@ -17,14 +17,15 @@ #include #include #include +#include #include #include namespace tensorwrapper::generate { +namespace { -Tensor random_orthogonal_matrix(std::size_t n, std::mt19937& gen) { - require_valid_n(n); - +Eigen::MatrixXd random_orthogonal_matrix_eigen(std::size_t n, + std::mt19937& gen) { Eigen::MatrixXd M(n, n); std::normal_distribution dist(0.0, 1.0); for(std::size_t i = 0; i < n; ++i) { @@ -34,16 +35,37 @@ Tensor random_orthogonal_matrix(std::size_t n, std::mt19937& gen) { } } Eigen::HouseholderQR qr(M); - const Eigen::MatrixXd Q = qr.householderQ(); + return Eigen::MatrixXd(qr.householderQ()); +} + +} // namespace + +template +Tensor random_orthogonal_matrix(std::size_t n, std::mt19937& gen) { + require_valid_n(n); - std::vector data(n * n); + const auto Q = random_orthogonal_matrix_eigen(n, gen); + + 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)); + data[i * n + j] = static_cast( + Q(static_cast(i), static_cast(j))); } } return utilities::make_tensor({n, n}, data); } +Tensor random_orthogonal_matrix(std::size_t n, std::mt19937& gen) { + return random_orthogonal_matrix(n, gen); +} + +#define DEFINE_RANDOM_ORTHOGONAL_MATRIX(TYPE) \ + template Tensor random_orthogonal_matrix(std::size_t n, \ + std::mt19937& gen); + +TW_APPLY_FLOATING_POINT_TYPES(DEFINE_RANDOM_ORTHOGONAL_MATRIX); + +#undef DEFINE_RANDOM_ORTHOGONAL_MATRIX + } // namespace tensorwrapper::generate diff --git a/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigenvalues.cpp b/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigenvalues.cpp index 98cc4e0f..ebac9fdb 100644 --- a/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigenvalues.cpp +++ b/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigenvalues.cpp @@ -14,42 +14,29 @@ * limitations under the License. */ -#include -#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})); +TEMPLATE_LIST_TEST_CASE("generate_eigenvalues", "", + types::floating_point_types) { + SECTION("n == 1") { + SymmetricMatrixSpec spec; + spec.n = 1; + spec.min_eigenvalue = 3.5; + auto gen = make_rng(1); + auto result = generate_eigenvalues(spec, gen); + auto corr = make_tensor({1}, std::vector{TestType{3.5}}); + REQUIRE(approximately_equal(result, corr)); } - return rv; -} -} // namespace -TEST_CASE("generate_eigenvalues") { SECTION("linear spacing") { SymmetricMatrixSpec spec; spec.n = 4; @@ -57,8 +44,10 @@ TEST_CASE("generate_eigenvalues") { 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}); + auto result = generate_eigenvalues(spec, gen); + auto corr = + make_tensor({4}, std::vector{TestType{1}, TestType{4}, + TestType{7}, TestType{10}}); REQUIRE(approximately_equal(result, corr)); } @@ -68,13 +57,20 @@ TEST_CASE("generate_eigenvalues") { 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())); + const auto lambda_min = static_cast(spec.min_eigenvalue); + const auto lambda_max = + static_cast(spec.min_eigenvalue * spec.condition_number); + const auto log_min = types::log(lambda_min); + const auto log_max = types::log(lambda_max); + const auto dlog = (log_max - log_min) / TestType{2}; + std::vector expected(3); + for(std::size_t i = 0; i < 3; ++i) { + expected[i] = types::exp(log_min + static_cast(i) * dlog); + } + auto gen = make_rng(1); + auto result = generate_eigenvalues(spec, gen); + auto corr = make_tensor({3}, expected); + REQUIRE(approximately_equal(result, corr)); } SECTION("degenerate spacing") { @@ -85,9 +81,11 @@ TEST_CASE("generate_eigenvalues") { 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)); } + auto result = generate_eigenvalues(spec, gen); + auto corr = + make_tensor({3}, std::vector{TestType{2.5}, TestType{2.5}, + TestType{2.5}}); + REQUIRE(approximately_equal(result, corr)); } SECTION("clustered spacing") { @@ -99,31 +97,27 @@ TEST_CASE("generate_eigenvalues") { 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); - } + auto gen = make_rng(spec.seed); + auto result = generate_eigenvalues(spec, gen); + + gen = make_rng(spec.seed); + auto result_copy = generate_eigenvalues(spec, gen); + REQUIRE(approximately_equal(result, result_copy)); + + SymmetricMatrixSpec plateau_spec = spec; + plateau_spec.spacing = EigenvalueSpacing::Degenerate; + auto plateau_gen = make_rng(1); + auto plateaus = + generate_eigenvalues(plateau_spec, plateau_gen); + REQUIRE(approximately_equal(result, plateaus, spec.cluster_width)); } SECTION("invalid n throws") { SymmetricMatrixSpec spec; spec.n = 0; auto gen = make_rng(1); - REQUIRE_THROWS_AS(generate_eigenvalues(spec, gen), + REQUIRE_THROWS_AS(generate_eigenvalues(spec, gen), std::invalid_argument); } } diff --git a/tests/cxx/unit_tests/tensorwrapper/generate/identity_matrix.cpp b/tests/cxx/unit_tests/tensorwrapper/generate/identity_matrix.cpp index c092a8de..d4cfa875 100644 --- a/tests/cxx/unit_tests/tensorwrapper/generate/identity_matrix.cpp +++ b/tests/cxx/unit_tests/tensorwrapper/generate/identity_matrix.cpp @@ -16,6 +16,7 @@ #include #include +#include #include #include @@ -24,21 +25,25 @@ using namespace tensorwrapper::generate; using namespace tensorwrapper::operations; using namespace tensorwrapper::utilities; -TEST_CASE("identity_matrix") { +TEMPLATE_LIST_TEST_CASE("identity_matrix", "", types::floating_point_types) { SECTION("2 by 2") { - auto result = identity_matrix(2); - auto corr = make_tensor({2, 2}, std::vector{1, 0, 0, 1}); + auto result = identity_matrix(2); + auto corr = + make_tensor({2, 2}, std::vector{TestType{1}, TestType{0}, + TestType{0}, TestType{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}); + auto result = identity_matrix(3); + auto corr = make_tensor( + {3, 3}, std::vector{TestType{1}, TestType{0}, TestType{0}, + TestType{0}, TestType{1}, TestType{0}, + TestType{0}, TestType{0}, TestType{1}}); REQUIRE(approximately_equal(result, corr)); } SECTION("invalid n throws") { - REQUIRE_THROWS_AS(identity_matrix(0), std::invalid_argument); + 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 index 6761c082..44b63936 100644 --- a/tests/cxx/unit_tests/tensorwrapper/generate/random_orthogonal_matrix.cpp +++ b/tests/cxx/unit_tests/tensorwrapper/generate/random_orthogonal_matrix.cpp @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include @@ -27,10 +28,11 @@ using namespace tensorwrapper::generate; using namespace tensorwrapper::operations; using namespace tensorwrapper::utilities; -TEST_CASE("random_orthogonal_matrix") { +TEMPLATE_LIST_TEST_CASE("random_orthogonal_matrix", "", + types::floating_point_types) { SECTION("shape") { auto gen = make_rng(7); - auto Q = random_orthogonal_matrix(3, gen); + auto Q = random_orthogonal_matrix(3, gen); auto buf = make_contiguous(Q.buffer()); REQUIRE(buf.shape().rank() == 2); REQUIRE(buf.shape().extent(0) == 3); @@ -39,20 +41,26 @@ TEST_CASE("random_orthogonal_matrix") { SECTION("orthogonality") { auto gen = make_rng(11); - auto Q = random_orthogonal_matrix(3, gen); + 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)); + auto ones = make_tensor( + {3}, std::vector{TestType{1}, TestType{1}, TestType{1}}); + auto ident = diagonal_matrix(ones); + constexpr double tol = std::is_same_v || + std::is_same_v || + std::is_same_v ? + 1e-5 : + 1e-12; + REQUIRE(approximately_equal(product, ident, tol)); } 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); + auto Q1 = random_orthogonal_matrix(4, gen1); + auto Q2 = random_orthogonal_matrix(4, gen2); REQUIRE(approximately_equal(Q1, Q2)); } } From b545cccc923a07466eeb61c9e9ef4472c1b53cb0 Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Sat, 30 May 2026 11:07:10 -0500 Subject: [PATCH 2/2] completes refactor of generate --- include/tensorwrapper/generate/add_noise.hpp | 36 ++++++++++ .../generate/generate_eigen_system.hpp | 28 ++++++++ src/tensorwrapper/generate/add_noise.cpp | 69 ++++++++++++++++-- .../generate/generate_eigen_system.cpp | 18 ++++- .../tensorwrapper/generate/add_noise.cpp | 72 ++++++++++++++----- .../generate/generate_eigen_system.cpp | 53 +++++++++----- 6 files changed, 233 insertions(+), 43 deletions(-) diff --git a/include/tensorwrapper/generate/add_noise.hpp b/include/tensorwrapper/generate/add_noise.hpp index ee10a448..f4215048 100644 --- a/include/tensorwrapper/generate/add_noise.hpp +++ b/include/tensorwrapper/generate/add_noise.hpp @@ -17,7 +17,9 @@ #pragma once #include #include +#include #include +#include namespace tensorwrapper::generate { @@ -26,6 +28,14 @@ namespace tensorwrapper::generate { * 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. * + * For sigma uncertain types the result has standard deviation @p t. For sigma + * interval types the result has radius @p t about the perturbed center. + * + * Explicit instantiations are provided only for types in + * @ref types::floating_point_types. + * + * @tparam T Element type of the returned tensor. + * * @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. @@ -34,9 +44,12 @@ namespace tensorwrapper::generate { * * @throw std::invalid_argument if @p t is negative. */ +template Tensor add_noise(const Tensor& matrix, double t, std::mt19937& gen); /** @brief Overload of add_noise that creates its own RNG from @p seed. + * + * @tparam T Element type of the returned tensor. * * @param[in] matrix The tensor to perturb. * @param[in] t Non-negative noise scale (standard deviation and clamp bound). @@ -47,6 +60,29 @@ Tensor add_noise(const Tensor& matrix, double t, std::mt19937& gen); * * @throw std::invalid_argument if @p t is negative. */ +template Tensor add_noise(const Tensor& matrix, double t, std::uint64_t seed = 42); +/** @brief Adds noise with element type `double`. + * + * Equivalent to `add_noise(matrix, t, gen)`. + */ +Tensor add_noise(const Tensor& matrix, double t, std::mt19937& gen); + +/** @brief Adds noise with element type `double` using @p seed. + * + * Equivalent to `add_noise(matrix, t, seed)`. + */ +Tensor add_noise(const Tensor& matrix, double t, std::uint64_t seed = 42); + +#define DECLARE_ADD_NOISE(TYPE) \ + extern template Tensor add_noise(const Tensor& matrix, double t, \ + std::mt19937& gen); \ + extern template Tensor add_noise(const Tensor& matrix, double t, \ + std::uint64_t seed); + +TW_APPLY_FLOATING_POINT_TYPES(DECLARE_ADD_NOISE); + +#undef DECLARE_ADD_NOISE + } // namespace tensorwrapper::generate diff --git a/include/tensorwrapper/generate/generate_eigen_system.hpp b/include/tensorwrapper/generate/generate_eigen_system.hpp index cbaeb653..c78306c4 100644 --- a/include/tensorwrapper/generate/generate_eigen_system.hpp +++ b/include/tensorwrapper/generate/generate_eigen_system.hpp @@ -16,8 +16,10 @@ #pragma once #include +#include #include #include +#include namespace tensorwrapper::generate { @@ -43,6 +45,11 @@ struct EigenSystem { * 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. * + * Explicit instantiations are provided only for types in + * @ref types::floating_point_types. + * + * @tparam T Element type of the returned tensors. + * * @param[in] spec Parameters controlling the matrix dimension, spectrum, and * random seed. * @@ -51,6 +58,27 @@ struct EigenSystem { * * @throw std::invalid_argument if @p spec.n is outside `[1, kMaxMatrixDim]`. */ +template EigenSystem generate_eigen_system(const SymmetricMatrixSpec& spec); +/** @brief Generates an eigen-system with element type `double`. + * + * Equivalent to `generate_eigen_system(spec)`. + * + * @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. + */ +EigenSystem generate_eigen_system(const SymmetricMatrixSpec& spec); + +#define DECLARE_GENERATE_EIGEN_SYSTEM(TYPE) \ + extern template EigenSystem generate_eigen_system( \ + const SymmetricMatrixSpec& spec); + +TW_APPLY_FLOATING_POINT_TYPES(DECLARE_GENERATE_EIGEN_SYSTEM); + +#undef DECLARE_GENERATE_EIGEN_SYSTEM + } // namespace tensorwrapper::generate diff --git a/src/tensorwrapper/generate/add_noise.cpp b/src/tensorwrapper/generate/add_noise.cpp index f6cf6b08..1c27808b 100644 --- a/src/tensorwrapper/generate/add_noise.cpp +++ b/src/tensorwrapper/generate/add_noise.cpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -32,27 +33,83 @@ std::vector shape_extents(const auto& shape) { return rv; } -} // namespace +template +double noise_center(const T& value) { + if constexpr(types::is_uncertain_v) { + return static_cast(value.mean()); + } else if constexpr(types::is_interval_v) { + return static_cast(value.median()); + } else { + return static_cast(value); + } +} -Tensor add_noise(const Tensor& matrix, double t, std::mt19937& gen) { +template +T apply_noise(double center, double delta, double t) { + if constexpr(types::is_uncertain_v) { + using value_type = typename T::value_t; + return T(static_cast(center + delta), + static_cast(t)); + } else if constexpr(types::is_interval_v) { + using value_type = typename T::value_t; + const value_type c = static_cast(center + delta); + const value_type w = static_cast(t); + return T(c - w, c + w); + } else { + return static_cast(center + delta); + } +} + +template +Tensor add_noise_impl(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()); + 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); } + for(auto& x : data) { + const auto center = noise_center(x); + const auto delta = std::clamp(dist(gen), -t, t); + x = apply_noise(center, delta, t); + } } return utilities::make_tensor(shape_extents(in_buf.shape()), data.begin(), data.end()); } +} // namespace + +template +Tensor add_noise(const Tensor& matrix, double t, std::mt19937& gen) { + return add_noise_impl(matrix, t, gen); +} + +template Tensor add_noise(const Tensor& matrix, double t, std::uint64_t seed) { auto gen = make_rng(seed); - return add_noise(matrix, t, gen); + return add_noise(matrix, t, gen); +} + +Tensor add_noise(const Tensor& matrix, double t, std::mt19937& gen) { + return add_noise(matrix, t, gen); +} + +Tensor add_noise(const Tensor& matrix, double t, std::uint64_t seed) { + return add_noise(matrix, t, seed); } +#define DEFINE_ADD_NOISE(TYPE) \ + template Tensor add_noise(const Tensor& matrix, double t, \ + std::mt19937& gen); \ + template Tensor add_noise(const Tensor& matrix, double t, \ + std::uint64_t seed); + +TW_APPLY_FLOATING_POINT_TYPES(DEFINE_ADD_NOISE); + +#undef DEFINE_ADD_NOISE + } // namespace tensorwrapper::generate diff --git a/src/tensorwrapper/generate/generate_eigen_system.cpp b/src/tensorwrapper/generate/generate_eigen_system.cpp index c3b67d7f..d308aafa 100644 --- a/src/tensorwrapper/generate/generate_eigen_system.cpp +++ b/src/tensorwrapper/generate/generate_eigen_system.cpp @@ -17,10 +17,12 @@ #include #include #include +#include #include namespace tensorwrapper::generate { +template EigenSystem generate_eigen_system(const SymmetricMatrixSpec& spec) { require_valid_n(spec.n); auto gen = make_rng(spec.seed); @@ -28,8 +30,8 @@ EigenSystem generate_eigen_system(const SymmetricMatrixSpec& spec) { EigenSystem rv; rv.n = n; - rv.eigenvalues = generate_eigenvalues(spec, gen); - rv.eigenvectors = random_orthogonal_matrix(n, gen); + rv.eigenvalues = generate_eigenvalues(spec, gen); + rv.eigenvectors = random_orthogonal_matrix(n, gen); const auto D = utilities::diagonal_matrix(rv.eigenvalues); @@ -42,4 +44,16 @@ EigenSystem generate_eigen_system(const SymmetricMatrixSpec& spec) { return rv; } +EigenSystem generate_eigen_system(const SymmetricMatrixSpec& spec) { + return generate_eigen_system(spec); +} + +#define DEFINE_GENERATE_EIGEN_SYSTEM(TYPE) \ + template EigenSystem generate_eigen_system( \ + const SymmetricMatrixSpec& spec); + +TW_APPLY_FLOATING_POINT_TYPES(DEFINE_GENERATE_EIGEN_SYSTEM); + +#undef DEFINE_GENERATE_EIGEN_SYSTEM + } // 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 index 9aaca48f..fd56794f 100644 --- a/tests/cxx/unit_tests/tensorwrapper/generate/add_noise.cpp +++ b/tests/cxx/unit_tests/tensorwrapper/generate/add_noise.cpp @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include @@ -28,51 +29,85 @@ 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)); +template +double elem_center(const T& value) { + if constexpr(types::is_uncertain_v) { + return static_cast(value.mean()); + } else if constexpr(types::is_interval_v) { + return static_cast(value.median()); + } else { + return static_cast(value); + } +} + +template +void require_within_noise(const T& in, const T& out, double t) { + if constexpr(types::is_uncertain_v) { + REQUIRE(out.sd() == Catch::Approx(t)); + REQUIRE(std::abs(out.mean() - in.mean()) <= t); + } else if constexpr(types::is_interval_v) { + REQUIRE(out.radius() == Catch::Approx(t)); + REQUIRE(std::abs(out.median() - in.median()) <= t); + } else { + REQUIRE(std::abs(elem_center(out) - elem_center(in)) <= t); } } } // namespace -TEST_CASE("add_noise") { - const auto matrix = make_tensor({2, 2}, std::vector{1, 2, 3, 4}); +TEMPLATE_LIST_TEST_CASE("add_noise", "", types::floating_point_types) { + const auto matrix = + make_tensor({2, 2}, std::vector{TestType{1}, TestType{2}, + TestType{3}, TestType{4}}); SECTION("t equals zero returns copy") { auto gen = make_rng(1); - auto out = add_noise(matrix, 0.0, gen); + 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 out = add_noise(matrix, t, gen); + using wtf::fp::float_cast; 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); + const auto in_val = + float_cast(in_buf.get_elem({i, j})); + const auto out_val = + float_cast(out_buf.get_elem({i, j})); + require_within_noise(in_val, out_val, 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)); + auto out1 = add_noise(matrix, t, 99); + auto out2 = add_noise(matrix, t, 99); + if constexpr(types::is_interval_v) { + using wtf::fp::float_cast; + auto b1 = make_contiguous(out1.buffer()); + auto b2 = make_contiguous(out2.buffer()); + for(std::size_t i = 0; i < 2; ++i) { + for(std::size_t j = 0; j < 2; ++j) { + const auto v1 = float_cast(b1.get_elem({i, j})); + const auto v2 = float_cast(b2.get_elem({i, j})); + REQUIRE(v1.lower() == Catch::Approx(v2.lower())); + REQUIRE(v1.upper() == Catch::Approx(v2.upper())); + } + } + } else { + REQUIRE(approximately_equal(out1, out2)); + } } SECTION("shape preserved") { auto gen = make_rng(3); - auto out = add_noise(matrix, 0.01, gen); + 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()); @@ -82,6 +117,7 @@ TEST_CASE("add_noise") { SECTION("invalid t throws") { auto gen = make_rng(1); - REQUIRE_THROWS_AS(add_noise(matrix, -1.0, gen), std::invalid_argument); + 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 index 3b72828b..805dad30 100644 --- a/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigen_system.cpp +++ b/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigen_system.cpp @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include @@ -27,14 +28,24 @@ using namespace tensorwrapper::generate; using namespace tensorwrapper::operations; using namespace tensorwrapper::utilities; -TEST_CASE("generate_eigen_system") { +namespace { +template +constexpr double eigen_system_tol = + std::is_same_v || std::is_same_v || + std::is_same_v ? + 1e-5 : + 1e-12; +} // namespace + +TEMPLATE_LIST_TEST_CASE("generate_eigen_system", "", + types::floating_point_types) { SECTION("shapes") { SymmetricMatrixSpec spec; spec.n = 4; spec.condition_number = 1e3; spec.spacing = EigenvalueSpacing::Linear; spec.seed = 11; - auto system = generate_eigen_system(spec); + auto system = generate_eigen_system(spec); REQUIRE(system.n == 4); REQUIRE( @@ -51,35 +62,37 @@ TEST_CASE("generate_eigen_system") { SymmetricMatrixSpec spec; spec.n = 3; spec.seed = 5; - auto system = generate_eigen_system(spec); + 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)); + REQUIRE(approximately_equal(sym, tran, eigen_system_tol)); } SECTION("orthonormal eigenvectors") { SymmetricMatrixSpec spec; spec.n = 3; spec.seed = 7; - auto system = generate_eigen_system(spec); + 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 ones = make_tensor( + {3}, std::vector{TestType{1}, TestType{1}, TestType{1}}); auto ident = diagonal_matrix(ones); - REQUIRE(approximately_equal(product, ident, 1e-12)); + REQUIRE( + approximately_equal(product, ident, eigen_system_tol)); } SECTION("exact eigenpair") { SymmetricMatrixSpec spec; spec.n = 4; spec.seed = 13; - auto system = generate_eigen_system(spec); + auto system = generate_eigen_system(spec); Tensor av; av("i,k") = system.matrix("i,j") * system.eigenvectors("j,k"); @@ -87,19 +100,20 @@ TEST_CASE("generate_eigen_system") { 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)); + REQUIRE(approximately_equal(av, scaled, eigen_system_tol)); } 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); + 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)); + REQUIRE(approximately_equal(system1.matrix, system2.matrix, + eigen_system_tol)); } SECTION("degenerate") { @@ -110,18 +124,23 @@ TEST_CASE("generate_eigen_system") { spec.spacing = EigenvalueSpacing::Degenerate; spec.n_clusters = 1; spec.seed = 7; - auto system = generate_eigen_system(spec); + auto system = generate_eigen_system(spec); - auto evals = make_tensor({2}, std::vector{1.0, 1.0}); + auto evals = + make_tensor({2}, std::vector{TestType{1.0}, TestType{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)); + auto ident = + make_tensor({2, 2}, std::vector{TestType{1}, TestType{0}, + TestType{0}, TestType{1}}); + REQUIRE(approximately_equal(system.matrix, ident, + eigen_system_tol)); } SECTION("invalid n throws") { SymmetricMatrixSpec spec; spec.n = 0; - REQUIRE_THROWS_AS(generate_eigen_system(spec), std::invalid_argument); + REQUIRE_THROWS_AS(generate_eigen_system(spec), + std::invalid_argument); } }