Affine

include/tensorwrapper/generate/identity_matrix.hpp

@@ -35,7 +35,7 @@ namespace tensorwrapper::generate {
 template<concepts::FloatingPoint T>
 Tensor identity_matrix(std::size_t n) {
     require_valid_n(n);
-    std::vector<T> values(n, 1.0);
+    std::vector<T> values(n, T{1});
     return utilities::diagonal_matrix(utilities::make_tensor({n}, values));
 }
 

include/tensorwrapper/types/floating_point.hpp

@@ -29,13 +29,25 @@ template<typename T>
 using uncertain_type = sigma::Uncertain<T>;
 using ufloat         = uncertain_type<float>;
 using udouble        = uncertain_type<double>;
+
 template<typename T>
 using interval_type = sigma::Interval<T>;
 using ifloat        = interval_type<float>;
 using idouble       = interval_type<double>;
 
+template<typename T>
+using affine_type = sigma::Affine<T>;
+using afloat      = affine_type<float>;
+using adouble     = affine_type<double>;
+
+template<typename T>
+using thresholded_affine_type = sigma::ThresholdedAffine<T>;
+using tafloat                 = thresholded_affine_type<float>;
+using tadouble                = thresholded_affine_type<double>;
+
 using floating_point_types =
-  std::tuple<float, double, ufloat, udouble, ifloat, idouble>;
+  std::tuple<float, double, ufloat, udouble, ifloat, idouble, afloat, adouble,
+             tafloat, tadouble>;
 
 template<typename T>
 constexpr bool is_uncertain_v =
@@ -46,44 +58,98 @@ constexpr bool is_interval_v =
   std::is_same_v<T, ifloat> || std::is_same_v<T, idouble>;
 
 template<typename T>
-T fabs(T value) {
+constexpr bool is_affine_v =
+  std::is_same_v<T, afloat> || std::is_same_v<T, adouble>;
+
+template<typename T>
+constexpr bool is_thresholded_affine_v =
+  std::is_same_v<T, tafloat> || std::is_same_v<T, tadouble>;
+
+template<typename T>
+constexpr bool is_uq_type_v = is_uncertain_v<T> || is_interval_v<T> ||
+                              is_affine_v<T> || is_thresholded_affine_v<T>;
+
+template<typename T, typename U, typename V>
+T construct_uq_type(const U& center, const V& radius) {
     if constexpr(is_uncertain_v<T>) {
-        return sigma::fabs(value);
+        return T(center, radius);
     } else if constexpr(is_interval_v<T>) {
-        return T(sigma::fabs(value));
+        return T(center - radius, center + radius);
+    } else if constexpr(is_affine_v<T> || is_thresholded_affine_v<T>) {
+        return T(center - radius, center + radius);
+    } else if constexpr(is_uq_type_v<T>) {
+        throw std::logic_error("UQ type not recognized in construct_uq_type.");
     } else {
-        return std::fabs(value);
+        return T(center + radius);
     }
 }
 
 template<typename T>
-T log(T value) {
+auto uq_center(const T& value) {
     if constexpr(is_uncertain_v<T>) {
-        return sigma::log(value);
+        return value.mean();
     } else if constexpr(is_interval_v<T>) {
-        return T(sigma::log(value));
+        return value.median();
+    } else if constexpr(is_affine_v<T> || is_thresholded_affine_v<T>) {
+        return value.center();
+    } else if constexpr(is_uq_type_v<T>) {
+        throw std::logic_error("UQ type not recognized in uq_center.");
     } else {
-        return std::log(value);
+        return value;
     }
 }
 
 template<typename T>
-T exp(T value) {
+auto uq_upper(const T& value) {
     if constexpr(is_uncertain_v<T>) {
-        return sigma::exp(value);
+        return value.mean() + value.sd();
     } else if constexpr(is_interval_v<T>) {
-        return T(sigma::exp(value));
+        return value.upper();
+    } else if constexpr(is_affine_v<T> || is_thresholded_affine_v<T>) {
+        return value.range().upper();
+    } else if constexpr(is_uq_type_v<T>) {
+        throw std::logic_error("UQ type not recognized in uq_upper.");
+    } else {
+        return value;
+    }
+}
+
+template<typename T, typename U>
+bool strictly_less(const T& lhs, const U& rhs) {
+    return uq_upper(lhs) < uq_upper(rhs);
+}
+
+template<typename T>
+T fabs(T value) {
+    if constexpr(is_uq_type_v<T>) {
+        return static_cast<T>(sigma::fabs(value));
+    } else {
+        return std::fabs(value);
+    }
+}
+
+template<typename T>
+T log(T value) {
+    if constexpr(is_uq_type_v<T>) {
+        return static_cast<T>(sigma::log(value));
+    } else {
+        return std::log(value);
+    }
+}
+
+template<typename T>
+T exp(T value) {
+    if constexpr(is_uq_type_v<T>) {
+        return static_cast<T>(sigma::exp(value));
     } else {
         return std::exp(value);
     }
 }
 
 template<typename T>
 T pow(T value, double pow) {
-    if constexpr(is_uncertain_v<T>) {
-        return sigma::pow(value, pow);
-    } else if constexpr(is_interval_v<T>) {
-        return T(sigma::pow(value, pow));
+    if constexpr(is_uq_type_v<T>) {
+        return static_cast<T>(sigma::pow(value, pow));
     } else {
         return std::pow(value, pow);
     }
@@ -95,13 +161,21 @@ T pow(T value, double pow) {
     MACRO_IN(tensorwrapper::types::ufloat);     \
     MACRO_IN(tensorwrapper::types::udouble);    \
     MACRO_IN(tensorwrapper::types::ifloat);     \
-    MACRO_IN(tensorwrapper::types::idouble);
+    MACRO_IN(tensorwrapper::types::idouble);    \
+    MACRO_IN(tensorwrapper::types::afloat);     \
+    MACRO_IN(tensorwrapper::types::adouble);    \
+    MACRO_IN(tensorwrapper::types::tafloat);    \
+    MACRO_IN(tensorwrapper::types::tadouble)
 } // namespace tensorwrapper::types
 
 WTF_REGISTER_FP_TYPE(tensorwrapper::types::ufloat);
 WTF_REGISTER_FP_TYPE(tensorwrapper::types::udouble);
 WTF_REGISTER_FP_TYPE(tensorwrapper::types::ifloat);
 WTF_REGISTER_FP_TYPE(tensorwrapper::types::idouble);
+WTF_REGISTER_FP_TYPE(tensorwrapper::types::afloat);
+WTF_REGISTER_FP_TYPE(tensorwrapper::types::adouble);
+WTF_REGISTER_FP_TYPE(tensorwrapper::types::tafloat);
+WTF_REGISTER_FP_TYPE(tensorwrapper::types::tadouble);
 
 #else
 template<typename T>
@@ -112,6 +186,14 @@ template<typename T>
 using interval_type = T;
 using ifloat        = float;
 using idouble       = double;
+template<typename T>
+using affine_type = T;
+using afloat      = float;
+using adouble     = double;
+template<typename T>
+using thresholded_affine_type = T;
+using tafloat                 = float;
+using tadouble                = double;
 
 using floating_point_types = std::tuple<float, double>;
 
@@ -121,6 +203,36 @@ constexpr bool is_uncertain_v = false;
 template<typename T>
 constexpr bool is_interval_v = false;
 
+template<typename T>
+constexpr bool is_affine_v = false;
+
+template<typename T>
+constexpr bool is_thresholded_affine_v = false;
+
+template<typename T>
+constexpr bool is_uq_type_v = false;
+
+template<typename T, typename U, typename V>
+T construct_uq_type(const U& center, const V&) {
+    return T(center);
+}
+
+template<typename T>
+T uq_center(const T& value) {
+    return value;
+}
+
+template<typename T>
+T uq_upper(const T& value) {
+    return value;
+}
+
+template<typename T, typename U>
+bool strictly_less(const T& lhs, const U& rhs) {
+    return lhs < rhs;
+    ;
+}
+
 template<typename T>
 T fabs(T value) {
     return std::fabs(value);

src/tensorwrapper/buffer/detail_/hash_utilities.hpp

@@ -60,10 +60,18 @@ template<typename T>
 void hash_input(hash_type& seed, const sigma::Uncertain<T>& value) {
     hash_input(seed, value.mean());
     hash_input(seed, value.sd());
+    // deps() is an unordered container, so its iteration order is not stable
+    // (e.g. across copies of the same value). Fold each dependency into an
+    // order-independent (commutative) accumulator before combining so that two
+    // equal values always produce the same hash.
+    hash_type deps_hash = 0;
     for(const auto& [dep, deriv] : value.deps()) {
-        hash_input(seed, dep);
-        hash_input(seed, deriv);
+        hash_type term = 0;
+        hash_input(term, dep);
+        hash_input(term, deriv);
+        deps_hash ^= term;
     }
+    hash_input(seed, deps_hash);
 }
 
 template<typename T>
@@ -72,6 +80,30 @@ void hash_input(hash_type& seed, const types::interval_type<T>& value) {
     hash_input(seed, value.upper());
 }
 
+template<typename T>
+void hash_input(hash_type& seed, const types::affine_type<T>& value) {
+    hash_input(seed, value.center());
+    // error_terms() is an unordered_map, so its iteration order is not stable
+    // (e.g. across copies of the same value). Fold each error term into an
+    // order-independent (commutative) accumulator before combining so that two
+    // equal affine forms always produce the same hash.
+    hash_type terms_hash = 0;
+    for(const auto& [label, coeff] : value.error_terms()) {
+        hash_type term = 0;
+        hash_input(term, label);
+        hash_input(term, coeff);
+        terms_hash ^= term;
+    }
+    hash_input(seed, terms_hash);
+}
+
+template<typename T>
+void hash_input(hash_type& seed,
+                const types::thresholded_affine_type<T>& value) {
+    hash_input(seed, value.affine());
+    hash_input(seed, value.threshold());
+}
+
 #endif
 
 class HashVisitor {

src/tensorwrapper/buffer/detail_/unary_operation_visitor.hpp

@@ -137,7 +137,7 @@ class ScalarMultiplicationVisitor : public UnaryOperationVisitor {
         // TODO: Change when public API changes to support other FP types
         auto scalar = wtf::fp::float_cast<double>(m_scalar_);
         pthis->scalar_multiplication(this->this_labels(), other_labels(),
-                                     scalar, *pother);
+                                     clean_t(scalar), *pother);
     }
 
 private:
@@ -150,18 +150,9 @@ class ApproximatelyEqualVisitor {
 
     template<typename FloatType>
     bool operator()(const std::span<FloatType> result) {
-        using clean_t = std::decay_t<FloatType>;
         for(std::size_t i = 0; i < result.size(); ++i) {
             auto abs_diff = types::fabs(result[i]);
-            double scalar_diff;
-            if constexpr(types::is_uncertain_v<clean_t>) {
-                scalar_diff = abs_diff.mean();
-            } else if constexpr(types::is_interval_v<clean_t>) {
-                scalar_diff = abs_diff.upper();
-            } else {
-                scalar_diff = abs_diff;
-            }
-            if(scalar_diff >= m_tol_) return false;
+            if(!types::strictly_less(abs_diff, m_tol_)) { return false; }
         }
         return true;
     }
@@ -176,11 +167,7 @@ struct InfinityNormVisitor {
         std::decay_t<FloatType> max_element{0.0};
         for(std::size_t i = 0; i < buffer.size(); ++i) {
             auto elem = types::fabs(buffer[i]);
-            if constexpr(types::is_interval_v<std::decay_t<FloatType>>) {
-                if(elem.upper() > max_element.upper()) max_element = elem;
-            } else {
-                if(elem > max_element) max_element = elem;
-            }
+            if(!types::strictly_less(elem, max_element)) { max_element = elem; }
         }
         return wtf::fp::make_float(max_element);
     }

src/tensorwrapper/diis/diis.cpp

@@ -24,17 +24,8 @@ namespace {
 
 struct Kernel {
     template<typename FloatType>
-    auto operator()(const std::span<FloatType>& t) {
-        using clean_type = std::decay_t<FloatType>;
-        double rv;
-        if constexpr(tensorwrapper::types::is_uncertain_v<clean_type>) {
-            rv = t[0].mean();
-        } else if constexpr(tensorwrapper::types::is_interval_v<clean_type>) {
-            rv = t[0].median();
-        } else {
-            rv = t[0];
-        }
-        return rv;
+    double operator()(const std::span<FloatType>& t) {
+        return static_cast<double>(tensorwrapper::types::uq_center(t[0]));
     }
 };
 

src/tensorwrapper/generate/add_noise.cpp

@@ -33,33 +33,6 @@ std::vector<std::size_t> shape_extents(const auto& shape) {
     return rv;
 }
 
-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);
-    }
-}
-
-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"); }
@@ -71,9 +44,9 @@ Tensor add_noise_impl(const Tensor& matrix, double t, std::mt19937& gen) {
     if(t > 0.0) {
         std::normal_distribution<double> dist(0.0, t);
         for(auto& x : data) {
-            const auto center = noise_center<T>(x);
+            const auto center = tensorwrapper::types::uq_center(x);
             const auto delta  = std::clamp(dist(gen), -t, t);
-            x                 = apply_noise<T>(center, delta, t);
+            x = tensorwrapper::types::construct_uq_type<T>(center, delta);
         }
     }