morsedec.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file morsedec.h
///
/// Morse decompositions.
/// This file contains the definition of a Morse decomposition class
/// for a discrete dynamical system represented by a combinatorial
/// cubical multivalued map. It also contains the interface to CHomP
/// for the purpose of computing the Conley index of each Morse set.
///
/// @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 January 3, 2007. Last revision: May 18, 2013.
 
 
#ifndef _CMGRAPHS_MORSEDEC_H_
#define _CMGRAPHS_MORSEDEC_H_
 
 
// include some standard C++ header files
#include <iostream>
#include <sstream>
#include <algorithm>
#include <cmath>
 
// include selected header files from the CHomP library
#include "chomp/system/textfile.h"
#include "chomp/homology/homology.h"
#include "chomp/struct/digraph.h"
 
// include local header files
#include "conindex.h"
 
 
// --------------------------------------------------
// --------------- MorseDecomposition ---------------
// --------------------------------------------------
 
/// The Morse decoposition class. This is a template whose parameters
/// define a class used to compute the combinatorial cubical multivalued map
/// on Morse sets, a class which corresponds to a single cube,
/// and a class which is used to keep sets of cubes
/// (the hashed set by default).
template <class mapcomp, class cubetype,
        class cubsettype = chomp::homology::hashedset<cubetype> >
class MorseDecomposition
{
        /// The type of a combinatorial cubical multivalued map.
        typedef chomp::homology::mvmap<cubetype,cubetype> cubmaptype;
 
public:
        /// The type of the map computation class.
        typedef mapcomp MapCompType;
 
        /// The type of a single cube.
        typedef cubetype CubeType;
 
        /// The type of a set of cubes used to store each Morse set.
        typedef cubsettype CubSetType;
 
        /// The type of the index pair.
        typedef IndexPair<mapcomp,cubetype,cubsettype> PairType;
 
        /// The type of the Conley index.
        typedef ConleyIndex<PairType> IndexType;
 
        /// The default constructor.
        MorseDecomposition (const mapcomp &_M = mapcomp ());
 
        /// The destructor.
        ~MorseDecomposition ();
 
        /// The copy constructor.
        MorseDecomposition (const MorseDecomposition &m);
 
        /// The assignment operator.
        MorseDecomposition<mapcomp,cubetype,cubsettype> &operator =
                (const MorseDecomposition<mapcomp,cubetype,cubsettype> &m);
 
        /// Swaps the internal data with another object (note: both must be
        /// initialized with the same map computation object reference).
        void swap (MorseDecomposition &m);
 
        // --- the combinatorial cubical multivalued map ---
        
        /// The map computation class.
        const mapcomp &M;
 
        // --- Morse Sets ---
        
        /// Changes (increases or decreases) the number of Morse sets.
        int setnumber (int n);
 
        /// Adds a cube to the given Morse set and increases the number
        /// of Morse sets if necessary.
        int add (int n, const cubetype &q);
 
        /// Returns the number of Morse sets.
        int count () const;
 
        /// Returns the n-th Morse set. Throws an exception in case
        /// of a wrong number.
        const cubsettype &operator [] (int n) const;
 
        /// Returns the n-th Morse set. Throws an exception in case
        /// of a wrong number.
        cubsettype &operator [] (int n);
 
        /// Returns the approximate distance between the two Morse sets.
        /// Warning: This computation may take some time for large sets.
        /// Note: The values returned by this function are tabulated, and
        /// are not updated when the Morse sets are modified with the 'add'
        /// function. Call with (-1,-1) to reset the tabulation.
        int distance (int n, int m) const;
 
        // --- Connections ---
 
        /// Adds a connection between the given Morse sets without adding
        /// any cubes to this connection (like in a direct connection).
        int addconn (int n, int m);
 
        /// Adds a cube to the connecting orbit between the given Morse sets.
        int addconn (int n, int m, const cubetype &q);
 
        /// Checks whether there exists a connection between given sets.
        /// Returns true if yes, false if no or wrong numbers.
        bool connected (int n, int m) const;
 
        /// Returns the connection between the given sets.
        /// Throws an exception if there is no connection or wrong numbers.
        const cubsettype &connection (int n, int m) const;
 
        /// Returns the approximate maximal distance of a cube in the orbit
        /// that connects the sets n and m (in either direction) from the
        /// k-th Morse set, i.e., how far the orbit is from the set.
        /// Note: This computation may take some time for large sets.
        int distance (int n, int m, int k) const;
 
        /// Returns the appropriate maximal distance between cubes that are
        /// contained in the Morse set n and the connecting orbit n-m or m-n.
        /// Note: The values returned by this function are tabulated, and
        /// are not updated when the Morse sets are modified with the 'add'
        /// function or when the connecting orbits are modified with the
        /// 'addconn' function. Call with (-1,-1) to reset the tabulation.
        int conndistance (int n, int m) const;
 
        /// Creates a graph representation of the Morse decomposition.
        /// Each Morse set corresponds to a vertex, and each connecting orbit
        /// gives rise to a corresponding edge.
        /// The graph given as an argument must be initially empty.
        int makegraph (chomp::homology::diGraph<> &g) const;
 
        // --- the Conley Indices of the Morse Components ---
 
        /// Computes the Conley index of the given Morse set.
        const IndexType &compute (int n);
 
        /// Computes the Conley indices of all the Morse sets.
        int compute ();
 
        /// Verifies if the Conley index for the given Morse set has been
        /// already computed.
        bool computed (int n) const;
 
        /// Retrieves the previously computed Conley index of the given
        /// Morse set. Throws an exception if the index has not yet
        /// been computed. Note that if the Morse set is chaned, then its
        /// index should be computed again using one of the functions above.
        const IndexType &index (int n) const;
 
        /// Sets the index of the given Morse set to the apriori known one.
        /// Does not compute it, but verifies if it is nontrivial.
        const IndexType &setindex (int n, const typename MorseDecomposition
                <mapcomp,cubetype,cubsettype>::IndexType &ind);
 
        /// Returns true if the Conley index of the given Morse set is
        /// trivial, and false if either nontrivial or not yet computed.
        bool trivial (int n) const;
 
        // --- operations on entire Morse decompositions ---
 
        /// Computes the set-wise intersection of two Morse decompositions.
        /// The Morse decomposition must be initially empty and different
        /// from the objects passed as arguments.
        void intersection
                (const MorseDecomposition<mapcomp,cubetype,cubsettype> &m,
                const MorseDecomposition<mapcomp,cubetype,cubsettype> &n);
 
        // --- joining sets ---
 
        /// Joins the two given Morse sets. Updates the connecting orbits.
        /// Resets the Conley index. The new Morse set replaces the n-th one,
        /// and the numbers of Morse sets beyond the m-th one are shifted
        /// downwards.
        int join (int n, int m);
 
        /// Make connecting orbits pass through the given Morse set,
        /// because its invariant part is in fact empty. The Morse set
        /// is deleted and its contents is included in connecting orbits.
        int passthru (int n);
 
        /// Make connecting orbits pass through the given Morse sets,
        /// because their invariant part is in fact empty. The Morse sets
        /// are deleted and their contents are included in the corresponding
        /// connecting orbits. Important: The array must contain the numbers
        /// of the Morse sets in the ascending order.
        int passthru (chomp::homology::hashedset<int> &setsToRemove);
 
        /// Joins small Morse sets with trivial indices with Morse sets
        /// whose indices are nontrivial. Only sets whose size, as well as
        /// the size of the connecting orbit is below the given limits are
        /// joined, provided their distance is also limited.
        int jointrivial (int_t maxsetsize, int_t maxconnsize,
                int maxdistance);
 
        // --- save to files / load from files ---
 
        /// Save the Morse sets and connecting orbits to files.
        /// Set numbers are attached to the given prefix.
        /// Connections are numbered N-M.
        int savetofiles (const char *prefix = 0);
 
        /// Load the Morse sets and connecting orbits from files.
        int loadfromfiles (const char *prefix = 0);
 
private:
        /// The number of Morse sets.
        int nsets;
 
        /// The Morse sets.
        chomp::homology::multitable<cubsettype *> sets;
 
        /// The Conley indices of the Morse sets.
        chomp::homology::multitable<IndexType *> indices;
        
        /// The table that stores the information on whether the indices
        /// are trivial or not: -1 = unknown, 0 = trivial, 1 = nontrivial.
        chomp::homology::multitable<int> nontrivial;
 
        /// Connecting orbits from one set to another.
        /// There is no connection iff the pointer is zero.
        chomp::homology::multitable<chomp::homology::multitable<cubsettype *> >
                conn;
 
        /// Were any connections defined?
        bool connDefined;
 
        /// Tabulated distances between the Morse sets.
        mutable chomp::homology::multitable<chomp::homology::multitable<int> >
                dist;
 
        /// Were any distances computed?
        mutable bool distDefined;
 
        /// Tabulated values of the function "conndistance".
        mutable chomp::homology::multitable<chomp::homology::multitable<int> >
                distconn;
 
        /// Were any values of conn distance defined?
        mutable bool distconnDefined;
 
        /// Initializes the connection array so that connections are defined
        /// from this point on.
        void setConnDefined (void);
 
        /// Release the memory (used in the destructor and copy constructor).
        /// Optionally, delete only Morse decompositions data for Morse sets
        /// whose numbers are greater or equal the given number.
        void deleteall (int firstnumber, int lastnumber);
 
        /// Copy all the data (used in the copy constructor and operator =).
        void copyall (const MorseDecomposition &m);
 
        /// Add one Morse set and initialize its data.
        void addone ();
 
        /// Verifies if two Morse sets are adjacent to each other.
        /// Note: If both sets are large, then this can be time-consuming.
        bool adjacent (int n, int m, int margin = 0) const;
 
        /// Returns the approximate distance between the two cubical sets.
        /// If minimum is set to true, then computes the minimal distance
        /// between any pair of cubes, one taken from X, the other from Y.
        /// If minimum is set to false, then computes the maximal distance
        /// of any cube in X from the set Y, that is, how far X is from Y
        /// (warning: the latter relation is not symmetric!)
        /// Note: This computation may take some time for large sets.
        static int cubsetdistance (const cubsettype &X, const cubsettype &Y,
                bool minimum);
 
}; /* class MorseDecomposition */
 
// --------------------------------------------------
 
template <class mapcomp, class cubetype, class cubsettype>
inline MorseDecomposition<mapcomp,cubetype,cubsettype>::MorseDecomposition
        (const mapcomp &_M): M (_M), nsets (0),
        sets (256), indices (256), nontrivial (256),
        connDefined (false), distDefined (false), distconnDefined (false)
{
        return;
} /* MorseDecomposition::MorseDecomposition */
 
template <class mapcomp, class cubetype, class cubsettype>
inline void MorseDecomposition<mapcomp,cubetype,cubsettype>::setConnDefined
        (void)
{
        // if connection arrays are already in use then do nothing
        if (connDefined)
                return;
 
        // initialize the connecting orbits' pointers if necessary
        for (int i = 0; i < nsets; ++ i)
        {
                for (int j = 0; j < nsets; ++ j)
                {
                        conn [i] [j] = 0;
                        conn [j] [i] = 0;
                }
        }
 
        // make a note that connection arrays are now in use
        connDefined = true;
 
        return;
} /* MorseDecomposition::setConnDefined */
 
template <class mapcomp, class cubetype, class cubsettype>
inline void MorseDecomposition<mapcomp,cubetype,cubsettype>::deleteall
        (int firstnumber, int lastnumber)
{
        int difference = lastnumber - firstnumber;
        if (difference <= 0)
                return;
 
        for (int i = firstnumber; i < lastnumber; ++ i)
        {
                // delete the sets
                if (sets [i])
                {
                        delete sets [i];
                        sets [i] = 0;
                }
 
                // delete the indices
                if (indices [i])
                {
                        delete indices [i];
                        indices [i] = 0;
                }
        }
 
        for (int i = lastnumber; i < nsets; ++ i)
        {
                // shift the sets
                sets [i - difference] = sets [i];
                sets [i] = 0;
 
                // shift the indices
                indices [i - difference] = indices [i];
                indices [i] = 0;
                
                // shift the nontriviality of Morse sets
                nontrivial [i - difference] = nontrivial [i];
        }
 
        // delete the connecting orbits
        if (connDefined)
        {
                for (int i = 0; i < nsets; ++ i)
                {
                        // delete connections
                        for (int j = firstnumber; j < lastnumber; ++ j)
                        {
                                if (conn [i] [j])
                                {
                                        delete conn [i] [j];
                                        conn [i] [j] = 0;
                                }
                        }
 
                        // shift connections
                        for (int j = lastnumber; j < nsets; ++ j)
                        {
                                conn [i] [j - difference] = conn [i] [j];
                                conn [i] [j] = 0;
                        }
                }
        }
 
        // shift tabulated distances
        if (distDefined)
        {
                for (int i = 0; i < nsets; ++ i)
                {
                        for (int j = lastnumber; j < nsets; ++ j)
                                dist [i] [j - difference] = dist [i] [j];
                }
        }
        if (distconnDefined)
        {
                for (int i = 0; i < nsets; ++ i)
                {
                        for (int j = lastnumber; j < nsets; ++ j)
                        {
                                distconn [i] [j - difference] =
                                        distconn [i] [j];
                        }
                }
        }
 
        if (connDefined)
        {
                for (int j = 0; j < nsets - difference; ++ j)
                {
                        // delete connections
                        for (int i = firstnumber; i < lastnumber; ++ i)
                        {
                                if (conn [i] [j])
                                {
                                        delete conn [i] [j];
                                        conn [i] [j] = 0;
                                }
                        }
 
                        // shift connections
                        for (int i = lastnumber; i < nsets; ++ i)
                        {
                                conn [i - difference] [j] = conn [i] [j];
                                conn [i] [j] = 0;
                        }
                }
        }
 
        // shift tabulated distances
        if (distDefined)
        {
                for (int j = 0; j < nsets - difference; ++ j)
                {
                        for (int i = lastnumber; i < nsets; ++ i)
                                dist [i - difference] [j] = dist [i] [j];
                }
        }
        if (distconnDefined)
        {
                for (int j = 0; j < nsets - difference; ++ j)
                {
                        for (int i = lastnumber; i < nsets; ++ i)
                        {
                                distconn [i - difference] [j] =
                                        distconn [i] [j];
                        }
                }
        }
 
        // update the number of Morse sets
        nsets -= difference;
        return;
} /* MorseDecomposition::deleteall */
 
template <class mapcomp, class cubetype, class cubsettype>
inline MorseDecomposition<mapcomp,cubetype,cubsettype>::~MorseDecomposition
        ()
{
        deleteall (0, nsets);
        return;
} /* MorseDecomposition::~MorseDecomposition */
 
template <class mapcomp, class cubetype, class cubsettype>
inline void MorseDecomposition<mapcomp,cubetype,cubsettype>::copyall
        (const MorseDecomposition<mapcomp,cubetype,cubsettype> &m)
{
        // copy the number of Morse sets
        nsets = m. nsets;
 
        for (int i = 0; i < nsets; ++ i)
        {
                // copy the Morse sets
                if (m. sets [i])
                        sets [i] = new cubsettype (*(m. sets [i]));
                else
                        sets [i] = 0;
 
                // copy the Conley indices
                if (m. indices [i])
                        indices [i] = new IndexType (*(m. indices [i]));
                else
                        indices [i] = 0;
                nontrivial [i] = m. nontrivial [i];
        }
 
        // copy the connecting orbits and tabulated distances
        connDefined = m. connDefined;
        if (m. connDefined)
        {
                for (int i = 0; i < nsets; ++ i)
                {
                        for (int j = 0; j < nsets; ++ j)
                        {
                                if (m. conn [i] [j])
                                {
                                        conn [i] [j] = new cubsettype
                                                (*(m. conn [i] [j]));
                                }
                                else
                                        conn [i] [j] = 0;
                        }
                }
        }
        distDefined = m. distDefined;
        if (m. distDefined)
        {
                for (int i = 0; i < nsets; ++ i)
                {
                        for (int j = 0; j < nsets; ++ j)
                        {
                                dist [i] [j] = m. dist [i] [j];
                        }
                }
        }
        distconnDefined = m. distconnDefined;
        if (m. distconnDefined)
        {
                for (int i = 0; i < nsets; ++ i)
                {
                        for (int j = 0; j < nsets; ++ j)
                        {
                                distconn [i] [j] = m. distconn [i] [j];
                        }
                }
        }
        return;
} /* MorseDecomposition::copyall */
 
template <class mapcomp, class cubetype, class cubsettype>
inline MorseDecomposition<mapcomp,cubetype,cubsettype>::MorseDecomposition
        (const MorseDecomposition<mapcomp,cubetype,cubsettype> &m): M (m. M)
{
        copyall (m);
        return;
} /* MorseDecomposition::MorseDecomposition */
 
template <class mapcomp, class cubetype, class cubsettype>
inline MorseDecomposition<mapcomp,cubetype,cubsettype> &
        MorseDecomposition<mapcomp,cubetype,cubsettype>::operator =
        (const MorseDecomposition<mapcomp,cubetype,cubsettype> &m)
{
        // protection from "M = M;"
        if (this == &m)
                return *this;
 
        // delete the currently allocated Morse sets and connecting orbits
        deleteall (0, nsets);
 
        // copy the Morse sets from the other decomposition
        copyall (m);
 
        return *this;
} /* MorseDecomposition::operator = */
 
template <class mapcomp, class cubetype, class cubsettype>
inline void MorseDecomposition<mapcomp,cubetype,cubsettype>::swap
        (MorseDecomposition<mapcomp,cubetype,cubsettype> &m)
{
        std::swap (this -> nsets, m. nsets);
        this -> sets. swap (m. sets);
        this -> indices. swap (m. indices);
        this -> nontrivial. swap (m. nontrivial);
        this -> conn. swap (m. conn);
        std::swap (this -> connDefined, m. connDefined);
        this -> dist. swap (m. dist);
        std::swap (this -> distDefined, m. distDefined);
        this -> distconn. swap (m. distconn);
        std::swap (this -> distconnDefined, m. distconnDefined);
        return;
} /* MorseDecomposition::swap */
 
template <class mapcomp, class cubetype, class cubsettype>
inline void MorseDecomposition<mapcomp,cubetype,cubsettype>::addone ()
{
        // initialize the new Morse set pointer
        sets [nsets] = new cubsettype;
 
        // initialize the new Conley index pointer
        indices [nsets] = 0;
        nontrivial [nsets] = -1;
 
        // initialize the connecting orbits' pointers if necessary
        if (connDefined)
        {
                for (int i = 0; i < nsets; ++ i)
                {
                        conn [nsets] [i] = 0;
                        conn [i] [nsets] = 0;
                }
                conn [nsets] [nsets] = 0;
        }
 
        // initialize the tabulated distances if necessary
        if (distDefined)
        {
                for (int i = 0; i < nsets; ++ i)
                {
                        dist [nsets] [i] = -1;
                        dist [i] [nsets] = -1;
                }
                dist [nsets] [nsets] = 0;
        }
        if (distconnDefined)
        {
                for (int i = 0; i < nsets; ++ i)
                {
                        distconn [nsets] [i] = -1;
                        distconn [i] [nsets] = -1;
                }
                distconn [nsets] [nsets] = 0;
        }
 
        // increase the stored number of Morse sets
        ++ nsets;
        return;
} /* MorseDecomposition::addone */
 
template <class mapcomp, class cubetype, class cubsettype>
inline bool MorseDecomposition<mapcomp,cubetype,cubsettype>::adjacent
        (int n, int m, int margin) const
{
        // define the right type of neighbors' iterator
        typedef chomp::homology::tNeighbors<typename cubetype::CoordType>
                neighborstype;
 
        // make sure the arguments are correct
        if ((n < 0) || (n >= nsets) || (m < 0) || (m >= nsets))
                return false;
        if (n == m)
                return true;
        if (margin < 0)
                margin = 0;
 
        // determine the smaller and the larger set
        bool n_m = (sets [n] -> size () < sets [m] -> size ());
        const cubsettype &smaller = n_m ? *(sets [n]) : *(sets [m]);
        const cubsettype &larger = n_m ? *(sets [m]) : *(sets [n]);
 
        // if both sets are empty, then the result is obvious
        if (larger. empty ())
                return false;
 
        // prepare an auxiliary buffer for coordinates of a cube
        int dim = larger [0]. dim ();
        typename cubetype::CoordType *c =
                new typename cubetype::CoordType [dim];
 
        // create an enhanced set if the margin is nonzero
        cubsettype enhanced;
        if (margin)
                enhanced = smaller;
        int start = 0;
        for (int i = 0; i < margin; ++ i)
        {
                int stop = enhanced. size ();
                for (int n = start; n < stop; ++ n)
                {
                        enhanced [n]. coord (c);
                        neighborstype nb (c, dim);
                        const typename cubetype::CoordType *nc;
                        while ((nc = nb. get ()) != 0)
                                enhanced. add (cubetype (nc, dim));
                }
                start = stop;
        }
 
        // choose the enhanced set or the original one,
        // depending on the margin
        const cubsettype &smallcheck = margin ? enhanced : smaller;
 
        // for each cube in the smaller set, check if any of its neighbors
        // is contained in the larger set
        for (int_t i = 0; i < smallcheck. size (); ++ i)
        {
                smallcheck [i]. coord (c);
                neighborstype n (c, dim);
                const typename cubetype::CoordType *nc;
                while ((nc = n. get ()) != 0)
                {
                        if (larger. check (cubetype (nc, dim)))
                        {
                                delete [] c;
                                return true;
                        }
                }
        }
        delete [] c;
        return false;
} /* MorseDecomposition::adjacent */
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MorseDecomposition<mapcomp,cubetype,cubsettype>::cubsetdistance
        (const cubsettype &X, const cubsettype &Y, bool minimum)
{
        // if one of the sets is empty, then say that the distance is zero
        if (X. empty () || Y. empty ())
                return 0;
 
        // prepare the arrays for the coordinates of the cubes to compare
        int dim = X [0]. dim ();
        typedef typename cubetype::CoordType coords;
        coords *c = new coords [dim + dim];
        coords *d = c + dim;
 
        // for each pair of cubes (one in the n-th set, another in the
        // m-th set), compute their distance and take the minimum
        int best_dist = -1;
        for (int_t i = 0; i < X. size (); ++ i)
        {
                X [i]. coord (c);
                int local_dist = -1;
                for (int_t j = 0; j < Y. size (); ++ j)
                {
                        Y [j]. coord (d);
                        int dist0 = 0;
                        for (int k = 0; k < dim; ++ k)
                                dist0 += (c [k] - d [k]) * (c [k] - d [k]);
                        dist0 = static_cast<int>
                                (std::sqrt (static_cast<double> (dist0)));
                        if ((local_dist < 0) || (local_dist > dist0))
                                local_dist = dist0;
                }
 
                if (minimum)
                {
                        if ((best_dist < 0) || (best_dist > local_dist))
                                best_dist = local_dist;
                }
                else
                {
                        if ((best_dist < 0) || (best_dist < local_dist))
                                best_dist = local_dist;
                }
        }
 
        delete [] c;
        return best_dist;
} /* MorseDecomposition::distance */
 
// --------------------------------------------------
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MorseDecomposition<mapcomp,cubetype,cubsettype>::setnumber
        (int n)
{
        if (n < 0)
                return 0;
        while (n > nsets)
                addone ();
        if (n < nsets)
                deleteall (n, nsets);
        return 0;
} /* MorseDecomposition::setnumber */
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MorseDecomposition<mapcomp,cubetype,cubsettype>::add
        (int n, const cubetype &q)
{
        if (n < 0)
                return 0;
        while (n > nsets)
                addone ();
        sets [n] -> add (q);
        return 0;
} /* MorseDecomposition::add */
 
// --------------------------------------------------
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MorseDecomposition<mapcomp,cubetype,cubsettype>::addconn
        (int n, int m)
{
        if ((n < 0) || (m < 0))
                return 0;
        setConnDefined ();
        while ((n > nsets) || (m > nsets))
                addone ();
        if (!conn [n] [m])
                conn [n] [m] = new cubsettype;
        return 0;
} /* MorseDecomposition::addconn */
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MorseDecomposition<mapcomp,cubetype,cubsettype>::addconn
        (int n, int m, const cubetype &q)
{
        if ((n < 0) || (m < 0))
                return 0;
        setConnDefined ();
        while ((n > nsets) || (m > nsets))
                addone ();
        if (!conn [n] [m])
                conn [n] [m] = new cubsettype;
        conn [n] [m] -> add (q);
        return 0;
} /* MorseDecomposition::addconn */
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MorseDecomposition<mapcomp,cubetype,cubsettype>::count () const
{
        return nsets;
} /* MorseDecomposition::count */
 
template <class mapcomp, class cubetype, class cubsettype>
inline const cubsettype &MorseDecomposition<mapcomp,cubetype,cubsettype>::
	operator [] (int n) const
{
        if ((n < 0) || (n >= nsets) || !sets [n])
                throw "Trying to access an undefined Morse set.";
        return *(sets [n]);
} /* MorseDecomposition::operator [] */
 
template <class mapcomp, class cubetype, class cubsettype>
inline cubsettype &MorseDecomposition<mapcomp,cubetype,cubsettype>::
	operator [] (int n)
{
        if ((n < 0) || (n >= nsets) || !sets [n])
                throw "Trying to access an undefined Morse set.";
        return *(sets [n]);
} /* MorseDecomposition::operator [] */
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MorseDecomposition<mapcomp,cubetype,cubsettype>::distance
        (int n, int m) const
{
        // make sure that n <= m
        if (n > m)
        {
                int tmp = n;
                n = m;
                m = tmp;
        }
 
        // if the number(s) are too large, return an undefined value
        if (m >= nsets)
                return -1;
 
        // if both numbers are negative, reset all the distances
        if (m < 0)
        {
                distDefined = false;
                return -1;
        }
 
        // if one of the numbers is negative and the array of distances
        // is not in use then no action is necessary
        if ((n < 0) && !distDefined)
        {
                return -1;
        }
 
        // if the sets are the same, their distance is zero
        if (n == m)
                return 0;
 
        // initialize the distance array if it was not done before
        if (!distDefined)
        {
                for (int i = 0; i < nsets; ++ i)
                {
                        for (int j = 0; j < nsets; ++ j)
                        {
                                dist [i] [j] = -1;
                        }
                }
                distDefined = true;
        }
 
        // if exactly one of the numbers is negative,
        // reset the other set of distances only
        if (n < 0)
        {
                for (int i = 0; i < nsets; ++ i)
                {
                        dist [i] [m] = -1;
                        dist [m] [i] = -1;
                }
                return -1;
        }
 
        // if the distance is known then return it
        if (dist [n] [m] >= 0)
                return dist [n] [m];
 
        // compute the distance and save it
        int d = cubsetdistance (*(sets [n]), *(sets [m]), true);
        dist [n] [m] = d;
        dist [m] [n] = d;
 
        // return the distance
        return d;
} /* MorseDecomposition::distance */
 
// --------------------------------------------------
 
template <class mapcomp, class cubetype, class cubsettype>
inline bool MorseDecomposition<mapcomp,cubetype,cubsettype>::connected
        (int n, int m) const
{
        return ((n >= 0) && (n < nsets) && (m >= 0) && (m < nsets) &&
                connDefined && conn [n] [m]);
} /* MorseDecomposition::connected */
 
template <class mapcomp, class cubetype, class cubsettype>
inline const cubsettype &MorseDecomposition<mapcomp,cubetype,cubsettype>::
	connection (int n, int m) const
{
        if ((n < 0) || (n >= nsets) || (m < 0) || (m >= nsets) ||
                !connDefined || !conn [n] [m])
                throw "Trying to access an undefined Morse set connection.";
        return *(conn [n] [m]);
} /* MorseDecomposition::connection */
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MorseDecomposition<mapcomp,cubetype,cubsettype>::distance
        (int n, int m, int k) const
{
        if ((n < 0) || (n >= nsets) || (m < 0) || (m >= nsets) ||
                (k < 0) || (k >= nsets))
                return -1;
        if ((n == m) || !connDefined || (!conn [n] [m] && !conn [m] [n]))
                return 0;
        return cubsetdistance (conn [n] [m] ? *(conn [n] [m]) :
                *(conn [m] [n]), *(sets [k]), false);
} /* MorseDecomposition::distance */
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MorseDecomposition<mapcomp,cubetype,cubsettype>::conndistance
        (int n, int m) const
{
        // if the number(s) are too large, return an undefined value
        if ((n >= nsets) || (m >= nsets))
                return -1;
 
        // if both numbers are negative, reset all the distances
        if ((n < 0) && (m < 0))
        {
                distconnDefined = false;
                return -1;
        }
 
        // if one of the numbers is negative and the tabulated distances
        // are not in use then no action is necessary
        if (((n < 0) || (m < 0)) && !distconnDefined)
        {
                return -1;
        }
 
        // initialize the distance array if it was not done before
        if (!distconnDefined)
        {
                for (int i = 0; i < nsets; ++ i)
                {
                        for (int j = 0; j < nsets; ++ j)
                        {
                                distconn [i] [j] = -1;
                        }
                }
                distconnDefined = true;
        }
 
        // if one of the numbers is negative, reset one set of distances only
        if ((n < 0) || (m < 0))
        {
                int num = (n >= 0) ? n : m;
                for (int i = 0; i < nsets; ++ i)
                {
                        distconn [i] [num] = -1;
                        distconn [num] [i] = -1;
                }
                return -1;
        }
 
        // if the distance is known then return it
        if (distconn [n] [m] >= 0)
                return distconn [n] [m];
 
        // compute the distance and save it
        int d = distance (n, m, n);
        distconn [n] [m] = d;
 
        // return the distance
        return d;
} /* MorseDecomposition::conndistance */
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MorseDecomposition<mapcomp,cubetype,cubsettype>::makegraph
        (chomp::homology::diGraph<> &g) const
{
        for (int i = 0; i < nsets; ++ i)
        {
                g. addVertex ();
                for (int j = 0; j < nsets; ++ j)
                {
                        if (connDefined && conn [i] [j])
                                g. addEdge (j);
                }
        }
        return 0;
} /* MorseDecomposition::makegraph */
 
// --------------------------------------------------
 
template <class mapcomp, class cubetype, class cubsettype>
inline const typename MorseDecomposition<mapcomp,cubetype,cubsettype>::
        IndexType &MorseDecomposition<mapcomp,cubetype,cubsettype>::compute
        (int n)
{
        if ((n < 0) || (n >= nsets) || !sets [n])
                throw "Trying to compute the Conley index of an undefined "
                        "Morse set.";
 
        // create an index pair and compute the map on the index pair
        PairType P (M);
        const cubsettype &cset = *(sets [n]);
        int_t csetsize = cset. size ();
        for (int_t i = 0; i < csetsize; ++ i)
                P. add (cset [i]);
        P. compute ();
 
        // create the Conley index object and compute the index
        if (indices [n])
                delete indices [n];
        indices [n] = new IndexType ();
        indices [n] -> setIndexPair (&P);
 
        indices [n] -> compute ();
        indices [n] -> setIndexPair (0);
 
        // remember if the index is nontrivial or not
        if (indices [n] -> trivial ())
                nontrivial [n] = 0;
        else
                nontrivial [n] = 1;
 
        return *(indices [n]);
} /* MorseDecomposition::compute */
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MorseDecomposition<mapcomp,cubetype,cubsettype>::compute ()
{
        for (int i = 0; i < nsets; ++ i)
                compute (i);
        return 0;
} /* MorseDecomposition::compute */
 
template <class mapcomp, class cubetype, class cubsettype>
inline bool MorseDecomposition<mapcomp,cubetype,cubsettype>::computed (int n)
        const
{
        return ((n >= 0) && (n < nsets) && indices [n]);
} /* MorseDecomposition::computed */
 
template <class mapcomp, class cubetype, class cubsettype>
inline const typename MorseDecomposition<mapcomp,cubetype,cubsettype>::
        IndexType &MorseDecomposition<mapcomp,cubetype,cubsettype>::index
        (int n) const
{
        if ((n < 0) || (n >= nsets) || !indices [n])
                throw "Trying to get an undefined Conley index.";
        return *(indices [n]);
} /* MorseDecomposition::index */
 
template <class mapcomp, class cubetype, class cubsettype>
inline const typename MorseDecomposition<mapcomp,cubetype,cubsettype>::
        IndexType &MorseDecomposition<mapcomp,cubetype,cubsettype>::setindex
        (int n, const typename
        MorseDecomposition<mapcomp,cubetype,cubsettype>::IndexType &ind)
{
        if ((n < 0) || (n >= nsets))
                throw "Trying to set a Conley index out of range.";
        if (indices [n])
                delete indices [n];
        indices [n] = new IndexType (ind);
        if (indices [n] -> trivial ())
                nontrivial [n] = 0;
        else
                nontrivial [n] = 1;
        return *(indices [n]);
} /* MorseDecomposition::setindex */
 
template <class mapcomp, class cubetype, class cubsettype>
inline bool MorseDecomposition<mapcomp,cubetype,cubsettype>::trivial
        (int n) const
{
        return ((n >= 0) && (n < nsets) && (nontrivial [n] == 0));
} /* MorseDecomposition::trivial */
 
// --------------------------------------------------
 
template <class mapcomp, class cubetype, class cubsettype>
inline void MorseDecomposition<mapcomp,cubetype,cubsettype>::intersection
        (const MorseDecomposition<mapcomp,cubetype,cubsettype> &m,
        const MorseDecomposition<mapcomp,cubetype,cubsettype> &n)
{
        for (int i = 0; i < m. nsets; ++ i)
        {
                for (int j = 0; j < n. nsets; ++ j)
                {
                        cubsettype X;
                        X. intersection (*(m. sets [i]), *(n. sets [j]));
                        if (X. empty ())
                                continue;
                        addone ();
                        sets [nsets - 1] -> swap (X);
                }
        }
        return;
} /* MorseDecomposition::intersection */
 
// --------------------------------------------------
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MorseDecomposition<mapcomp,cubetype,cubsettype>::join
        (int n, int m)
{
        if ((n < 0) || (n >= nsets) || (m < 0) || (m >= nsets) || (n == m))
                throw "Incorrect numbers of Morse sets to join.";
 
        // show a message
        chomp::homology::sbug << "* Joining Morse sets: " <<
                n << " and " << m << "(" << sets [n] -> size () <<
                " - " << ((conn [n] [m] ? conn [n] [m] -> size () :
                static_cast<int_t> (0)) +
                (conn [m] [n] ? conn [m] [n] -> size () :
                static_cast<int_t> (0))) << " - " <<
                sets [m] -> size () << ").\n";
 
        // reset tabulated distances
        if (distDefined)
        {
                for (int i = 0; i < nsets; ++ i)
                {
                        dist [i] [n] = -1;
                        dist [n] [i] = -1;
                }
        }
        if (distconnDefined)
        {
                for (int i = 0; i < nsets; ++ i)
                {
                        distconn [i] [n] = -1;
                        distconn [n] [i] = -1;
                }
        }
 
        // join the sets
        sets [n] -> add (*(sets [m]));
        if (connDefined && conn [n] [m])
                sets [n] -> add (*(conn [n] [m]));
        if (connDefined && conn [m] [n])
                sets [n] -> add (*(conn [m] [n]));
 
        // adjust the connecting orbits
        if (connDefined)
        {
                for (int k = 0; k < nsets; ++ k)
                {
                        if ((k == n) || (k == m))
                                continue;
 
                        // adjust the connection (n+m) --> k
                        if (conn [n] [k])
                        {
                                if (conn [m] [k])
                                {
                                        conn [n] [k] ->
                                                add (*(conn [m] [k]));
                                }
                                conn [n] [k] -> remove (*(sets [n]));
                        }
                        else if (conn [m] [k])
                        {
                                conn [n] [k] = conn [m] [k];
                                conn [m] [k] = 0;
                                conn [n] [k] -> remove (*(sets [n]));
                        }
 
                        // adjust the connection k --> (n+m)
                        if (conn [k] [n])
                        {
                                if (conn [k] [m])
                                {
                                        conn [k] [n] ->
                                                add (*(conn [k] [m]));
                                }
                                conn [k] [n] -> remove (*(sets [n]));
                        }
                        else if (conn [k] [m])
                        {
                                conn [k] [n] = conn [k] [m];
                                conn [k] [m] = 0;
                                conn [k] [n] -> remove (*(sets [n]));
                        }
                }
        }
 
        // reset the Conley index of the n-th Morse set if necessary
        if (indices [n] && nontrivial [m])
        {
                delete indices [n];
                indices [n] = 0;
                nontrivial [n] = -1;
        }
 
        // delete the m-th Morse set which has been joined to the n-th one
        deleteall (m, m + 1);
        return 0;
} /* MorseDecomposition::join */
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MorseDecomposition<mapcomp,cubetype,cubsettype>::passthru (int n)
{
        using chomp::homology::sbug;
 
        if ((n < 0) || (n >= nsets))
                throw "Incorrect numbers of Morse set to pass through.";
 
        // show a message
//      sbug << "* Passing through Morse set " << n << " (" <<
//              sets [n] -> size () << " cubes).\n";
 
        // join connecting orbits and reset tabulated distances
        for (int i = 0; i < nsets; ++ i)
        {
                if (i == n)
                        continue;
                for (int j = 0; j < nsets; ++ j)
                {
                        if (j == n)
                                continue;
                //      if (conn [i] [n] && !conn [n] [j])
                //      {
                //              sbug << "WARNIING: Passing through "
                //                      "an attractor: " << i << " -> " <<
                //                      n << " | " << j << ".\n";
                //      }
                //      else if (!conn [i] [n] && conn [n] [j])
                //      {
                //              sbug << "WARNIING: Passing through "
                //                      "a repeller: " << i << " | " <<
                //                      n << " -> " << j << ".\n";
                //      }
                //      else if (!conn [i] [n] && !conn [n] [j])
                //              continue;
                        if (connDefined && conn [i] [n] && conn [n] [j])
                        {
                        //      sbug << "(orbits joined: " << i <<
                        //              " --[" << conn [i] [n] -> size () <<
                        //              "]-> " << n << "[" <<
                        //              sets [n] -> size () << "] --[" <<
                        //              conn [n] [j] -> size () << "]-> " <<
                        //              j << ")\n";
                                // modify the connecting orbit
                                if (!conn [i] [j])
                                        conn [i] [j] = new cubsettype;
                                conn [i] [j] -> add (*(conn [i] [n]));
                                conn [i] [j] -> add (*(conn [n] [j]));
                                conn [i] [j] -> add (*(sets [n]));
 
                                // reset the corresponding distances
                                if (distconnDefined)
                                {
                                        distconn [i] [j] = -1;
                                        distconn [j] [i] = -1;
                                }
                        }
                }
        }
 
        // delete the n-th Morse set
        deleteall (n, n + 1);
        return 0;
} /* MorseDecomposition::passthru */
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MorseDecomposition<mapcomp,cubetype,cubsettype>::passthru
        (chomp::homology::hashedset<int> &setsToRemove)
{
        using chomp::homology::sbug;
 
        int nSetsToRemove = setsToRemove. size ();
//      if (nSetsToRemove < 100)
        {
                for (int index = nSetsToRemove - 1; index >= 0; -- index)
                {
                        this -> passthru (setsToRemove [index]);
                }
        }
 
        // TODO: A more efficient version for higher numbers of sets.
 
        return 0;
} /* MorseDecomposition::passthru */
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MorseDecomposition<mapcomp,cubetype,cubsettype>::jointrivial
        (int_t maxsetsize, int_t maxconnsize, int maxdistance)
{
        using chomp::homology::sbug;
        sbug << "- Join trivial: maxsetsize = " << maxsetsize <<
                ", maxconnsize = " << maxconnsize <<
                ", maxdistance = " << maxdistance << ".\n";
 
        // compute the Conley indices of all the Morse sets if necessary
        for (int n = 0; n < nsets; ++ n)
        {
                if (nontrivial [n] < 0)
                        compute (n);
        }
 
        // join pairs of Morse sets
        bool allow_trivial = false;
        while (1)
        {
                // output information for debugging
                if (false) // DEBUG
                {
                        using std::setw;
                        sbug << "Morse decomposition connection set sizes:\n";
                        sbug << "               ";
                        for (int i = 0; i < nsets; ++ i)
                                sbug << setw (7) << i;
                        sbug << "\n";
                        for (int i = 0; i < nsets; ++ i)
                        {
                                sbug << setw (3) << i << " (" <<
                                        setw (7) << sets [i] -> size ();
                                switch (nontrivial [i])
                                {
                                case -1: sbug << ")- "; break;
                                case 1: sbug << ")+ "; break;
                                case 0: sbug << "): "; break;
                                default: sbug << ")[" << nontrivial [i] <<
                                        "] "; break;
                                }
                                for (int j = 0; j < nsets; ++ j)
                                {
                                        if (connDefined && conn [i] [j])
                                        {
                                                sbug << setw (7) <<
                                                        conn [i] [j] ->
                                                        size ();
                                        }
                                        else
                                                sbug << "    ---";
                                }
                                sbug << "\n";
                        }
                }
 
                // the numbers of best choice in the loop below
                int best_n = -1, best_m = -1;
                // the size of what is added to the Morse set
        //      int_t best_size = 0;
                // the size of the nontrivial/large Morse set
        //      int_t best_nsize = 0;
                // the best distance between the Morse sets found so far
                int best_dist = 0;
                // the best distance between the connecting orbit and the set
                int best_conndist = 0;
 
                // check all pairs of sets and find the optimal ones to join
                for (int n = 0; n < nsets; ++ n)
                {
                        for (int m = 0; m < nsets; ++ m)
                        {
                                // --- check the necessary conditions ---
 
                                // make sure these are two different sets
                                if (n != m);
                                else
                                        continue;
 
                                // remember the sizes of the two sets
                                int_t nsize = sets [n] -> size ();
                                int_t msize = sets [m] -> size ();
 
                                // check if the n-th set is nontrivial
                                // or very large, and the m-th set is trivial
                                if (!nontrivial [m] && (nontrivial [n] ||
                                        (nsize > maxsetsize) ||
                                        allow_trivial));
                                else
                                        continue;
 
                                // count the size of connecting orbits
                                int_t conn_nm =
                                        (connDefined && conn [n] [m]) ?
                                        conn [n] [m] -> size () :
                                        static_cast<int_t> (0);
                                int_t conn_mn =
                                        (connDefined && conn [m] [n]) ?
                                        conn [m] [n] -> size () :
                                        static_cast<int_t> (0);
                        //      int_t add_size = msize + conn_nm + conn_mn;
 
                                // make sure the size limits are observed
                                if ((msize < maxsetsize) &&
                                        (conn_nm < maxconnsize) &&
                                        (conn_mn < maxconnsize));
                                else
                                        continue;
 
                                // make sure that no loop will be created
                                bool looping = false;
                                if (connDefined)
                                {
                                        for (int k = 0; k < nsets; ++ k)
                                        {
                                                if ((conn [n] [k] &&
                                                        conn [k] [m]) ||
                                                        (conn [m] [k] &&
                                                        conn [k] [n]))
                                                {
                                                        looping = true;
                                                        break;
                                                }
                                        }
                                }
                                if (looping)
                                        continue;
 
                                // compute the dist. between the Morse sets
                                int dist0 = distance (n, m);
                        //      sbug << "- Dist (" << n << "," << m <<
                        //              ") = " << dist0 << ".\n";
                                if (dist0 <= maxdistance);
                                else
                                        continue;
 
                                // compute the maximal distance between cubes
                                // in the connecting orbit and the Morse set
                                int conndist0 = conndistance (n, m);
                        //      sbug << "- ConnDist (" << n << "," << m <<
                        //              ") = " << conndist << ".\n";
                                if (conndist0 <= maxdistance);
                                else
                                        continue;
 
                                // --- check if this is a better choice ---
 
                                // make sure that the distance is not going
                                // to be increased
                                if ((best_n < 0) || (dist0 + conndist0 <=
                                        best_dist + best_conndist));
                                else
                                        continue;
 
                                // check if the size of what is added
                                // to the nontrivial set does not exceed
                                // the size found in the previous case
                        //      if ((best_n < 0) || (best_size < add_size));
                        //      else
                        //              continue;
 
                                // choose this pair of Morse sets
                                best_n = n;
                                best_m = m;
                        //      best_size = add_size;
                        //      best_nsize = nsize;
                                best_dist = dist0;
                        }
                }
 
                // if no good pair has been found...
                if ((best_n < 0) || (best_m < 0))
                {
                        // if loose criteria was already set, break the loop
                        if (allow_trivial)
                                break;
 
                        // loosen the criteria of checking the sets to join
                        allow_trivial = true;
                        continue;
                }
 
                // join the two Morse sets (keep the n-th index if necessary)
                join (best_n, best_m);
        }
 
        // compute the missing Conley indices
        for (int n = 0; n < nsets; ++ n)
        {
                if (nontrivial [n] < 0)
                        compute (n);
        }
 
        return 0;
} /* MorseDecomposition::jointrivial */
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MorseDecomposition<mapcomp,cubetype,cubsettype>::savetofiles
        (const char *prefix)
{
        for (int n = 0; n < nsets; ++ n)
        {
                // create a file for this Morse set
                std::ostringstream ofname;
                if (prefix)
                        ofname << prefix;
                ofname << n;
                std::string fname = ofname. str ();
                ofname << ".cub";
                std::ofstream out (ofname. str (). c_str ());
 
                // save the nontriviality of the index information
                out << "; This is the Morse set no. " << n << ".\n";
                out << "; Its Conley index ";
                switch (nontrivial [n])
                {
                case 0: out << "is trivial.\n"; break;
                case 1: out << "is nontrivial.\n"; break;
                default: out << "has not been computed.\n"; break;
                }
 
                // save the Conley index itself
                if (indices [n])
                {
                        out << "; The Conley index: " <<
                                *(indices [n]) << "\n";
                }
 
                // save the Morse set to a file
                out << *(sets [n]);
                out. close ();
 
                // save the connecting orbits to files
                for (int m = 0; m < nsets; ++ m)
                {
                        // if there is no connection, skip the file
                        if (!connDefined || !conn [n] [m])
                                continue;
 
                        // create a file name for the connection
                        std::ostringstream connname;
                        connname << fname << '-' << m << ".cub";
 
                        // save the connection
                        std::ofstream outc (connname. str (). c_str ());
                        outc << "; This is a covering of all the computed "
                                "connecting orbits\n"
                                "; that run from the Morse set no. " << n <<
                                " to " << m << ".\n";
                        if (conn [n] [m] -> empty ())
                        {
                                outc << "\n; Note: This set is empty "
                                        "which means that the orbit\n"
                                        "; was most likely not computed.\n";
                        }
                        else
                        {
                                outc << *(conn [n] [m]);
                        }
                }
        }
        return 0;
} /* MorseDecomposition::savetofiles */
 
template <class mapcomp, class cubetype, class cubsettype>
inline int MorseDecomposition<mapcomp,cubetype,cubsettype>::loadfromfiles
        (const char *prefix)
{
        for (int n = 0; ; ++ n)
        {
                // create a file name for this Morse set
                std::ostringstream fname;
                if (prefix)
                        fname << prefix;
                fname << n;
 
                // open a file if possible
                std::ifstream in (fname. str (). c_str ());
 
                // if the file does not exist, then quit the reading process
                if (!in)
                        break;
 
                // read the Conley index information if any is available
                int nontrivialindex = -1;
                IndexType *theIndex = 0;
                while (in. peek () == ';')
                {
                        // read the comment line from the input file
                        std::string line;
                        getline (in, line);
 
                        // look for the info that the index is nontrivial
                        if (line. find ("Conley index is nontrivial") !=
                                std::string::npos)
                        {
                                nontrivialindex = 1;
                                continue;
                        }
 
                        // look for the information that the index is trivial
                        if (line. find ("Conley index is trivial") !=
                                std::string::npos)
                        {
                                nontrivialindex = 0;
                                continue;
                        }
 
                        // look for the actual Conley index to read
                        std::string indHeader ("Conley index: ");
                        std::string::size_type pos = line. find (indHeader);
                        if (pos != std::string::npos)
                        {
                                std::istringstream str (line);
                                std::string::size_type len = pos +
                                        indHeader. size ();
                                for (unsigned int i = 0; i < len; ++ i)
                                        str. get ();
                                theIndex = new IndexType ();
                                str >> *theIndex;
                                continue;
                        }
                }
 
                // add a Morse set and read it
                this -> addone ();
                in >> *(sets [n]);
                in. close ();
                indices [n] = theIndex;
                nontrivial [n] = nontrivialindex;
        }
 
        // load the connections from files
        for (int n = 0; n < nsets; ++ n)
        {
                for (int m = 0; m < nsets; ++ m)
                {
                        // create a file name for this connection
                        std::ostringstream fname;
                        if (prefix)
                                fname << prefix;
                        fname << n << '-' << m;
 
                        // try opening the file for this connection
                        std::ifstream in (fname. str (). c_str ());
                        if (!in)
                                continue;
 
                        // create the connection and read it
                        addconn (n, m);
                        in >> *(conn [n] [m]);
                }
        }
        return 0;
} /* MorseDecomposition::loadfromfiles */
 
 
#endif // _CMGRAPHS_MORSEDEC_H_