Updates "Generate" to support multiple types.

include/tensorwrapper/generate/add_noise.hpp

@@ -17,7 +17,9 @@
 #pragma once
 #include <cstdint>
 #include <random>
+#include <tensorwrapper/concepts/floating_point.hpp>
 #include <tensorwrapper/tensor/tensor.hpp>
+#include <tensorwrapper/types/floating_point.hpp>
 
 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<concepts::FloatingPoint T>
 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<concepts::FloatingPoint T>
 Tensor add_noise(const Tensor& matrix, double t, std::uint64_t seed = 42);
 
+/** @brief Adds noise with element type `double`.
+ *
+ *  Equivalent to `add_noise<double>(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<double>(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<TYPE>(const Tensor& matrix, double t, \
+                                           std::mt19937& gen);             \
+    extern template Tensor add_noise<TYPE>(const Tensor& matrix, double t, \
+                                           std::uint64_t seed);
+
+TW_APPLY_FLOATING_POINT_TYPES(DECLARE_ADD_NOISE);
+
+#undef DECLARE_ADD_NOISE
+
 } // namespace tensorwrapper::generate

include/tensorwrapper/generate/generate_eigen_system.hpp

@@ -16,8 +16,10 @@
 
 #pragma once
 #include <cstddef>
+#include <tensorwrapper/concepts/floating_point.hpp>
 #include <tensorwrapper/generate/generate_utils.hpp>
 #include <tensorwrapper/tensor/tensor.hpp>
+#include <tensorwrapper/types/floating_point.hpp>
 
 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<concepts::FloatingPoint T>
 EigenSystem generate_eigen_system(const SymmetricMatrixSpec& spec);
 
+/** @brief Generates an eigen-system with element type `double`.
+ *
+ *  Equivalent to `generate_eigen_system<double>(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<TYPE>( \
+      const SymmetricMatrixSpec& spec);
+
+TW_APPLY_FLOATING_POINT_TYPES(DECLARE_GENERATE_EIGEN_SYSTEM);
+
+#undef DECLARE_GENERATE_EIGEN_SYSTEM
+
 } // namespace tensorwrapper::generate

include/tensorwrapper/generate/generate_eigenvalues.hpp

@@ -16,8 +16,10 @@
 
 #pragma once
 #include <random>
+#include <tensorwrapper/concepts/floating_point.hpp>
 #include <tensorwrapper/generate/generate_utils.hpp>
 #include <tensorwrapper/tensor/tensor.hpp>
+#include <tensorwrapper/types/floating_point.hpp>
 
 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<concepts::FloatingPoint T>
+Tensor generate_eigenvalues(const SymmetricMatrixSpec& spec, std::mt19937& gen);
+
+/** @brief Generates eigenvalues with element type `double`.
+ *
+ *  Equivalent to `generate_eigenvalues<double>(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<TYPE>( \
+      const SymmetricMatrixSpec& spec, std::mt19937& gen);
+
+TW_APPLY_FLOATING_POINT_TYPES(DECLARE_GENERATE_EIGENVALUES);
+
+#undef DECLARE_GENERATE_EIGENVALUES
+
 } // namespace tensorwrapper::generate

include/tensorwrapper/generate/identity_matrix.hpp

@@ -16,6 +16,7 @@
 
 #pragma once
 #include <cstddef>
+#include <tensorwrapper/concepts/floating_point.hpp>
 #include <tensorwrapper/generate/generate_utils.hpp>
 #include <tensorwrapper/utilities/diagonal_matrix.hpp>
 #include <tensorwrapper/utilities/make_tensor.hpp>
@@ -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<concepts::FloatingPoint T>
+Tensor identity_matrix(std::size_t n) {
     require_valid_n(n);
-    std::vector<double> values(n, 1.0);
-    auto values_tensor = utilities::make_tensor({n}, values);
-    return utilities::diagonal_matrix(values_tensor);
+    std::vector<T> 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<double>(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<double>(n);
 }
 
 } // namespace tensorwrapper::generate

include/tensorwrapper/generate/random_orthogonal_matrix.hpp

@@ -16,7 +16,9 @@
 
 #pragma once
 #include <random>
+#include <tensorwrapper/concepts/floating_point.hpp>
 #include <tensorwrapper/tensor/tensor.hpp>
+#include <tensorwrapper/types/floating_point.hpp>
 
 namespace tensorwrapper::generate {
 
@@ -25,13 +27,38 @@ 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.
  *
  *  @return A rank-2 tensor whose columns form an orthonormal basis.
  *
  *  @throw std::invalid_argument if @p n is outside the allowed range.
  */
+template<concepts::FloatingPoint T>
+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<double>(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<TYPE>(std::size_t n, \
+                                                          std::mt19937& gen);
+
+TW_APPLY_FLOATING_POINT_TYPES(DECLARE_RANDOM_ORTHOGONAL_MATRIX);
+
+#undef DECLARE_RANDOM_ORTHOGONAL_MATRIX
+
 } // namespace tensorwrapper::generate

src/tensorwrapper/generate/add_noise.cpp

@@ -20,6 +20,7 @@
 #include <tensorwrapper/buffer/contiguous.hpp>
 #include <tensorwrapper/generate/add_noise.hpp>
 #include <tensorwrapper/generate/generate_utils.hpp>
+#include <tensorwrapper/types/floating_point.hpp>
 #include <tensorwrapper/utilities/make_tensor.hpp>
 #include <vector>
 
@@ -32,27 +33,83 @@ std::vector<std::size_t> shape_extents(const auto& shape) {
     return rv;
 }
 
-} // namespace
+template<concepts::FloatingPoint T>
+double noise_center(const T& value) {
+    if constexpr(types::is_uncertain_v<T>) {
+        return static_cast<double>(value.mean());
+    } else if constexpr(types::is_interval_v<T>) {
+        return static_cast<double>(value.median());
+    } else {
+        return static_cast<double>(value);
+    }
+}
 
-Tensor add_noise(const Tensor& matrix, double t, std::mt19937& gen) {
+template<concepts::FloatingPoint T>
+T apply_noise(double center, double delta, double t) {
+    if constexpr(types::is_uncertain_v<T>) {
+        using value_type = typename T::value_t;
+        return T(static_cast<value_type>(center + delta),
+                 static_cast<value_type>(t));
+    } else if constexpr(types::is_interval_v<T>) {
+        using value_type   = typename T::value_t;
+        const value_type c = static_cast<value_type>(center + delta);
+        const value_type w = static_cast<value_type>(t);
+        return T(c - w, c + w);
+    } else {
+        return static_cast<T>(center + delta);
+    }
+}
+
+template<concepts::FloatingPoint T>
+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<const double>(in_buf);
-    std::vector<double> data(in_data.begin(), in_data.end());
+    auto in_data = buffer::get_raw_data<const T>(in_buf);
+    std::vector<T> data(in_data.begin(), in_data.end());
 
     if(t > 0.0) {
         std::normal_distribution<double> dist(0.0, t);
-        for(auto& x : data) { x += std::clamp(dist(gen), -t, t); }
+        for(auto& x : data) {
+            const auto center = noise_center<T>(x);
+            const auto delta  = std::clamp(dist(gen), -t, t);
+            x                 = apply_noise<T>(center, delta, t);
+        }
     }
 
     return utilities::make_tensor(shape_extents(in_buf.shape()), data.begin(),
                                   data.end());
 }
 
+} // namespace
+
+template<concepts::FloatingPoint T>
+Tensor add_noise(const Tensor& matrix, double t, std::mt19937& gen) {
+    return add_noise_impl<T>(matrix, t, gen);
+}
+
+template<concepts::FloatingPoint T>
 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<T>(matrix, t, gen);
+}
+
+Tensor add_noise(const Tensor& matrix, double t, std::mt19937& gen) {
+    return add_noise<double>(matrix, t, gen);
+}
+
+Tensor add_noise(const Tensor& matrix, double t, std::uint64_t seed) {
+    return add_noise<double>(matrix, t, seed);
 }
 
+#define DEFINE_ADD_NOISE(TYPE)                                      \
+    template Tensor add_noise<TYPE>(const Tensor& matrix, double t, \
+                                    std::mt19937& gen);             \
+    template Tensor add_noise<TYPE>(const Tensor& matrix, double t, \
+                                    std::uint64_t seed);
+
+TW_APPLY_FLOATING_POINT_TYPES(DEFINE_ADD_NOISE);
+
+#undef DEFINE_ADD_NOISE
+
 } // namespace tensorwrapper::generate

src/tensorwrapper/generate/generate_eigen_system.cpp

@@ -17,19 +17,21 @@
 #include <tensorwrapper/generate/generate_eigen_system.hpp>
 #include <tensorwrapper/generate/generate_eigenvalues.hpp>
 #include <tensorwrapper/generate/random_orthogonal_matrix.hpp>
+#include <tensorwrapper/types/floating_point.hpp>
 #include <tensorwrapper/utilities/diagonal_matrix.hpp>
 
 namespace tensorwrapper::generate {
 
+template<concepts::FloatingPoint T>
 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);
+    rv.eigenvalues  = generate_eigenvalues<T>(spec, gen);
+    rv.eigenvectors = random_orthogonal_matrix<T>(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<double>(spec);
+}
+
+#define DEFINE_GENERATE_EIGEN_SYSTEM(TYPE)            \
+    template EigenSystem generate_eigen_system<TYPE>( \
+      const SymmetricMatrixSpec& spec);
+
+TW_APPLY_FLOATING_POINT_TYPES(DEFINE_GENERATE_EIGEN_SYSTEM);
+
+#undef DEFINE_GENERATE_EIGEN_SYSTEM
+
 } // namespace tensorwrapper::generate