mapcomp_t.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file mapcomp_t.h
///
/// Map computation routines using tight enclosures (not necessarily
/// products of intervals).
/// This file contains the definition of routines for the computation
/// of combinatorial cubical multivalued maps.
///
/// @author Pawel Pilarczyk
///
/////////////////////////////////////////////////////////////////////////////
 
// Copyright (C) 1997-2014 by Pawel Pilarczyk.
//
// 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 19, 2012. Last revision: May 18, 2013.
 
 
#ifndef _CMGRAPHS_MAPCOMP_T_H_
#define _CMGRAPHS_MAPCOMP_T_H_
 
 
// include some standard C++ header files
#include <iostream>
#include <fstream>
#include <algorithm>
#include <memory>
#include <new>
 
// include selected header files from the CHomP library
#include "chomp/system/textfile.h"
#include "chomp/system/timeused.h"
#include "chomp/struct/autoarray.h"
#include "chomp/struct/digraph.h"
 
// include local header files
#include "config.h"
#include "typedefs.h"
#include "typeintv.h"
#include "mapdist.h"
#include "mapcomp.h"
#include "coverrect.h"
#include "utils.h"
 
 
// --------------------------------------------------
// -------------- MapComputationTight ---------------
// --------------------------------------------------
 
/// The generic map computation routine that computes a rigorous
/// cubical multivalued map based on a function that computes the image
/// of an interval vector as a rectangular set. The "mapcomp" class must
/// have an appropriate method called "compute" which is used to compute
/// the map in interval arithmetic.
template <class mapcomp, class cubetype,
        class cubsettype = chomp::homology::hashedset<cubetype> >
class MapComputationTight:
        public MapComputation<mapcomp,cubetype,cubsettype>
{
public:
        /// The default constructor.
        MapComputationTight (const double *_offset, const double *_width,
                int _intwidth, int _subdivdepth,
                const mapcomp &_M = mapcomp ());
 
        /// The operator for computing the image of a box as a set of boxes,
        /// as it is in a combinatorial cubical multivalued map.
        /// The image is restricted to the given codomain (if provided).
        /// The integral coefficients of cubes are transformed to real
        /// numbers according to the rectangular area defined by offset
        /// from the origin and its width, as well as the width
        /// of this area in terms of integral coordinates.
        int operator () (const cubetype &q, cubsettype *img,
                chomp::homology::diGraph<> *g, const cubsettype *codomain,
                const cubsettype *disjoint, bool throwIfCropped) const;
 
}; /* class MapComputationTight */
 
// --------------------------------------------------
 
template <class mapcomp, class cubetype, class cubsettype>
inline MapComputationTight<mapcomp,cubetype,cubsettype>::MapComputationTight
        (const double *_offset, const double *_width,
        int _intwidth, int _subdivdepth, const mapcomp &_M):
        MapComputation<mapcomp,cubetype,cubsettype>
        (_offset, _width, _intwidth, _subdivdepth, _M)
{
        return;
} /* MapComputationTight::MapComputationTight */
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MapComputationTight<mapcomp,cubetype,cubsettype>::operator ()
        (const cubetype &q, cubsettype *img,
        chomp::homology::diGraph<> *g, const cubsettype *codomain,
        const cubsettype *disjoint, bool throwIfCropped) const
{
        using chomp::homology::auto_array;
 
        typedef typename cubetype::CoordType CoordType;
 
        // prepare the coordinates of the input cube
        auto_array<CoordType> coordPtr (new CoordType [spaceDim]);
        CoordType *coord (coordPtr. get ());
        q. coord (coord);
 
        // prepare an interval vector that represents the input cube
        IntervalVectorType x (spaceDim);
        for (int i = 0; i < spaceDim; ++ i)
        {
                x [i] = IntervalType (coord [i], coord [i] + 1);
                x [i] /= this -> intwidth;
                x [i] *= this -> width [i];
                x [i] += this -> offset [i];
        }
        resetRounding ();
 
        // compute the image of the cube in terms of a rectangular set
        RectSetType rectSet (x, x); // dummy rect2set
        this -> M. compute (x, rectSet, spaceDim, coord,
                this -> subdivdepth);
 
        // compute a product of intervals that bounds the rectangular set
        IntervalVectorType boundingBox (rectSet);
 
        // rescale the bounding box (towards the integer coordinates)
        for (int i = 0; i < spaceDim; ++ i)
        {
                boundingBox [i] -= this -> offset [i];
                boundingBox [i] /= this -> width [i];
                boundingBox [i] *= this -> intwidth;
        }
        resetRounding ();
 
        // prepare arrays for integer ranges of the image cubes
        auto_array<CoordType> leftPtr (new CoordType [spaceDim]);
        CoordType *left = leftPtr. get ();
        auto_array<CoordType> rightPtr (new CoordType [spaceDim]);
        CoordType *right = rightPtr. get ();
 
        // compute integer ranges of the image cubes
        this -> encloseIntervalInt (boundingBox, spaceDim, left, right);
 
        // crop the image of this box to the given area if necessary
        if (this -> cropping)
        {
                cropRanges (left, right, spaceDim, this -> intwidth,
                        throwIfCropped, this -> cropped);
        }
 
        // make sure that the size of the image is reasonably small
        checkImageSize (left, right, spaceDim,
                this -> maxImgDiam, this -> maxImgVol,
                this -> maxImgDiamAllowed, this -> maxImgVolAllowed);
 
        // prepare an array for integer widths of the phase space
        auto_array<CoordType> intWidthPtr (new CoordType [spaceDim]);
        CoordType *intWidth (intWidthPtr. get ());
        for (int i = 0; i < spaceDim; ++ i)
                intWidth [i] = this -> intwidth;
 
        // prepare a variable for saving image cubes
        bool usingImg = img && img -> empty ();
        cubsettype localImg;
        cubsettype &image (usingImg ? *img : localImg);
 
        // cover the computed image in a tight way
        coverRectSet (rectSet, this -> offset, this -> width, intWidth,
                spaceDim, image, left, right, codomain, disjoint);
 
        // add corresponding edges to the graph if necessary
        if (g)
        {
                int_t imgSize = image. size ();
                for (int_t i = 0; i < imgSize; ++ i)
                {
                        const cubetype &qImg (image [i]);
                        int n = codomain -> getnumber (qImg);
                        if (n < 0)
                                throw "An edge to outside the codomain found.";
                        g -> addEdge (n);
                }
        }
 
        // switch rounding mode back to the default one if necessary
        resetRounding ();
 
        return 0;
} /* MapComputationTight::operator () */
 
 
#endif // _CMGRAPHS_MAPCOMP_T_H_