mapdist.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file mapdist.h
///
/// Map computation distance checker.
/// This file contains the definition of routines for the verification
/// of whether the interval images of cubes computed while computing
/// the combinatorial map get closer to or farther from some subspaces.
/// Please, note that the routines defined here should only be used
/// from "cmsingle", not from "cmgraphs" itself. This is only a temporary
/// hack introduced into the code to gather the data in question,
/// don't rely on it to be efficient or user-friendly. Sorry about this!
///
/// @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 September 16, 2010. Last revision: January 12, 2011.
 
 
#ifndef _CMGRAPHS_MAPDIST_H_
#define _CMGRAPHS_MAPDIST_H_
 
 
// include some standard C++ header files
#include <iostream>
#include <fstream>
#include <sstream>
#include <algorithm>
#include <memory>
#include <new>
#include <cmath>
 
// include selected header files from the CHomP library
#include "chomp/system/textfile.h"
#include "chomp/system/timeused.h"
#include "chomp/homology/homology.h"
 
// include local header files
#include "config.h"
#include "typedefs.h"
#include "typeintv.h"
 
 
// --------------------------------------------------
// ------------------ MapDistance -------------------
// --------------------------------------------------
 
/// The map distance change tracker.
/// Note: This class is not thread-safe.
template <class cubetype,
        class cubsettype = chomp::homology::hashedset<cubetype> >
class MapDistance
{
public:
        /// The default constructor.
        MapDistance ();
 
        /// The destructor that saves all the collected information to files.
        ~MapDistance ();
 
        /// The operator for gathering the information on a given cube.
        void operator () (const typename cubetype::CoordType *coord, int dim,
                const IntervalType *x, const IntervalType *y) const;
 
private:
        /// Should the distances be really checked? Please, set this
        /// to 'true' only if you know what you are doing. Otherwise,
        /// leave it as 'false'.
        static const bool checkDistanceChanges = false;
 
        /// Are the X and Y coordinates rotated? If so then the distance
        /// in this direction will be measured vertically,
        /// just towards the X axis.
        static const bool rotatedXY = false;
 
        /// The sets of boxes with the distance checked:
        /// closer, farther, undetermined,
        /// for each hyperplane separately.
        mutable chomp::homology::multitable<cubsettype> distCubes;
 
        /// The dimension of the phase space.
        mutable int dim;
 
        /// The consecutive number of sets for file naming purposes.
        static int number;
 
        /// Computes the number of subspaces, depending on the dimension.
        static int countSubspaces (int dim);
 
        /// Computes the lower and upper bounds for the distance
        /// of the interval box from the given subspace.
        /// Note: It can actually compute a multiple of the distance
        /// or the squared distance, but this is OK for the purpose
        /// of comparing the distances if they get smaller or larger.
        static void distance (int subspace, const IntervalType *box, int dim,
                double &distMin, double &distMax);
 
}; /* class MapDistance */
 
// --------------------------------------------------
 
template <class cubetype, class cubsettype>
int MapDistance<cubetype,cubsettype>::number = 0;
 
// --------------------------------------------------
 
template <class cubetype, class cubsettype>
inline MapDistance<cubetype,cubsettype>::MapDistance (): dim (0)
{
        return;
} /* MapDistance::MapDistance */
 
template <class cubetype, class cubsettype>
inline MapDistance<cubetype,cubsettype>::~MapDistance ()
{
        if (!checkDistanceChanges)
                return;
 
        int_t countCubes = 0;
        int nSubspaces = countSubspaces (dim);
        for (int s = 0; s < nSubspaces; ++ s)
        {
                const char *changeFile [] =
                        {"c", "f", "u"};
                const char *changeDesc [] =
                        {"is closer to ", "is farther from ",
                                "is neither closer to nor farther from "};
                for (int n = 0; n < 3; ++ n)
                {
                        const cubsettype &cubes = distCubes [3 * s + n];
                        if (cubes. empty ())
                                continue;
                        countCubes += cubes. size ();
                        std::ostringstream filename;
                        filename << "_" << number << "_" << s <<
                                changeFile [n] << ".cub";
                        std::ofstream f (filename. str (). c_str ());
                        chomp::homology::sbug << "* Saving " <<
                                cubes. size () << " cubes to " <<
                                filename. str () << "... ";
                        f << "; Cubes whose image " << changeDesc [n] <<
                                " the subspace no. " << s << "\n"
                                "; in comparison to the distance of the "
                                " box itself from that subspace.\n";
                        f << cubes;
                        chomp::homology::sbug << "Done.\n";
                }
        }
        if (countCubes)
                ++ number;
        return;
} /* MapDistance::~MapDistance */
 
template <class cubetype, class cubsettype>
inline int MapDistance<cubetype,cubsettype>::countSubspaces (int dim)
{
        // return 2^d - d - 1, which is the number of possible subsets
        // that have at least 2 elements each, and which correspond
        // to the subspaces of interest
        int subspaces = (1 << dim) - dim - 1;
 
        // add one more computation in dim > 2: the minimal distance
        // from the codimension-1 hyperplanes
        if (dim > 2)
                ++ subspaces;
        return subspaces;
} /* MapDistance::countSubspaces */
 
template <class cubetype, class cubsettype>
void MapDistance<cubetype,cubsettype>::distance (int subspace,
        const IntervalType *box, int dim, double &distMin, double &distMax)
{
        // distance from the diagonal or from the X axis in the 2D case
        if ((subspace == 0) && (dim == 2))
        {
                double xLeft;
                double xRight;
                if (rotatedXY)
                {
                        xLeft = box [1]. leftBound ();
                        xRight = box [1]. rightBound ();
                }
                else
                {
                        IntervalType boxDist = box [0] - box [1];
                        xLeft = boxDist. leftBound ();
                        xRight = boxDist. rightBound ();
                }
                double xLeftAbs = std::abs (xLeft);
                double xRightAbs = std::abs (xRight);
                distMax = std::max (xLeftAbs, xRightAbs);
                if ((xLeft <= 0) && (xRight >= 0))
                        distMin = 0;
                else
                        distMin = std::min (xLeftAbs, xRightAbs);
        }
        // distance from the x=y=z line in the 3D case
        else if ((subspace == 3) && (dim == 3))
        {
                // the line has the direction vector (a,b,c)=(1,1,1);
                // the formula for the projected point is:
                // \frac{(a,b,c) \dot (x,y,z)}{a^2+b^2+c^2} (a,b,c),
                // so we need to compute the distance from this point
                // to (x,y,z); this formula was given by Hiroe;
                // in order to speed up the computations, let's compute
                // 9 times the square of the distance
                const IntervalType &x = box [0];
                const IntervalType &y = box [1];
                const IntervalType &z = box [2];
                IntervalType dist = sqr (-2 * x + y + z) +
                        sqr (x - 2 * y + z) + sqr (x + y - 2 * z);
                distMin = dist. leftBound ();
                if (distMin < 0)
                        distMin = 0;
                distMax = dist. rightBound ();
        }
        // the minimal distance from the three codim-1 subspaces
        else if ((subspace == 4) && (dim == 3))
        {
                for (int sub = 0; sub < dim; ++ sub)
                {
                        double subMin = 0;
                        double subMax = 0;
                        distance (sub, box, dim, subMin, subMax);
                        if (!sub || (distMin > subMin))
                                distMin = subMin;
                        if (!sub || (distMax < subMax))
                                distMax = subMax;
                }
        }
        // distance from the 2D subspaces in the 3D case
        else if (dim == 3)
        {
                // determine the coordinates: 0 = yz, 1 = xz, 2 = xy
                int coord0 = (subspace > 0) ? 0 : 1;
                int coord1 = (subspace < 2) ? 2 : 1;
                IntervalType boxDist = box [coord0] - box [coord1];
                double xLeft = boxDist. leftBound ();
                double xRight = boxDist. rightBound ();
                double xLeftAbs = std::abs (xLeft);
                double xRightAbs = std::abs (xRight);
                distMax = std::max (xLeftAbs, xRightAbs);
                if ((xLeft <= 0) && (xRight >= 0))
                        distMin = 0;
                else
                        distMin = std::min (xLeftAbs, xRightAbs);
        }
        // the remaining cases are not yet implemented
        else
        {
                throw "Distance measuring not implemented for the "
                        "requested case.\n"
                        "Please, set 'checkDistanceChanges' in 'mapdist.h' "
                        "to 'false'\nand re-compile your program.";
        }
        return;
} /* MapDistance::distance */
 
// --------------------------------------------------
 
template <class cubetype, class cubsettype>
inline void MapDistance<cubetype,cubsettype>::operator ()
        (const typename cubetype::CoordType *coord, int dim,
        const IntervalType *x, const IntervalType *y) const
{
        if (!checkDistanceChanges)
                return;
 
        // update the dimension
        if (!(this -> dim))
                this -> dim = dim;
 
        // check the distance changes from all the subspaces of interest
        int nSubspaces = countSubspaces (dim);
        for (int s = 0; s < nSubspaces; ++ s)
        {
                // determine the distance bounds for the original box
                double xDistMin = 0;
                double xDistMax = 0;
                distance (s, x, dim, xDistMin, xDistMax);
 
                // determine the distance bounds for its image
                double yDistMin = 0;
                double yDistMax = 0;
                distance (s, y, dim, yDistMin, yDistMax);
 
                // if the image is closer
                if (xDistMin >= yDistMax)
                {
                        distCubes [3 * s + 0]. add (cubetype (coord, dim));
                }
                // if the image is farther
                else if (xDistMax <= yDistMin)
                {
                        distCubes [3 * s + 1]. add (cubetype (coord, dim));
                }
                // if the distance change is undertermined
                else
                {
                        distCubes [3 * s + 2]. add (cubetype (coord, dim));
                }
        }
        // swich the rounding direction to the neutral one
        resetRounding ();
        return;
} /* MapDistance::operator () */
 
 
#endif // _CMGRAPHS_MAPDIST_H_