Zono.hpp

Go to the documentation of this file.

#ifndef ZONOOPT_ZONO_HPP_
#define ZONOOPT_ZONO_HPP_
 
#include "ConZono.hpp"
 
namespace ZonoOpt
{
 
using namespace detail;
 
class Zono : public ConZono
{
    public:
 
        // constructors
        Zono() { sharp = true; }
 
        Zono(const Eigen::SparseMatrix<zono_float>& G, const Eigen::Vector<zono_float, -1>& c,
            const bool zero_one_form=false)
        {
            set(G, c, zero_one_form);
            sharp = true;
        }
 
        // virtual destructor
        ~Zono() override = default;
 
        // set method
        void set(const Eigen::SparseMatrix<zono_float>& G, const Eigen::Vector<zono_float, -1>& c,
            bool zero_one_form=false);
 
        HybZono* clone() const override
        {
            return new Zono(*this);
        }
 
        std::unique_ptr<Zono> reduce_order(int n_o);
 
        // generator conversion between [-1,1] and [0,1]
        void convert_form() override;
 
        // display methods
        std::string print() const override;
 
    protected:
 
        bool do_is_empty(const OptSettings &settings, OptSolution* solution) const override;
 
        Box do_bounding_box(const OptSettings &settings, OptSolution* solution) override;
 
        zono_float do_support(const Eigen::Vector<zono_float, -1>& d, const OptSettings &settings,
            OptSolution* solution) override;
};
 
// forward declarations
std::unique_ptr<Zono> interval_2_zono(const Box& box);
 
std::unique_ptr<Zono> make_regular_zono_2D(zono_float radius, int n_sides, bool outer_approx=false, const Eigen::Vector<zono_float, 2>& c=Eigen::Vector<zono_float, 2>::Zero());
 
// implementation
inline void Zono::set(const Eigen::SparseMatrix<zono_float>& G, const Eigen::Vector<zono_float, -1>& c,
    const bool zero_one_form)
{
    // check dimensions
    if (G.rows() != c.size())
    {
        throw std::invalid_argument("Zono: inconsistent dimensions.");
    }
 
    // zonotope parameters
    this->G = G;
    this->c = c;
    this->nG = static_cast<int>(this->G.cols());
    this->n = static_cast<int>(this->G.rows());
    this->zero_one_form = zero_one_form;
 
    // abstract zono parameters
    this->nGc = this->nG;
    this->nGb = 0;
    this->nC = 0;
    this->Gc = this->G;
    this->Gb.resize(this->n, 0);
    this->A.resize(0, this->nG);
    this->Ac = this->A;
    this->Ab.resize(0, 0);
    this->b.resize(0);
}
 
inline void Zono::convert_form()
{
    Eigen::Vector<zono_float, -1> c;
    Eigen::SparseMatrix<zono_float> G;
 
    if (!this->zero_one_form) // convert to [0,1] generators
    {
        c = this->c - this->G*Eigen::Vector<zono_float, -1>::Ones(this->nG);
        G = 2.0*this->G;
 
        set(G, c, true);
    }
    else // convert to [-1,1] generators
    {
        c = this->c + 0.5*this->G*Eigen::Vector<zono_float, -1>::Ones(this->nG);
        G = 0.5*this->G;
 
        set(G, c, false);
    }
}
 
inline std::string Zono::print() const
{
    std::stringstream ss;
    ss << "Zono: " << std::endl;
    ss << "n: " << this->n << std::endl;
    ss << "nG: " << this->nG << std::endl;
    ss << "G: " << Eigen::Matrix<zono_float, -1, -1>(this->G) << std::endl;
    ss << "c: " << this->c << std::endl;
    ss << "zero_one_form: " << this->zero_one_form;
    return ss.str();
}
 
inline bool Zono::do_is_empty(const OptSettings &settings, OptSolution* solution) const
{
    if (this->n == 0)
        return true;
    else
        return false;
}
 
inline Box Zono::do_bounding_box(const OptSettings &settings, OptSolution* solution)
{
    // convert to [-1,1] form
    if (this->zero_one_form) this->convert_form();
 
    // init bounds
    Eigen::Vector<zono_float, -1> l = this->c;
    Eigen::Vector<zono_float, -1> u = this->c;
 
    // compute bounds
    for (int i=0; i<this->nG; ++i)
    {
        l -= this->G.col(i).cwiseAbs();
        u += this->G.col(i).cwiseAbs();
    }
 
    // return zonotope bounding box
    return {l, u};
}
 
inline std::unique_ptr<Zono> Zono::reduce_order(const int n_o)
{
    // check validity
    if (n_o < this->n)
        throw std::invalid_argument("Zono reduce_order: desired order is less than dimension of set");
 
    // trivial case
    if (this->nG <= n_o)
        return std::make_unique<Zono>(*this);
 
    // convert to [-1,1] form
    if (this->zero_one_form) this->convert_form();
 
    // sort columns by decreasing 2-norm
 
    // vector of 2-norms and indices
    std::vector<std::pair<int, zono_float>> sort_vec;
    for (int i=0; i<this->nG; ++i)
    {
        sort_vec.emplace_back(i, this->G.col(i).norm());
    }
 
    // sort
    auto comp = [](const std::pair<int, zono_float>& a, const std::pair<int, zono_float>& b) -> bool
    {
        return a.second > b.second;
    };
    std::sort(sort_vec.begin(), sort_vec.end(), comp);
 
    // zonotope to keep
    const int n_K = n_o - this->n;
    Eigen::SparseMatrix<zono_float> G_K (this->n, n_K);
    std::vector<Eigen::Triplet<zono_float>> triplets;
    for (int i=0; i < n_K; ++i)
    {
        const int k = sort_vec[i].first; // column
        for (Eigen::SparseMatrix<zono_float>::InnerIterator it(this->G, k); it; ++it)
        {
            triplets.emplace_back(it.row(), i, it.value());
        }
    }
    G_K.setFromSortedTriplets(triplets.begin(), triplets.end());
    const Zono K (G_K, this->c);
 
    // zonotope to over-approximate
    Eigen::SparseMatrix<zono_float> G_L (this->n, this->nG - n_K);
    triplets.clear();
    for (int i=n_K; i < this->nG; ++i)
    {
        const int k = sort_vec[i].first; // column
        for (Eigen::SparseMatrix<zono_float>::InnerIterator it(this->G, k); it; ++it)
        {
            triplets.emplace_back(it.row(), i-n_K, it.value());
        }
    }
    G_L.setFromSortedTriplets(triplets.begin(), triplets.end());
    Zono L (G_L, Eigen::Vector<zono_float, -1>::Zero(this->n));
 
    // get bounding box
    const auto L_R = interval_2_zono(L.bounding_box());
 
    // minkowski sum
    auto Z = minkowski_sum(K, *L_R);
 
    // check that dynamic cast is valid
    if (!Z->is_zono())
        throw std::runtime_error("In Zono::ReduceOrder, return type is not a zonotope?");
 
    // cast to zono and return
    return std::unique_ptr<Zono>(dynamic_cast<Zono*>(Z.release()));
}
 
inline zono_float Zono::do_support(const Eigen::Vector<zono_float, -1>& d, const OptSettings &settings, OptSolution* solution)
{
    if (this->zero_one_form) this->convert_form();
 
    zono_float h = d.dot(this->c);
    const Eigen::Matrix<zono_float, -1, -1> Gd = this->G.toDense();
    for (int i=0; i<this->nG; ++i)
    {
        h += std::abs(d.dot(Gd.col(i)));
    }
    return h;
}
 
} // namespace ZonoOpt
 
#endif