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-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 January 3, 2007. Last revision: February 29, 2008.


#ifndef _CMGRAPHS_MORSEDEC_H_
#define _CMGRAPHS_MORSEDEC_H_


// include some standard C++ header files
#include <iostream>
#include <sstream>
#include <algorithm>

// include selected header files from the CHomP library
#include "chomp/system/textfile.h"
#include "chomp/homology/homology.h"
#include "chomp/struct/digraph.h"

// include selected header files from the CAPD library
#include "interval/doubleInterval.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);

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

        /// 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 maxsetsize, int 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;

        /// Tabulated distances between the Morse sets.
        mutable chomp::homology::multitable<chomp::homology::multitable<int> >
                dist;

        /// Tabulated values of the function "conndistance".
        mutable chomp::homology::multitable<chomp::homology::multitable<int> >
                distconn;

        /// 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)
{
        return;
} /* MorseDecomposition::MorseDecomposition */

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 and tabulated distances
        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 and tabulated distances
                for (int j = lastnumber; j < nsets; ++ j)
                {
                        conn [i] [j - difference] = conn [i] [j];
                        conn [i] [j] = 0;
                        dist [i] [j - difference] = dist [i] [j];
                        distconn [i] [j - difference] = distconn [i] [j];
                }
        }
        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 and tabulated distances
                for (int i = lastnumber; i < nsets; ++ i)
                {
                        conn [i - difference] [j] = conn [i] [j];
                        conn [i] [j] = 0;
                        dist [i - difference] [j] = dist [i] [j];
                        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
                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;
                        dist [i] [j] = m. dist [i] [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>::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
        for (int i = 0; i < nsets; ++ i)
        {
                conn [nsets] [i] = 0;
                conn [i] [nsets] = 0;
        }
        conn [nsets] [nsets] = 0;

        // initialize the tabulated distances
        for (int i = 0; i < nsets; ++ i)
        {
                dist [nsets] [i] = -1;
                dist [i] [nsets] = -1;
                distconn [nsets] [i] = -1;
                distconn [i] [nsets] = -1;
        }
        dist [nsets] [nsets] = 0;
        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 n (c, dim);
                        const typename cubetype::CoordType *nc;
                        while ((nc = n. 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 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 i = 0; i < X. size (); ++ i)
        {
                X [i]. coord (c);
                int local_dist = -1;
                for (int j = 0; j < Y. size (); ++ j)
                {
                        Y [j]. coord (d);
                        int dist = 0;
                        for (int k = 0; k < dim; ++ k)
                                dist += (c [k] - d [k]) * (c [k] - d [k]);
                        dist = static_cast<int> (std::sqrt (dist));
                        if ((local_dist < 0) || (local_dist > dist))
                                local_dist = dist;
                }

                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;
        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;
        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 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)
        {
                for (int i = 0; i < nsets; ++ i)
                {
                        for (int j = 0; j < nsets; ++ j)
                                dist [i] [j] = -1;
                }
                return -1;
        }

        // if one of the numbers is negative, reset one 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 sets are the same, their distance is zero
        if (n == m)
                return 0;

        // 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) &&
                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) ||
                !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) || (!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))
        {
                for (int i = 0; i < nsets; ++ i)
                {
                        for (int j = 0; j < nsets; ++ j)
                                distconn [i] [j] = -1;
                }
                return -1;
        }

        // 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 (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 csetsize = cset. size ();
        for (int 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 () : 0) +
                (conn [m] [n] ? conn [m] [n] -> size () : 0)) << " - " <<
                sets [m] -> size () << ").\n";

        // reset tabulated distances
        for (int i = 0; i < nsets; ++ i)
        {
                dist [i] [n] = -1;
                dist [n] [i] = -1;
                distconn [i] [n] = -1;
                distconn [n] [i] = -1;
        }

        // join the sets
        sets [n] -> add (*(sets [m]));
        if (conn [n] [m])
                sets [n] -> add (*(conn [n] [m]));
        if (conn [m] [n])
                sets [n] -> add (*(conn [m] [n]));

        // adjust the connecting orbits
        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 (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
                                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>::jointrivial
        (int maxsetsize, int 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 (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 best_size = 0;
                // the size of the nontrivial/large Morse set
                int 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 nsize = sets [n] -> size ();
                                int 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 conn_nm = conn [n] [m] ?
                                        conn [n] [m] -> size () : 0;
                                int conn_mn = conn [m] [n] ?
                                        conn [m] [n] -> size () : 0;
                                int 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;
                                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 dist = distance (n, m);
                        //      sbug << "- Dist (" << n << "," << m <<
                        //              ") = " << dist << ".\n";
                                if (dist <= maxdistance);
                                else
                                        continue;

                                // compute the maximal distance between cubes
                                // in the connecting orbit and the Morse set
                                int conndist = conndistance (n, m);
                        //      sbug << "- ConnDist (" << n << "," << m <<
                        //              ") = " << conndist << ".\n";
                                if (conndist <= 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) || (dist + conndist <=
                                        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 = dist;
                        }
                }

                // 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 ();
                std::ofstream out (fname. 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 (!conn [n] [m])
                                continue;

                        // create a file name for the connection
                        std::ostringstream connname;
                        connname << fname << '-' << m;

                        // 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
                        conn [n] [m] = new cubsettype;
                        in >> *(conn [n] [m]);
                }
        }
        return 0;
} /* MorseDecomposition::loadfromfiles */


#endif // _CMGRAPHS_MORSEDEC_H_

---