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-2008 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 2 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, write to the
// Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
// MA 02111-1307, USA.

// Started on February 11, 2008. Last revision: April 13, 2008.


#ifndef _CMGRAPHS_UTILS_H_
#define _CMGRAPHS_UTILS_H_


// include some standard C++ header files
#include <string>
#include <sstream>
#include <fstream>
#include <cmath>
#include <algorithm>
#include <vector>

// 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 local header files
#include "config.h"
#include "typedefs.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 */


// --------------------------------------------------
// ----------- 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;
                        double paramWidth = paramRight [i] - paramLeft [i];
                        if (curCoord [j] > 0)
                        {
                                leftMapParam [i] = paramLeft [i] +
                                        curCoord [j] * paramWidth /
                                        paramSubdiv [j];
                        }
                        if (curCoord [j] < paramSubdiv [j] - 1)
                        {
                                rightMapParam [i] = paramLeft [i] +
                                        (curCoord [j] + 1) * paramWidth /
                                        paramSubdiv [j];
                        }
                        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;

        // make a correction to the minimal requested number of parts
        if (minCount < 1)
                minCount = 1;

        // 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 = 1;
        for (int i = 0; i < paramDim; ++ i)
                volume *= paramSubdiv [i];
        volume /= minCount;

        // 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];
                double side = std::exp (std::log (volume) / (paramDim - i));
                double delta = (i == paramDim - 1) ? 0.99999 : 0.5;
                maxIter [index] = static_cast<int>
                        (paramSubdiv [index] / side + delta);
        //      chomp::homology::sbug << "Vol = " << volume <<
        //              ", side = " << side <<
        //              ", maxIter = " << maxIter [index] << "\n";
                if (maxIter [index] > paramSubdiv [index])
                        maxIter [index] = paramSubdiv [index];
                if (maxIter [index] <= 0)
                        maxIter [index] = 1;

                // divide the volume by the approximate side length
                volume /= paramSubdiv [index];
                volume *= maxIter [index];
        //      chomp::homology::sbug << " x " << maxIter [index] << "\n";
        }

        // fill in the array with the subdivision coordinates
        for (int i = 0; i < paramDim; ++ i)
        {
                subCoords [i]. push_back (0);
                for (int j = 1; j < maxIter [i]; ++ j)
                {
                        subCoords [i]. push_back ((paramSubdiv [i] * j +
                                maxIter [i] / 2) / maxIter [i]);
                }
                subCoords [i]. push_back (paramSubdiv [i]);
        }

        // 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 */


// --------------------------------------------------
// ------------- ranges of coordinates --------------
// --------------------------------------------------

/// Updates the minimal and maximal coordinates of the box
/// which contains all the Morse sets.
inline void coordMinMax (const theMorseDecompositionType &morseDec,
        spcCoord *coordLeftMin, spcCoord *coordRightMax)
{
        spcCoord coord [spaceDim];
        int nSets = morseDec. count ();
        for (int n = 0; n < nSets; ++ n)
        {
                const spcCubes &theSet = morseDec [n];
                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 */


// --------------------------------------------------
// --------- 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 */


#endif // _CMGRAPHS_UTILS_H_

---