intconvex.cpp

Go to the documentation of this file.

/*****************************************************************************
    
    intconvex.cpp -- See declarations in intconvex.hpp.

    This file is a part of the Arageli library.

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

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

#include "config.hpp"

#if !defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE) || \
    defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE_intconvex)

#include <limits>

#include "type_traits.hpp"
#include "exception.hpp"
#include "big_int.hpp"
#include "vector.hpp"
#include "skeleton.hpp"

#include "intconvex.hpp"


namespace Arageli
{

template <typename Gen, typename IntGen>
void intconvex_simple (const Gen& gen, IntGen& intgen)
{
    // Note, there is a difference between notation used here and
    // in the source of the idea (Emelichev V.A., Kovalev M.M...).
    // In the book a polyhedron is extended to a cone by adding the last
    // component in coordinate. But here we extend by the first component.

    std::cout << "\nDebug output for intconvex_simple.";
    
    ARAGELI_ASSERT_0(gen.ncols() > 0);
    
    typedef typename Gen::element_type T;
    typedef typename Gen::size_type size_type;
    typedef vector<T, false> Vec;
    typedef vector<size_type, false> Idxvec;

    ARAGELI_ASSERT_0(type_traits<T>::is_integer_number);    // temporary limitation

    size_type d = gen.ncols()-1;    // dimention of the space with polyhedron

    std::cout << "\nd = " << d;

    
    // Separate vertices and directions of recession in gen.

    Gen verts;  // firstly, here the vertices from gen will be stored
    Idxvec reces;   // indeces of the directions of recession in gen
    intgen.assign(0, gen.ncols());

    for(size_type i = 0; i < gen.nrows(); ++i)
    {
        ARAGELI_ASSERT_0(!is_negative(gen(i, 0)));
        if(!is_null(gen(i, 0))) // if i-th generatrix is a vertex
            verts.insert_row(verts.nrows(), gen.copy_row(i));
        else    // if i-th generatirx is a direction of recession
        {
            reces.push_back(i);
            intgen.insert_row(intgen.nrows(), gen.copy_row(i));
        }
    }

    std::cout
        << "\nverts.nrows() = " << verts.nrows()
        << "\nreces.size() = " << reces.size()
        << "\nreces = " << reces;


    // Build in verts parametric representaion (points) of an intersection
    // of set Q' = closing(union(Q, q^1, ..., q^t)) (p. 107, bottom)
    // and hyperplane x0 = 1.

    ARAGELI_ASSERT_0
    (
        verts.nrows()*pow(big_int(2), reces.size()) <=
        big_int(std::numeric_limits<size_type>::max())
    );

    vector<bool, false> reces_mask(reces.size());
    size_type nreccomb = power(size_type(2), reces_mask.size());    // WARNING! Way not pow?
    size_type norigverts = verts.nrows();
    //verts.reserve(norigverts*nreccomb, verts.ncols());
    
    // pass the first combination because already done:
    if(!reces_mask.is_empty())reces_mask[0] = true;

    std::cout << "\nnreccomb = " << nreccomb;

    
    for(size_type i = 1; i < nreccomb; ++i) // along to all possible reces_mask
    {
        Vec reccomb(verts.ncols());
        for(size_type j = 0; j < reces_mask.size(); ++j)
            if(reces_mask[j])reccomb += gen.copy_row(reces[j]);

        for(size_type j = 0; j < norigverts; ++j)
            verts.insert_row
            (
                verts.nrows(),
                verts.copy_row(j) + verts(j, 0)*reccomb
            );

        // move to the next combination
        for(size_type j = 0; j < reces_mask.size(); ++j)
            if(reces_mask[j])reces_mask[j] = false;
            else
            {
                reces_mask[j] = true;
                break;
            }
    }


    // Build minimal implicit representation of the polytope
    // with vertices verts.
    // At the same time system verts is being partially reduced.

    std::cout
        << "\nbefore skeleton: verts.nrows() = " << verts.nrows();

    Gen ineqs, q, eqs;
    skeleton(verts, ineqs, q, eqs);

    std::cout
        << "\nafter skeleton: verts.nrows() = " << verts.nrows()
        << "\nineqs.nrows() = " << ineqs.nrows()
        << "\neqs.nrows() = " << eqs.nrows();


    // Build bounding box for set G (1.6).
    
    Vec minlim(d), maxlim(d);   // limits of the bounding box

    for(size_type i = 0; i < verts.nrows(); ++i)
    {
        Vec v = verts.copy_row(i);
        v /= safe_reference(v.front()); // note, round to zero

        for(size_type j = 0; j < minlim.size(); ++j)
        {
            if(v[j+1] < minlim[j])minlim[j] = v[j+1];
            if(v[j+1] > maxlim[j])maxlim[j] = v[j+1];
        }
    }

    std::cout
        << "\nminlim = " << minlim
        << "\nmaxlim = " << maxlim
        << "\nsize box = " << maxlim - minlim + 1;


    // Look at any integer point in the bounding box (BB) and determine
    // whether G includes it with the help of the system of inequation ineqs.
    // Add all such points to intgen.

    // the number of integer points in the bounding box
    T nbb = product(maxlim - minlim + unit<T>());

    std::cout << "\nnbb = " << nbb;

    Vec point = minlim;
    point.push_front(unit<T>());    // the first component is always 1

    for(T i = 0; i < nbb; ++i)  // through all integer points in the BB
    {
        // determine if G includes this point
        if
        (
            all_greater_equal(ineqs*point, null<T>()) &&
            is_null(eqs*point)
        )
            intgen.insert_row(intgen.nrows(), point);

        // move to the next point in BB
        for(size_type j = 1; j < point.size(); ++j)
            if(point[j] == maxlim[j-1])point[j] = minlim[j-1];
            else
            {
                ++point[j];
                break;
            }
    }

    std::cout << "\nintgen.nrows() = " << intgen.nrows();

}

}


#if ARAGELI_DEBUG_LEVEL > 3

#include "matrix.hpp"

namespace Arageli
{

template void intconvex_simple
(
    const matrix<big_int, true>& gen,
    matrix<big_int, true>& intgen
);


template void intconvex_simple
(
    const matrix<big_int, false>& gen,
    matrix<big_int, false>& intgen
);


template void intconvex_simple
(
    const matrix<int, true>& gen,
    matrix<int, true>& intgen
);


template void intconvex_simple
(
    const matrix<int, false>& gen,
    matrix<int, false>& intgen
);


}

#endif


#endif  // #if !defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE) || ...

---