coverrect.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file coverrect.h
///
/// Covering a rotated rectangular set or a paralellepiped with cubes
/// with respect to a fixed grid in R^n.
///
/// @author Pawel Pilarczyk, Tomasz Kapela
///
/////////////////////////////////////////////////////////////////////////////
 
// Copyright (C) 1997-2014 by Pawel Pilarczyk and Tomasz Kapela.
//
// This file is part of my research software package. This is free software:
// you can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// This software is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this software; see the file "license.txt". If not,
// please, see <https://www.gnu.org/licenses/>.
 
// Started on July 16, 2012. Last revision: July 8, 2014.
 
 
#ifndef _CMGRAPHS_COVERRECT_H_
#define _CMGRAPHS_COVERRECT_H_
 
 
// include some standard C++ header files
#include <iostream>
#include <algorithm>
#include <string>
#include <new>
 
// include selected header files from the CHomP library
#include "chomp/system/config.h"
#include "chomp/system/textfile.h"
#include "chomp/struct/autoarray.h"
 
// include all the necessary header files from the CAPD library
#ifdef _CAPD29_
#include "capd/vectalg/lib.h"
#include "capd/dynset/lib.h"
#include "capd/dynsys/lib.h"
#else
//#include "capd/capdlib.h"
#define CAPD_USER_NAMESPACE capd
#define CAPD_DEFAULT_DIMENSION spaceDim
#include "capd/fdcapdlib.h"
#undef CAPD_USER_NAMESPACE
#undef CAPD_DEFAULT_DIMENSION
#endif
 
// include local header files
#include "typeintv.h"
#include "typerect.h"
 
 
// --------------------------------------------------
// ----------------- COVER RECT SET -----------------
// --------------------------------------------------
 
/// Covers a rectangular set by cubes with respect
/// to a uniform subdivision of a rectangular area
/// defined by the offset of the lower left corner
/// and the width of the area.
/// \param rectSet - the rectangular set to be covered
/// \param offset - the real offset of the central box in the phase space
/// \param width - the real width of the central box in the phase space
/// \param intWidth - the width of the grid on the central box
/// in each direction
/// \param dim - the dimension of the phase space
/// \param image - the set of cubes in which to store the cubes
/// that cover the given rectangular set
/// \param left - the left bound on grid elements in the image
/// \param right - the right bound on grid elements in the image
/// \param codomain - a codomain (possibly very large) to which
/// the image should be restricted
/// \param disjoint - a set of cubes from which the image
/// should be disjoint (these cubes are not added even if
/// they might actually intersect the rectangular set)
 
template <class DoubleArray1, class DoubleArray2,
        class CoordType, class CubSetType>
void coverRectSet (const RectSetType &rectSet,
        const DoubleArray1 &offset, const DoubleArray2 &width,
        const CoordType *intWidth, int dim,
        CubSetType &image,
        const CoordType *left, const CoordType *right,
        const CubSetType *codomain, const CubSetType *disjoint)
{
        // extract the type of a single cube
        typedef typename CubSetType::value_type CubeType;
 
        // extract the types related to the rectangular sets
        typedef capd::C0Rect2Set C0SetType;
        typedef typename C0SetType::MatrixType MatrixType;
        typedef typename C0SetType::VectorType VectorType;
 
        // iterate all the cubes in the given rectangular box;
        // there is no need to check any cubes outside of this area,
        // e.g. because of phase space cropping;
        // each cube is represented by its minimal (lower left etc.) corner,
        // e.g., the cube [0,1]x[3,4]x[-2,-1] is represented by (0,3,-2)
        chomp::homology::tRectangle<CoordType> r (left, right, dim);
        const CoordType *coord;
 
        // C0RectSet has representation x + C*r0 + B*r;
        // first we transform it into two representations
        // x + C * w1   where w1 = r0 + (C_1*B)*r,
        // x + B * w2   where w2 = r + (B_1*C)*r0;
        // then for given box Z we transform it
        // to those representations computing
        // z1 = C_1 * (Z - x),
        // z2 = B_1 * (Z - x);
        // finally, if w1 is disjoint from z1 or w2 is disjoint from z2
        // then box Z is disjoint from rectSet.
        using capd::matrixAlgorithms::gaussInverseMatrix;
        MatrixType C_1 = gaussInverseMatrix (
#ifdef PREDICTOR_CORRECTOR
                rectSet. corrector. get_C ()
#else
                rectSet. get_C ()
#endif
        );
        MatrixType B_1 = gaussInverseMatrix (
#ifdef PREDICTOR_CORRECTOR
                rectSet. corrector. get_B ()
#else
                rectSet. get_B ()
#endif
        );
 
        VectorType w1 =
#ifdef PREDICTOR_CORRECTOR
                rectSet. corrector. get_r0 () +
                (C_1 * rectSet. corrector. get_B ()) *
                rectSet. corrector. get_r ();
#else
                rectSet. get_r0 () +
                (C_1 * rectSet. get_B ()) * rectSet. get_r ();
#endif
        VectorType w2 =
#ifdef PREDICTOR_CORRECTOR
                rectSet. corrector. get_r () +
                (B_1 * rectSet. corrector. get_C ()) *
                rectSet. corrector. get_r0 ();
#else
                rectSet. get_r () +
                (B_1 * rectSet. get_C ()) * rectSet. get_r0 ();
#endif
 
        // get the center of the rectangular set
        VectorType x =
#ifdef PREDICTOR_CORRECTOR
                rectSet. corrector. get_x ();
#else
                rectSet. get_x ();
#endif
 
        // prepare non-initialized vectors to be used across the while loop
        VectorType Z (dim, false);
        VectorType Z_x (x. dimension (), false);
        VectorType z1 (x. dimension (), false);
        VectorType z2 (x. dimension (), false);
 
        while ((coord = r. get ()) != 0)
        {
                // create a cube that will be considered
                CubeType q (coord, dim);
 
                // if the cube is not in the codomain then skip it
                if (codomain && !codomain -> check (q))
                        continue;
 
                // if the cube is in the set from which the image should be
                // disjoint then also skip it
                if (disjoint && disjoint -> check (q))
                        continue;
 
                // compute the actual box corresponding to the cube
                for (int i = 0; i < dim; ++ i)
                {
                        Z [i] = IntervalType (coord [i], coord [i] + 1);
                        Z [i] /= intWidth [i];
                        Z [i] *= width [i];
                        Z [i] += offset [i];
                }
                Z_x = Z - x;
 
                // check if the box is disjoint from the rectangular set
                z1 = C_1 * Z_x;
                if (capd::vectalg::intersectionIsEmpty (z1, w1))
                        continue;
                z2 = B_1 * Z_x;
                if(capd::vectalg::intersectionIsEmpty(z2, w2))
                        continue;
 
                // add the cube to the image if the box is not proved
                // to be disjoint from the rectangular set
                image. add (q);
        }
 
        // switch rounding mode back to the default one if necessary
        resetRounding ();
 
        return;
} /* coverRectSet */
 
 
#endif // _CMGRAPHS_COVERRECT_H_