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-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 12, 2007. Last revision: April 14, 2008.


#ifndef _CMGRAPHS_COMPMDEC_H_
#define _CMGRAPHS_COMPMDEC_H_


// include some standard C++ header files
#include <algorithm>
#include <new>
#include <iostream>
#include <iomanip>
#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/system/timeused.h"

// include local header files
#include "config.h"
#include "typedefs.h"
#include "conindex.h"
#include "morsedec.h"
#include "indcache.h"
#include "invpart.h"
#include "utils.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 _n, chomp::homology::multitable
                <chomp::homology::multitable<int> > &_connections,
                chomp::homology::multitable<int> &_conncount,
                const int *_numTranslate);

        /// Adds an element v at the connection from source to target.
        void operator () (int source, int target, int v);

private:
        /// The number of vertices between the connections are recorded.
        int 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> >
                &connections;

        /// The numbers of elements in each connecting orbit.
        chomp::homology::multitable<int> &conncount;

        /// The translation table for vertex numbers in the connections.
        const int *numTranslate;

}; /* class TConnTable */

inline TConnTable::TConnTable (int _n,
        chomp::homology::multitable<chomp::homology::multitable<int> > &
        _connections, chomp::homology::multitable<int> &_conncount,
        const int *_numTranslate):
        n (_n), connections (_connections), conncount (_conncount),
        numTranslate (_numTranslate)
{
        int maxcount = n * n;
        for (int i = 0; i < maxcount; ++ i)
                conncount [i] = 0;
        return;
} /* TConnTable::TConnTable */

inline void TConnTable::operator () (int source, int target, int v)
{
        int 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 spcCubes &_X, const int *_numTranslate);

        /// Adds an element v at the connection from source to target.
        void operator () (int source, int target, int v);

private:
        /// The Morse decomposition in which to record the connections.
        theMorseDecompositionType &morseDec;

        /// The set which contains cubes to add.
        const spcCubes &X;

        /// The translation table for vertex numbers in the connections.
        const int *numTranslate;

}; /* class TConnMorse */

inline TConnMorse::TConnMorse (theMorseDecompositionType &_morseDec,
        const spcCubes &_X, const int *_numTranslate):
        morseDec (_morseDec), X (_X), numTranslate (_numTranslate)
{
        return;
} /* TConnMorse::TConnMorse */

inline void TConnMorse::operator () (int source, int target, int v)
{
        morseDec. addconn (source, 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), theMap);
                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)
                                {
                                        subdivNonemptyInv = 1;
                                        return -1;
                                }
                                if (X. empty ())
                                {
                                        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
                {
                        // 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 = true;
                        theCubMap. cropped = false;
                        P. compute ();
                        theCubMap. cropping = crop;

                        // if no cropping took place then compute the index
                        if (ignoreIsolationForConleyIndex ||
                                !theCubMap. cropped)
                        {
                                // 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? maxRefineSize0 : maxRefineSize1;

                // if the set is too large to subdivide then cancel
                if (theSet. size () > maxRefineSize)
                {
                        sbug << "Failure: The set size " <<
                                theSet. size () << " is beyond the max "
                                "allowed " << maxRefineSize << ".\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 */


// --------------------------------------------------
// ----------- 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 areguments "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 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)
{
        using chomp::homology::sbug;
        using chomp::homology::diGraph;
        using chomp::homology::multitable;

        // do nothing if the set X is empty
        if (X. empty ())
                return;

        // determine if it is necessary to compute the connecting orbits
        // in terms of cubical sets
        bool computeconn = connOrbits || (maxJoinSize > 0);

        // compute the components of the chain recurrent set of the graph
        multitable<int> compVertices, compEnds;
        diGraph<> *sccAddr = 0;
        diGraph<> *transpAddr = /*(savebasmin || savebasmax) ? &gT :*/ 0;
        int nSets = chomp::homology::SCC (g, compVertices, compEnds,
                sccAddr, transpAddr);

        // create the Morse decomposition data structure
        // and count the number of cubes in the SCC
        morseDec. setnumber (nSets);
        int countSCCcubes = 0;
        for (int n = 0; n < nSets; ++ n)
        {
                // determine the vertex range of this set
                int beginVertex = n ? compEnds [n - 1] : 0;
                int endVertex = compEnds [n];
                countSCCcubes += endVertex - beginVertex;

                // add the cubes to the Morse set
                for (int i = beginVertex; i < endVertex; ++ i)
                        morseDec. add (n, X [compVertices [i]]);
        }

        // compute the connecting orbits for the Morse decomposition
        {
                // collapse Morse sets to single vertices of a graph
                int *newNumbers = computeconn ? new int [X. size ()] : 0;
                diGraph<> collapsed;
                chomp::homology::collapseVertices (g, nSets, compVertices,
                        compEnds, collapsed, newNumbers);

                // prepare a translation table for vertex numbers
                int *oldNumbers = computeconn ?
                        new int [X. size () - countSCCcubes + nSets] : 0;
                if (oldNumbers)
                {
                        for (int i = 0; i < X. size (); ++ i)
                                oldNumbers [newNumbers [i]] = i;
                }
                if (newNumbers)
                        delete [] newNumbers;

                // compute the connecting orbits and the connection graph
                diGraph<> connGraph;
                TConnMorse conn (morseDec, X, oldNumbers);
                if (computeconn)
                {
                        chomp::homology::inclusionGraph (collapsed, nSets,
                                connGraph, conn);
                }
                else
                {
                        chomp::homology::inclusionGraph (collapsed, nSets,
                                connGraph);
                }
                if (oldNumbers)
                        delete [] oldNumbers;

                // store the connections in the Morse decomposition structure
                for (int v = 0; v < nSets; ++ v)
                {
                        int nEdges = connGraph. countEdges (v);
                        for (int e = 0; e < nEdges; ++ e)
                        {
                                int w = connGraph. getEdge (v, e);
                                morseDec. addconn (v, w);
                        }
                }
        }

        // 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 ();

        // retrieve the Conley indices from cache if possible
        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;

        // compute the Conley indices of the Morse sets
        spcCubes wrongCubes;
        spcCubes attrCubes;
        spcCubes skipCubes;
        std::vector<int> subdivNonemptyInv (nSets, -1);
        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 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 the Conley index "
                                "of the Morse set no. " << n << " (" <<
                                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;
                                wrongCubes. 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";
                        }
                }
        }

        // make a note of the Morse sets whose invariant part is empty
        std::vector<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. 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";
                }
        }

        // 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
        for (int i = passThru. size () - 1; i >= 0; -- i)
        {
                morseDec. passthru (passThru [i]);
        }

        // 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 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 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 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";
                }
        }

        // 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";

        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
        (const theMapType &theMap,
        const theCubMapType theCubMap0, const theCubMapType theCubMap1,
        const double *offset, const double *width, int 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)
{
        using chomp::homology::sbug;
        using chomp::homology::diGraph;
        using chomp::homology::multitable;

        // create an empty Morse decomposition
        theMorseDecompositionType *morseDecPtr =
                new theMorseDecompositionType (theCubMap0);
        theMorseDecompositionType &morseDec = *morseDecPtr;

        // 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));

        // make the initial subdivision(s)
        for (int subdivDepth = 1; subdivDepth <= finalDepth; ++ subdivDepth)
        {
                // subdivide all the cubes in X
                {
                        spcCubes Y;
                        subdivideCubes (X, Y);
                        X. swap (Y);
                }

                // if the initial subdivision depth has not been reached yet
                // then skip the computation of the invariant part
                if (subdivDepth < initialDepth)
                        continue;

                // indicate the number of cubes after the subdivision
                sbug << "Depth " << subdivDepth << ": " <<
                        X. size () << ". ";

                // prepare a local multivalued cubical map
                // at this subdivision depth
                theCubMapType localCubMap (offset, width,
                        1 << subdivDepth, theMap);

                // extract the multivalued cubical map
                // from the Morse decomposition if this is the final depth
                const theCubMapType &theCubMap =
                        (subdivDepth == finalDepth) ?
                        theCubMap0 : localCubMap;

                // compute the graph of the multivalued cubical map on X
                sbug << "F... ";
                diGraph<> g;
                computeMapGraph (X, g, theCubMap);

                // show the average size of the image of a cube
                int avgImgSize = g. countEdges () /
                        (X. size () ? X. size () : 1);
                sbug << "[avg " << avgImgSize << "] ";

                // compute the requested Morse decomposition
                // if this is the final subdivision level
                if (subdivDepth == finalDepth)
                {
                        sbug << "Morse dec...\n";
                        invMorseDec (X, g, morseDec,
                                offset, width, finalDepth,
                                theMap, theCubMap0, theCubMap1,
                                skipIndices, cacheFileName, morseDecComputed,
                                wrongIndices, skippedIndices, attractors,
                                connOrbits, leftMapParam, rightMapParam,
                                paramCount);
                        break;
                }

                // compute the invariant part of the map on X
                sbug << "Inv... ";
                invariantPart (X, g);
                sbug << X. size () << " left.\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 level " << subdivDepth << ".\n";
                        f << X << "\n";
                        sbug << "Inv(X) saved to " << filename << ".\n";
                }
        }

        return morseDecPtr;
} /* computeMorseDecomposition */


#endif // _CMGRAPHS_COMPMDEC_H_

---