utils.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file utils.h
///
/// Utilites and helper functions.
/// This file contains the definitions of various little useful utility
/// procedures for the Conley-Morse graphs computation software.
///
/// @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 February 11, 2008. Last revision: March 16, 2016.
 
 
#ifndef _CMGRAPHS_UTILS_H_
#define _CMGRAPHS_UTILS_H_
 
 
// include some standard C++ header files
#include <string>
#include <sstream>
#include <fstream>
#include <iostream>
#include <ios>
#include <iomanip>
#include <cmath>
#include <algorithm>
#include <vector>
#include <new>
 
// include selected header files from the CHomP library
#include "chomp/system/config.h"
#include "chomp/system/textfile.h"
#include "chomp/cubes/cube.h"
#include "chomp/multiwork/mw.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 "typedyns.h"
//#include "eigenval.h"
 
 
// --------------------------------------------------
// ----------------- file existence -----------------
// --------------------------------------------------
 
/// Returns 'true' iff the given file exists and it is allowed to read it.
inline bool fileExists (const char *filename)
{
        std::ifstream in (filename);
        return !!in;
} /* fileExists */
 
 
// --------------------------------------------------
// ---------- file name without extension -----------
// --------------------------------------------------
 
/// Returns the file name without its extension if any.
/// If the file didn't have any extension, returns the original file name.
inline std::string fileNameWithoutExt (const std::string &fileName)
{
        unsigned int pos = fileName. size ();
        if (!pos)
                return fileName;
        while (1)
        {
                int ch = fileName [-- pos];
                if ((pos == 0) || (ch == '/') || (ch == '\\') || (ch == ':'))
                        return fileName;
                if (ch == '.')
                        return std::string (fileName, 0, pos);
        }
} /* fileNameWithoutExt */
 
 
// --------------------------------------------------
// ----------- local value of a variable ------------
// --------------------------------------------------
 
/// Initializes the given variable with the value provided
/// and restores the previous value at the end of the current block.
/// The restoration of the previous value takes place upon destruction
/// of this object. This action effects in a local value of a given variable
/// valid in the current block only (e.g., until 'return').
template <class T>
class local_value
{
public:
        /// The only allowed constructor.
        local_value (T &_variable, const T &_value):
                variable (_variable), previous_value (_variable)
                {variable = _value; return;}
 
        /// The destructor which restores the original value of the variable.
        ~local_value () {variable = previous_value; return;}
 
private:
        /// A reference of the variable.
        T &variable;
 
        /// The original value of the variable.
        T previous_value;
 
}; /* class local_value */
 
 
// --------------------------------------------------
// ---------------- parameter cubes -----------------
// --------------------------------------------------
 
/// Verifies whether the given box is located at the boundary
/// of the small region which is either not at the boundary
/// of the large region, or is adjacent to some other similar region.
template <class intType1, class intType2, class intType3>
inline bool internalBoundaryPoint (const intType1 *coord,
        const intType2 *left, const intType2 *right,
        const intType3 *limit, int dim)
{
        for (int i = 0; i < dim; ++ i)
        {
                if ((coord [i] == left [i]) && (left [i] != 0))
                {
                        return true;
                }
                else if ((coord [i] == right [i] - 1) &&
                        (right [i] != limit [i]))
                {
                        return true;
                }
        }
        return false;
} /* internalBoundaryPoint */
 
/// Computes the real coordinates of the parameter cube
/// which corresponds to the given box.
/// Uses the globally defined ranges.
template <class intType>
void computeParam (const intType *curCoord,
        double *leftMapParam, double *rightMapParam)
{
        for (int i = 0; i < paramCount; ++ i)
        {
                leftMapParam [i] = paramLeft [i];
                rightMapParam [i] = paramRight [i];
                for (int j = 0; j < paramDim; ++ j)
                {
                        if (i != paramSelect [j])
                                continue;
#ifndef NONLINEAR_PARAMETERS
                        double paramWidth = paramRight [i] - paramLeft [i];
#endif
                        if (curCoord [j] > 0)
                        {
#ifdef NONLINEAR_PARAMETERS
                                leftMapParam [i] = nonlinearParameter
                                        (j, curCoord [j]);
#else
                                leftMapParam [i] = paramLeft [i] +
                                        curCoord [j] * paramWidth /
                                        paramSubdiv [j];
#endif
                        }
#ifdef CONFIG_PARAMPOINTS
                        rightMapParam [i] = leftMapParam [i];
#else
                        if (curCoord [j] < paramSubdiv [j] - 1)
                        {
#ifdef NONLINEAR_PARAMETERS
                                rightMapParam [i] = nonlinearParameter
                                        (j, curCoord [j] + 1);
#else
                                rightMapParam [i] = paramLeft [i] +
                                        (curCoord [j] + 1) * paramWidth /
                                        paramSubdiv [j];
#endif
                        }
#endif
                        break;
                }
        }
        return;
} /* computeParam */
 
/// Determines heuristically an optimal subdivision of a large parameter
/// region into the given minimal number of smaller regions.
/// Fills in the provided array of vectors with the division coordinates.
/// Allocates a rectangle for iterating this region and returns its pointer.
/// Note: The rectangle is built upon the two arrays provided; they should
/// not be deleted or otherwise the pointers stored within the rectangle
/// become invalid.
inline parRect *subRectangles (parCoord *zeroIter, parCoord *maxIter,
        std::vector<parCoord> *subCoords, int minCount)
{
        using chomp::homology::sbug;
        const bool debug = true;
 
        // use a fixed distance between boxes if requested to
        int boxStep = 0;
 
        // make a correction to the minimal requested number of parts
        if (minCount == 0)
                minCount = -4;
        if (minCount < 0)
        {
                boxStep = -minCount;
                minCount = 0;
        }
 
        // prepare the increasing order of all the sides
        // of the full parameter rectangular area
        bool taken [paramDim];
        for (int i = 0; i < paramDim; ++ i)
                taken [i] = false;
        int increasing [paramDim];
        for (int i = 0; i < paramDim; ++ i)
        {
                // find an index with the minimal value but not yet taken
                int minIndex = -1;
                for (int j = 0; j < paramDim; ++ j)
                {
                        if (taken [j])
                                continue;
                        if ((minIndex < 0) ||
                                (paramSubdiv [j] < paramSubdiv [minIndex]))
                        {
                                minIndex = j;
                        }
                }
 
                // take this index
                increasing [i] = minIndex;
                taken [minIndex] = true;
        }
 
        // compute the recommended maximal volume of each rectangular region
        double volume = 0;
        if (minCount > 0)
        {
                volume = 1;
                for (int i = 0; i < paramDim; ++ i)
                        volume *= paramSubdiv [i];
                volume /= minCount;
        }
 
        // show debugging information if necessary
        if (debug)
        {
                sbug << "Computing subdivisions - detailed information:\n";
                sbug << " boxStep = " << boxStep << ", minCount = " <<
                        minCount << ", volume = " << volume << "\n";
        }
 
        // prepare the number of parts in each direction
        // and fill in the maximal corner of the rectangle for iterating
        for (int i = 0; i < paramDim; ++ i)
        {
                // determine the recommended side length
                int index = increasing [i];
                int localBoxStep = (boxStep < paramSubdiv [index]) ?
                        boxStep : paramSubdiv [index];
                double side = volume ?
                        (std::exp (std::log (volume) / (paramDim - i))) :
                        (static_cast<double> (paramSubdiv [index]) /
                        ((paramSubdiv [index] + localBoxStep - 1) /
                        localBoxStep));
                double delta = (volume && (i == paramDim - 1)) ?
                        0.99999 : 0.5;
                maxIter [index] = static_cast<parCoord>
                        (paramSubdiv [index] / side + delta);
                if (maxIter [index] > paramSubdiv [index])
                        maxIter [index] = paramSubdiv [index];
                if (maxIter [index] <= 0)
                        maxIter [index] = 1;
                if (debug)
                {
                        sbug << " index " << index << ": remaining vol = " <<
                        volume << ", side = " << side << ", maxIter = " <<
                        maxIter [index] << "\n";
                }
 
                // divide the volume by the approximate side length
                if (volume)
                {
                        volume /= paramSubdiv [index];
                        volume *= maxIter [index];
                }
        //      chomp::homology::sbug << " x " << maxIter [index] << "\n";
        }
 
        // fill in the array with the subdivision coordinates
        int minWidth = 0;
        int maxWidth = 0;
        for (int i = 0; i < paramDim; ++ i)
        {
                if (debug)
                {
                        sbug << " dir " << i << ":";
                }
                subCoords [i]. push_back (0);
                for (int j = 1; j < maxIter [i]; ++ j)
                {
                        parCoord coord = static_cast<parCoord>
                                ((static_cast<long> (paramSubdiv [i] +
                                maxIter [i] - 1) * j + maxIter [i] - 1) /
                                maxIter [i] - j);
                        subCoords [i]. push_back (coord);
                        if (debug)
                        {
                                parCoord width = subCoords [i] [j] -
                                        subCoords [i] [j - 1] + 1;
                                if ((j == 1) || (width < minWidth))
                                        minWidth = width;
                                if ((j == 1) || (maxWidth < width))
                                        maxWidth = width;
                                sbug << " " << coord;
                        }
                }
                subCoords [i]. push_back (paramSubdiv [i]);
                if (debug)
                {
                        parCoord lastWidth = subCoords [i] [maxIter [i]] -
                                subCoords [i] [maxIter [i] - 1];
                        sbug << " " << paramSubdiv [i] << " (width " <<
                                minWidth << "-" << maxWidth << "; " <<
                                lastWidth << ")\n";
                }
        }
 
        // prepare the minimal corner of the rectangle for iterating
        for (int i = 0; i < paramDim; ++ i)
                zeroIter [i] = 0;
 
        // return the rectangle in question
        return new parRect (zeroIter, maxIter, paramDim);
} /* subRectangles */
 
 
// --------------------------------------------------
// ------ string representation of coordinates ------
// --------------------------------------------------
 
/// Returns the number of digits of the given non-negative integer number.
template <class intType>
inline int countDigits (intType x)
{
        intType tenPower = 10;
        int digits = 1;
        while (1)
        {
                if (x < tenPower)
                        return digits;
                intType nextPower = tenPower * 10;
                if (nextPower < tenPower)
                        return digits;
                tenPower = nextPower;
                ++ digits;
        }
} /* countDigits */
 
/// Generates an underscore-separated list of non-negative coordinates
/// padded with zeros to the same width assuming that their values
/// are below the given limits.
template <class intType1, class intType2>
inline std::string coord2str (const intType1 *coords,
        const intType2 *maxCoords, int dim)
{
        std::ostringstream str;
        for (int i = 0; i < dim; ++ i)
        {
                if (i)
                        str << '_';
                int maxDigits = countDigits (maxCoords [i] - 1);
                int digits = countDigits (coords [i]);
                for (int d = digits; d < maxDigits; ++ d)
                        str << '0';
                str << coords [i];
        }
        return str. str ();
} /* coord2str */
 
/// Generates an underscore-separated list of non-negative coordinates.
template <class intType>
inline std::string coord2str (const intType *coords, int dim)
{
        std::ostringstream str;
        for (int i = 0; i < dim; ++ i)
        {
                if (i)
                        str << '_';
                str << coords [i];
        }
        return str. str ();
} /* coord2str */
 
/// Generates an underscore-separated list of non-negative coordinates.
template <class cubeType>
inline std::string coord2str (const cubeType &q)
{
        using chomp::homology::auto_array;
 
        // determine the number of coordinates
        int dim = q. dim ();
        if (dim <= 0)
                return std::string ("");
 
        // prepare an array of coordinates (to be deallocated automatically)
        typedef typename cubeType::CoordType coordType;
        auto_array<coordType> coords (new coordType [dim]);
        q. coord (coords. get ());
 
        // compute the string representation of the coordinates
        return coord2str (coords. get (), dim);
} /* coord2str */
 
 
// --------------------------------------------------
// ------------- ranges of coordinates --------------
// --------------------------------------------------
 
/// A class whose objects store, update and show coordinate ranges.
class CoordMinMax
{
public:
        /// The default constructor.
        CoordMinMax ();
 
        /// The constructor from precomputed values.
        CoordMinMax (const spcCoord *leftMin, const spcCoord *rightMax);
 
        /// Returns true if at least one cube was taken for the ranges.
        bool defined () const;
 
        /// Updates the coordinate ranges for the given set of cubes.
        void operator () (const spcCubes &theSet);
 
//      /// Updates the coordinate ranges for the entire Morse decomposition.
//      void operator () (const theMorseDecompositionType &morseDec);
 
        /// The minimal coordinates.
        spcCoord coordLeftMin [spaceDim];
 
        /// The strict upper bound for the maximal coordinates.
        spcCoord coordRightMax [spaceDim];
 
}; /* class CoordMinMax */
 
inline CoordMinMax::CoordMinMax ()
{
        for (int i = 0; i < spaceDim; ++ i)
        {
                coordLeftMin [i] = 1 << finalDepth;
                coordRightMax [i] = 0;
        }
        return;
} /* CoordMinMax::CoordMinMax */
 
inline CoordMinMax::CoordMinMax (const spcCoord *leftMin,
        const spcCoord *rightMax)
{
        for (int i = 0; i < spaceDim; ++ i)
        {
                coordLeftMin [i] = leftMin [i];
                coordRightMax [i] = rightMax [i];
        }
        return;
} /* CoordMinMax::CoordMinMax */
 
inline bool CoordMinMax::defined () const
{
        return (coordLeftMin [0] < coordRightMax [0]);
} /* CoordMinMax::defined */
 
inline void CoordMinMax::operator () (const spcCubes &theSet)
{
        spcCoord coord [spaceDim];
        int nCubes = theSet. size ();
        for (int i = 0; i < nCubes; ++ i)
        {
                theSet [i]. coord (coord);
                for (int j = 0; j < spaceDim; ++ j)
                {
                        if (coordLeftMin [j] > coord [j])
                                coordLeftMin [j] = coord [j];
                        if (coordRightMax [j] <= coord [j])
                                coordRightMax [j] = coord [j] + 1;
                }
        }
        return;
} /* CoordMinMax::operator () */
 
//inline void CoordMinMax::operator ()
//      (const theMorseDecompositionType &morseDec)
//{
//      int nSets = morseDec. count ();
//      for (int n = 0; n < nSets; ++ n)
//              this -> operator () (morseDec [n]);
//      return;
//} /* CoordMinMax::operator () */
 
/// Outputs the ranges of coordinates to an output stream.
/// In case the ranges are undefined, outputs "[unknown]".
inline std::ostream &operator << (std::ostream &out,
        const CoordMinMax &range)
{
        if (!range. defined ())
        {
                out << "[unknown]";
                return out;
        }
        for (int i = 0; i < spaceDim; ++ i)
        {
                out << (i ? " x [" : "[") <<
                        range. coordLeftMin [i] << "," <<
                        range. coordRightMax [i] << "]";
        }
        return out;
} /* operator << */
 
/// Shows the real coordinates of the coordinate range.
template <class OutStream>
inline OutStream &showRealCoords (OutStream &out, const CoordMinMax &range)
{
        if (!range. defined ())
        {
                out << "[unknown]";
                return out;
        }
        const spcCoord spaceResolution (1 << finalDepth);
        for (int i = 0; i < spaceDim; ++ i)
        {
                double rangeMin = range. coordLeftMin [i] * spaceWidth [i] /
                        spaceResolution + spaceOffset [i];
                double rangeMax = range. coordRightMax [i] * spaceWidth [i] /
                        spaceResolution + spaceOffset [i];
                out << (i ? " x [" : "[") << rangeMin << "," <<
                        rangeMax << "]";
        }
        return out;
} /* showRealCoords */
 
/// Verifies if the given set of cubes is contained in the interior
/// of the phase space box.
inline bool checkIsolation (const spcCubes &theSet, int depth)
{
        spcCoord coord [spaceDim];
        int nCubes = theSet. size ();
        spcCoord rightBound = (1 << depth) - 1;
        for (int i = 0; i < nCubes; ++ i)
        {
                theSet [i]. coord (coord);
                for (int j = 0; j < spaceDim; ++ j)
                {
                        if ((coord [j] <= 0) || (coord [j] >= rightBound))
                                return false;
                }
        }
        return true;
} /* checkIsolation */
 
 
// --------------------------------------------------
// --------- subdividing phase space cubes ----------
// --------------------------------------------------
 
/// Subdivides one set of cubes to another set of cubes
/// with respect to twice finer grid.
template <class typeCubes>
void subdivideCubes (const typeCubes &X, typeCubes &result)
{
        typedef typename typeCubes::value_type tCube;
        typedef typename typeCubes::value_type::CoordType tCoord;
 
        if (X. empty ())
                return;
        int dim = X [0]. dim ();
 
        tCoord cmin [tCube::MaxDim];
        tCoord cmax [tCube::MaxDim];
        int n = X. size ();
        for (int i = 0; i < n; ++ i)
        {
                X [i]. coord (cmin);
                for (int j = 0; j < dim; ++ j)
                {
                        cmin [j] *= 2;
                        cmax [j] = cmin [j] + 2;
                }
                chomp::homology::tRectangle<tCoord> r (cmin, cmax, dim);
                const tCoord *c;
                while ((c = r. get ()) != 0)
                        result. add (tCube (c, dim));
        }
        return;
} /* subdivideCubes */
 
/// Subdivides a cube to a set of cubes with respect to twice finer grid.
template <class typeCube, class typeCubes>
void subdivideCube (const typeCube &q, typeCubes &result)
{
        typedef typename typeCubes::value_type::CoordType tCoord;
 
        int dim = q. dim ();
 
        tCoord cmin [typeCube::MaxDim];
        tCoord cmax [typeCube::MaxDim];
        q. coord (cmin);
        for (int j = 0; j < dim; ++ j)
        {
                cmin [j] *= 2;
                cmax [j] = cmin [j] + 2;
        }
        chomp::homology::tRectangle<tCoord> r (cmin, cmax, dim);
        const tCoord *c;
        while ((c = r. get ()) != 0)
                result. add (typeCube (c, dim));
        return;
} /* subdivideCube */
 
/// Embeds a set of cubes to another set of cubes
/// with respect to twice coarser grid.
template <class typeCubes>
void embedCubes (const typeCubes &X, typeCubes &result)
{
        typedef typename typeCubes::value_type tCube;
        typedef typename typeCubes::value_type::CoordType tCoord;
 
        if (X. empty ())
                return;
        int dim = X [0]. dim ();
 
        tCoord coord [tCube::MaxDim];
        int n = X. size ();
        for (int i = 0; i < n; ++ i)
        {
                X [i]. coord (coord);
                for (int j = 0; j < dim; ++ j)
                        coord [j] /= 2;
                result. add (tCube (coord, dim));
        }
        return;
} /* embedCubes */
 
/// Enhances the given set of cubes provided number of times.
inline void enhanceCubes (spcCubes &X, int enhanceTimes, int subdivDepth)
{
        if (X. empty ())
                return;
        spcCoord coordBound = 1 << subdivDepth;
        int dim = X [0]. dim ();
        spcCoord *c = new spcCoord [dim];
 
        int cur = 0;
        while (enhanceTimes --)
        {
                int curBound = X. size ();
                for (; cur < curBound; ++ cur)
                {
                        X [cur]. coord (c);
                        chomp::homology::tNeighbors<spcCoord> n (c, dim);
                        const spcCoord *nc = 0;
                        while ((nc = n. get ()) != 0)
                        {
                                bool outOfSpace = false;
                                for (int i = 0; i < dim; ++ i)
                                {
                                        if ((nc [i] < 0) ||
                                                (nc [i] >= coordBound))
                                        {
                                                outOfSpace = true;
                                                break;
                                        }
                                }
                                if (!outOfSpace)
                                        X. add (spcCube (nc, dim));
                        }
                }
        }
 
        delete [] c;
        return;
} /* enhanceCubes */
 
 
// --------------------------------------------------
// --------------- phase space cubes ----------------
// --------------------------------------------------
 
/// Returns the volume of the given number of cubes in the phase space
/// at the provided subdivision level.
inline double cubesVolume (int_t nCubes, int subdiv)
{
        const spcCoord intWidth (static_cast<spcCoord> (1) << subdiv);
        double cubeVolume (1.0);
        for (int i = 0; i < spaceDim; ++ i)
                cubeVolume *= spaceWidth [i] / intWidth;
        return nCubes * cubeVolume;
} /* cubesVolume */
 
/// Crops the ranges, marks this fact if it occurred,
/// and throws an exception if requested to.
template <class CoordType>
inline void cropRanges (CoordType *left, CoordType *right, int dim,
        int intwidth, bool throwIfCropped, bool &cropped)
{
        for (int i = 0; i < dim; ++ i)
        {
                if (left [i] < 0)
                {
                        left [i] = (right [i] > 0) ?
                                0 : (right [i] - 1);
                        if (throwIfCropped)
                                throw "Image sticks outside: left.";
                        cropped = true;
                }
                if (right [i] > intwidth)
                {
                        right [i] = (left [i] < intwidth) ?
                                intwidth : (left [i] + 1);
                        if (throwIfCropped)
                                throw "Image sticks outside: right.";
                        cropped = true;
                }
        }
        return;
} /* cropRanges */
 
/// Checks if the size of the image does not exceed the globally defined
/// maximal diameter and maximal volume (throws exceptions should this
/// be the case), and updates variables that store maximal values
/// of these quantities.
/// Returns the volume of the entire rectangular box.
template <class CoordType>
inline long checkImageSize (const CoordType *left, const CoordType *right,
        int dim, int &maxImgDiam, int &maxImgVol,
        int maxImgDiamAllowed, int maxImgVolAllowed)
{
        long volume = 1;
        for (int i = 0; i < dim; ++ i)
        {
                int size = right [i] - left [i];
                if (maxImgDiam < size)
                        maxImgDiam = size;
                if (size > maxImgDiamAllowed)
                        throw "Excessive diameter of box image.";
                volume *= size;
                if (volume > maxImgVolAllowed)
                        throw "Excessive volume of box image.";
        }
        if (maxImgVol < volume)
                maxImgVol = volume;
        return volume;
} /* checkImageSize */
 
/// Checks if the given object is in the array of the given size.
/// Returns true if yes, false otherwise.
template <class ObjectType, class ArrayType>
inline bool objectInArray (const ObjectType &theObject,
        const ArrayType &theArray, int_t arraySize)
{
        for (int_t n = 0; n < arraySize; ++ n)
        {
                if (theArray [n] == theObject)
                        return true;
        }
        return false;
} /* objectInArray */
 
/// Replaces the hashed set of objects with another one
/// that contains the same objects, but some of the objects
/// are moved selectively to the front, taking the suggested
/// step size in going through the initial hashed set.
/// The hashed set must have the property that adding an element to it
/// several times causes this element to be added only once.
/// Note: The step should be a (small) number N such that
/// groups of N consecutive objects are clustered together,
/// and the purpose of this procedure is to take representatives
/// of these clusters to the front of the hashed set.
template <class HashSetType>
inline void moveSomeToFront (HashSetType &X, int_t step)
{
        HashSetType Y;
        int_t size = X. size ();
 
        // take elements at positions that are multiples of the provided step
        const int_t step1 = 25013 * step;
        for (int_t i = step1 / 7; i < size; i += step1)
                Y. add (X [i]);
        const int_t step2 = 503 * step;
        for (int_t i = step2 / 13; i < size; i += step2)
                Y. add (X [i]);
        for (int_t i = step >> 1; i < size; i += step)
                Y. add (X [i]);
 
        // take elements scattered in the set
        const int_t step3 = 751;
        for (int_t i = step3 >> 1; i < size; i += step3)
                Y. add (X [i]);
        const int_t step4 = 73;
        for (int_t i = step4 >> 1; i < size; i += step4)
                Y. add (X [i]);
        const int_t step5 = 7;
        for (int_t i = step5 >> 1; i < size; i += step5)
                Y. add (X [i]);
 
        // take all the elements
        for (int_t i = 0; i < size; ++ i)
                Y. add (X [i]);
 
        // check that the transformation was correct
        if (Y. size () != size)
                throw "Wrong size of the new set in 'moveSomeToFront'.";
 
        // replace the result with the new set
        Y. swap (X);
 
        return;
} /* moveSomeToFront */
 
/// Computes the volume of the intersection of the unitary cube
/// with the provided coordinates with the given product of intervals.
template <class CoordType, class IntervalType>
inline double intersectionVolume (const CoordType *c,
        const IntervalType *intervalBox, int dim)
{
        double vol = 1;
        for (int i = 0; i < dim; ++ i)
        {
                // if the interval is disjoint (to the left of the cube)
                if (c [i] >= intervalBox [i]. rightBound ())
                        return 0;
                // if the interval is disjoint (to the right of the cube)
                if (c [i] + 1 <= intervalBox [i]. leftBound ())
                        return 0;
                // if the left endpoint of the interval is inside the cube
                if (c [i] < intervalBox [i]. leftBound ())
                {
                        // if the right endpoint of the interval is inside
                        if (c [i] + 1 > intervalBox [i]. rightBound ())
                        {
                                vol *= intervalBox [i]. rightBound () -
                                        intervalBox [i]. leftBound ();
                        }
                        // if the right endpoint is outside of the cube
                        else
                        {
                                vol *= c [i] + 1 - intervalBox [i]. leftBound ();
                        }
                }
                // if the left endpoint of the interval is outside the cube
                else
                {
                        // if the right endpoint of the interval is inside
                        if (c [i] + 1 > intervalBox [i]. rightBound ())
                        {
                                vol *= intervalBox [i]. rightBound () - c [i];
                        }
                        // if the right endpoint is outside of the cube
                        else
                        {
                        //      vol *= 1;
                        }
                }
        }
        return vol;
} /* intersectionVolume */
 
 
// --------------------------------------------------
// ------------------ set wrapper -------------------
// --------------------------------------------------
 
/// An imitation of an array whose entries are 1 iff the index
/// belongs to the given set.
/// This class may be useful if one needs to compute the subgraph
/// of a directed graph restricted to a given set,
/// using the subgraph method of the class diGraph.
template <class elemType, class setType>
class set2arrayWrapper
{
public:
        /// The constructor of a set wrapper.
        set2arrayWrapper (const setType &_theSet):
                theSet (_theSet)
        {
                return;
        }
 
        /// The operator [] which returns 1
        /// iff the element belongs to the set.
        int operator [] (const elemType &n) const
        {
                return theSet. check (n) ? 1 : 0;
        }
 
private:
        /// A reference to the wrapped set.
        const setType &theSet;
}; /* class set2arrayWrapper */
 
 
// --------------------------------------------------
// ------------------ graph utils -------------------
// --------------------------------------------------
 
/// Computes the graph corresponding to the restriction of the
/// combinatorial map to the provided subset.
template <class SetType, class GraphType>
inline void computeRestrictedGraph (GraphType &restr, const GraphType &g,
        const SetType &full, const SetType &subset)
{
        int_t subsetCount = subset. size ();
        for (int_t i = 0; i < subsetCount; ++ i)
        {
                restr. addVertex ();
                int_t vertex = full. getnumber (subset [i]);
                int_t nEdges = g. countEdges (vertex);
                for (int_t j = 0; j < nEdges; ++ j)
                {
                        int_t target = g. getEdge (vertex, j);
                        int_t number = subset. getnumber
                                (full [target]);
                        if (number >= 0)
                                restr. addEdge (number);
                }
        }
        return;
} /* computeRestrictedGraph */
 
 
// --------------------------------------------------
// ----------------- output numbers -----------------
// --------------------------------------------------
 
/// Creates a text representation of a floating-point number
/// at the prescribed precision.
inline std::string double2string (const double &x, int precision)
{
        std::ostringstream out;
        out. setf (std::ios::fixed);
        out. precision (precision);
        out << x;
        return out. str ();
} /* double2string */
 
 
// --------------------------------------------------
// --------------- progress indicator ---------------
// --------------------------------------------------
 
/// Shows a progress indicator to the screen as an integer number
/// followed by an additional character (e.g. 'k').
inline void showProgress (int n, char c)
{
        chomp::homology::scon << std::setw (9) << n << c <<
                "\b\b\b\b\b\b\b\b\b\b";
        return;
} /* showProgress */
 
 
#endif // _CMGRAPHS_UTILS_H_