Intervals.hpp

Go to the documentation of this file.

#ifndef ZONOOPT_INTERVAL_UTILITIES_HPP_
#define ZONOOPT_INTERVAL_UTILITIES_HPP_
 
#include <limits>
#include <string>
#include <vector>
#include <set>
#include <cmath>
#include <cassert>
#include <Eigen/Dense>
#include <Eigen/Sparse>
 
/*
Reference: 
"Applied Interval Analysis"
Luc Jaulin, Michel Kieffer, Olivier Didrit, Eric Walter
*/
 
namespace ZonoOpt {
 
    using namespace detail;
 
    // forward declarations
    struct Interval;
    struct IntervalView;
 
 
    template <typename Derived>
    struct IntervalBase
    {
        zono_float& y_min() { return static_cast<Derived*>(this)->get_y_min(); }
 
        zono_float& y_max() { return static_cast<Derived*>(this)->get_y_max(); }
 
        const zono_float& y_min() const { return static_cast<const Derived*>(this)->get_y_min(); }
 
        const zono_float& y_max() const { return static_cast<const Derived*>(this)->get_y_max(); }
 
        void set(const zono_float min, const zono_float max)
        {
            y_min() = min;
            y_max() = max;
        }
 
        void add_assign(const Derived& x1, const Derived& x2)
        {
            y_min() = x1.y_min() + x2.y_min();
            y_max() = x1.y_max() + x2.y_max();
        }
 
        void subtract_assign(const Derived& x1, const Derived& x2)
        {
            y_min() = x1.y_min() - x2.y_max();
            y_max() = x1.y_max() - x2.y_min();
        }
 
        void multiply_assign(const Derived& x1, const Derived& x2)
        {
            zono_float a = x1.y_min() * x2.y_min();
            zono_float b = x1.y_min() * x2.y_max();
            zono_float c = x1.y_max() * x2.y_min();
            zono_float d = x1.y_max() * x2.y_max();
 
            y_min() = std::min({a, b, c, d});
            y_max() = std::max({a, b, c, d});
        }
 
        void multiply_assign(const Derived& x1, zono_float alpha)
        {
            if (alpha >= 0)
            {
                y_min() = x1.y_min() * alpha;
                y_max() = x1.y_max() * alpha;
            }
            else
            {
                y_min() = x1.y_max() * alpha;
                y_max() = x1.y_min() * alpha;
            }
        }
 
        void inverse()
        {
            auto& min = y_min();
            auto& max = y_max();
 
            if (std::abs(min) < zono_eps && std::abs(max) < zono_eps)
            {
                min = std::numeric_limits<zono_float>::infinity();
                max = -std::numeric_limits<zono_float>::infinity();
            }
            else if (min > 0 || max < 0)
            {
                min = one / max;
                max = one / min;
            }
            else if (std::abs(min) < zono_eps && max > 0)
            {
                min = one / max;
                max = std::numeric_limits<zono_float>::infinity();
            }
            else if (min < 0 && std::abs(max) < zono_eps)
            {
                max = one / min;
                min = -std::numeric_limits<zono_float>::infinity();
            }
            else
            {
                min = -std::numeric_limits<zono_float>::infinity();
                max = std::numeric_limits<zono_float>::infinity();
            }
        }
 
        void divide_assign(const Derived& x1, const Derived& x2)
        {
            Derived inv = x2;
            inv.inverse();
            multiply_assign(x1, inv);
        }
 
        void intersect_assign(const Derived& x1, const Derived& x2)
        {
            y_min() = std::max(x1.y_min(), x2.y_min());
            y_max() = std::min(x1.y_max(), x2.y_max());
        }
 
        bool is_empty() const
        {
            return y_min() - y_max() > zono_eps;
        }
 
        bool contains(zono_float y) const
        {
            return y >= y_min() - zono_eps && y <= y_max() + zono_eps;
        }
 
        bool is_single_valued() const
        {
            return std::abs(y_max() - y_min()) < zono_eps;
        }
 
        zono_float width() const
        {
            return std::abs(y_max() - y_min());
        }
 
        void sin_assign(const Derived& x)
        {
            if (x.y_max() - x.y_min() >= two*pi)
            {
                y_max() = one;
                y_min() = -one;
            }
            else
            {
                // shift domain to [-pi, pi]
                zono_float u = x.y_max(); // init
                zono_float l = x.y_min(); // init
                while (u > pi)
                {
                    u -= two*pi;
                    l -= two*pi;
                }
                while (l < -pi)
                {
                    u += two*pi;
                    l += two*pi;
                }
 
                // get bounds
                y_max() = ((l < pi/two && pi/two < u) || (l < -3*pi/two && -3*pi/two < u)) ? one : std::max(std::sin(u), std::sin(l));
                y_min() = ((l < -pi/two && -pi/two < u) || (l < 3*pi/two && 3*pi/two < u)) ? -one : std::min(std::sin(u), std::sin(l));
            }
        }
 
        void cos_assign(const Derived& x)
        {
            Derived x_sin;
            x_sin.set(x.y_min() + pi/two, x.y_max() + pi/two);
            sin_assign(x_sin);
        }
 
        void tan_assign(const Derived& x)
        {
            if (x.y_max() - x.y_min() >= pi)
            {
                y_max() = std::numeric_limits<zono_float>::infinity();
                y_min() = -std::numeric_limits<zono_float>::infinity();
            }
            else
            {
                // shift domain to [-pi, pi]
                zono_float u = x.y_max(); // init
                zono_float l = x.y_min(); // init
                while (u > pi)
                {
                    u -= two*pi;
                    l -= two*pi;
                }
                while (l < -pi)
                {
                    u += two*pi;
                    l += two*pi;
                }
 
                // get bounds
                if ((l < -pi/two && -pi/two < u) || (l < pi/two && pi/two < u))
                {
                    y_max() = std::numeric_limits<zono_float>::infinity();
                    y_min() = -std::numeric_limits<zono_float>::infinity();
                }
                else
                {
                    y_max() = std::tan(u);
                    y_min() = std::tan(l);
                }
            }
        }
 
        void arcsin_assign(const Derived& x)
        {
            assert(x.y_min() >= -one && x.y_min() <= one);
            assert(x.y_max() >= -one && x.y_max() <= one);
            y_min() = std::asin(x.y_min());
            y_max() = std::asin(x.y_max());
        }
 
        void arccos_assign(const Derived& x)
        {
            assert(x.y_min() >= -one && x.y_min() <= one);
            assert(x.y_max() >= -one && x.y_max() <= one);
            y_min() = std::acos(x.y_max());
            y_max() = std::acos(x.y_min());
        }
 
        void arctan_assign(const Derived& x)
        {
            y_min() = std::atan(x.y_min());
            y_max() = std::atan(x.y_max());
        }
 
        void exp_assign(const Derived& x)
        {
            y_min() = std::exp(x.y_min());
            y_max() = std::exp(x.y_max());
        }
    };
 
    struct Interval : IntervalBase<Interval>
    {
        // members
        zono_float lb;
 
        zono_float ub;
 
        // constructor
        Interval() : lb(0), ub(0) {}
 
        Interval(zono_float y_min, zono_float y_max) : lb(y_min), ub(y_max) {}
 
        Interval* clone() const
        {
            return new Interval(*this);
        }
 
        // get methods
 
        zono_float& get_y_min() { return lb; }
 
        zono_float& get_y_max() { return ub; }
 
        const zono_float& get_y_min() const { return lb; }
 
        const zono_float& get_y_max() const { return ub; }
 
        // operators
 
        Interval operator+(const Interval& other) const
        {
            Interval result;
            result.add_assign(*this, other);
            return result;
        }
 
        Interval operator-(const Interval& other) const
        {
            Interval result;
            result.subtract_assign(*this, other);
            return result;
        }
 
        Interval operator*(const Interval& other) const
        {
            Interval result;
            result.multiply_assign(*this, other);
            return result;
        }
 
        Interval operator*(zono_float alpha) const
        {
            Interval result;
            result.multiply_assign(*this, alpha);
            return result;
        }
 
        Interval inv() const
        {
            Interval result = *this;
            result.inverse();
            return result;
        }
 
        Interval operator/(const Interval& other) const
        {
            Interval result;
            result.divide_assign(*this, other);
            return result;
        }
 
        Interval intersect(const Interval& other) const
        {
            Interval result;
            result.intersect_assign(*this, other);
            return result;
        }
 
        zono_float center() const
        {
            return (ub + lb) / two;
        }
 
        // as interval view
        IntervalView as_view();
 
        // print methods
        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;
        }
 
        Interval sin() const
        {
            Interval result;
            result.sin_assign(*this);
            return result;
        }
 
        Interval cos() const
        {
            Interval result;
            result.cos_assign(*this);
            return result;
        }
 
        Interval tan() const
        {
            Interval result;
            result.tan_assign(*this);
            return result;
        }
 
        Interval arcsin() const
        {
            Interval result;
            result.arcsin_assign(*this);
            return result;
        }
 
        Interval arccos() const
        {
            Interval result;
            result.arccos_assign(*this);
            return result;
        }
 
        Interval arctan() const
        {
            Interval result;
            result.arctan_assign(*this);
            return result;
        }
 
        Interval exp() const
        {
            Interval result;
            result.exp_assign(*this);
            return result;
        }
    };
 
    struct IntervalView : IntervalBase<IntervalView>
    {
        // members
 
        zono_float* lb_ptr = nullptr;
 
        zono_float* ub_ptr = nullptr;
 
        // constructor
 
        IntervalView(zono_float* y_min, zono_float* y_max) : lb_ptr(y_min), ub_ptr(y_max) {}
 
        // assignment
 
        template <typename Derived>
        IntervalView& operator=(const IntervalBase<Derived>& other)
        {
            y_min() = other.y_min();
            y_max() = other.y_max();
            return *this;
        }
 
        // to Interval
        Interval to_interval() const;
 
        // methods
 
        zono_float& get_y_min() { return *lb_ptr; }
 
        zono_float& get_y_max() { return *ub_ptr; }
 
        const zono_float& get_y_min() const { return *lb_ptr; }
 
        const zono_float& get_y_max() const { return *ub_ptr; }
    };
 
    // implementations
    inline Interval IntervalView::to_interval() const
    {
        return Interval(*lb_ptr, *ub_ptr);
    }
 
    inline IntervalView Interval::as_view()
    {
        return IntervalView(&lb, &ub);
    }
 
    class Box
    {
    public:
 
        // constructors
 
        Box() = default;
 
        explicit Box(const size_t size);
 
        explicit Box(const std::vector<Interval>& 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
 
        IntervalView operator[](size_t i);
 
        Interval operator[](size_t i) const;
        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;
 
        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
 
        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);
 
        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);
    };
 
} // namespace ZonoOpt
 
 
#endif