Intervals.hpp

Go to the documentation of this file.

#ifndef ZONOOPT_INTERVAL_UTILITIES_HPP_
#define ZONOOPT_INTERVAL_UTILITIES_HPP_
 
#include <string>
#include <vector>
#include <set>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <cmath>
#include <boost/numeric/interval.hpp>
 
namespace ZonoOpt
{
    namespace detail
    {
        // need this to get inverse hyperbolics to compile with MSVC
        template<class T>
        struct rounded_transc_fixed : boost::numeric::interval_lib::rounded_transc_std<T> {
            static T asinh_down(const T& x) { using std::asinh; return asinh(x); }
            static T asinh_up  (const T& x) { using std::asinh; return asinh(x); }
            static T acosh_down(const T& x) { using std::acosh; return acosh(x); }
            static T acosh_up  (const T& x) { using std::acosh; return acosh(x); }
            static T atanh_down(const T& x) { using std::atanh; return atanh(x); }
            static T atanh_up  (const T& x) { using std::atanh; return atanh(x); }
        };
    }
 
    using namespace detail;
 
 
    class Interval
    {
    public:
        // constructor
 
        Interval() : _val(zero, zero) {}
 
        Interval(const zono_float y_min, const zono_float y_max) : _val(y_min, y_max) {}
 
        Interval* clone() const
        {
            return new Interval(*this);
        }
 
        // get methods
 
        zono_float lb() const { return _val.lower(); }
 
        zono_float ub() const { return _val.upper(); }
 
        // operators
 
        Interval operator+(const Interval& other) const
        {
            return Interval(this->_val + other._val);
        }
 
        Interval operator+(const zono_float alpha) const
        {
            return Interval(this->_val + alpha);
        }
 
        Interval operator-(const Interval& other) const
        {
            return Interval(this->_val - other._val);
        }
 
        Interval operator-(const zono_float alpha) const
        {
            return Interval(this->_val - alpha);
        }
 
        Interval operator*(const Interval& other) const
        {
            return Interval(this->_val * other._val);
        }
 
        Interval operator*(const zono_float alpha) const
        {
            return Interval(this->_val * alpha);
        }
 
        Interval operator/(const Interval& other) const
        {
            return Interval(this->_val / other._val);
        }
 
        Interval operator/(const zono_float alpha) const
        {
            return Interval(this->_val / alpha);
        }
 
        Interval inv() const
        {
            return Interval(boost::numeric::interval_lib::multiplicative_inverse(_val));
        }
 
        Interval intersect(const Interval& other) const
        {
            return Interval(boost::numeric::intersect(_val, other._val));
        }
 
        bool contains(const zono_float x) const
        {
            return x >= this->lb() - zono_eps && x <= this->ub() + zono_eps;
        }
 
        bool is_single_valued() const
        {
            return this->width() < zono_eps;
        }
 
        bool is_empty() const
        {
            return boost::numeric::empty(_val);
        }
 
        zono_float center() const
        {
            return boost::numeric::median(_val);
        }
 
        zono_float width() const
        {
            return _val.upper() - _val.lower();
        }
 
        Interval radius() const
        {
            const zono_float r = this->width() / two;
            return Interval(-r, r);
        }
 
 
        Interval abs() const
        {
            return Interval(boost::numeric::abs(_val));
        }
 
        Interval sqrt() const
        {
            return Interval(boost::numeric::sqrt(_val));
        }
 
        Interval pow(const int n) const
        {
            return Interval(boost::numeric::pow(_val, n));
        }
 
        Interval nth_root(const int n) const
        {
            return Interval(boost::numeric::nth_root(_val, n));
        }
 
        Interval exp() const
        {
            return Interval(boost::numeric::exp(_val));
        }
 
        Interval log() const
        {
            return Interval(boost::numeric::log(_val));
        }
 
        Interval sin() const
        {
            return Interval(boost::numeric::sin(_val));
        }
 
        Interval cos() const
        {
            return Interval(boost::numeric::cos(_val));
        }
 
        Interval tan() const
        {
            return Interval(boost::numeric::tan(_val));
        }
 
        Interval arcsin() const
        {
            return Interval(boost::numeric::asin(_val));
        }
 
        Interval arccos() const
        {
            return Interval(boost::numeric::acos(_val));
        }
 
        Interval arctan() const
        {
            return Interval(boost::numeric::atan(_val));
        }
 
        Interval sinh() const
        {
            return Interval(boost::numeric::sinh(_val));
        }
 
        Interval cosh() const
        {
            return Interval(boost::numeric::cosh(_val));
        }
 
        Interval tanh() const
        {
            return Interval(boost::numeric::tanh(_val));
        }
 
        Interval arcsinh() const
        {
            return Interval(boost::numeric::asinh(_val));
        }
 
        Interval arccosh() const
        {
            return Interval(boost::numeric::acosh(_val));
        }
 
        Interval arctanh() const
        {
            return Interval(boost::numeric::atanh(_val));
        }
 
        std::string print() const
        {
            return "Interval: [" + std::to_string(lb()) + ", " + std::to_string(ub()) + "]";
        }
 
        friend std::ostream& operator<<(std::ostream& os, const Interval& interval)
        {
            os << interval.print();
            return os;
        }
 
 
    private:
 
        typedef boost::numeric::interval_lib::policies<
            boost::numeric::interval_lib::save_state<rounded_transc_fixed<zono_float>>,
            boost::numeric::interval_lib::checking_base<zono_float>
        > interval_policy;
 
        boost::numeric::interval<zono_float, interval_policy> _val;
 
        explicit Interval(const boost::numeric::interval<zono_float, interval_policy>& val) : _val(val) {}
    };
 
 
    class Box
    {
    public:
        // constructors
 
        Box() = default;
 
        explicit Box(const size_t size);
 
        explicit Box(const std::vector<Interval>& vals);
 
        explicit Box(const Eigen::Vector<Interval, -1>& vals);
 
        Box(const Eigen::Vector<zono_float, -1>& x_lb, const Eigen::Vector<zono_float, -1>& x_ub);
 
        // virtual destructor
        virtual ~Box() = default;
 
        // copy assignment
        Box& operator=(const Box& other);
 
        // copy constructor
        Box(const Box& other);
 
        // element-wise assignment, access
 
        Interval get_element(int i) const;
 
 
        void set_element(int i, const Interval& val);
 
        size_t size() const;
 
        // project onto box
 
        virtual void project(Eigen::Ref<Eigen::Vector<zono_float, -1>> x) const;
 
        virtual Box* clone() const;
 
        // access bounds
        const Eigen::Vector<zono_float, -1>& lower() const { return x_lb; }
 
        const Eigen::Vector<zono_float, -1>& upper() const { return x_ub; }
 
        // width of box
        zono_float width() const;
 
        Box radius() const;
 
        Eigen::Vector<zono_float, -1> center() const;
 
        // operator overloading
 
        Box operator+(const Box& other) const;
 
        Box operator-(const Box& other) const;
 
        Box operator*(const Box& other) const;
 
        Box operator*(zono_float alpha) const;
 
        Box operator/(const Box& other) const;
 
        // interval contractors
 
        bool contract(const Eigen::SparseMatrix<zono_float, Eigen::RowMajor>& A, const Eigen::Vector<zono_float, -1>& b,
                      int iter);
 
        bool contract_subset(const Eigen::SparseMatrix<zono_float, Eigen::RowMajor>& A_rm,
                             const Eigen::Vector<zono_float, -1>& b, int iter,
                             const Eigen::SparseMatrix<zono_float>& A, const std::set<int>& inds,
                             int tree_search_depth);
 
        Box linear_map(const Eigen::Matrix<zono_float, -1, -1>& A) const;
 
        Box linear_map(const Eigen::SparseMatrix<zono_float, Eigen::RowMajor>& A) const;
 
        Interval dot(const Eigen::Vector<zono_float, -1>& x) const;
 
        void permute(const Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic>& P);
 
        std::string print() const;
 
        friend std::ostream& operator<<(std::ostream& os, const Box& box);
 
    protected:
        // members
        Eigen::Vector<zono_float, -1> x_lb;
 
        Eigen::Vector<zono_float, -1> x_ub;
 
    private:
        // back end for contraction operator
 
        virtual bool contract_helper(const Eigen::SparseMatrix<zono_float, Eigen::RowMajor>& A,
                                     const Eigen::Vector<zono_float, -1>& b, const int iter,
                                     const std::set<int>& constraints);
 
        bool contract_forward(const Eigen::SparseMatrix<zono_float, Eigen::RowMajor>& A,
                              const Eigen::Vector<zono_float, -1>& b, const std::set<int>& constraints);
 
        void contract_backward(const Eigen::SparseMatrix<zono_float, Eigen::RowMajor>& A,
                               const Eigen::Vector<zono_float, -1>& b, const std::set<int>& constraints);
 
        static void get_vars_cons(const Eigen::SparseMatrix<zono_float>& A,
                                  const Eigen::SparseMatrix<zono_float, Eigen::RowMajor>& A_rm,
                                  std::set<int>& constraints, std::set<int>& vars, const std::set<int>& new_vars,
                                  int depth, int max_depth);
 
    private:
        friend class MI_Box;
    };
 
    // mixed-integer box (some bounds are fixed)
    class MI_Box final : public Box
    {
    public:
        // constructors
 
        MI_Box() = default;
 
        MI_Box(const Eigen::Vector<zono_float, -1>& x_lb, const Eigen::Vector<zono_float, -1>& x_ub,
               const std::pair<int, int>& idx_b, bool zero_one_form);
 
        // clone operation
        Box* clone() const override;
 
        // project
        void project(Eigen::Ref<Eigen::Vector<zono_float, -1>> x) const override;
 
        // get binary indices
        const std::pair<int, int>& binary_indices() const { return idx_b; }
 
    protected:
        bool contract_helper(const Eigen::SparseMatrix<zono_float, Eigen::RowMajor>& A,
                             const Eigen::Vector<zono_float, -1>& b, const int iter,
                             const std::set<int>& constraints) override;
 
    private:
        std::pair<int, int> idx_b;
        zono_float bin_low = zero, bin_high = one;
    };
 
 
    class IntervalMatrix
    {
    public:
        IntervalMatrix() : rows_(0), cols_(0)
        {
        }
 
        IntervalMatrix(size_t rows, size_t cols, const std::vector<Eigen::Triplet<Interval>>& triplets);
 
        IntervalMatrix(const Eigen::SparseMatrix<zono_float>& mat_lb,
                       const Eigen::SparseMatrix<zono_float>& mat_ub);
 
        IntervalMatrix(const Eigen::Matrix<zono_float, -1, -1>& mat_lb,
                       const Eigen::Matrix<zono_float, -1, -1>& mat_ub);
 
        IntervalMatrix(const Eigen::Matrix<Interval, -1, -1>& mat);
 
        Eigen::SparseMatrix<zono_float> center() const;
 
        Eigen::SparseMatrix<zono_float> diam() const;
 
        IntervalMatrix radius() const;
 
        zono_float width() const;
 
        Box operator*(const Eigen::Vector<zono_float, -1>& v) const;
 
        Box operator*(const Box& box) const;
 
        IntervalMatrix operator*(const Eigen::SparseMatrix<zono_float, Eigen::RowMajor>& A) const;
 
        IntervalMatrix operator*(const IntervalMatrix& other) const;
 
        IntervalMatrix operator+(const IntervalMatrix& other) const;
 
        IntervalMatrix operator-(const IntervalMatrix& other) const;
 
        size_t rows() const { return rows_; }
 
        size_t cols() const { return cols_; }
 
        std::string print() const;
 
        friend std::ostream& operator<<(std::ostream& os, const IntervalMatrix& interval_matrix);
 
    private:
        size_t rows_, cols_;
        std::vector<std::vector<std::pair<size_t, Interval>>> mat_; // rows->cols->vals
    };
} // namespace ZonoOpt
 
#endif