thresholded_affine.hpp

Go to the documentation of this file.

#pragma once
#include <cmath>
#include <sigma/affine/affine.hpp>

 
namespace sigma {

template<typename ValueType>
class ThresholdedAffine {
public:
    using size_type = std::size_t;

    using value_t = ValueType;

    using error_term_t = typename Affine<ValueType>::error_term_t;

    using error_terms_t = typename Affine<ValueType>::error_terms_t;

    using interval_t = typename Affine<ValueType>::interval_t;

    using affine_t = Affine<ValueType>;

    struct Threshold {
        value_t value;
        explicit Threshold(value_t v) : value(v) {}
    };
 
    // --- Constructors and Assignment ----------------------------------------

    ThresholdedAffine() : m_threshold_(default_threshold().value) {}

    explicit ThresholdedAffine(value_t center,
                               Threshold t = default_threshold()) :
      ThresholdedAffine(interval_t(center), t) {}

    ThresholdedAffine(value_t lo, value_t hi,
                      Threshold t = default_threshold()) :
      ThresholdedAffine(interval_t(lo, hi), t) {
        apply_threshold_();
    }

    explicit ThresholdedAffine(const interval_t& interval,
                               Threshold t = default_threshold()) :
      m_affine_(interval), m_threshold_(t.value) {
        apply_threshold_();
    }

    ThresholdedAffine(value_t center, error_terms_t radii,
                      Threshold t = default_threshold()) :
      m_affine_(center, std::move(radii)), m_threshold_(t.value) {
        apply_threshold_();
    }

    ThresholdedAffine(affine_t a, value_t threshold) :
      m_affine_(std::move(a)), m_threshold_(threshold) {
        apply_threshold_();
    }

    ThresholdedAffine(const ThresholdedAffine& other) = default;

    ThresholdedAffine(ThresholdedAffine&& other) noexcept = default;

    ThresholdedAffine& operator=(const ThresholdedAffine& other) = default;

    ThresholdedAffine& operator=(ThresholdedAffine&& other) noexcept = default;
 
    // -- State Accessors -----------------------------------------------------

    static constexpr Threshold default_threshold() { return Threshold(0.001); }

    interval_t range() const {
        if(m_affine_.empty()) return interval_t{};
        auto r = m_affine_.radius();
        return interval_t(m_affine_.center() - r, m_affine_.center() + r);
    }

    value_t center() const { return m_affine_.center(); }

    const error_terms_t& error_terms() const { return m_affine_.error_terms(); }

    value_t radius() const { return m_affine_.radius(); }

    void set_center(value_t c) { m_affine_.set_center(c); }

    void add_error_term(error_term_t error_term, value_t r) {
        m_affine_.add_error_term(error_term, r);
        apply_threshold_();
    }

    bool contains(value_t v) const { return range().contains(v); }

    bool contains(const interval_t& i) const { return range().contains(i); }

    bool contains(const ThresholdedAffine& other) const {
        return range().contains(other.range());
    }

    bool empty() const noexcept { return m_affine_.empty(); }

    std::string print_affine_form() const {
        return m_affine_.print_affine_form();
    }

    std::string print_interval_form() const {
        return range().print_interval_form();
    }

    const affine_t& affine() const { return m_affine_; }

    value_t threshold() const { return m_threshold_; }
 
    // -- Arithmetic Operators ------------------------------------------------

    ThresholdedAffine operator-() const {
        auto result      = *this;
        result.m_affine_ = -m_affine_;
        return result;
    }

    ThresholdedAffine& operator+=(value_t value) {
        m_affine_ += value;
        return *this;
    }

    ThresholdedAffine& operator+=(const ThresholdedAffine& other) {
        m_affine_ += other.m_affine_;
        apply_threshold_();
        return *this;
    }

    ThresholdedAffine operator+(value_t value) const {
        return ThresholdedAffine(*this) += value;
    }

    ThresholdedAffine operator+(const ThresholdedAffine& other) const {
        return ThresholdedAffine(*this) += other;
    }

    ThresholdedAffine& operator-=(value_t value) {
        m_affine_ -= value;
        return *this;
    }

    ThresholdedAffine& operator-=(const ThresholdedAffine& other) {
        m_affine_ -= other.m_affine_;
        apply_threshold_();
        return *this;
    }

    ThresholdedAffine operator-(value_t value) const {
        return ThresholdedAffine(*this) -= value;
    }

    ThresholdedAffine operator-(const ThresholdedAffine& other) const {
        return ThresholdedAffine(*this) -= other;
    }

    ThresholdedAffine& operator*=(value_t value) {
        m_affine_ *= value;
        apply_threshold_();
        return *this;
    }

    ThresholdedAffine& operator*=(const ThresholdedAffine& other) {
        m_affine_ *= other.m_affine_;
        apply_threshold_();
        return *this;
    }

    ThresholdedAffine operator*(value_t value) const {
        return ThresholdedAffine(*this) *= value;
    }

    ThresholdedAffine operator*(const ThresholdedAffine& other) const {
        return ThresholdedAffine(*this) *= other;
    }

    ThresholdedAffine& operator/=(value_t value) {
        m_affine_ /= value;
        apply_threshold_();
        return *this;
    }

    ThresholdedAffine& operator/=(const ThresholdedAffine& other) {
        // Minimum absolute value of other over its interval.
        value_t other_total = other.m_affine_.radius();
        value_t b_min       = std::fabs(other.m_affine_.center()) - other_total;
        if(b_min <= value_t(0)) {
            throw std::domain_error(
              "ThresholdedAffine division by interval containing zero");
        }
 
        m_affine_ /= other.m_affine_;
 
        apply_threshold_();
        return *this;
    }

    ThresholdedAffine operator/(value_t value) const {
        return ThresholdedAffine(*this) /= value;
    }

    ThresholdedAffine operator/(const ThresholdedAffine& other) const {
        return ThresholdedAffine(*this) /= other;
    }
 
    // -- Comparison Operators ------------------------------------------------

    bool operator==(const ThresholdedAffine& other) const {
        return m_affine_ == other.m_affine_ &&
               m_threshold_ == other.m_threshold_;
    }

    bool operator!=(const ThresholdedAffine& other) const {
        return !(*this == other);
    }
 
private:
    void apply_threshold_() {
        if(m_affine_.empty()) return;
        auto total_r = m_affine_.radius();
 
        auto terms = m_affine_.error_terms();
        if(terms.empty()) return;
 
        error_terms_t new_terms;
        for(auto&& [sym, coeff] : terms) {
            if(std::fabs(coeff) < std::numeric_limits<value_t>::epsilon()) {
                continue;
            } else if(total_r != value_t(0) &&
                      std::fabs(coeff) / total_r < m_threshold_) {
                continue;
            } else {
                new_terms[sym] = coeff;
            }
        }
        m_affine_ = affine_t(m_affine_.center(), std::move(new_terms));
    }

    affine_t m_affine_;

    value_t m_threshold_;
};
 
// -- Non-member functions ----------------------------------------------------

template<typename ValueType>
std::ostream& operator<<(std::ostream& os,
                         const ThresholdedAffine<ValueType>& a) {
    os << a.range();
    return os;
}

template<typename ValueType>
ThresholdedAffine<ValueType> operator*(ValueType value,
                                       const ThresholdedAffine<ValueType>& a) {
    return a * value;
}
 
// -- Operations --------------------------------------------------------------

template<typename T>
ThresholdedAffine<T> abs(const ThresholdedAffine<T>& a) {
    return ThresholdedAffine<T>(abs(a.affine()), a.threshold());
}

template<typename T>
ThresholdedAffine<T> fabs(const ThresholdedAffine<T>& a) {
    return ThresholdedAffine<T>(fabs(a.affine()), a.threshold());
}

template<typename T>
ThresholdedAffine<T> sqrt(const ThresholdedAffine<T>& a) {
    return ThresholdedAffine<T>(sqrt(a.affine()), a.threshold());
}

template<typename T>
ThresholdedAffine<T> exp(const ThresholdedAffine<T>& a) {
    return ThresholdedAffine<T>(exp(a.affine()), a.threshold());
}

template<typename T>
ThresholdedAffine<T> log(const ThresholdedAffine<T>& a) {
    return ThresholdedAffine<T>(log(a.affine()), a.threshold());
}

template<typename T, typename U>
ThresholdedAffine<T> pow(const ThresholdedAffine<T>& a, const U& exp) {
    return ThresholdedAffine<T>(pow(a.affine(), exp), a.threshold());
}

using TAFloat = ThresholdedAffine<float>;

using TADouble = ThresholdedAffine<double>;
 
} // namespace sigma