cone.hpp

Go to the documentation of this file.

/*****************************************************************************
    
    cone.hpp -- a representation of a polyhedral cone and operations on it.

    This file is a part of the Arageli library.

    Copyright (C) Sergey S. Lyalin, 2006
    University of Nizhni Novgorod, Russia, 2005--2006

*****************************************************************************/

#ifndef _ARAGELI_cone_hpp_
#define _ARAGELI_cone_hpp_

#include "config.hpp"

#include "factory.hpp"
#include "frwrddecl.hpp"
#include "big_int.hpp"
#include "matrix.hpp"


namespace Arageli
{


// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


struct fromnull_t {};

const fromnull_t fromnull = fromnull_t();

struct fromspace_t {};

const fromspace_t fromspace = fromspace_t();

struct fromineq_t {};

const fromineq_t fromineq = fromineq_t();

struct fromeq_t {};

const fromeq_t fromeq = fromeq_t();

struct fromgen_t {};

const fromgen_t fromgen = fromgen_t();

struct frombasis_t {};

const frombasis_t frombasis = frombasis_t();

struct fromdual_t {};

const fromdual_t fromdual = fromdual_t();



class cannot_represent_cone : public exception {};


template <typename Out, typename T, typename M, typename CFG>
void output_list (Out& out, const cone<T, M, CFG>& x);


template <typename T, typename M> class cone_default_config {};



template <typename T, typename M, typename CFG>
class cone
{
public:

    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


    typedef T inequation_element_type;
    typedef T equation_element_type;
    typedef T generatrix_element_type;
    typedef T basis_element_type;

    typedef M inequation_type;
    typedef M equation_type;
    typedef M generatrix_type;
    typedef M basis_type;

    typedef typename M::difference_type dim_type;
    typedef typename M::difference_type difference_type;
    typedef typename M::size_type size_type;



    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++




    cone () : state(STATE_ALL) {}



    template <typename D>
    cone (const D& dim, const fromnull_t&)
    :   basis_m(0, dim),
        generatrix_m(0, dim),
        state(STATE_PARAMETRIC)
    { if(Arageli::is_null(dim))set_flag(STATE_IMPLICIT); }
         


    template <typename D>
    cone (const D& dim, const fromspace_t&)
    :   inequation_m(0, dim),
        equation_m(0, dim),
        state(STATE_IMPLICIT)
    { if(Arageli::is_null(dim))set_flag(STATE_PARAMETRIC); }
         


    template <typename M1>
    cone (const M1& ineqmat, const fromineq_t&)
    :   inequation_m(ineqmat),
        equation_m(0, ineqmat.ncols()),
        state(STATE_IMPLICIT_VALID)
    {}



    template <typename M1>
    cone (const M1& eqmat, const fromeq_t&)
    :   inequation_m(0, eqmat.ncols()),
        equation_m(eqmat),
        state(STATE_IMPLICIT_VALID)
    {}



    template <typename M1>
    cone (const M1& genmat, const fromgen_t&)
    :   generatrix_m(genmat),
        basis_m(0, genmat.ncols()),
        state(STATE_PARAMETRIC_VALID)
    {}



    template <typename M1>
    cone (const M1& basismat, const frombasis_t&)
    :   generatrix_m(0, basismat.ncols()),
        basis_m(basismat),
        state(STATE_PARAMETRIC_VALID)
    {}



    template <typename Cone1>
    cone (const Cone1& c, const fromdual_t&);




    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



    bool is_implicit_valid () const
    {
        ARAGELI_ASSERT_1(is_correct_flags());
        return get_flag(STATE_IMPLICIT_VALID);
    }

    bool is_parametric_valid () const
    {
        ARAGELI_ASSERT_1(is_correct_flags());
        return get_flag(STATE_PARAMETRIC_VALID);
    }

    bool is_all_valid () const
    {
        ARAGELI_ASSERT_1(is_correct_flags());
        return get_flag(STATE_VALID);
    }



    bool is_implicit_normal () const
    {
        ARAGELI_ASSERT_1(is_correct_flags());
        return get_flag(STATE_IMPLICIT);
    }



    bool is_parametric_normal () const
    {
        ARAGELI_ASSERT_1(is_correct_flags());
        return get_flag(STATE_PARAMETRIC);
    }


    bool is_all_normal () const
    {
        ARAGELI_ASSERT_1(is_correct_flags());
        return get_flag(STATE_ALL);
    }




    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++




    void validate_implicit () const
    {
        ARAGELI_ASSERT_1(is_correct_flags());
        if(!get_flag(STATE_IMPLICIT_VALID))force_validate_implicit();
    }


    void validate_parametric () const
    {
        ARAGELI_ASSERT_1(is_correct_flags());
        if(!get_flag(STATE_PARAMETRIC_VALID))force_validate_parametric();
    }

    void validate_all () const
    {
        ARAGELI_ASSERT_1(is_correct_flags());

        if(!get_flag(STATE_IMPLICIT_VALID))force_validate_implicit();
        if(!get_flag(STATE_PARAMETRIC_VALID))force_validate_parametric();
    }
    

    void normalize_implicit () const
    {
        ARAGELI_ASSERT_1(is_correct_flags());
        if(!get_flag(STATE_IMPLICIT))force_normalize_implicit();
    }


    void normalize_parametric () const
    {
        ARAGELI_ASSERT_1(is_correct_flags());
        if(!get_flag(STATE_PARAMETRIC))force_normalize_parametric();
    }
    
    void normalize_all () const
    {
        ARAGELI_ASSERT_1(is_correct_flags());
        if(!get_flag(STATE_ALL))force_normalize_all();
    }




    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


    

    void set_valid_implicit ()
    {
        ARAGELI_ASSERT_1(is_correct_flags());
        set_flag(STATE_IMPLICIT_VALID);
    }


    void set_valid_parametric ()
    {
        ARAGELI_ASSERT_1(is_correct_flags());
        set_flag(STATE_PARAMETRIC_VALID);
    }


    void set_valid_all ()
    {
        ARAGELI_ASSERT_1(is_correct_flags());
        set_flag(STATE_VALID);
    }


    void set_normal_implicit ()
    {
        ARAGELI_ASSERT_1(is_correct_flags());
        set_flag(STATE_IMPLICIT);
    }


    void set_normal_parametric ()
    {
        ARAGELI_ASSERT_1(is_correct_flags());
        set_flag(STATE_PARAMETRIC);
    }


    void set_normal_all ()
    {
        ARAGELI_ASSERT_1(is_correct_flags());
        set_flag(STATE_ALL);
    }



    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


    

    bool is_pointed () const;

    bool is_subspace () const;

    bool is_space () const;

    bool is_null () const;

    bool is_point () const { return is_null(); }

    dim_type space_dim () const;

    dim_type dim () const;

    bool is_bodily () const
    { return dim() == space_dim(); }

    bool is_simplicial () const;



    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


    template <typename V1>
    bool strict_inside (const V1& x) const
    {
        normalize_implicit();
        return
            all_greater(inequation_m*x, null<typename V1::element_type>()) &&
            Arageli::is_null(equation_m*x);
    }

    template <typename V1>
    bool strict_outside (const V1& x) const
    {
        normalize_implicit();
        return
            !all_greater_equal(inequation_m*x, null<typename V1::element_type>()) ||
            !Arageli::is_null(equation_m*x);
    }

    template <typename V1>
    bool unstrict_inside (const V1& x) const
    {
        normalize_implicit();
        return
            all_greater_equal(inequation_m*x, null<typename V1::element_type>()) &&
            Arageli::is_null(equation_m*x);
    }

    template <typename V1>
    bool unstrict_outside (const V1& x) const
    {
        normalize_implicit();
        return
            !all_greater(inequation_m*x, null<typename V1::element_type>()) ||
            !Arageli::is_null(equation_m*x);
    }


    
    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



    inequation_type inequation () const;


    const equation_type& equation () const;


    generatrix_type generatrix () const;


    const basis_type& basis () const;

    const equation_type& min_ambient_equation () const;
    
    basis_type min_ambient_basis () const;
    //{ return inverse(min_ambient_equation()); }   // It's incorrect.
    
    equation_type max_embedded_equation () const;
    //{ return inverse(max_embedded_basis()); } // It's incorrect.

    const basis_type& max_embedded_basis () const;

    inequation_type min_inequation () const;

    generatrix_type min_generatrix () const;



    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


    const inequation_type& inequation_matrix () const
    { return inequation_m; }

    inequation_type& inequation_matrix ()
    {
        clear_flag(STATE_PARAMETRIC | STATE_IMPLICIT_NORMAL);
        return inequation_m;
    }

    const equation_type& equation_matrix () const
    { return equation_m; }

    equation_type& equation_matrix ()
    {
        clear_flag(STATE_PARAMETRIC | STATE_IMPLICIT_NORMAL);
        return equation_m;
    }

    const generatrix_type& generatrix_matrix () const
    { return generatrix_m; }
    
    generatrix_type& generatrix_matrix ()
    {
        clear_flag(STATE_IMPLICIT | STATE_PARAMETRIC_NORMAL);
        return generatrix_m;
    }

    const basis_type& basis_matrix () const
    { return basis_m; }
    
    basis_type& basis_matrix ()
    {
        clear_flag(STATE_IMPLICIT | STATE_PARAMETRIC_NORMAL);
        return basis_m;
    }



    void pack ();

    void keep_only_implicit ();

    void keep_only_parametric ();
    
    void assign_null ();

    matrix<bool> relation () const
    {
        validate_all();
        return generatrix_m*transpose(inequation_m);
    }

    
    // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


    template <typename Cone1>
    cone& intersection (const Cone1& x);

    template <typename Cone1>
    cone& conic_union (const Cone1& x);


private:

    void force_validate_implicit () const;
    void force_validate_parametric () const;
    void force_normalize_implicit () const;
    void force_normalize_parametric () const;
    void force_normalize_all () const;

    template <typename M1, typename M2>
    inline void extend_rep_matrix (M1& res, const M2& x) const
    {
        // WARNING! Here we add 2*x.nrows() rows in res,
        // but only (x.nrows()+1) is enough.
        res.insert_matrix_bottom(x);
        res.insert_matrix_bottom(-x);
    }

    template <typename M1, typename M2>
    inline void fast_extend_rep_matrix (M1& res, const M2& x) const
    {
        res.insert_matrix_bottom(x);
        res.insert_matrix_bottom(-x);
    }

    enum State
    {
        STATE_IMPLICIT_VALID    = 0x01,
        STATE_IMPLICIT_NORMAL   = 0x01 << 1,
        STATE_PARAMETRIC_VALID  = 0x01 << 2,
        STATE_PARAMETRIC_NORMAL = 0x01 << 3,

        STATE_VALID = STATE_IMPLICIT_VALID | STATE_PARAMETRIC_VALID,
        STATE_NORMAL = STATE_IMPLICIT_NORMAL | STATE_PARAMETRIC_NORMAL,
        STATE_IMPLICIT = STATE_IMPLICIT_VALID | STATE_IMPLICIT_NORMAL,
        STATE_PARAMETRIC = STATE_PARAMETRIC_VALID | STATE_PARAMETRIC_NORMAL,

        STATE_ALL = STATE_IMPLICIT | STATE_PARAMETRIC
    };

    bool get_flag (int st) const { return (state & st) == st; }
    void set_flag (int st) const { state = State(state | st); }
    void clear_flag (int st) const { state = State(state & ~st); }


    mutable inequation_type inequation_m;
    mutable equation_type equation_m;
    mutable generatrix_type generatrix_m;
    mutable basis_type basis_m;

    mutable State state;

    
    // We suppose that the following functions will be called
    // only in the debug purposes.
    
    bool is_correct_flags () const
    {
        return
            (get_flag(STATE_IMPLICIT_NORMAL) ?
                get_flag(STATE_IMPLICIT_VALID) : true) &&
            (get_flag(STATE_PARAMETRIC_NORMAL) ?
                get_flag(STATE_PARAMETRIC_VALID) : true);
    }

    bool is_there_valid () const
    {
        return
            get_flag(STATE_IMPLICIT_VALID) &&
                inequation_m.ncols() == equation_m.ncols() ||
            get_flag(STATE_PARAMETRIC_VALID) &&
                generatrix_m.ncols() == basis_m.ncols();
    }

    
    template <typename Out, typename T2, typename M2, typename CFG2>
    friend void output_list (Out& out, const cone<T2, M2, CFG2>& x);

};


template <typename TT, typename M, typename CFG>
struct factory<cone<TT, M, CFG> >
{
private:

    typedef cone<TT, M, CFG> T;

public:

    static const bool is_specialized = true;

    //static const T& unit ();
    
    //static const T& unit (const T& x);

    //static const T& opposite_unit ();
    
    //static const T& opposite_unit (const T& x);

    static const T& null ()
    {
        static const T null_s = T();
        return null_s;
    }
    
    static const T& null (const T& x)
    { return null(); }

};



template <typename Out, typename T, typename M, typename CFG>
inline Out& operator<< (Out& out, const cone<T, M, CFG>& x)
{ output_list(out, x); return out; }



template
<
    typename T1, typename M1, typename CFG1,
    typename T2, typename M2, typename CFG2
>
inline cone<T1, M1, CFG1> intersection
(
    cone<T1, M1, CFG1> a,
    const cone<T2, M2, CFG2>& b
)
{ return a.intersection(b); }


template
<
    typename T1, typename M1, typename CFG1,
    typename T2, typename M2, typename CFG2
>
inline cone<T1, M1, CFG1> cone_union
(
    cone<T1, M1, CFG1> a,
    const cone<T2, M2, CFG2>& b
)
{ return a.cone_union(b); }


template
<
    typename T1, typename M1, typename CFG1,
    typename T2, typename M2, typename CFG2
>
int cmp
(
    const cone<T1, M1, CFG1>& a,
    const cone<T2, M2, CFG2>& b
)
{
    // WARNING! Temporary implimentation. Maybe it doesn't work always.
    
    if(int c = cmp(a.space_dim(), b.space_dim()))return c;
    if(int c = cmp(a.dim(), b.dim()))return c;
    if(is_null(a.dim()) && is_null(b.dim()))return 0;

    ARAGELI_ASSERT_1
    (
        a.min_ambient_equation().nrows() ==
        b.min_ambient_equation().nrows()
    );

    a.normalize_implicit();
    b.normalize_implicit();

    if(int c = cmp(a.inequation_matrix().nrows(), b.inequation_matrix().nrows()))
        return -c;

    // WARNING! Too expensive!

    vector<vector<T1>, false> arows(a.inequation_matrix().nrows());
    vector<vector<T2>, false> brows(b.inequation_matrix().nrows());
    ARAGELI_ASSERT_1(arows.size() == brows.size());

    for(typename cone<T1, M1, CFG1>::size_type i = 0; i < arows.size(); ++i)
    {
        arows[i] = a.inequation_matrix().copy_row(i);
        arows[i] /= _Internal::content(arows[i].begin(), arows[i].end());
        brows[i] = b.inequation_matrix().copy_row(i);
        brows[i] /= _Internal::content(brows[i].begin(), brows[i].end());
    }

    std::sort(arows);
    std::sort(brows);

    if(int c = cmp(arows, brows))return c;
    
    if(int c = cmp(a.equation_matrix().nrows(), b.equation_matrix().nrows()))
        return -c;

    if(a.equation_matrix().nrows() == 0)return 0;
    
    if
    (
        int c = cmp
        (
            rref(matrix<rational<T1>, false>(a.equation_matrix())),
            rref(matrix<rational<T1>, false>(b.equation_matrix()))
        )
    )
        return c;

    return 0;
}


} // namesapce Arageli


#ifdef ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE
    #define ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE_cone
    #include "cone.cpp"
    #undef  ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE_cone
#endif

#endif  // #ifndef _ARAGELI_cone_hpp_

---