affine.hpp

Go to the documentation of this file.

#pragma once
#include <atomic>
#include <optional>
#include <sigma/interval/interval.hpp>
#include <sstream>
#include <unordered_map>

 
namespace sigma {

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

    using value_t = ValueType;

    using error_term_t = size_type;

    using error_terms_t = std::unordered_map<error_term_t, value_t>;

    using interval_t = Interval<value_t>;
 
    // --- Constructors and Assignment ----------------------------------------

    Affine() = default;

    Affine(value_t center) : Affine(interval_t(center, center)) {}

    Affine(value_t lo, value_t hi) : Affine(interval_t(lo, hi)) {}

    explicit Affine(const interval_t& interval) {
        if(interval.empty()) { return; }
        m_center_ = interval.median();
        if(interval.radius() > 0) {
            m_error_terms_[make_error_term()] = interval.radius();
        }
    }

    Affine(value_t center, error_terms_t radii) :
      m_center_(center), m_error_terms_(std::move(radii)) {}

    Affine(const Affine& other) = default;

    Affine(Affine&& other) noexcept = default;

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

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

    interval_t range() const;

    value_t center() const {
        assert_not_empty_();
        return *m_center_;
    }

    const error_terms_t& error_terms() const { return m_error_terms_; }

    value_t radius() const;

    void set_center(value_t center) { m_center_ = center; }

    void add_error_term(error_term_t error_term, value_t radius) {
        if(empty()) set_center(value_t{0});
        m_error_terms_[error_term] = radius;
    }

    bool contains(value_t value) const;

    bool contains(const interval_t& interval) const;

    bool contains(const Affine& affine) const {
        return contains(affine.range());
    }

    bool empty() const noexcept { return !m_center_; }

    std::string print_affine_form() const;

    std::string print_interval_form() const {
        return range().print_interval_form();
    }
 
    // -- Arithmetic Operators ------------------------------------------------

    Affine operator-() const;

    Affine& operator+=(value_t value) {
        if(empty()) { return *this = Affine(value); }
        (*m_center_) += value;
        return *this;
    }

    Affine& operator+=(const Affine& other);

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

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

    Affine& operator-=(value_t value) {
        if(empty()) { return *this = Affine(-value); }
        (*m_center_) -= value;
        return *this;
    }

    Affine& operator-=(const Affine& other);

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

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

    Affine& operator*=(value_t value) {
        if(empty()) { return *this; }
        (*m_center_) *= value;
        for(auto&& [error_symbol, error_term_i] : m_error_terms_) {
            error_term_i *= value;
        }
        return *this;
    }

    Affine& operator*=(const Affine& other);

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

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

    Affine& operator/=(value_t value) {
        if(value == 0) { throw std::domain_error("Division by zero"); }
        return *this *= value_t(1.0 / value);
    }

    Affine& operator/=(const Affine& other);

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

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

    Affine apply_affine_transform(value_t alpha, value_t zeta,
                                  value_t delta) const {
        value_t new_center = empty() ? zeta : alpha * center() + zeta;
        error_terms_t new_error_terms;
        for(auto&& [error_symbol, error_term_i] : m_error_terms_) {
            new_error_terms[error_symbol] = alpha * error_term_i;
        }
        new_error_terms[make_error_term()] = delta;
        return Affine(new_center, std::move(new_error_terms));
    }

    Affine multiplicative_inverse() const;
 
    // -- Comparison Operators ------------------------------------------------

    bool operator==(const Affine& other) const {
        if(empty() != other.empty()) { return false; }
        if(empty()) { return true; }
        if(m_center_ != other.m_center_) { return false; }
        return m_error_terms_ == other.m_error_terms_;
    }

    bool operator!=(const Affine& other) const { return !(*this == other); }

    static error_term_t make_error_term() {
        static std::atomic<size_type> s_next_id{0};
        return s_next_id.fetch_add(1);
    }
 
private:
    void assert_not_empty_() const {
        if(empty()) { throw std::domain_error("Affine form is empty"); }
    }

    std::optional<value_t> m_center_;

    error_terms_t m_error_terms_;
};
 
// -- Non-member functions
// ---------------------------------------------------

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

template<typename ValueType>
Affine<ValueType> operator*(ValueType value, const Affine<ValueType>& a) {
    return a * value;
}
 
// -- Out-of-line definitions
// ------------------------------------------------
 
template<typename ValueType>
auto Affine<ValueType>::range() const -> interval_t {
    if(empty()) { return interval_t(); }
    auto r = radius();
    return interval_t(center() - r, center() + r);
}
 
template<typename ValueType>
auto Affine<ValueType>::radius() const -> value_t {
    assert_not_empty_();
    value_t r = 0;
    for(auto&& [error_symbol, error_term_i] : m_error_terms_) {
        r += std::fabs(error_term_i);
    }
    return r;
}
 
template<typename ValueType>
auto Affine<ValueType>::contains(value_t value) const -> bool {
    if(empty()) { return false; }
    return range().contains(value);
}
 
template<typename ValueType>
auto Affine<ValueType>::contains(const interval_t& interval) const -> bool {
    if(interval.empty()) { return true; }
    if(empty()) { return false; }
    return range().contains(interval);
}
 
template<typename ValueType>
std::string Affine<ValueType>::print_affine_form() const {
    if(empty()) { return "∅"; }
    std::stringstream ss;
    ss << center();
    for(auto&& [error_symbol, error_term_i] : m_error_terms_) {
        ss << " +/- " << error_term_i;
    }
    return ss.str();
}
 
template<typename ValueType>
auto Affine<ValueType>::operator-() const -> Affine {
    if(empty()) { return *this; }
    value_t new_center = -center();
    error_terms_t new_error_terms;
    for(auto&& [error_symbol, error_term_i] : m_error_terms_) {
        new_error_terms[error_symbol] = -error_term_i;
    }
    return Affine(new_center, new_error_terms);
}
 
template<typename ValueType>
auto Affine<ValueType>::operator+=(const Affine& other) -> Affine& {
    if(empty()) { return *this = other; }
    if(other.empty()) { return *this; }
    value_t new_center            = center() + other.center();
    error_terms_t new_error_terms = m_error_terms_;
    for(auto&& [error_symbol, error_term_i] : other.m_error_terms_) {
        new_error_terms[error_symbol] += error_term_i;
    }
    return *this = Affine(new_center, new_error_terms);
}
 
template<typename ValueType>
auto Affine<ValueType>::operator-=(const Affine& other) -> Affine& {
    return *this += -other;
}
 
template<typename ValueType>
auto Affine<ValueType>::operator*=(const Affine& other) -> Affine& {
    if(empty() || other.empty()) { return *this = Affine(); }
    value_t new_center = center() * other.center();
    error_terms_t new_error_terms;
    value_t new_radius = 0;
    for(auto&& [error_symbol, error_term_i] : m_error_terms_) {
        new_error_terms[error_symbol] = error_term_i * other.center();
        new_radius += std::fabs(new_error_terms[error_symbol]);
    }
    for(auto&& [error_symbol, error_term_j] : other.m_error_terms_) {
        new_error_terms[error_symbol] += error_term_j * center();
        new_radius += std::fabs(new_error_terms[error_symbol]);
    }
    auto correction                    = radius() * other.radius();
    new_error_terms[make_error_term()] = correction;
    return *this                       = Affine(new_center, new_error_terms);
}
 
template<typename ValueType>
auto Affine<ValueType>::operator/=(const Affine& other) -> Affine& {
    // Multiply by this by 1 / other
    return *this *= other.multiplicative_inverse();
}
 
template<typename ValueType>
auto Affine<ValueType>::multiplicative_inverse() const -> Affine {
    assert_not_empty_();
    if(contains(0)) { throw std::domain_error("Division by zero"); }
    // Compute the affine transformation which transforms other to 1 / other
    auto other_range = range();
    auto abs_inf     = std::fabs(other_range.lower());
    auto abs_sup     = std::fabs(other_range.upper());
    auto a           = std::min(abs_inf, abs_sup);
    auto b           = std::max(abs_inf, abs_sup);
    auto alpha       = value_t(-1.0) / (b * b);
    auto lo          = value_t(1.0) / a - alpha * a;
    auto hi          = value_t(2.0) / b;
    interval_t interval(std::min(lo, hi), std::max(lo, hi));
    auto zeta  = std::fabs(interval.median());
    auto delta = interval.radius();
    return apply_affine_transform(alpha, zeta, delta);
}

using AFloat = Affine<float>;

using ADouble = Affine<double>;
 
} // namespace sigma
 
#include <sigma/affine/eigen_compat.hpp>
#include <sigma/affine/operations/operations.hpp>