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
*/
 
#ifndef M_PI
    #define M_PI 3.14159265358979323846
#endif
 
namespace ZonoOpt {
 
    // 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 = 1.0 / max;
                max = 1.0 / min;
            }
            else if (std::abs(min) < zono_eps && max > 0)
            {
                min = 1.0 / max;
                max = std::numeric_limits<zono_float>::infinity();
            }
            else if (min < 0 && std::abs(max) < zono_eps)
            {
                max = 1.0 / 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() >= 2*M_PI)
            {
                y_max() = 1;
                y_min() = -1;
            }
            else
            {
                // shift domain to [-pi, pi]
                zono_float u = x.y_max(); // init
                zono_float l = x.y_min(); // init
                while (u > M_PI)
                {
                    u -= 2*M_PI;
                    l -= 2*M_PI;
                }
                while (l < -M_PI)
                {
                    u += 2*M_PI;
                    l += 2*M_PI;
                }
 
                // get bounds
                y_max() = ((l < M_PI/2 && M_PI/2 < u) || (l < -3*M_PI/2 && -3*M_PI/2 < u)) ? 1 : std::max(std::sin(u), std::sin(l));
                y_min() = ((l < -M_PI/2 && -M_PI/2 < u) || (l < 3*M_PI/2 && 3*M_PI/2 < u)) ? -1 : std::min(std::sin(u), std::sin(l));
            }
        }
 
        void cos_assign(const Derived& x)
        {
            Derived x_sin;
            x_sin.set(x.y_min() + M_PI/2, x.y_max() + M_PI/2);
            sin_assign(x_sin);
        }
 
        void tan_assign(const Derived& x)
        {
            if (x.y_max() - x.y_min() >= M_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 > M_PI)
                {
                    u -= 2*M_PI;
                    l -= 2*M_PI;
                }
                while (l < -M_PI)
                {
                    u += 2*M_PI;
                    l += 2*M_PI;
                }
 
                // get bounds
                if ((l < -M_PI/2 && -M_PI/2 < u) || (l < M_PI/2 && M_PI/2 < 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() >= -1 && x.y_min() <= 1);
            assert(x.y_max() >= -1 && x.y_max() <= 1);
            y_min() = std::asin(x.y_min());
            y_max() = std::asin(x.y_max());
        }
 
        void arccos_assign(const Derived& x)
        {
            assert(x.y_min() >= -1 && x.y_min() <= 1);
            assert(x.y_max() >= -1 && x.y_max() <= 1);
            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) / 2;
        }
 
        // 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)
        {
            x_lb.resize(static_cast<Eigen::Index>(size));
            x_ub.resize(static_cast<Eigen::Index>(size));
        }
 
        explicit Box(const std::vector<Interval>& vals)
        {
            x_lb.resize(static_cast<Eigen::Index>(vals.size()));
            x_ub.resize(static_cast<Eigen::Index>(vals.size()));
 
            for (Eigen::Index i=0; i<vals.size(); i++)
            {
                this->x_lb(i) = vals[i].y_min();
                this->x_ub(i) = vals[i].y_max();
            }
        }
 
        Box(const Eigen::Vector<zono_float, -1>& x_lb, const Eigen::Vector<zono_float, -1>& x_ub)
        {
            if (x_lb.size() != x_ub.size())
                throw std::invalid_argument("x_l and x_u must have the same size");
            this->x_lb = x_lb;
            this->x_ub = x_ub;
        }
 
        // virtual destructor
        virtual ~Box() = default;
 
        // copy assignment
        Box& operator=(const Box& other)
        {
            if (this != &other)
            {
                this->x_lb = other.x_lb;
                this->x_ub = other.x_ub;
            }
            return *this;
        }
 
        // copy constructor
        Box(const Box& other)
        {
            this->x_lb = other.x_lb;
            this->x_ub = other.x_ub;
        }
 
        // element-wise assignment, access
 
        IntervalView operator[](const size_t i)
        {
            if (i >= x_lb.size())
                throw std::out_of_range("Index out of range");
            return IntervalView(&x_lb(static_cast<Eigen::Index>(i)), &x_ub(static_cast<Eigen::Index>(i)));
        }
 
        Interval operator[](const size_t i) const
        {
            if (i >= x_lb.size())
                throw std::out_of_range("Index out of range");
            return Interval(x_lb(static_cast<Eigen::Index>(i)), x_ub(static_cast<Eigen::Index>(i)));
        }
 
        size_t size() const
        {
            return x_lb.size();
        }
 
        // project onto box
 
        virtual void project(Eigen::Ref<Eigen::Vector<zono_float, -1>> x) const
        {
            if (x.size() != x_lb.size())
                throw std::invalid_argument("x must have the same size as the Box");
            x = x.cwiseMax(x_lb).cwiseMin(x_ub);
        }
 
        virtual Box* clone() const
        {
            return new Box(*this);
        }
 
        // 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
        {
            zono_float w = 0;
            for (Eigen::Index i=0; i<x_lb.size(); i++)
            {
                w += x_ub(i) - x_lb(i);
            }
            return w;
        }
 
        Eigen::Vector<zono_float, -1> center() const
        {
            Eigen::Vector<zono_float, -1> c (this->size());
            for (Eigen::Index i=0; i<this->size(); i++)
            {
                c(i) = (*this)[i].center();
            }
            return c;
        }
 
        // operator overloading
 
        Box operator+(const Box& other) const
        {
            if (this->size() != other.size())
                throw std::invalid_argument("Box addition: inconsistent dimensions");
            Box out =  *this;
            for (size_t i=0; i<this->size(); ++i)
            {
                out[i] = (*this)[i] + other[i];
            }
            return out;
        }
 
        Box operator-(const Box& other) const
        {
            if (this->size() != other.size())
                throw std::invalid_argument("Box subtraction: inconsistent dimensions");
            Box out =  *this;
            for (size_t i=0; i<this->size(); ++i)
            {
                out[i] = (*this)[i] - other[i];
            }
            return out;
        }
 
        Box operator*(const Box& other) const
        {
            if (this->size() != other.size())
                throw std::invalid_argument("Box multiplication: inconsistent dimensions");
            Box out =  *this;
            for (size_t i=0; i<this->size(); ++i)
            {
                out[i] = (*this)[i] * other[i];
            }
            return out;
        }
 
        Box operator*(zono_float alpha) const
        {
            Box out = *this;
            for (size_t i=0; i<this->size(); ++i)
            {
                out[i] = (*this)[i]*alpha;
            }
            return out;
        }
 
        Box operator/(const Box& other) const
        {
            if (this->size() != other.size())
                throw std::invalid_argument("Box division: inconsistent dimensions");
            Box out = *this;
            for (size_t i=0; i<this->size(); ++i)
            {
                out[i] = (*this)[i] / other[i];
            }
            return out;
        }
 
        // interval contractors
 
        bool contract(const Eigen::SparseMatrix<zono_float, Eigen::RowMajor>& A, const Eigen::Vector<zono_float, -1>& b, const int iter)
        {
            if (iter <= 0)
                throw std::invalid_argument("iter must be positive");
 
            // contract over all constraints
            std::set<int> constraints;
            for (int i=0; i<A.rows(); i++)
            {
                constraints.insert(i);
            }
 
            // run contractor
            return contract_helper(A, b, iter, constraints);
        }
 
        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)
        {
            if (iter <= 0)
                throw std::invalid_argument("iter must be positive");
            if (A_rm.rows() != A.rows() || A_rm.cols() != A.cols())
                throw std::invalid_argument("A_rm must equal A");
 
            // get affected constraints
            std::set<int> all_constraints, all_vars;
 
            get_vars_cons(A, A_rm, all_constraints, all_vars, inds, 0, tree_search_depth);
 
            // run contractor
            return contract_helper(A_rm, b, iter, all_constraints);
        }
 
        Box linear_map(const Eigen::Matrix<zono_float, -1, -1>& A) const
        {
            // input handling
            if (A.cols() != x_lb.size())
                throw std::invalid_argument("Matrix A must have the same number of columns as the size of the Box");
 
            // declare
            Box y(A.rows());
 
            // linear map
            for (int i=0; i<A.rows(); i++)
            {
                y[i] = Interval(0, 0);
                for (int j=0; j<A.cols(); j++)
                {
                    y[i].add_assign(y[i], ((*this)[j]*A(i, j)).as_view());
                }
            }
            return y;
        }
 
        Box linear_map(const Eigen::SparseMatrix<zono_float, Eigen::RowMajor>& A) const
        {
            // input handling
            if (A.cols() != x_lb.size())
                throw std::invalid_argument("Matrix A must have the same number of columns as the size of the Box");
 
            // declare
            Box y (A.rows());
 
            // linear map
            for (int i=0; i<A.rows(); i++)
            {
                y[i] = Interval(0, 0);
                for (Eigen::SparseMatrix<zono_float, Eigen::RowMajor>::InnerIterator it(A, i); it; ++it)
                {
                    y[i].add_assign(y[i], ((*this)[it.col()]*it.value()).as_view());
                }
            }
            return y;
        }
 
        Interval dot(const Eigen::Vector<zono_float, -1>& x) const
        {
            // input handling
            if (x.size() != x_lb.size())
                throw std::invalid_argument("Vector x must have the same size as the Box");
 
            // declare
            Interval y(0, 0);
 
            // linear map
            for (int i=0; i<this->x_lb.size(); i++)
                y.add_assign(y, ((*this)[i]*x(i)));
            return y;
        }
 
        void permute(const Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic>& P)
        {
            this->x_lb = P*this->x_lb;
            this->x_ub = P*this->x_ub;
        }
 
        std::string print() const
        {
            std::stringstream ss;
            ss << "Box: " << std::endl;
            for (size_t i=0; i<x_lb.size(); i++)
            {
                ss << "  " << (*this)[i] << std::endl;
            }
            return ss.str();
        }
 
        friend std::ostream& operator<<(std::ostream& os, const Box& box)
        {
            os << box.print();
            return os;
        }
 
    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)
        {
            for (int i=0; i<iter; i++)
            {
                // loop through constraints
                for (const int k : constraints)
                {
                    // forward propagate
                    Interval y(0, 0);
                    std::vector<int> cols; // keeping track of columns
                    for (Eigen::SparseMatrix<zono_float, Eigen::RowMajor>::InnerIterator it(A, k); it; ++it)
                    {
                        y.add_assign(y, (*this)[it.col()].to_interval()*it.value());
                        cols.push_back(static_cast<int>(it.col()));
                    }
 
                    // check validity
                    if (!y.contains(b(k)))
                        return false;
                    else
                        y = Interval(b(k), b(k));
 
                    // backward propagate
                    for (const int col : cols)
                    {
                        Interval x = y; // init
                        zono_float a_col=1;
                        for (Eigen::SparseMatrix<zono_float, Eigen::RowMajor>::InnerIterator it(A, k); it; ++it)
                        {
                            if (it.col() != col)
                            {
                                x.add_assign(x, (*this)[it.col()].to_interval()*(-it.value()));
                            }
                            else
                            {
                                a_col = it.value();
                            }
                        }
 
                        // update interval
                        (*this)[col].intersect_assign((*this)[col], (x * (1/a_col)).as_view());
                    }
                }
            }
 
            return true; // constraints valid
        }
 
        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)
        {
            // immediately copy over constraints and vars
            vars.insert(new_vars.begin(), new_vars.end());
 
            // find new constraints
            std::set<int> new_constraints;
            for (const int i : new_vars)
            {
                for (Eigen::SparseMatrix<zono_float>::InnerIterator it(A, i); it; ++it)
                {
                    if (!constraints.count(static_cast<int>(it.row())))
                    {
                        new_constraints.insert(static_cast<int>(it.row()));
                    }
                }
            }
 
            // find new vars
            std::set<int> new_new_vars;
            for (const int i : new_constraints)
            {
                for (Eigen::SparseMatrix<zono_float, Eigen::RowMajor>::InnerIterator it(A_rm, i); it; ++it)
                {
                    if (!vars.count(static_cast<int>(it.col())))
                    {
                        new_new_vars.insert(static_cast<int>(it.col()));
                    }
                }
            }
 
            // add new constraints to set
            constraints.insert(new_constraints.begin(), new_constraints.end());
 
            // recurse if able
            depth++;
            if (depth < max_depth && new_new_vars.empty())
                get_vars_cons(A, A_rm, constraints, vars, new_new_vars, depth, max_depth);
        }
    };
 
} // namespace ZonoOpt::detail
 
 
#endif