mapcomp.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file mapcomp.h
///
/// Map computation routines.
/// 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, 2006. Last revision: March 16, 2016.
 
 
#ifndef _CMGRAPHS_MAPCOMP_H_
#define _CMGRAPHS_MAPCOMP_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 "utils.h"
 
 
// --------------------------------------------------
// ----------------- MapComputation -----------------
// --------------------------------------------------
 
/// A generic map computation routine that computes a rigorous
/// cubical multivalued map based on a function that computes the image
/// of intervals using interval arithmetic. The "mapcomp" class must have
/// a 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 MapComputation
{
public:
        /// The type of the underlying interval map.
        typedef mapcomp MapCompType;
 
        /// The type of the graph used in the computation of the image.
        typedef chomp::homology::diGraph<> GraphType;
 
        /// The map used to compute an interval enclosure of the image
        /// of an interval box.
        const mapcomp &M;
 
        /// The default constructor.
        /// Remembers the offset and width of the phase space,
        /// and the width (the same in each direction) in integers.
        /// Note: 'intwidth' does not have to be a power of 2.
        /// Note: The arrays of the offset and width of the phase space
        /// must exist until the object is deleted; the same for the map.
        MapComputation (const double *_offset, const double *_width,
                int _intwidth, int _subdivdepth,
                const mapcomp &_M);
 
        /// The destructor.
        ~MapComputation ();
 
        /// 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;
 
        /// Using cache for the map.
        bool cache;
 
        /// Is cropping of images to the designated box in effect?
        mutable bool cropping;
 
        /// Was the image cropped at least once? Reset this variable
        /// to "false" in order to detect image cropping again.
        mutable bool cropped;
 
        /// The maximal image diameter encountered so far.
        /// Set this variable to zero to gather meaningful information.
        static int maxImgDiam;
 
        /// The maximal image volume encountered so far.
        /// Set this variable to zero to gather meaningful information.
        static int maxImgVol;
 
        /// Returns the width of the phase space in terms of
        /// the number of boxes.
        int intWidth () const;
 
        /// Returns the subdivision depth at which the computations are done.
        int subdivDepth () const;
 
        /// A helper function that encloses an interval box
        /// in the interior of a box with integer coordinates.
        template <class IntervalArray>
        static void encloseIntervalInt (const IntervalArray &box, int dim,
                typename cubetype::CoordType *left,
                typename cubetype::CoordType *right);
 
        /// Cleans the cache of the cubical map, e.g., if the underlying
        /// interval map has changed.
        void cleanCache () const;
 
        /// Sets new (more or less restrictive) limits for the maximal
        /// allowed image diameter and volume.
        static void setAllowedImgSize (int diameter, int volume);
 
protected:
        /// The offset of the cubical rectangle.
        const double *offset;
 
        /// The width of the rectangle in each direction.
        const double *width;
 
        /// The width of the rectangle in terms of the number of cubes.
        int intwidth;
 
        /// The binary subdivision depth.
        int subdivdepth;
 
        /// Use interval arithmetic to compute the coordinate scope.
        /// The image scaled to the interval grid may be saved
        /// if an array of intervals is provided.
        /// Note: 'coord' may point to the same location
        /// as 'left' or 'right'.
        int compute (const typename cubetype::CoordType *coord, int dim,
                typename cubetype::CoordType *left,
                typename cubetype::CoordType *right,
                IntervalType *scaledImage) const;
 
        /// The maximal allowed image diameter.
        static int maxImgDiamAllowed;
 
        /// The maximal allowed image volume.
        static int maxImgVolAllowed;
 
private:
        /// The number of times the cache was used successfully.
        mutable int_t cacheused;
 
        /// The cache domain.
        mutable cubsettype computed;
 
        /// The cache left vertex.
        mutable chomp::homology::multitable<cubetype> leftcache;
 
        /// The cache right vertex.
        mutable chomp::homology::multitable<cubetype> rightcache;
 
        /// The cached interval vectors for images.
        mutable chomp::homology::multitable<IntervalType *> imagecache;
 
        /// Map distance measuring object.
        /// This is an auxiliary object for gathering information
        /// on the distance of a cube from certain subspaces
        /// in comparison to the distance of its image.
        /// Please, see the definition of this class
        /// for more details.
        MapDistance<cubetype,cubsettype> checkDistance;
 
}; /* class MapComputation */
 
// --------------------------------------------------
 
template <class mapcomp, class cubetype, class cubsettype>
int MapComputation<mapcomp,cubetype,cubsettype>::maxImgDiam (0);
 
template <class mapcomp, class cubetype, class cubsettype>
int MapComputation<mapcomp,cubetype,cubsettype>::maxImgVol (0);
 
template <class mapcomp, class cubetype, class cubsettype>
int MapComputation<mapcomp,cubetype,cubsettype>::maxImgDiamAllowed
        (maxImageDiameter);
 
template <class mapcomp, class cubetype, class cubsettype>
int MapComputation<mapcomp,cubetype,cubsettype>::maxImgVolAllowed
        (maxImageVolume);
 
// --------------------------------------------------
 
template <class mapcomp, class cubetype, class cubsettype>
inline MapComputation<mapcomp,cubetype,cubsettype>::MapComputation
        (const double *_offset, const double *_width,
        int _intwidth, int _subdivdepth, const mapcomp &_M):
        M (_M), cache (false), cropping (false), cropped (false),
        offset (_offset), width (_width),
        intwidth (_intwidth), subdivdepth (_subdivdepth),
        cacheused (0), leftcache (8192), rightcache (8192)
{
        testIntervals ();
        return;
} /* MapComputation::MapComputation */
 
template <class mapcomp, class cubetype, class cubsettype>
inline MapComputation<mapcomp,cubetype,cubsettype>::~MapComputation ()
{
        // DEBUG
//      if (cache && !computed. empty ())
//              chomp::homology::sbug << "MAP: " << computed. size () <<
//                      " images cached, used " << cacheused << " times.\n";
 
        // DEBUG: Save the map.
        if (cache && false)
        {
                std::ofstream f ("f.map");
                for (int_t i = 0; i < computed. size (); ++ i)
                {
                        typename cubetype::CoordType c [cubetype::MaxDim];
                        f << "[" << computed [i];
                        computed [i]. coord (c);
                        for (int j = 0; j < computed [i]. dim (); ++ j)
                                ++ c [j];
                        f << cubetype (c, computed [i]. dim ()) << "] ";
                        f << "[" << leftcache [i] << rightcache [i] << "]\n";
                }
        }
 
        // delete the cached scaled images
#ifdef CONFIG_IFRACTION
        for (int_t i = 0; i < computed. size (); ++ i)
                delete [] (imagecache [i]);
#endif
 
        return;
} /* MapComputation::~MapComputation */
 
// --------------------------------------------------
 
template <class mapcomp, class cubetype, class cubsettype>
template <class IntervalArray>
inline void MapComputation<mapcomp,cubetype,cubsettype>::encloseIntervalInt
        (const IntervalArray &box, int dim,
        typename cubetype::CoordType *left,
        typename cubetype::CoordType *right)
{
        typedef typename cubetype::CoordType CoordType;
 
        // make sure that the actual image is contained
        // in the INTERIOR of the box with integer coordinates
        for (int i = 0; i < dim; ++ i)
        {
                double left_b = box [i]. leftBound ();
                left [i] = static_cast<CoordType> (left_b);
                for (int j = 0; (j < 4) && (left [i] >= left_b); ++ j)
                        -- left [i];
                if ((left [i] >= left_b) || (left_b - left [i] > 2))
                {
                        throw "Interval enclosure: Left bound out of range.";
                }
 
                double right_b = box [i]. rightBound ();
                right [i] = static_cast<CoordType> (right_b);
                for (int j = 0; (j < 4) && (right [i] <= right_b); ++ j)
                        ++ right [i];
                if ((right [i] <= right_b) || (right [i] - right_b > 2))
                {
                        throw "Interval enclosure: "
                                "Right bound out of range.";
                }
        }
        return;
} /* MapComputation::encloseIntervalInt */
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MapComputation<mapcomp,cubetype,cubsettype>::compute
        (const typename cubetype::CoordType *coord, int dim,
        typename cubetype::CoordType *left,
        typename cubetype::CoordType *right,
        IntervalType *scaledImage) const
{
        using chomp::homology::auto_array;
 
        // compute the actual box corresponding to the cube
        auto_array<IntervalType> xPtr (new IntervalType [dim]);
        IntervalType *x = xPtr. get ();
        for (int i = 0; i < dim; ++ i)
        {
                x [i] = IntervalType (coord [i], coord [i] + 1);
                x [i] /= intwidth;
                x [i] *= width [i];
                x [i] += offset [i];
        }
        // swich the rounding direction to the neutral one
        resetRounding ();
 
        // compute the image of this box
        auto_array<IntervalType> yPtr
                (scaledImage ? 0 : new IntervalType [dim]);
        IntervalType *y = scaledImage ? scaledImage : yPtr. get ();
        M (x, y, dim, coord, subdivdepth);
        checkDistance (coord, dim, x, y);
 
        // measure the expansion of the map and output it to the screen
        // (temporary code for the "plasma" model, Dec 7, 2010)
        if (false)
        {
                using chomp::homology::slog;
                slog << "Exp " << cubetype (coord, dim) << ": (";
                for (int i = 0; i < dim; ++ i)
                {
                        double expansion = (y [i]. rightBound () -
                                y [i]. leftBound ()) /
                                (x [i]. rightBound () - x [i]. leftBound ());
                        slog << (i ? "," : "") << expansion;
                }
                slog << ")\n";
        }
 
        // rescale the image box (towards the integer coordinates)
        for (int i = 0; i < dim; ++ i)
        {
                y [i] -= offset [i];
                y [i] /= width [i];
                y [i] *= intwidth;
        }
        // swich the rounding direction to the neutral one
        resetRounding ();
 
        encloseIntervalInt (y, dim, left, right);
 
        return 0;
} /* MapComputation::compute */
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MapComputation<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 arrays for the range of the image box
        int dim = q. dim ();
        auto_array<coordType> leftPtr (new coordType [dim]);
        coordType *left = leftPtr. get ();
        auto_array<coordType> rightPtr (new coordType [dim]);
        coordType *right = rightPtr. get ();
 
        // check if the image of this cube has already been computed
        int_t number = cache ? computed. getnumber (q) :
                static_cast<int_t> (-1);
 
        // prepare an interval vector for scaled image (if necessary)
#ifdef CONFIG_IFRACTION
        auto_array<IntervalType> scaledImagePtr (new IntervalType [dim]);
        IntervalType *scaledImage (scaledImagePtr. get ());
#else
        IntervalType *scaledImage (0);
#endif
 
        // compute the image of this box if necessary
        if (number < 0)
        {
                // determine the cordinates of the lower left vertex
                q. coord (left);
 
                // compute the image of the box
                compute (left, dim, left, right, scaledImage);
 
                // store this image in the cache if requested to
                if (cache)
                {
                        leftcache [computed. size ()] = cubetype (left, dim);
                        rightcache [computed. size ()] =
                                cubetype (right, dim);
#ifdef CONFIG_IFRACTION
                        imagecache [computed. size ()] =
                                scaledImagePtr. release ();
#endif
                        computed. add (q);
                }
        }
 
        // use the cached values for the image of this box
        else
        {
                // retrieve the coordinates of the left and right ends
                leftcache [number]. coord (left);
                rightcache [number]. coord (right);
#ifdef CONFIG_IFRACTION
                scaledImage = imagecache [number];
#endif
 
                // increase the right coorinates in case of wrapping
                for (int i = 0; i < dim; ++ i)
                {
                        if (right [i] <= left [i])
                                right [i] += intwidth;
                }
 
                // make a note that the cached value has been used
                ++ cacheused;
        }
 
        // crop the image of this box to the given area if necessary
        if (cropping)
        {
                cropRanges (left, right, dim,
                        intwidth, throwIfCropped, cropped);
        }
 
        // make sure that the size of the image is reasonably small
        long volume = checkImageSize (left, right, dim,
                maxImgDiam, maxImgVol, maxImgDiamAllowed, maxImgVolAllowed);
 
        // compute the set of cubes that is the image of the given box;
        // if codomain is defined then compute the restriction of the box
        if (codomain)
        {
                // if the size of codomain is smaller than the volume of the
                // image box then check all the cubes if they are in the box
                if (codomain && (codomain -> size () < volume + 8))
                {
                        // prepare an array for the coordinates of cubes
                        auto_array<coordType> coordPtr (new coordType [dim]);
                        coordType *c = coordPtr. get ();
                        int_t codomSize = codomain -> size ();
 
                        // check all the cubes in the codomain
                        for (int_t i = 0; i < codomSize; ++ i)
                        {
                                // retrieve the coordinates of the cube
                                const cubetype &q = (*codomain) [i];
                                if (disjoint && disjoint -> check (q))
                                        continue;
                                q. coord (c);
 
#ifdef CONFIG_IFRACTION
                                // check the percentage of intersection
                                if (intersectionVolume (c, scaledImage,
                                        dim) < CONFIG_IFRACTION)
                                {
                                        continue;
                                }
#else
                                // check if the cube is inside the rectangle
                                bool inside = true;
                                for (int j = 0; j < dim; ++ j)
                                {
                                        if ((c [j] < left [j]) ||
                                                (c [j] >= right [j]))
                                        {
                                                inside = false;
                                                break;
                                        }
                                }
                                if (!inside)
                                        continue;
#endif
                                // add the cube to the image being computed
                                if (img)
                                        img -> add (q);
                                if (g)
                                        g -> addEdge (i);
                        }
                }
 
                // otherwise, process all the cubes in the image rectangle
                else
                {
                        chomp::homology::tRectangle<typename
                                cubetype::CoordType> r (left, right, dim);
                        const typename cubetype::CoordType *c;
                        while ((c = r. get ()) != 0)
                        {
                                cubetype q (c, dim);
                                if (disjoint && disjoint -> check (q))
                                        continue;
                                int_t n = codomain ?
                                        codomain -> getnumber (q) : -1;
                                if (n < 0)
                                        continue;
#ifdef CONFIG_IFRACTION
                                // check the percentage of intersection
                                if (intersectionVolume (c, scaledImage,
                                        dim) < CONFIG_IFRACTION)
                                {
                                        continue;
                                }
#endif
                                if (img)
                                        img -> add (q);
                                if (g)
                                        g -> addEdge (n);
                        }
                }
        }
 
        // otherwise, add all the cubes in the entire rectangular box
        else if (img)
        {
                chomp::homology::tRectangle<typename cubetype::CoordType>
                        r (left, right, dim);
                const typename cubetype::CoordType *c;
                while ((c = r. get ()) != 0)
                {
                        cubetype q (c, dim);
                        if (disjoint && disjoint -> check (q))
                                continue;
#ifdef CONFIG_IFRACTION
                        // check the percentage of intersection
                        if (intersectionVolume (c, scaledImage,
                                        dim) < CONFIG_IFRACTION)
                        {
                                continue;
                        }
#endif
                        img -> add (q);
                }
        }
 
        return 0;
} /* MapComputation::operator () */
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MapComputation<mapcomp,cubetype,cubsettype>::intWidth ()
        const
{
        return intwidth;
} /* MapComputation::intWidth */
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MapComputation<mapcomp,cubetype,cubsettype>::subdivDepth ()
        const
{
        return subdivdepth;
} /* MapComputation::subdivDepth */
 
template <class mapcomp, class cubetype, class cubsettype>
inline void MapComputation<mapcomp,cubetype,cubsettype>::cleanCache () const
{
        cubsettype empty;
        computed = empty;
        return;
} /* MapComputation::cleanCache */
 
template <class mapcomp, class cubetype, class cubsettype>
inline void MapComputation<mapcomp,cubetype,cubsettype>::setAllowedImgSize
        (int diameter, int volume)
{
//      chomp::homology::sbug << "Setting max img diam & vol to " <<
//              diameter << " and " << volume << ".\n";
        maxImgDiamAllowed = diameter;
        maxImgVolAllowed = volume;
        return;
} /* MapComputation::setAllowedImgSize */
 
 
#endif // _CMGRAPHS_MAPCOMP_H_