compmdec.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file compmdec.h
///
/// Computation of Morse decompositions.
/// This file contains the definition of a function for the computation
/// of a Morse decomposition for a given map, as well as several auxiliary
/// procedures necessary to compute this Morse decomposition in an efficient
/// way by restricting the phase space to the invariant part computed
/// with respect to a coarser grid. The computed Morse decomposition
/// is processed further in order to eliminate spurious Morse sets
/// wherever possible.
///
/// @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: August 22, 2014.
 
 
#ifndef _CMGRAPHS_COMPMDEC_H_
#define _CMGRAPHS_COMPMDEC_H_
 
 
// include some standard C++ header files
#include <algorithm>
#include <new>
#include <iostream>
#include <iomanip>
#include <ios>
#include <string>
#include <fstream>
#include <sstream>
#include <cstdio>
 
// include selected header files from the CHomP library
#include "chomp/system/textfile.h"
#include "chomp/cubes/pointset.h"
#include "chomp/struct/digraph.h"
#include "chomp/struct/multitab.h"
#include "chomp/struct/hashsets.h"
#include "chomp/system/timeused.h"
 
// include local header files
#include "config.h"
#include "typedefs.h"
#include "conindex.h"
#include "morsedec.h"
#include "typedyns.h"
#include "indcache.h"
#include "invpart.h"
#include "utils.h"
#include "spacewrap.h"
#include "procmdec.h"
#include "morsecache.h"
#include "eigenval.h"
#include "dotgraph.h"
#include "datatext.h"
 
 
// --------------------------------------------------
// ---- connection recorders for inclusionGraph -----
// --------------------------------------------------
 
/// A simple class for storing connections in an array that uses
/// the "multitable" class.
class TConnTable
{
public:
        /// The constructor.
        TConnTable (int_t _n, chomp::homology::multitable
                <chomp::homology::multitable<int_t> > &_connections,
                chomp::homology::multitable<int_t> &_conncount,
                const int_t *_numTranslate);
 
        /// Adds an element v at the connection from source to target.
        void operator () (int_t source, int_t target, int_t v);
 
private:
        /// The number of vertices between the connections are recorded.
        int_t n;
 
        /// The tables of connectiong orbits; those for v -> w are stored
        /// at the index n * v + w.
        chomp::homology::multitable<chomp::homology::multitable<int_t> >
                &connections;
 
        /// The numbers of elements in each connecting orbit.
        chomp::homology::multitable<int_t> &conncount;
 
        /// The translation table for vertex numbers in the connections.
        const int_t *numTranslate;
 
}; /* class TConnTable */
 
inline TConnTable::TConnTable (int_t _n,
        chomp::homology::multitable<chomp::homology::multitable<int_t> > &
        _connections, chomp::homology::multitable<int_t> &_conncount,
        const int_t *_numTranslate):
        n (_n), connections (_connections), conncount (_conncount),
        numTranslate (_numTranslate)
{
        int_t maxcount = n * n;
        for (int_t i = 0; i < maxcount; ++ i)
                conncount [i] = 0;
        return;
} /* TConnTable::TConnTable */
 
inline void TConnTable::operator () (int_t source, int_t target, int_t v)
{
        int_t pos = n * source + target;
        connections [pos] [conncount [pos] ++] = numTranslate [v];
        return;
} /* TConnTable::operator () */
 
// --------------------------------------------------
 
/// A simple class for storing connections in an array that uses
/// the "multitable" class.
class TConnMorse
{
public:
        /// The constructor.
        TConnMorse (theMorseDecompositionType &_morseDec,
                const std::vector<int> &_setNumbers,
                const spcCubes &spcX, const int_t *_numTranslate);
 
        /// Adds an element v at the connection from source to target.
        void operator () (int_t source, int_t target, int_t v);
 
private:
        /// The Morse decomposition in which to record the connections.
        theMorseDecompositionType &morseDec;
 
        /// The numbers of Morse sets in the Morse decomposition
        /// that correspond to their indices given by source and target.
        const std::vector<int> &setNumbers;
 
        /// The set which contains cubes to add.
        const spcCubes &X;
 
        /// The translation table for vertex numbers in the connections.
        const int_t *numTranslate;
 
}; /* class TConnMorse */
 
inline TConnMorse::TConnMorse (theMorseDecompositionType &_morseDec,
        const std::vector<int> &_setNumbers, const spcCubes &spcX,
        const int_t *_numTranslate): morseDec (_morseDec),
        setNumbers (_setNumbers), X (spcX), numTranslate (_numTranslate)
{
        return;
} /* TConnMorse::TConnMorse */
 
inline void TConnMorse::operator () (int_t source, int_t target, int_t v)
{
        morseDec. addconn (setNumbers [source], setNumbers [target],
                X [numTranslate [v]]);
        return;
} /* TConnMorse::operator () */
 
 
// --------------------------------------------------
// ------------ Conley index computation ------------
// --------------------------------------------------
 
/// This small auxiliary procedure sets a dummy Conley index to the given
/// Morse set. The Conley index is set to be trivial, unless it is requested
/// that a number of a Betti number is provided to be set to 1.
inline void setDummyIndex (theMorseDecompositionType &morseDec, int n,
        int nonzero = -1)
{
        int betti [spaceDim + 1];
        for (int i = 0; i <= spaceDim; ++ i)
                betti [i] = (i == nonzero) ? 1 : 0;
        theConleyIndexType ind;
        ind. define (betti, spaceDim);
        morseDec. setindex (n, ind);
        return;
} /* setDummyIndex */
 
/// Verifies that the provided Morse set can be used as a basis
/// for a Conley index pair and computes the Conley index.
/// If this verification fails then subdivides the Morse set
/// and repeats the procedure restricted to the invariant part of the set.
/// Returns 0 if succeeds or -1 in the case of failure.
/// Makes a note of having reached the subdivision limit or having come
/// across empty invariant part in the provided variable
/// which is set to 1 in the former case and 0 in the latter
/// (otherwise it is not modified).
/// If a problem with isolation is encountered then sets "noIsolation" to 1.
/// If the index pair is computed then sets the value of the variable
/// which remembers if this is an attractor (empty exit set) or not.
inline int computeConleyIndex (theMorseDecompositionType &morseDec, int n,
        const double *offset, const double *width,
        int subdivDepth, const theMapType &theMap,
        const theCubMapType &theCubMap0, const theCubMapType &theCubMap1,
        int &subdivNonemptyInv, int &noIsolation, int &attractor)
{
        using chomp::homology::sbug;
 
        // check the images to make sure that all of them are contained
        // in the rectangular region of the given depth
        spcCubes X;
        for (int subdivisions = 0;; ++ subdivisions)
        {
                // determine the combinatorial map to use
                theCubMapType subdivCubMap (offset, width,
                        1 << (subdivDepth + subdivisions),
                        subdivDepth + subdivisions, theMap);
                if (subdivisions > 1)
                {
                        subdivCubMap. setAllowedImgSize
                                (maxImageDiameter << 3, maxImageVolume << 3);
                }
                subdivCubMap. cache = false;
                const theCubMapType &theCubMap =
                        (subdivisions > 1) ? subdivCubMap :
                        (subdivisions ? theCubMap1 : theCubMap0);
                bool crop = theCubMap. cropping;
 
                // restrict the set to its invariant part if necessary
                if (subdivisions)
                {
                        try
                        {
                                int result = invariantPart (X, theCubMap);
                                if (result < 0)
                                {
                                        sbug << "An error occurred while "
                                                "computing Inv(X).\n";
                                        subdivNonemptyInv = 1;
                                        return -1;
                                }
                                if (X. empty ())
                                {
                                        sbug << "Inv(X) turned out to be "
                                                "empty.\n";
                                        setDummyIndex (morseDec, n);
                                        subdivNonemptyInv = 0;
                                        return 0;
                                }
                        }
                        catch (const char *msg)
                        {
                                subdivNonemptyInv = 1;
                                sbug << "Failed: " << msg << "\n";
                                sbug << "Failure: Unable to compute Inv.\n";
                                return -1;
                        }
                }
 
                // determine the set of cubes for which to compute the index
                const spcCubes &theSet = subdivisions ? X : morseDec [n];
 
                // try computing the Conley index using this set
                // as a basis for the index pair
                try
                {
                        // set appropriate space wrapping
                        localSpaceWrapping wrapping
                                (1 << (subdivDepth + subdivisions));
 
                        // create an index pair object
                        theIndexPairType P (theCubMap);
 
                        // add cubes to the core part of the index pair
                        int size = theSet. size ();
                        for (int i = 0; i < size; ++ i)
                                P. add (theSet [i]);
 
                        // compute the map on the index pair
                        theCubMap. cropping = wrapping ? false : true;
                        theCubMap. cropped = false;
                        P. compute ();
                        theCubMap. cropping = crop;
 
                        // correctness check: if cropping is not allowed
                        // then an attempt of cropping should have caused
                        // an exception thrown, so this condition should
                        // never happen
                        if (!ignoreIsolationForConleyIndex &&
                                theCubMap. cropped)
                        {
                                throw "Image cropped without permission!!!";
                        }
 
                        // verify if this is an attractor
                        attractor = P. countExit () ? 0 : 1;
 
                        // create the Conley index object
                        theConleyIndexType ind;
                        ind. setIndexPair (&P);
 
                        // compute the index
                        ind. compute ();
                        ind. setIndexPair (0);
 
                        // set the index of the Morse set
                        morseDec. setindex (n, ind);
                        return 0;
                }
                catch (const char *msg)
                {
                        theCubMap. cropping = crop;
                        sbug << "Failed: " << msg << "\n";
                }
 
                // if the maximal subdivision level was reached then cancel
                if (subdivisions == refineDepth)
                {
                        sbug << "Failure: Max refinement depth " <<
                                refineDepth << " reached.\n";
                        subdivNonemptyInv = 1;
                        if (theCubMap. cropped)
                                noIsolation = 1;
                        return -1;
                }
 
                // determine the maximal allowed size of a set to refine
                int maxRefineSize =
                        subdivisions ? maxRefineSize1 : maxRefineSize0;
 
                // if the set is too large to subdivide then cancel
                if (theSet. size () > maxRefineSize)
                {
                        sbug << "Failure: Not subdividing" <<
                                (subdivisions ? " anymore" : "") <<
                                ", because the set size " <<
                                theSet. size () << "\nis beyond the max "
                                "allowed " << maxRefineSize << ".\n";
                //      sbug << "DEBUG: cropped = " << theCubMap. cropped << ".\n";
                //      sbug << "DEBUG: noIsolation = " << noIsolation << ".\n";
                        subdivNonemptyInv = 1;
                        if (theCubMap. cropped)
                                noIsolation = 1;
                        return -1;
                }
 
                // subdivide the set
                sbug << "- subdiv " << (subdivDepth + subdivisions + 1) <<
                        ": " << theSet. size () << " -> ";
                spcCubes Y;
                subdivideCubes (theSet, Y);
                spcCubes emptySet;
                X = emptySet;
                X. swap (Y);
                sbug << X. size () << " ";
        }
 
        return -1;
} /* computeConleyIndex */
 
 
// --------------------------------------------------
// -------------- chain recurrent set ---------------
// --------------------------------------------------
 
/// Computes the chain recurrent set of X.
/// Note: The resulting set may be taken the same as X
/// if one wants to replace the set with its chain recurrent set.
template <class typeCubes>
void chainRecurrentSet (const typeCubes &X, typeCubes &result,
        const chomp::homology::diGraph<> &g)
{
        using chomp::homology::diGraph;
        using chomp::homology::multitable;
 
        typeCubes Y;
 
        // make the result an empty set if the set X is empty
        if (X. empty ())
        {
                result. swap (Y);
                return;
        }
 
        // compute the components of the chain recurrent set of the graph
        multitable<int_t> compVertices (8192);
        multitable<int_t> compEnds (8192);
        diGraph<> *sccAddr = 0;
        diGraph<> *transpAddr = 0;
        int_t nSets = chomp::homology::SCC (g, compVertices, compEnds,
                sccAddr, transpAddr);
 
        // create the union of the boxes corresponding to the SCCs
        for (int_t n = 0; n < nSets; ++ n)
        {
                // determine the vertex range of this set
                int_t beginVertex = n ? compEnds [n - 1] :
                        static_cast<int_t> (0);
                int_t endVertex = compEnds [n];
 
                // add the cubes to the Morse set
                for (int_t i = beginVertex; i < endVertex; ++ i)
                        Y. add (X [compVertices [i]]);
        }
 
        result. swap (Y);
        return;
} /* chainRecurrentSet */
 
 
// --------------------------------------------------
// ----------- inv & Morse decomposition ------------
// --------------------------------------------------
 
/// Computes a Conley-Morse decomposition of ChainRecSet (X).
/// The combinatorial cubical multivalued map F is encoded
/// in the directed graph provided.
/// The Morse decomposition structure must be initially empty.
/// The arguments "offset", "width", "subdivDepth" and "theMap"
/// are only necessary for the purpose of refining trivial Morse sets.
/// The global constants "refineDepth", "maxRefineSize0" and "maxRefineSize1"
/// determine, respectively, the maximal allowed refinement depth,
/// and the maximal size of a cubical set which is refined if necessary.
/// Further simplification of the computed Morse decomposition is done
/// by joining Morse sets (see the description in "morsedec.h" for details).
/// If requested, the Conley indices of Morse sets which have at least
/// the given number of boxes are not computed (they are set to be trivial).
/// If a file name prefix is nonempty then an attempt is being made
/// to read cached Morse decomposition information from a file.
template <class typeCubes, class typeMorseDec>
void invMorseDec (const typeCubes &X, const chomp::homology::diGraph<> &g,
        typeMorseDec &morseDec,
        const double *offset, const double *width, int subdivDepth,
        const theMapType &theMap,
        const theCubMapType &theCubMap0, const theCubMapType &theCubMap1,
        int_t skipIndices,
        const std::string &cacheFileName, bool &morseDecComputed,
        std::vector<int> &wrongIndices, std::vector<int> &skippedIndices,
        std::vector<int> &attractors, bool connOrbits,
        const double *leftMapParam, const double *rightMapParam,
        int paramCount, chomp::homology::timeused &timeSubdiv)
{
        using chomp::homology::sbug;
        using chomp::homology::diGraph;
        using chomp::homology::multitable;
        using chomp::homology::hashedset;
        using chomp::homology::timeused;
 
        // do nothing if the set X is empty
        if (X. empty ())
        {
                sbug << "0.\n";
                return;
        }
 
        // compute the components of the chain recurrent set of the graph
        multitable<int_t> compVertices (8192);
        multitable<int_t> compEnds (8192);
        diGraph<> *sccAddr = 0;
        diGraph<> *transpAddr = /*(savebasmin || savebasmax) ? &gT :*/ 0;
        int_t nSets_big = chomp::homology::SCC (g, compVertices, compEnds,
                sccAddr, transpAddr);
        int_t nCubes = (nSets_big > 0) ? compEnds [nSets_big - 1] : 0;
        sbug << nCubes << " (" << nSets_big << " sets).\n";
        sbug << "Volume of the union of the Morse sets: " <<
                cubesVolume (nCubes, subdivDepth) << ".\n";
        sbug << "Morse decomp computed in " << timeSubdiv << ".\n";
        int nSets = static_cast<int> (nSets_big);
        if (static_cast<int_t> (nSets) != nSets_big)
                throw "The number of Morse sets exceeds reasonable limits.";
 
        // stop here if there are no recurrent components
        // (this correction was suggested by Marcio Gameiro; thanks!)
        if (!nSets)
                return;
 
        // prepare file names if necessary
        bool cacheIndices = !cacheFileName. empty ();
        std::string indFileNameStr = cacheFileName;
        const char *indFileName = indFileNameStr. c_str ();
        std::string lockFileNameStr = cacheFileName + ".lock";
        const char *lockFileName = lockFileNameStr. c_str ();
 
        // create the Morse decomposition data structure
        // and count the number of cubes in the SCCs
        timeused timeSetMorseDec;
        timeSetMorseDec = 0;
        sbug << "Initializing the MorseDec data structure... ";
        morseDec. setnumber (nSets);
        int_t countSCCcubes = 0;
        for (int n = 0; n < nSets; ++ n)
        {
                // determine the vertex range of this set
                int_t beginVertex = n ? compEnds [n - 1] :
                        static_cast<int_t> (0);
                int_t endVertex = compEnds [n];
                countSCCcubes += endVertex - beginVertex;
 
                // add the cubes to the Morse set
                for (int_t i = beginVertex; i < endVertex; ++ i)
                        morseDec. add (n, X [compVertices [i]]);
        }
        sbug << timeSetMorseDec << ".\n";
 
        // verify if all the Morse sets are isolated
        // from the outside of the box by at least one cube layer
        sbug << "Checking isolation of the Morse sets... ";
        bool isolated = true;
        for (int n = 0; n < nSets; ++ n)
        {
                if (checkIsolation (morseDec [n], subdivDepth))
                        continue;
                isolated = false;
                break;
        }
        sbug << (isolated ? "[isolated]\n" : "[no isol]\n");
 
        // retrieve the Conley indices from cache if possible
        timeused timeConleyIndices;
        timeConleyIndices = 0;
        bool cachedAvailable = false;
        IndexCache cached (nSets);
        if (cacheIndices && fileExists (indFileName))
        {
                if (fileExists (lockFileName))
                        sbug << "! Skipping cached indices (locked).\n";
                else if (cached. read (indFileName, morseDec) != 0)
                        sbug << "! Error while reading cached indices.\n";
                else
                {
                        sbug << "[using cached Conley indices]\n";
                        cachedAvailable = true;
                }
        }
        morseDecComputed = !cachedAvailable;
 
        // find small Morse sets whose invariant part is empty
        std::vector<int> subdivNonemptyInv (nSets, -1);
        if (!cachedAvailable)
        {
                timeused timeEmptyInv;
                timeEmptyInv = 0;
                sbug << "Checking small Morse sets for empty inv...\n";
                for (int n = 0; n < nSets; ++ n)
                {
                        // if the size of the set is too large then skip it
                        if (morseDec [n]. size () > maxRefineSize0)
                                continue;
 
                        // refine the set and check if its inv is nonempty
                        if (checkEmptyInv (morseDec [n], finalDepth,
                                offset, width, theMap,
                                theCubMap0, theCubMap1))
                        {
                                subdivNonemptyInv [n] = 0;
                                cached. emptyInv [n] = true;
                        }
                }
                sbug << "Invariant part of small Morse sets checked in " <<
                        timeEmptyInv << ".\n";
        }
 
#ifdef CONFIG_NOCONNECTIONS
        sbug << "Note: Skippnig the computation of connecting orbits.\n";
//      bool computeconn = false;
#else
        {
                // determine if it is necessary to compute the connecting
                // orbits in terms of cubical sets
                bool computeconn = connOrbits || (maxJoinSize > 0);
 
                // compute the connecting orbits for the Morse decomposition
                timeused timeConnOrb;
                timeConnOrb = 0;
 
                // prepare numbers of Morse sets with possibly nonempty inv
                std::vector<int> setNumbers1;
                for (int n = 0; n < nSets; ++ n)
                {
                        if (cachedAvailable && cached. emptyInv [n])
                                continue;
                        if (subdivNonemptyInv [n] == 0)
                                continue;
                        setNumbers1. push_back (n);
                }
                int nSets1 (setNumbers1. size ());
 
                sbug << "Computing connecting orbits between " << nSets1 <<
                        " sets... ";
 
                // prepare tables with the corresponding components
                multitable<int_t> compVertices1 (8192);
                multitable<int_t> compEnds1 (8192);
                int cur (0);
                for (int i = 0; i < nSets1; ++ i)
                {
                        int n (setNumbers1 [i]);
                        for (int j = n ? compEnds [n - 1] : 0;
                                j < compEnds [n]; ++ j)
                        {
                                compVertices1 [cur ++] = compVertices [j];
                        }
                        compEnds1 [i] = cur;
                }
 
                // create a graph in which possibly nontrivial Morse sets
                // are collapsed to single vertices of the graph
                int_t *newNumbers = computeconn ? new int_t [X. size ()] : 0;
                diGraph<> collapsed;
                chomp::homology::collapseVertices (g, nSets1, compVertices1,
                        compEnds1, collapsed, newNumbers);
 
                // prepare a translation table for vertex numbers
                int_t *oldNumbers = computeconn ? new int_t [X. size ()] : 0;
                if (oldNumbers)
                {
                        for (int_t i = 0; i < X. size (); ++ i)
                                oldNumbers [newNumbers [i]] = i;
                }
                if (newNumbers)
                        delete [] newNumbers;
 
                // compute the connecting orbits and the connection graph
                diGraph<> connGraph;
                if (computeconn)
                {
                        TConnMorse conn (morseDec, setNumbers1,
                                X, oldNumbers);
                        chomp::homology::inclusionGraph (collapsed, nSets1,
                                connGraph, conn);
                }
                else
                {
                        chomp::homology::inclusionGraph (collapsed, nSets1,
                                connGraph);
                }
                if (oldNumbers)
                        delete [] oldNumbers;
 
                // store the connections in the Morse decomposition structure
                for (int v = 0; v < nSets1; ++ v)
                {
                        int_t nEdges = connGraph. countEdges (v);
                        for (int_t e = 0; e < nEdges; ++ e)
                        {
                                int_t w = connGraph. getEdge (v, e);
                                morseDec. addconn (setNumbers1 [v],
                                        setNumbers1 [w]);
                        }
                }
 
                // show the time used for computing connecting orbits
                sbug << timeConnOrb << ".\n";
        }
#endif
 
        // compute the Conley indices of the Morse sets
        spcCubes wrongCubes;
        spcCubes attrCubes;
        spcCubes skipCubes;
        for (int n = 0; n < nSets; ++ n)
        {
                if (cachedAvailable)
                {
                        if (cached. wrongIndex [n])
                                wrongCubes. add (morseDec [n] [0]);
                        if (cached. attractor [n])
                                attrCubes. add (morseDec [n] [0]);
                        if (cached. skipped [n])
                                skipCubes. add (morseDec [n] [0]);
                        continue;
                }
 
                // skip the computation of the Conley index
                // if it is known that its invariant part is empty
                if (subdivNonemptyInv [n] == 0)
                {
                        setDummyIndex (morseDec, n);
                        continue;
                }
 
                // skip the computation of the Conley index if requested to
                if ((skipIndices > 0) &&
                        (morseDec [n]. size () >= skipIndices))
                {
                        sbug << "* Skipping ConIndex (" << n << ").\n";
                        setDummyIndex (morseDec, n);
                //???   subdivNonemptyInv [n] = 1;
                        cached. skipped [n] = true;
                        skipCubes. add (morseDec [n] [0]);
                }
                // or compute the Conley index
                else
                {
                        sbug << "* Computing Conley index of Morse Set " <<
                                n << " of " << nSets << " (" <<
                                morseDec [n]. size () << " cubes)...\n";
                        int noIsolation = -1;
                        int attractor = -1;
                        int result = computeConleyIndex (morseDec, n,
                                offset, width, subdivDepth, theMap,
                                theCubMap0, theCubMap1,
                                subdivNonemptyInv [n], noIsolation,
                                attractor);
                        if (attractor > 0)
                        {
                                cached. attractor [n] = true;
                                attrCubes. add (morseDec [n] [0]);
                        }
                        if (((result < 0) || (noIsolation > 0)) &&
                                setConleyIndex (morseDec, n, subdivDepth,
                                leftMapParam, rightMapParam, paramCount))
                        {
                                sbug << "* The index was set to an apriori "
                                        "known one.\n";
                        }
                        else if (result < 0)
                        {
                                sbug << "* Failed to compute the index.\n";
                                setDummyIndex (morseDec, n);
                                cached. wrongIndex [n] = true;
                                cached. skipped [n] = true;
                                skipCubes. add (morseDec [n] [0]);
                        }
                        else if (noIsolation > 0)
                        {
                                sbug << "* No isolation detected. The index "
                                        "could not be computed reliably.\n";
                                cached. wrongIndex [n] = true;
                                wrongCubes. add (morseDec [n] [0]);
                        }
                        else
                        {
                                sbug << "* The computed index is " <<
                                        (morseDec. trivial (n) ?
                                        "" : "non") << "trivial.\n";
                        }
                }
        }
        sbug << "Conley indices determined in " <<
                timeConleyIndices << ".\n";
 
        // make a note of the Morse sets whose invariant part is empty
        timeused timePassThru;
        timePassThru = 0;
        hashedset<int> passThru;
        for (int n = 0; n < nSets; ++ n)
        {
                // if the cached result is available then use it
                if (cachedAvailable)
                {
                        if (cached. emptyInv [n])
                                passThru. push_back (n);
                        continue;
                }
 
                // if refinements were already tried and failed then skip it
                if (subdivNonemptyInv [n] == 1)
                        continue;
 
                // if refinements have already proved that inv is empty
                // then make a note of it
                if (subdivNonemptyInv [n] == 0)
                {
                        passThru. push_back (n);
                        cached. emptyInv [n] = true;
                        continue;
                }
 
                // if the Conley index is nontrivial then inv is nonempty
                if (morseDec. computed (n) && !morseDec. trivial (n))
                        continue;
 
                // if the size of the set is too large then skip it
                if (morseDec [n]. size () > maxRefineSize0)
                        continue;
 
                // refine the set and check if its inv is nonempty
                sbug << "Refining the set no. " << n << " (" <<
                        morseDec [n]. size () << " cubes)...\n";
                if (checkEmptyInv (morseDec [n], finalDepth,
                        offset, width, theMap, theCubMap0, theCubMap1))
                {
                        sbug << "Empty inv. The set will be removed.\n";
                        passThru. push_back (n);
                        cached. emptyInv [n] = true;
                }
                else
                {
                        sbug << "Invariant part may be nonempty.\n";
                }
        }
        sbug << "Morse sets with trivial indices analysed in " <<
                timePassThru << ".\n";
 
        // save the Conley indices and the information on empty inv
        // to a file if they were computed and a file name prefix is provided
        if (cacheIndices && !cachedAvailable && !fileExists (indFileName) &&
                !fileExists (lockFileName))
        {
                // create the lock file
                std::ofstream lockFile (lockFileName);
                lockFile. close ();
 
                // write the cached index information
                cached. write (indFileName, morseDec);
 
                // remove the lock file
                std::remove (lockFileName);
        }
 
        // remove Morse sets whose invariant part is proven to be empty
        if (!passThru. empty ())
        {
                sbug << "Removing " << passThru. size () <<
                        " spurious Morse sets... ";
                timeused timeRemoving;
                timeRemoving = 0;
#ifdef CONFIG_NOCONNECTIONS
                if ((passThru. size () > 20) &&
                        (passThru. size () > 0.8 * morseDec. count ()))
                {
                        sbug << "[Q] ";
 
                        // determine the sizes of the Morse decompositions
                        int_t oldSize = morseDec. count ();
                        int_t newSize = morseDec. count () -
                                passThru. size ();
 
                        // create a new Morse decomposition
                        typeMorseDec newMorseDec (theCubMap0);
                        newMorseDec. setnumber (newSize);
 
                        // swap in the data from the old Morse decomposition
                        int_t newNumber = 0;
                        for (int_t oldNumber = 0; oldNumber < oldSize;
                                ++ oldNumber)
                        {
                                // if this set is to be skipped
                                if (passThru. check (oldNumber))
                                        continue;
 
                                // if this set is to remain
                                newMorseDec [newNumber]. swap
                                        (morseDec [oldNumber]);
                                if (morseDec. computed (oldNumber))
                                {
                                        newMorseDec. setindex (newNumber,
                                                morseDec. index (oldNumber));
                                }
                                ++ newNumber;
                        }
                        if (newNumber != newSize)
                                throw "A bug in the pass-through loop.";
 
                        // swap the Morse decompositions
                        morseDec. swap (newMorseDec);
                }
                else
                {
                        for (int i = passThru. size () - 1; i >= 0; -- i)
                        {
                                morseDec. passthru (passThru [i]);
                        }
                }
#else
                for (int i = passThru. size () - 1; i >= 0; -- i)
                {
                        morseDec. passthru (passThru [i]);
                }
#endif
                sbug << timeRemoving << ".\n";
        }
 
        // manipulate the Morse decomposition to make it a coarser one
        if (maxJoinSize > 0)
        {
                sbug << "* " << morseDec. count () << " Morse sets; sizes: ";
                for (int i = 0; i < morseDec. count (); ++ i)
                        sbug << (i ? " " : "") << morseDec [i]. size ();
                sbug << "\n! Joining trivial Morse sets (size: " <<
                        maxJoinSize << ", conn: " << maxJoinConnection <<
                        ", dist: " << maxJoinDistance << ")...\n";
                morseDec. jointrivial (maxJoinSize, maxJoinConnection,
                        maxJoinDistance);
        }
 
        // remember the new number of Morse sets in the Morse decomposition
        nSets = morseDec. count ();
 
        // determine the numbers of new Morse sets whose indices could not
        // be computed
        for (int_t i = 0; i < wrongCubes. size (); ++ i)
        {
                int n = 0;
                while ((n < nSets) && !morseDec [n]. check (wrongCubes [i]))
                        ++ n;
                if (n < nSets)
                        wrongIndices. push_back (n);
                else
                {
                        wrongIndices. push_back (-1);
                        sbug << "Note: A Morse set with a wrong index "
                                "could not be identified.\n";
                }
        }
 
        // determine the numbers of new Morse sets whose indices were skipped
        for (int_t i = 0; i < skipCubes. size (); ++ i)
        {
                int n = 0;
                while ((n < nSets) && !morseDec [n]. check (skipCubes [i]))
                        ++ n;
                if (n < nSets)
                        skippedIndices. push_back (n);
                else
                {
                        skippedIndices. push_back (-1);
                        sbug << "Note: A Morse set whose index was skipped "
                                "could not be identified.\n";
                }
        }
 
        // determine the numbers of new Morse sets which are attractors
        for (int_t i = 0; i < attrCubes. size (); ++ i)
        {
                int n = 0;
                while ((n < nSets) && !morseDec [n]. check (attrCubes [i]))
                        ++ n;
                if (n < nSets)
                        attractors. push_back (n);
                else
                {
                        attractors. push_back (-1);
                        sbug << "Note: A Morse set which is an attractor "
                                "could not be identified.\n";
                }
        }
 
        // verify if all the new Morse sets are isolated
        // from the outside of the box by at least one cube layer
        sbug << "Checking isolation of the final Morse sets... ";
        bool finalIsolated = true;
        for (int n = 0; n < nSets; ++ n)
        {
                if (checkIsolation (morseDec [n], subdivDepth))
                        continue;
                finalIsolated = false;
                break;
        }
        sbug << (finalIsolated ? "[isolated]\n" : "[no isol]\n");
 
        // show a summary information
        sbug << "* " << morseDec. count () << " Morse sets [" <<
                attrCubes. size () << " attractor(s)]; sizes: ";
        for (int i = 0; i < morseDec. count (); ++ i)
                sbug << (i ? " " : "") << morseDec [i]. size ();
        sbug << ".\n";
 
        // compute the min and max coordinates of Morse sets
        CoordMinMax coordMinMax;
        int_t nMorseSets = morseDec. count ();
        for (int_t nSet = 0; nSet < nMorseSets; ++ nSet)
                coordMinMax (morseDec [nSet]);
        if (coordMinMax. defined ())
        {
                sbug << "The computed Morse sets are in " <<
                        coordMinMax << ".\n";
                sbug << "Real coordinates: ";
                showRealCoords (sbug, coordMinMax) << ".\n";
        }
 
        return;
} /* invMorseDec */
 
 
// --------------------------------------------------
// -------------- Morse Decomposition ---------------
// --------------------------------------------------
 
/// Computes the Morse decomposition using all the pre- and postprocessing.
/// The target Morse decomposition object must be already initialized
/// with the right map which is to be used to compute the Conley indices
/// of the Morse sets.
/// If a file name prefix is nonempty then an attempt is being made
/// to read cached Morse decomposition information from a file.
inline theMorseDecompositionType *computeMorseDecomposition
        (theMapType &theMap,
        const theCubMapType theCubMap0, const theCubMapType theCubMap1,
        const double *offset, const double *width, int_t skipIndices,
        const std::string &cacheFileName,
        const std::string &pictureFileName,
        const std::string &cubesFilePrefix,
        const std::string &morseDecFileName,
        const std::string &graphFileName,
        const std::string &procFilePrefix,
        const std::string &mapOptFileName,
        bool &morseDecComputed,
        std::vector<int> &wrongIndices, std::vector<int> &skippedIndices,
        std::vector<int> &attractors, bool connOrbits,
        const double *leftMapParam, const double *rightMapParam,
        int paramCount)
{
        using chomp::homology::sbug;
        using chomp::homology::sout;
        using chomp::homology::diGraph;
        using chomp::homology::timeused;
 
        // start measuring time of the subdivisions
        timeused timeSubdiv;
        timeSubdiv = 0;
 
        // check for the presence of local map optimization hints
        // and use them if they are available
        bool mapOptRead = false;
        if (!mapOptFileName. empty () &&
                fileExists (mapOptFileName. c_str ()))
        {
                std::ifstream in (mapOptFileName. c_str ());
                std::string str;
                in >> str;
                chomp::homology::ignorecomments (in);
                if ((str. size () > 0) && (in. peek () == 'D'))
                {
                        chomp::multiwork::mwData data;
                        text2data (str, data);
                        theMap. loadInternals (data);
                        mapOptRead = true;
                        sbug << "[using previous map optimization data]\n";
                }
        }
 
        // create an empty Morse decomposition
        theMorseDecompositionType *morseDecPtr =
                new theMorseDecompositionType (theCubMap0);
        theMorseDecompositionType &morseDec = *morseDecPtr;
 
        // retrieve the Morse decomposition and the auxiliary data
        // from a temporary binary compressed file if the file is available
        if (!morseDecFileName. empty () &&
                fileExists (morseDecFileName. c_str ()))
        {
                sbug << "Loading Morse dec cache... ";
                using chomp::homology::program_time;
                double start_load_time (program_time);
                loadMorseDecCache (morseDecFileName. c_str (), morseDec,
                        wrongIndices, skippedIndices, attractors);
                double end_load_time (program_time);
                sbug << "(" << static_cast<long> ((end_load_time -
                        start_load_time) * 100) / 100 << " sec.)\n";
        }
 
        // otherwise compute the Morse decomposition
        else
        {
 
        // create a cubical set with one cube
        spcCoord zero [spaceDim];
        for (int i = 0; i < spaceDim; ++ i)
                zero [i] = 0;
        spcCubes X;
        X. add (spcCube (zero, spaceDim));
 
        // TEMPORARY CODE for GLOBCLOG (discussion on 2010-05-06)
        spcCubes restriction;
        bool restrictionRead = false;
        {
                std::ifstream restr ("restriction.cub");
                if (restr)
                {
                        sout << "\n*** Restricting to a subset: ";
                        restr >> restriction;
                        sout << restriction. size () << " cubes. ***\n";
                        restrictionRead = true;
                }
        }
 
        // make the subdivision(s)
        for (int subdivDepth = 1; subdivDepth <= finalDepth; ++ subdivDepth)
        {
                // use the restriction cubes for X
                // or subdivide all the cubes in X
                {
                        if (subdivDepth >= initialDepth)
                                sbug << "Depth " << subdivDepth << ": ";
                        if (restrictionRead && (subdivDepth <= initialDepth))
                        {
                                if (subdivDepth == initialDepth)
                                        X. swap (restriction);
                        }
                        else
                        {
                                spcCubes Y;
                                subdivideCubes (X, Y);
                                X. swap (Y);
                        }
                        if (subdivDepth >= initialDepth)
                                sbug << X. size () << ". ";
                }
 
                // if the initial subdivision depth has not been reached yet
                // then skip the computation of the invariant part
                if (subdivDepth < initialDepth)
                        continue;
 
                // shuffle the set of cubes so that the ODE integration
                // sapmles the entire phase space more or less uniformly
                moveSomeToFront (X, 1 << spaceDim);
 
                // prepare a local multivalued cubical map
                // at this subdivision depth
                theCubMapType localCubMap (offset, width,
                        1 << subdivDepth, subdivDepth, theMap);
 
                // extract the multivalued cubical map
                // from the Morse decomposition if this is the final depth
                const theCubMapType &theCubMap =
                        (subdivDepth == finalDepth) ?
                        theCubMap0 : localCubMap;
 
                // set the default limitations on image size
                theCubMap. setAllowedImgSize (maxImageDiameter,
                        maxImageVolume);
 
                // measure time of the construction of the graph of the map
                timeused timeMap;
                timeMap = 0;
 
#ifdef CONFIG_ODEVART
                // prepare a set of cubes to be added to X
                // to make a necessary correction for isolation
                spcCubes Xadd;
#endif
                // compute the graph of the multivalued cubical map on X
                sbug << "F... ";
                diGraph<> g;
                while (1)
                {
                        // try the computation with the current map
                        try
                        {
#ifdef CONFIG_ODEVART
                                // compute the graph of the map and compute
                                // a set that should be added to X
                                // to preserve the isolation condition
                                {
                                        spcCubes empty;
                                        Xadd. swap (empty);
                                }
                                int_t maxAdd (std::max
                                        (static_cast<int_t> (10000),
                                        X. size () / 100));
                                computeMapGraphIsol (X, g, theCubMap, Xadd,
                                        maxAdd);
#else
                                // compute the graph of the map
                                computeMapGraph (X, g, theCubMap);
#endif
 
                                // show the avg size of the image of a cube
                                int_t avgImgSize = g. countVertices () ?
                                        (g. countEdges () /
                                        g. countVertices ()) : 0;
                                sbug << "[avg " << avgImgSize << "] ";
 
                                // show the time of computation of the map
                                sbug << "(" << double2string
                                        (static_cast<double> (timeMap), 1) <<
                                        " sec) ";
 
                                // if there are no edges in the graph then OK
                                if (!g. countEdges ())
                                        break;
 
                                // if it is not necessary to repeat then OK
                                if (!theMap. adjust (true, subdivDepth))
                                        break;
                        }
                        // in case of failure, try adjusting the map;
                        // if it is not possible to repeat then quit
                        catch (const char *msg)
                        {
                                resetRounding ();
                                sbug << "Failed (" << double2string
                                        (static_cast<double> (timeMap), 1) <<
                                        " sec).\nFailure reason: " <<
                                        msg << ".\n";
                                if (!theMap. adjust (false, subdivDepth))
                                        throw;
                        }
                        catch (const std::exception &e)
                        {
                                resetRounding ();
                                sbug << "Failed (" << double2string
                                        (static_cast<double> (timeMap), 1) <<
                                        " sec).\nFailure reason: " <<
                                        e. what () << "\n";
                                if (!theMap. adjust (false, subdivDepth))
                                        throw;
                        }
 
                        // forget the map cache and the graph
                        // before trying again
                        theCubMap. cleanCache ();
                        diGraph<> dummy;
                        g. swap (dummy);
                        timeMap. reset ();
                }
 
#ifdef CONFIG_ODEVART
                // add the computed correction to the set of cubes X
                X. add (Xadd);
#endif
 
                // compute the requested Morse decomposition
                // if this is the final subdivision level
                if (subdivDepth == finalDepth)
                {
                        // save the local map optimization data to a file
                        if (!mapOptRead && !mapOptFileName. empty () &&
                                !fileExists (mapOptFileName. c_str ()))
                        {
                                chomp::multiwork::mwData data;
                                theMap. saveInternals (data);
                                std::string str (data2text (data));
                                std::ofstream out (mapOptFileName. c_str ());
                                out << str << " D\n";
                                sbug << "[map optimization data saved]\n";
                        }
 
                        // set new limitations on image size to make sure
                        // that the computations are not interrupted
                        theCubMap0. setAllowedImgSize
                                (maxImageDiameter << 3, maxImageVolume << 3);
                        theCubMap1. setAllowedImgSize
                                (maxImageDiameter << 3, maxImageVolume << 3);
 
                        // compute the Morse decomposition
                        sbug << "Morse dec... ";
                        invMorseDec (X, g, morseDec,
                                offset, width, subdivDepth,
                                theMap, theCubMap0, theCubMap1,
                                skipIndices, cacheFileName, morseDecComputed,
                                wrongIndices, skippedIndices, attractors,
                                connOrbits, leftMapParam, rightMapParam,
                                paramCount, timeSubdiv);
 
                        // turn on space wrapping if necessary
                        localSpaceWrapping wrapping (1 << subdivDepth);
 
                        // run any extra analysis of the Morse decomposition,
                        // unless already done before
                        bool alreadyProcessed = false;
                        if (!procFilePrefix. empty ())
                        {
                                std::string procMarkerName =
                                        procFilePrefix + ".mark";
                                if (fileExists (procMarkerName. c_str ()))
                                        alreadyProcessed = true;
                                else
                                {
                                        std::ofstream markerFile
                                                (procMarkerName. c_str ());
                                        markerFile. close ();
                                }
                        }
                        if (!alreadyProcessed)
                        {
                                processMorseDec (morseDec, X, g,
                                        theCubMap0, theCubMap1,
                                        cubesFilePrefix, procFilePrefix);
                        }
                        break;
                }
 
                // compute the invariant part of the map on X
#ifdef CONFIG_NOCONNECTIONS
                sbug << "ChRec... ";
                chainRecurrentSet (X, X, g);
#else
                sbug << "Inv... ";
                invariantPart (X, g);
#endif
                sbug << X. size () << " ";
#ifdef CONFIG_ENHANCEINV
                sbug << "E" << CONFIG_ENHANCEINV << ": ";
                enhanceCubes (X, CONFIG_ENHANCEINV, subdivDepth);
                sbug << X. size () << " ";
#endif
                sbug << (checkIsolation (X, subdivDepth) ?
                        "[isolated]\n" : "[no isol]\n");
 
                // save the result to a file (for pictures for the paper)
                if (false)
                {
                        std::ostringstream fstr;
                        fstr << "_inv" << subdivDepth;
                        std::string filename = fstr. str ();
                        std::ofstream f (filename. c_str ());
                        f << "; Inv(X) at depth " << subdivDepth << ".\n";
                        f << X << "\n";
                        sbug << "Inv(X) saved to " << filename << ".\n";
                }
        }
 
        // end of the computation of the Morse decomposition
        }
 
        // save the computed Morse sets to separate files
        if (!cubesFilePrefix. empty ())
        {
                sout << "Saving the Morse decomposition to '" <<
                        cubesFilePrefix << "*'...\n";
                morseDec. savetofiles (cubesFilePrefix. c_str ());
        }
 
        // save the computed Conley-Morse graph for the "dot" program
        if (!graphFileName. empty () &&
                !fileExists (graphFileName. c_str ()))
        {
                // create the Morse graph from the Morse decomposition
                chomp::homology::diGraph<> connGraph;
                morseDec. makegraph (connGraph);
                chomp::homology::diGraph<> morseGraph;
                transitiveReduction (connGraph, morseGraph);
 
                // gather the sizes of each Morse set into a vector
                int nSets = morseGraph. countVertices ();
                std::vector<int_t> sizes (nSets);
                for (int n = 0; n < nSets; ++ n)
                        sizes [n] = morseDec [n]. size ();
 
                // gather the Conley indices and the eigenvalues
                std::vector<theConleyIndexType> indices (nSets);
                std::vector<IndexEigenValues> eigenValues (nSets);
                for (int n = 0; n < nSets; ++ n)
                {
                        indices [n] = morseDec. index (n);
                        eigenValues [n] = IndexEigenValues (indices [n]);
                }
 
                sout << "Saving the Conley-Morse graph to '" <<
                        graphFileName << "'...\n";
                std::ofstream graphFile (graphFileName. c_str ());
                writeDotGraph (graphFile, morseGraph, sizes,
                        indices, eigenValues,
                        wrongIndices, skippedIndices, attractors);
                graphFile << "\n";
                graphFile. close ();
        }
 
        // save the Morse decomposition and the auxiliary data
        // to a temporary binary compressed file if requested to
        if (!morseDecFileName. empty () &&
                !fileExists (morseDecFileName. c_str ()))
        {
                sbug << "Saving Morse dec cache... ";
                using chomp::homology::program_time;
                double start_save_time (program_time);
                saveMorseDecCache (morseDecFileName. c_str (), morseDec,
                        wrongIndices, skippedIndices, attractors);
                double end_save_time (program_time);
                sbug << "(" << static_cast<long> ((end_save_time -
                        start_save_time) * 100) / 100 << " sec.)\n";
 
                // DEBUG: verify if the Morse decomposition can be retrieved
                // correctly from the saved file, to make sure that this works
                if (false)
                {
                        double start_verif_time (program_time);
                        sbug << "DEBUG: Verifying Morse Dec cache... ";
                        theMorseDecompositionType morseDec1 (theCubMap0);
                        std::vector<int> wrongIndices1;
                        std::vector<int> skippedIndices1;
                        std::vector<int> attractors1;
                        loadMorseDecCache (morseDecFileName. c_str (),
                                morseDec1,
                                wrongIndices1, skippedIndices1, attractors1);
                        if (morseDec. count () != morseDec1. count ())
                                throw "Wrong number of sets.";
                        for (int n = 0; n < morseDec. count (); ++ n)
                        {
                                if (!(morseDec [n] == morseDec1 [n]))
                                        throw "Corrupted sets.";
                        }
                        if (!(wrongIndices == wrongIndices1))
                                throw "Corrupted wrongIndices.";
                        if (!(skippedIndices == skippedIndices1))
                                throw "Corrupted skippedIndices.";
                        if (!(attractors == attractors1))
                                throw "Corrupted attractors.";
                        double end_verif_time (program_time);
                        sbug << "(" << static_cast<long> ((end_verif_time -
                                start_verif_time) * 100) / 100 <<
                                " sec.) All OK.\n";
                }
        }
 
        return morseDecPtr;
} /* computeMorseDecomposition */
 
 
#endif // _CMGRAPHS_COMPMDEC_H_