mapgraph.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file mapgraph.h
///
/// Computation of the graph representation of a combinatorial map.
/// This file contains the definition of a function for the computation
/// of a directed graph which represents a combinatorial map
/// restricted to a given set of cubes.
///
/// @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 12, 2007. Last revision: October 9, 2022.
 
 
#ifndef _CMGRAPHS_MAPGRAPH_H_
#define _CMGRAPHS_MAPGRAPH_H_
 
 
// include selected header files from the CHomP library
#include "chomp/system/textfile.h"
#include "chomp/cubes/pointset.h"
#include "chomp/cubes/cube.h"
 
// include local header files
#include "config.h"
#include "typedefs.h"
#include "typedyns.h"
#include "utils.h"
#include "spacewrap.h"
 
 
// --------------------------------------------------
// -------------- compute cubical map ---------------
// --------------------------------------------------
 
/// Computes the cubical map on 'X' and fills out the graph 'g'.
template <class typeCubes, class typeGraph, class typeCubMap>
inline void computeMapGraph (const typeCubes &X, typeGraph &g,
        const typeCubMap &theCubMap, bool cropping = true)
{
        // set space wrapping if necessary
        localSpaceWrapping wrapping (theCubMap. intWidth ());
 
        // remember the previous cropping settings and set the new ones
        local_value<bool> local_cropping (theCubMap. cropping,
                cropping ? true : theCubMap. cropping);
        local_value<bool> local_cropped (theCubMap. cropped,
                theCubMap. cropped);
 
        // compute the images of all the cubes
        int_t nX = X. size ();
        for (int_t cur = 0; cur < nX; ++ cur)
        {
                // add a vertex to the graph that corresponds to the cube
                g. addVertex ();
 
                // compute the image of the cube
                theCubMap (X [cur], 0, &g, &X, 0, false);
 
                // show progress indicator if necessary
                if (cur && !(cur % ((spaceDim >= 3) ? 1000 : 10000)))
                {
                        chomp::homology::scon << std::setw (8) <<
                                (cur / 1000) << "k\b\b\b\b\b\b\b\b\b";
                }
        }
 
        return;
} /* computeMapGraph */
 
#ifdef CONFIG_ODEVART
/// Computes the cubical map on 'X' and fills out the graph 'g'.
/// Increases 'X' to an isolating neighborhood as much as possible
/// within the phase space.
/// A positive upper bound on the number of elements that can be added
/// must be given, and should not be extremely large (e.g. 10000 might be
/// a good choice).
template <class typeCubes, class typeGraph, class typeCubMap>
inline void computeMapGraphIsol (const typeCubes &X, typeGraph &g,
        const typeCubMap &theCubMap, typeCubes &Xadd, int_t maxAdd)
{
        using chomp::homology::sbug;
 
        // make sure that the max add limit is non-negative
        if (maxAdd < Xadd. size ())
                maxAdd = Xadd. size ();
 
        // set space wrapping if necessary
        localSpaceWrapping wrapping (theCubMap. intWidth ());
 
        // remember the previous cropping settings and set the new ones
        local_value<bool> local_cropping (theCubMap. cropping, true);
        local_value<bool> local_cropped (theCubMap. cropped, false);
 
        // compute the images of all the cubes in X and construct Xexit
        spcCubes Xexit;
        int_t nX = X. size ();
        for (int_t cur = 0; cur < nX; ++ cur)
        {
                // add a vertex to the graph that corresponds to the cube
                g. addVertex ();
 
                // compute the image of the cube (without restrictions,
                // with just cropping to the phase space)
                spcCubes img;
                theCubMap (X [cur], &img, 0, 0, 0, false);
 
                // add cubes to Xexit if necessary
                // and add corresponding edges to the graph
                int_t nImg (img. size ());
                for (int_t i = 0; i < nImg; ++ i)
                {
                        const spcCube &q (img [i]);
                        int_t n (X. getnumber (q));
                        if (n >= 0)
                        {
                                g. addEdge (n);
                        }
                        else
                        {
                                n = Xexit. add (q);
                                g. addEdge (nX + maxAdd + n);
                        }
                }
 
                // show progress indicator if necessary
                if (cur && !(cur % ((spaceDim >= 3) ? 1000 : 10000)))
                {
                        chomp::homology::scon << std::setw (8) <<
                                (cur / 1000) << "k\b\b\b\b\b\b\b\b\b";
                }
        }
        if (theCubMap. cropped)
        {
                sbug << "(cropped) ";
                theCubMap. cropped = false;
        }
        sbug << "(XexitBegin " << Xexit. size () << ") ";
 
        // compute the images of all the cubes in Xexit and move those cubes
        // to X whose images intersect X
        chomp::homology::hashedset<int_t> addNumbers;
        bool modified = false;
        int_t processed = nX;
        do
        {
                modified = false;
                for (int_t cur = 0; cur < Xexit. size (); ++ cur)
                {
                        // skip the cube if it is already in Xadd
                        if (addNumbers. check (nX + maxAdd + cur))
                                continue;
 
                        // show progress indicator if necessary
                        ++ processed;
                        if ((spaceDim >= 3) && !(processed % 1000))
                        {
                                chomp::homology::scon << std::setw (8) <<
                                        (processed / 1000) <<
                                        "k\b\b\b\b\b\b\b\b\b";
                        }
 
                        // compute the image of the cube
                        // (cropped to the phase space)
                        spcCubes img;
                        theCubMap (Xexit [cur], &img, 0, 0, 0, false);
 
                        // check if the image intersects X or Xadd
                        bool intersects = false;
                        int_t nImg (img. size ());
                        for (int_t i = 0; i < nImg; ++ i)
                        {
                                const spcCube &q (img [i]);
                                if (X. check (q) || Xadd. check (q))
                                {
                                        intersects = true;
                                        break;
                                }
                        }
 
                        // if the image does not intersect X nor Xadd
                        // then nothing more needs to be done
                        if (!intersects)
                                continue;
 
                        // if intersets then schedule the cube to be
                        // added to X, add the cubes in its image to Xexit,
                        // and add the corresponding edges to the graph
                        if (Xadd. size () >= maxAdd)
                                throw "Too many cubes must be added to X.";
                        modified = true;
                        Xadd. add (Xexit [cur]);
                        addNumbers. add (nX + maxAdd + cur);
                        g. addVertex ();
                        for (int_t i = 0; i < nImg; ++ i)
                        {
                                const spcCube &q (img [i]);
                                int_t n (X. getnumber (q));
                                if (n >= 0)
                                {
                                        g. addEdge (n);
                                        continue;
                                }
                                n = Xadd. getnumber (q);
                                if (n >= 0)
                                {
                                        g. addEdge (nX + n);
                                        continue;
                                }
                                n = Xexit. add (q);
                                g. addEdge (nX + maxAdd + n);
                        }
                }
        } while (modified);
 
        sbug << "(XexitEnd " << Xexit. size () << ") ";
        sbug << "(" << g. countVertices () << "v," <<
                g. countEdges () << "e) ";
        if (!Xadd. empty ())
                sbug << "(Xadd " << Xadd. size () << ") ";
 
        // correct the outgoing edges in the graph 'g':
        // change the numbers corresponding to Xadd and eliminate
        // those corresponding to the remainder of Xexit
        typeGraph h;
        int_t XsizeNew (X. size () + Xadd. size ());
        for (int_t v = 0; v < XsizeNew; ++ v)
        {
                h. addVertex ();
                for (int_t e = 0; e < g. countEdges (v); ++ e)
                {
                        int_t w (g. getEdge (v, e));
                        if (w < XsizeNew)
                        {
                                h. addEdge (w);
                        }
                        else
                        {
                                int_t pos (addNumbers. getnumber (w));
                                if (pos >= 0)
                                        h. addEdge (nX + pos);
                        }
                }
        }
        g. swap (h);
 
        sbug << "(" << g. countVertices () << "v," <<
                g. countEdges () << "e) ";
 
        return;
} /* computeMapGraphIsol */
#endif // CONFIG_ODEVART
 
 
#endif // _CMGRAPHS_MAPGRAPH_H_