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);
 
        zono_float get_volume();
 
        // 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&, OptSolution*) const override;
 
        Box do_bounding_box(const OptSettings&, OptSolution*) override;
 
        zono_float do_support(const Eigen::Vector<zono_float, -1>& d, const OptSettings&,
            OptSolution*) 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());
 
// utilities
namespace detail
{
 
    // generated by Gemini
    template<typename T>
    void combinations_util(const std::vector<T>& elements, const size_t k, const size_t start_index,
                           std::vector<T>& current_combination,
                           std::vector<std::vector<T>>& result)
    {
        // Base case: Combination is complete (size k reached)
        if (current_combination.size() == k)
        {
            result.push_back(current_combination);
            return;
        }
 
        // Recursive step: Iterate over available elements
        for (size_t i = start_index; i < elements.size(); ++i)
        {
            // 🌟 Pruning optimization: stop if there aren't enough elements left to form a k-combination
            if (elements.size() - i < k - current_combination.size())
            {
                return;
            }
 
            // 1. Include elements[i]
            current_combination.push_back(elements[i]);
 
            // 2. Recurse (start from i + 1 to prevent repeats)
            combinations_util(elements, k, i + 1, current_combination, result);
 
            // 3. Backtrack (remove element to try the next possibility)
            current_combination.pop_back();
        }
    }
 
    // generated by Gemini
    template<typename T>
    std::vector<std::vector<T>> get_combinations(const std::vector<T>& input_set, const size_t k)
    {
        if (k > input_set.size())
        {
            return {}; // Cannot choose k elements from a smaller set
        }
 
        // Copy set elements to a vector. std::set ensures they are already sorted.
        std::vector<T> elements(input_set.begin(), input_set.end());
 
        std::vector<std::vector<T>> result;
        std::vector<T> current_combination;
 
        // Start the recursive process from the first element (index 0)
        combinations_util(elements, k, 0, current_combination, result);
        return result;
    }
 
} // namespace detail
 
 
// 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&, OptSolution*) const
{
    if (this->n == 0)
        return true;
    else
        return false;
}
 
inline Box Zono::do_bounding_box(const OptSettings&, OptSolution*)
{
    // 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[static_cast<size_t>(i)].first; // column
        for (Eigen::SparseMatrix<zono_float>::InnerIterator it(this->G, k); it; ++it)
        {
            triplets.emplace_back(static_cast<int>(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[static_cast<size_t>(i)].first; // column
        for (Eigen::SparseMatrix<zono_float>::InnerIterator it(this->G, k); it; ++it)
        {
            triplets.emplace_back(static_cast<int>(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&, OptSolution*)
{
    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;
}
 
 
    inline zono_float Zono::get_volume()
{
    // require [0,1] form
    if (!this->zero_one_form) this->convert_form();
 
    // get dense form of generator matrix
    const Eigen::Matrix<zono_float, -1, -1> Gd = this->G.toDense();
 
    // get cols of generator matrix
    std::vector<int> cols;
    cols.reserve(static_cast<size_t>(this->nG));
    for (int i=0; i<this->nG; ++i)
    {
        cols.push_back(i);
    }
 
    // get combinations
    const auto combs = get_combinations<int>(cols, static_cast<size_t>(this->n));
 
    // get determinants
    zono_float det_sum = zero;
    Eigen::LDLT<Eigen::Matrix<zono_float, -1, -1>> ldlt;
 
    for (const auto& comb : combs)
    {
        const Eigen::Matrix<zono_float, -1, -1> G_comb = Gd(Eigen::placeholders::all, comb);
        const Eigen::Matrix<zono_float, -1, -1> GT_G = G_comb.transpose()*G_comb;
        ldlt.compute(GT_G);
        const Eigen::Vector<zono_float, -1> D = ldlt.vectorD();
        det_sum += std::sqrt(D.prod());
    }
 
    // volume is det_sum
    return det_sum;
}
 
 
 
} // namespace ZonoOpt
 
#endif