conindex.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file conindex.h
///
/// Conley index computation routines.
/// This file contains the definition of routines for the computation
/// of the Conley index using the CHomP library.
///
/// @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 July 19, 2006. Last revision: June 22, 2008.


#ifndef _CMGRAPHS_CONINDEX_H_
#define _CMGRAPHS_CONINDEX_H_


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

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

// include selected header files from the CAPD library
#include "interval/doubleInterval.h"

// include local header files
#include "config.h"


// --------------------------------------------------
// ------- the CLAPACK eigenvalues subroutine -------
// --------------------------------------------------

extern "C" {
int dgeev_ (char *jobvl, char *jobvr, long int *n,
        double * a, long int *lda, double *wr, double *wi, double *vl, 
        long int *ldvl, double *vr, long int *ldvr, double *work, 
        long int *lwork, long int *info);
}


// --------------------------------------------------
// ----------------- MapComputation -----------------
// --------------------------------------------------

/// The identity map class. This is a sample class that shows how one should
/// define a function class for the "MapComputation" class below.
class IdentityMap
{
public:
        /// Computes the image of a given box in terms of a rectangular box.
        int operator () (const interval *x, interval *y, int dim) const;

}; /* class IdentityMap */

inline int IdentityMap::operator () (const interval *x, interval *y, int dim)
        const
{
        for (int i = 0; i < dim; ++ i)
                y [i] = x [i];
        return 0;
} /* IdentityMap::operator () */

// --------------------------------------------------

/// The undefined map class.
class UndefinedMap
{
public:
        /// Computes the image of a given box in terms of a rectangular box.
        int operator () (const interval *x, interval *y, int dim) const;

}; /* class UndefinedMap */

inline int UndefinedMap::operator () (const interval *, interval *, int)
        const
{
        throw "Trying to use an undefined map.";
        return 0;
} /* UndefinedMap::operator () */

// --------------------------------------------------

/// The generic map computation routine that computes a rigorous
/// cubical multivalued map based on a function that computes the image
/// of intervals using interval arithmetic. The "mapcomp" class must have
/// a static method called "compute" which is used to compute the map
/// in interval arithmetic (see the sample class IdentityMap for details).
template <class mapcomp, class cubetype,
        class cubsettype = chomp::homology::hashedset<cubetype> >
class MapComputation
{
public:
        /// The default constructor.
        MapComputation (const double *_offset, const double *_width,
                int _intwidth, const mapcomp &_M = mapcomp ());

        /// The destructor.
        ~MapComputation ();

        /// The operator for computing the image of a box in terms of a set
        /// of boxes (a multivalued cubical map). The integral coefficients
        /// of cubes are transformed to real numbers according to the
        /// rectangular area defined by offset from the origin and its width,
        /// as well as the width of this area in terms of integral
        /// coefficients. Note: 'intwidth' does not have to be a power of 2.
        int operator () (const cubetype &q, cubsettype &img) const;

        /// Using cache for the map.
        bool cache;

        /// Is cropping of images to the designated box in effect?
        mutable bool cropping;

        /// Was the image cropped at least once? Reset this variable
        /// to "false" in order to detect image cropping again.
        mutable bool cropped;

        /// The maximal image diameter encountered so far.
        /// Set this variable to zero to gather meaningful information.
        static int maxImgDiam;

        /// The maximal image volume encountered so far.
        /// Set this variable to zero to gather meaningful information.
        static int maxImgVol;

private:
        /// The number of times the cache was used successfully.
        mutable int cacheused;

        /// The offset of the cubical rectangle.
        const double *offset;

        /// The width of the rectangle in each direction.
        const double *width;

        /// The width of the rectangle in terms of the number of cubes.
        int intwidth;

        /// The map used to compute an interval enclosure of the image
        /// of an interval box.
        const mapcomp &M;

        /// The cache domain.
        mutable cubsettype computed;

        /// The cache left vertex.
        mutable chomp::homology::multitable<cubetype> leftcache;

        /// The cache right vertex.
        mutable chomp::homology::multitable<cubetype> rightcache;

        /// Use interval arithmetic to compute the coordinate scope.
        /// Note: 'coord' may point to the same location
        /// as 'left' or 'right'.
        int compute (const typename cubetype::CoordType *coord, int dim,
                typename cubetype::CoordType *left,
                typename cubetype::CoordType *right) const;

}; /* class MapComputation */

// --------------------------------------------------

template <class mapcomp, class cubetype, class cubsettype>
int MapComputation<mapcomp,cubetype,cubsettype>::maxImgDiam;

template <class mapcomp, class cubetype, class cubsettype>
int MapComputation<mapcomp,cubetype,cubsettype>::maxImgVol;

// --------------------------------------------------

template <class mapcomp, class cubetype, class cubsettype>
inline MapComputation<mapcomp,cubetype,cubsettype>::MapComputation
        (const double *_offset, const double *_width,
        int _intwidth, const mapcomp &_M):
        cache (false), cropping (false), cropped (false), cacheused (0),
        offset (_offset), width (_width), intwidth (_intwidth), M (_M)
{
        return;
} /* MapComputation::MapComputation */

template <class mapcomp, class cubetype, class cubsettype>
inline MapComputation<mapcomp,cubetype,cubsettype>::~MapComputation ()
{
        // DEBUG
//      if (cache && !computed. empty ())
//              chomp::homology::sbug << "MAP: " << computed. size () <<
//                      " images cached, used " << cacheused << " times.\n";

        // DEBUG: Save the map.
        if (cache && false)
        {
                std::ofstream f ("f.map");
                for (int i = 0; i < computed. size (); ++ i)
                {
                        typename cubetype::CoordType c [cubetype::MaxDim];
                        f << "[" << computed [i];
                        computed [i]. coord (c);
                        for (int j = 0; j < computed [i]. dim (); ++ j)
                                ++ c [j];
                        f << cubetype (c, computed [i]. dim ()) << "] ";
                        f << "[" << leftcache [i] << rightcache [i] << "]\n";
                }
        }
        return;
} /* MapComputation::~MapComputation */

template <class mapcomp, class cubetype, class cubsettype>
inline int MapComputation<mapcomp,cubetype,cubsettype>::compute
        (const typename cubetype::CoordType *coord, int dim,
        typename cubetype::CoordType *left,
        typename cubetype::CoordType *right) const
{
        // compute the actual box corresponding to the cube
        std::auto_ptr<interval> xPtr (new interval [dim]);
        interval *x = xPtr. get ();
        for (int i = 0; i < dim; ++ i)
        {
                x [i] = interval (coord [i], coord [i] + 1);
                x [i] /= intwidth;
                x [i] *= width [i];
                x [i] += offset [i];
        }

        // compute the image of this box
        std::auto_ptr<interval> yPtr (new interval [dim]);
        interval *y = yPtr. get ();
        M (x, y, dim);

        // rescale the image box to the integer coordinates and make sure
        // that the actual image is contained in the INTERIOR of the box
        for (int i = 0; i < dim; ++ i)
        {
                y [i] -= offset [i];
                y [i] /= width [i];
                y [i] *= intwidth;

                double left_b = y [i]. leftBound ();
                left [i] =
                        static_cast<typename cubetype::CoordType> (left_b);
                for (int j = 0; (j < 4) && (left [i] >= left_b); ++ j)
                        -- left [i];
                if ((left [i] >= left_b) || (left_b - left [i] > 2))
                        throw "Map computation: Left bound out of range.";

                double right_b = y [i]. rightBound ();
                right [i] =
                        static_cast<typename cubetype::CoordType> (right_b);
                for (int j = 0; (j < 4) && (right [i] <= right_b); ++ j)
                        ++ right [i];
                if ((right [i] <= right_b) || (right [i] - right_b > 2))
                        throw "Map computation: Right bound out of range.";
        }

        return 0;
} /* MapComputation::compute */

template <class mapcomp, class cubetype, class cubsettype>
inline int MapComputation<mapcomp,cubetype,cubsettype>::operator ()
        (const cubetype &q, cubsettype &img) const
{
        typedef typename cubetype::CoordType coordType;

        // prepare arrays for the range of the image box
        int dim = q. dim ();
        std::auto_ptr<coordType> leftPtr (new coordType [dim]);
        coordType *left = leftPtr. get ();
        std::auto_ptr<coordType> rightPtr (new coordType [dim]);
        coordType *right = rightPtr. get ();

        // check if the image of this cube has already been computed
        int number = cache ? computed. getnumber (q) : -1;

        // compute the image of this box if necessary
        if (number < 0)
        {
                q. coord (left);
                compute (left, dim, left, right);

                // store this image in the cache if requested to
                if (cache)
                {
                        leftcache [computed. size ()] = cubetype (left, dim);
                        rightcache [computed. size ()] =
                                cubetype (right, dim);
                        computed. add (q);
                }
        }

        // use the cached values for the image of this box
        else
        {
                ++ cacheused;
                leftcache [number]. coord (left);
                rightcache [number]. coord (right);
        }

        // crop the image of this box to the given area if necessary
        if (cropping)
        {
                for (int i = 0; i < dim; ++ i)
                {
                        if (left [i] < 0)
                        {
                                left [i] = 0;
                                if (right [i] <= left [i])
                                        right [i] = left [i] + 1;
                                cropped = true;
                        }
                        if (right [i] > intwidth)
                        {
                                right [i] = intwidth;
                                if (left [i] >= right [i])
                                        left [i] = right [i] - 1;
                                cropped = true;
                        }
                }
        }

        // make sure that the size of the image is reasonably small
        int volume = 1;
        for (int i = 0; i < dim; ++ i)
        {
                int size = right [i] - left [i];
                if (maxImgDiam < size)
                        maxImgDiam = size;
                if (size > maxImageDiameter)
                        throw "Excessive diameter of box image.";
                volume *= size;
                if (volume > maxImageVolume)
                        throw "Excessive volume of box image.";
        }
        if (maxImgVol < volume)
                maxImgVol = volume;

        // cover the image box with cubes
        chomp::homology::tRectangle<typename cubetype::CoordType>
                r (left, right, dim);
        const typename cubetype::CoordType *c;
        while ((c = r. get ()) != 0)
        {
                img. add (cubetype (c, dim));
        }

        return 0;
} /* MapComputation::operator () */

// --------------------------------------------------

/// The undefined map computation class.
template<class cubetype, class cubsettype>
class UndefinedComputation
{
public:
        int operator () (const cubetype &q, cubsettype &img) const;

}; /* class UndefinedComputation */

template<class cubetype, class cubsettype>
int UndefinedComputation<cubetype, cubsettype>::operator ()
        (const cubetype &q, cubsettype &img) const
{
        throw "Trying to use undefined map computation.";
        return 0;
} /* UndefinedComputation::operator () */


// --------------------------------------------------
// ------------------- INDEX PAIR -------------------
// --------------------------------------------------
// The interface of a class that feeds cubes to the
// Conley index computation procedure.

/// The empty index pair class simulates an empty index pair.
/// This is a general template for you to build an appropriate index pair.
/// This might just be an interface to access an array of cubes in which
/// the isolating neighborhood occupies the front of the array, and then
/// cubes that form the exit set follow. Afterwards, all the other cubes
/// that appear in the images of those two groups of cubes are listed.
class EmptyIndexPair
{
public:
        /// Returns the dimension of the underlying space.
        int dim () const {return 1;}

        /// Returns the number of cubes in the isolating neighborhood.
        int countInv () const {return 0;}

        /// Returns the number of cubes in the exit set.
        int countExit () const {return 0;}

        /// Returns the coordinates of the n-th cube
        /// (first Inv, then Exit, followed by outside).
        int getcube (int n, int *coord) const {return 0;}

        /// Returns the number of cubes in the image of the n-th cube.
        int imgcount (int n) const {return 0;}

        /// Returns the number of the n-th cube in the image.
        int getimgcube (int n, int i) const {return 0;}

}; /* class EmptyIndexPair */

// --------------------------------------------------

/// A generic class that computes an index pair given the isolating
/// neighborhood and a means to compute the multivalued cubical map.
template <class mapcomp, class cubetype = chomp::homology::cube,
        class cubsettype = chomp::homology::hashedset<cubetype> >
class IndexPair
{
public:
        /// The only allowed constructor.
        IndexPair (const mapcomp &_M = mapcomp ());

        /// Adds a cube to S.
        int add (const cubetype &q);

        /// Returns the dimension of the phase space.
        int dim () const;

        /// Returns the number of cubes in the isolating neighborhood.
        int countInv () const;
        
        /// Returns the invariant part as a set of cubes.
        const cubsettype &getInv () const;

        /// Returns the number of cubes in the exit set.
        int countExit () const;

        /// Returns the exit set as a set of cubes.
        const cubsettype &getExit () const;

        /// Retrieves the coordinates of the n-th cube.
        int getcube (int n, int *coord) const;

        /// Returns the number of cubes in the image of the n-th cube
        /// (which is either in S, or in the exit set).
        int imgcount (int n) const;

        /// Returns the number of the i-th cube in the image of the n-th cube.
        int getimgcube (int n, int i) const;

        /// Computes the map on one cube. If the cube is in S, then its image
        /// that sticks out of S is added to the exit set.
        int compute (int n);

        /// Computes the map on the entire set S, and also on the exit set,
        /// so that the whole information necessary to compute the index map
        /// on this index pair is obatined and stored within this object.
        int compute ();

        /// Clears the sets and the stored map.
        int clear ();

private:
        /// The isolating neighborhood whose Conley index is computed.
        cubsettype S;

        /// The image of S without S.
        cubsettype exitSet;

        /// Additional images of cubes fromthe exit set.
        cubsettype outside;

        /// The multivalued cubical map that has been computed so far.
        chomp::homology::mvmap<cubetype,cubetype> F;

        /// The map object.
        const mapcomp &M;

}; /* class IndexPair */

// --------------------------------------------------

template <class mapcomp, class cubetype, class cubsettype>
inline IndexPair<mapcomp,cubetype,cubsettype>::IndexPair
        (const mapcomp &_M): M (_M)
{
        return;
} /* IndexPair::IndexPair */

template <class mapcomp, class cubetype, class cubsettype>
inline int IndexPair<mapcomp,cubetype,cubsettype>::dim () const
{
        if (!S. empty ())
                return S [0]. dim ();
        else if (!exitSet. empty ())
                return exitSet [0]. dim ();
        else
                return -1;
} /* IndexPair::dim */

template <class mapcomp, class cubetype, class cubsettype>
inline int IndexPair<mapcomp,cubetype,cubsettype>::countInv () const 
{
        return S. size ();
} /* IndexPair::countInv */

template <class mapcomp, class cubetype, class cubsettype>
inline const cubsettype &IndexPair<mapcomp,cubetype,cubsettype>::getInv ()
        const
{
        return S;
} /* IndexPair::getInv */

template <class mapcomp, class cubetype, class cubsettype>
inline int IndexPair<mapcomp,cubetype,cubsettype>::countExit () const 
{
        return exitSet. size ();
} /* IndexPair::countExit */

template <class mapcomp, class cubetype, class cubsettype>
inline const cubsettype &IndexPair<mapcomp,cubetype,cubsettype>::getExit ()
        const
{
        return exitSet;
} /* IndexPair::getExit */

template <class mapcomp, class cubetype, class cubsettype>
inline int IndexPair<mapcomp,cubetype,cubsettype>::getcube
        (int n, int *coord) const
{
        int countS = S. size ();
        if (n < countS)
        {
                S [n]. coord (coord);
                return 0;
        }
        int countExit = exitSet. size ();
        if (n < countS + countExit)
        {
                exitSet [n - countS]. coord (coord);
                return 0;
        }
        outside [n - countS - countExit]. coord (coord);
        return 0;
} /* IndexPair::getcube */

template <class mapcomp, class cubetype, class cubsettype>
inline int IndexPair<mapcomp,cubetype,cubsettype>::imgcount (int n) const
{
        int countS = S. size ();
        if (n < countS)
                return F (S [n]). size ();
        else
                return F (exitSet [n - countS]). size ();
} /* IndexPair::imgcount */

template <class mapcomp, class cubetype, class cubsettype>
inline int IndexPair<mapcomp,cubetype,cubsettype>::getimgcube
        (int n, int i) const
{
        int countS = S. size ();
        const cubetype &q = (n < countS) ? (F (S [n])) [i] :
                (F (exitSet [n - countS])) [i];
        int num = S. getnumber (q);
        if (num >= 0)
                return num;
        num = exitSet. getnumber (q);
        if (num >= 0)
                return (countS + num);
        int countExit = exitSet. size ();
        num = outside. getnumber (q);
        if (num >= 0)
                return (countS + countExit + num);
        throw "Index pair not prepared.";
//      num = countS + exitSet. size () + outside. size ();
//      outside. add (q);
//      return num;
} /* IndexPair::getimgcube */

// --------------------------------------------------

template <class mapcomp, class cubetype, class cubsettype>
inline int IndexPair<mapcomp,cubetype,cubsettype>::compute (int n)
{
        // determine the cube whose image is to be computed
        int countS = S. size ();
        const cubetype &q = (n < countS) ? S [n] : exitSet [n - countS];

        // compute the image of this cube
        cubsettype img;
        M (q, img);
        F [q]. add (img);

        // add the cubes in the image to the exit set or to the outside set
        for (int i = 0; i < img. size (); ++ i)
        {
                if ((n < countS) && !S. check (img [i]))
                        exitSet. add (img [i]);
                if (n >= countS)
                {
                        if (S. check (img [i]))
                                throw "Wrong index pair detected.";
                        else if (!exitSet. check (img [i]))
                                outside. add (img [i]);
                }
        }

        // verify if the size of the index pair is within the given limit
        if (countS + exitSet. size () + outside. size () > maxIndexPairSize)
                throw "Excessive index pair size.";

        return 0;
} /* IndexPair::compute */

template <class mapcomp, class cubetype, class cubsettype>
inline int IndexPair<mapcomp,cubetype,cubsettype>::compute ()
{
        chomp::homology::sbug << "Computing the index pair... ";

        int countS = S. size ();
        for (int i = 0; i < countS; ++ i)
                compute (i);
        int countSX = countS + exitSet. size ();
        for (int i = countS; i < countSX; ++ i)
                compute (i);

        chomp::homology::sbug << S. size () << " + " << exitSet. size () <<
                " + " << outside. size () << " cubes.\n";
        return 0;
} /* IndexPair::compute */

template <class mapcomp, class cubetype, class cubsettype>
inline int IndexPair<mapcomp,cubetype,cubsettype>::add (const cubetype &q)
{
        typename cubetype::CoordType coord [cubetype::MaxDim];
        q. coord (coord);
        cubetype r (coord, q. dim ());
        S. add (r);
        return 0;
} /* IndexPair::add */

template <class mapcomp, class cubetype, class cubsettype>
inline int IndexPair<mapcomp,cubetype,cubsettype>::clear ()
{
        cubsettype empty;
        S = empty;
        exitSet = empty;
        outside = empty;
        chomp::homology::mvmap<cubetype,cubetype> nomap;
        F = nomap;
        return 0;
} /* IndexPair::clear */


// --------------------------------------------------
// ------------------ CONLEY INDEX ------------------
// --------------------------------------------------

template <class IndexPair, class euclidom>
class ConleyIndex;

template <class IndexPair, class euclidom>
std::istream &operator >> (std::istream &in,
        ConleyIndex<IndexPair,euclidom> &c);

template <class IndexPair, class euclidom>
std::ostream &operator << (std::ostream &out,
        const ConleyIndex<IndexPair,euclidom> &c);

/// The class that computes and returns properties of the Conley index.
/// The index pair object is used to get the cubical index pair.
template <class IndexPair = EmptyIndexPair,
        class euclidom = chomp::homology::integer>
class ConleyIndex
{
public:
        /// The default constructor.
        ConleyIndex (const IndexPair *c = 0);

        /// The copy constructor.
        ConleyIndex (const ConleyIndex<IndexPair,euclidom> &c);

        /// The destructor.
        ~ConleyIndex ();

        /// The assignment operator.
        ConleyIndex<IndexPair,euclidom> &operator =
                (const ConleyIndex<IndexPair,euclidom> &c);

        /// Sets the index pair to be used to compute the Conley index.
        void setIndexPair (const IndexPair *f);

        /// Computes the Conley index. If requested, uses the two-layer
        /// construction. Otherwise, first tries the flat model, and switches
        /// to the two-layer construction only in the case of failure.
        int compute (bool twolayer = false);

        /// Reduces the index map using some shift equivalences to certain
        /// extent. Shrinks the matrices and generators if requested to.
        int reduce (bool shrink = true);
        
        /// Truncates the Conley index to the new top level.
        void truncate (int newlevel);

        /// Computes the eigenvalues of the index map.
        /// The eigenvalues are appended to the given vectors of doubles
        /// with the "push_back" method.
        /// Returns the number of eigenvalues or -1 in case of failure.
        int eigenvalues (int level, std::vector<double> &Re,
                std::vector<double> &Im, bool nonzero = true,
                bool sorted = true) const;

        /// Verifies if the index is trivial. Returns true if yes,
        /// false if it is non-trivial. Note: This information is tabulated.
        bool trivial () const;

        /// Retrieves the nth Betti number.
        int BettiNumber (int n) const;

        /// Retrieves the 'i'th [torsion] coefficient at the nth level.
        /// Returns 0 if 'i' is beyond the available scope.
        euclidom Coefficient (int n, int i) const;

        /// Retrieves the top level, or the maximal dimension of the index.
        /// The numbers of available levels are from 0 to dim inclusive.
        /// Returns -1 if the index is undefined.
        int dim () const;

        /// Retrieves the matrix of the index map at the given level.
        /// Returns 0 if the map at this level is zero.
        const chomp::homology::mmatrix<euclidom> *Map (int n) const;

        /// Retrieves a human-readable text representation of the map.
        std::string MapString (int n, const char *newline = "\n") const;

        /// Retrieves a human-readable text showing the eigenvalues.
        std::string EigenString (int n) const;

        /// Retrieves a human-readable text representation of the homology.
        std::string HomString () const;

        /// Defines the Conley index by the given Betti numbers
        /// and the given maps in homology.
        void define (const int *betti,
                const chomp::homology::mmatrix<euclidom> *matr,
                int _toplevel);

        /// Defines the Conley index by the given Betti numbers
        /// and adds the identity map.
        void define (const int *betti, int _toplevel);

private:
        /// An object that feeds this class with cubes.
        const IndexPair *cf;

        /// The number of the top homology level.
        int toplevel;

        /// The torsion coefficients of the corresponding homology
        /// generators. If delta () == 1, then no torsion.
        std::vector<std::vector<euclidom> > coef;

        /// The matrices of the index map at each level.
        chomp::homology::mmatrix<euclidom> *indmap;

        /// Converts the Conley index to a string (without the index pair).
        friend std::ostream &operator << <> (std::ostream &out,
                const ConleyIndex<IndexPair,euclidom> &c);

        /// Retrieves the Conley index from a string
        /// (without the index pair).
        friend std::istream &operator >> <> (std::istream &in,
                ConleyIndex<IndexPair,euclidom> &c);

        /// Rounds a real number to the given precision.
        /// If the absolute value of the number is below the given threshold
        /// then the number is truncated to 0.
        static double round (double x, double epsilon, double threshold);

        /// Variable that stores the tabulated information on whether
        /// the index is nontrivial. Values used: -1 = unknown, 0 = trivial,
        /// 1 = nontrivial.
        mutable int nontrivialindex;

}; /* class ConleyIndex */

// --------------------------------------------------

/// An auxiliary class that represents complex numbers.
/// It is used in the procedure for computing eigenvalues.
struct ppComplex
{
        double re, im;
        ppComplex (double _re, double _im): re (_re), im (_im) {}

}; /* struct ppComplex */

inline bool operator < (const ppComplex &x, const ppComplex &y)
{
        if (x. re < y. re)
                return true;
        else if (x. re > y. re)
                return false;
        return (x. im < y. im);
} /* operator < */

// --------------------------------------------------

template <class IndexPair, class euclidom>
inline ConleyIndex<IndexPair,euclidom>::ConleyIndex (const IndexPair *c)
{
        cf = c;
        toplevel = -1;
        indmap = 0;
        nontrivialindex = -1;
        return;
} /* ConleyIndex::ConleyIndex */

template <class IndexPair, class euclidom>
inline ConleyIndex<IndexPair,euclidom>::ConleyIndex
        (const ConleyIndex<IndexPair,euclidom> &c)
{
        cf = c. cf;
        toplevel = c. toplevel;
        coef = c. coef;
        if (toplevel >= 0)
                indmap = new chomp::homology::mmatrix<euclidom> [toplevel + 1];
        else
                indmap = 0;
        for (int i = 0; i <= toplevel; ++ i)
                indmap [i] = c. indmap [i];
        nontrivialindex = c. nontrivialindex;
        return;
} /* ConleyIndex::ConleyIndex */

template <class IndexPair, class euclidom>
inline ConleyIndex<IndexPair,euclidom> &
        ConleyIndex<IndexPair,euclidom>::operator =
        (const ConleyIndex<IndexPair,euclidom> &c)
{
        cf = c. cf;
        toplevel = c. toplevel;
        coef = c. coef;
        if (indmap)
                delete [] indmap;
        if (toplevel >= 0)
                indmap = new chomp::homology::mmatrix<euclidom> [toplevel + 1];
        else
                indmap = 0;
        for (int i = 0; i <= toplevel; ++ i)
                indmap [i] = c. indmap [i];
        nontrivialindex = c. nontrivialindex;
        return *this;
} /* ConleyIndex::operator = */

template <class IndexPair, class euclidom>
inline ConleyIndex<IndexPair,euclidom>::~ConleyIndex ()
{
        if (indmap)
                delete [] indmap;
        return;
} /* ConleyIndex::~ConleyIndex */

template <class IndexPair, class euclidom>
inline void ConleyIndex<IndexPair,euclidom>::setIndexPair
        (const IndexPair *f)
{
        cf = f;
        return;
} /* ConleyIndex::setIndexPair */

template <class IndexPair, class euclidom>
double ConleyIndex<IndexPair,euclidom>::round (double x, double epsilon,
        double threshold)
{
        if ((x < 0) && (x > -threshold))
                return 0;
        if ((x > 0) && (x < threshold))
                return 0;
        x /= epsilon;
        long xl = static_cast<long> (x);
        if (xl - x > 0.5)
                -- xl;
        else if (x - xl > 0.5)
                ++ xl;
        return (xl * epsilon);
} /* ConleyIndex::round */

template <class IndexPair, class euclidom>
inline int ConleyIndex<IndexPair,euclidom>::BettiNumber (int n) const
{
        if ((n < 0) || (n > toplevel))
                return 0;
        int size = coef [n]. size ();
        int count = 0;
        for (int i = 0; i < size; ++ i)
                if (coef [n] [i]. delta () == 1)
                        ++ count;
        return count;
} /* ConleyIndex::BettiNumber */

template <class IndexPair, class euclidom>
inline const chomp::homology::mmatrix<euclidom> *
        ConleyIndex<IndexPair,euclidom>::Map (int n) const
{
        if (!indmap || (n < 0) || (n > toplevel))
                return 0;
        return (indmap + n);
} /* ConleyIndex::Map */

template <class IndexPair, class euclidom>
inline euclidom ConleyIndex<IndexPair,euclidom>::Coefficient
        (int n, int i) const
{
        euclidom zero;
        zero = 0;
        if ((n < 0) || (n > toplevel))
                return zero;
        int size = coef [n]. size ();
        if ((i < 0) || (i >= size))
                return zero;
        return coef [n] [i];
} /* ConleyIndex::Coefficient */

template <class IndexPair, class euclidom>
inline int ConleyIndex<IndexPair,euclidom>::dim () const
{
        return toplevel;
} /* ConleyIndex::dim */

template <class IndexPair, class euclidom>
inline std::string ConleyIndex<IndexPair,euclidom>::MapString (int n,
        const char *newline) const
{
        if (!indmap || (n < 0) || (n > toplevel))
                return std::string ();
        std::ostringstream out;
        const chomp::homology::mmatrix<euclidom> &m = indmap [n];
        bool nonzero = false;
        for (int j = 0; j < m. getncols (); j ++)
                if (m. getcol (j). size ())
                {
                        if (!nonzero)
                        {
                                out << "Map " << n << ":" << newline;
                                nonzero = true;
                        }
                        out << "#" << (j + 1) << " = " << m. getcol (j) <<
                                newline;
                }
        if (!nonzero)
                out << "Map" << n << " = 0" << newline;
        return out. str ();
} /* ConleyIndex::MapString */

template <class IndexPair, class euclidom>
inline std::string ConleyIndex<IndexPair,euclidom>::HomString () const
{
        if (coef. size () <= 0)
                return std::string ("0");
        std::ostringstream out;
        int coefsize = coef. size ();
        for (int i = 0; i < coefsize; ++ i)
        {
                bool nonzero = false;
                if (i)
                        out << ", ";
                int counter = 0;
                int size = coef [i]. size ();
                for (int j = 0; j < size; ++ j)
                {
                        if (coef [i] [j]. delta () > 1)
                        {
                                if (nonzero)
                                        out << "+";
                                if (counter)
                                {
                                        out << "Z";
                                        if (counter > 1)
                                                out << "^" << counter;
                                        out << "+";
                                        counter = 0;
                                }
                                out << "Z_" << coef [i] [j];
                                nonzero = true;
                        }
                        else
                                ++ counter;
                }
                if (counter)
                {
                        if (nonzero)
                                out << "+";
                        out << "Z";
                        if (counter > 1)
                                out << "^" << counter;
                        nonzero = true;
                }
                if (!nonzero)
                        out << "0";
        }
        return out. str ();
} /* ConleyIndex::HomString */

template <class IndexPair, class euclidom>
inline int ConleyIndex<IndexPair,euclidom>::compute (bool twolayer)
{
        using namespace chomp::homology;

        // measure the time of creating data structures
        sbug << "Preparing to compute the Conley index...\n";
        timeused prepTime;

        // create the cubical map with the shifted dimension or not
        int shifted = twolayer ? 1 : 0; // raise dimension (see below)
        SetOfCubes X, A, Y, B;
        SetOfCubes Xflat, Aflat;
        int dim = cf -> dim ();
        int countS = cf -> countInv ();
        int c [Cube::MaxDim];
        typename Cube::CoordType coord [Cube::MaxDim];
        coord [dim] = 2;
        for (int i = 0; i < countS; ++ i)
        {
                cf -> getcube (i, c);
                for (int j = 0; j < dim; ++ j)
                {
                        coord [j] = static_cast<typename Cube::CoordType>
                                (c [j]);
                        if (coord [j] != c [j])
                                throw "Cube coordinatess out of range.";
                }
                Cube q (coord, dim + shifted);
                X. add (q);
                Y. add (q);
                Cube qflat (coord, dim);
                Xflat. add (qflat);
        }
        int countExit = cf -> countExit ();
        coord [dim] = 1;
        for (int i = 0; i < countExit; ++ i)
        {
                cf -> getcube (countS + i, c);
                for (int j = 0; j < dim; ++ j)
                {
                        coord [j] = static_cast<typename Cube::CoordType>
                                (c [j]);
                        if (coord [j] != c [j])
                                throw "Cube coordinatess out of range.";
                }
                Cube q (coord, dim + shifted);
                A. add (q);
                B. add (q);
                Cube qflat (coord, dim);
                Aflat. add (qflat);
        }
        coord [dim] = 0;
        for (int i = 0; i < countExit; ++ i)
        {
                int countImg = cf -> imgcount (countS + i);
                for (int j = 0; j < countImg; ++ j)
                {
                        int n = cf -> getimgcube (countS + i, j);
                        cf -> getcube (n, c);
                        for (int j = 0; j < dim; ++ j)
                        {
                                coord [j] = static_cast<typename
                                        Cube::CoordType> (c [j]);
                                if (coord [j] != c [j])
                                        throw "Coordinatess out of range.";
                        }
                        Cube q (coord, dim + shifted);
                        B. add (q);
                }
        }

        cubicalmap F;
        for (int i = 0; i < countS + countExit; ++ i)
        {
                if (i < countS)
                        coord [dim] = 2;
                else
                        coord [dim] = 1;
                cf -> getcube (i, c);
                for (int j = 0; j < dim; ++ j)
                        coord [j] = (typename Cube::CoordType) (c [j]);
                const Cube q (coord, dim + shifted);
                F [q];
                int countImg = cf -> imgcount (i);
                for (int j = 0; j < countImg; ++ j)
                {
                        int n = cf -> getimgcube (i, j);
                        cf -> getcube (n, c);
                        for (int j = 0; j < dim; ++ j)
                                coord [j] = (typename Cube::CoordType)
                                        (c [j]);
                        Cube qtest (coord, dim);
                        if (Xflat. check (qtest))
                                coord [dim] = 2;
                        else if (Aflat. check (qtest))
                                coord [dim] = 1;
                        else
                                coord [dim] = 0;
                        Cube r (coord, dim + shifted);
                        F [q]. add (r);
                }
        }

        // release the flat sets of cubes
        SetOfCubes emptyset;
        Xflat = emptyset;
        Aflat = emptyset;

        // DEBUG !!!
        static char debugfilename [] = "_/_x.cub";
        if (debugfilename [1] == '9')
                debugfilename [1] = 'A';
        if (debugfilename [1] < 'Z')
                ++ (debugfilename [1]);

        // DEBUG !!!
        if (false)
        {
                { debugfilename [3] = 'x'; std::ofstream f (debugfilename); f << X; }
                { debugfilename [3] = 'a'; std::ofstream f (debugfilename); f << A; }
                { debugfilename [3] = 'y'; std::ofstream f (debugfilename); f << Y; }
                { debugfilename [3] = 'b'; std::ofstream f (debugfilename); f << B; }
                { debugfilename [3] = 'f'; std::ofstream f (debugfilename); f << F; }
        }

        // indicate the time used for the preparation
        sbug << "Preparation completed in " << prepTime << ".\n";
        prepTime = 0;

        // measure the time of creating data structures
        sbug << "Computing the Conley index map in homology...\n";
        timeused compTime;

        // compute the map in homology
        chain<euclidom> *homXA = 0, *homYB = 0;
        chainmap<euclidom> *homF = 0;
        int maxlevelXA = -1, maxlevelYB = -1, maxlevelF = -1;
        bool error = false;
        try
        {
                chomp::homology::outputstream::mute mcon (scon);
                chomp::homology::outputstream::mute mout (sout);
                chomp::homology::outputstream::mute mbug (sbug);
                // since "Homology2l" is not well tested yet,
                // it is safer to raise the dimension (see above)
                // and use the alternative solution at the moment
        //      if (twolayer)
        //      {
        //              maxlevelF = Homology2l (F, "F", X, "X", A, "A",
        //                      homXA, maxlevelXA, homF, 0xFF);
        //              maxlevelYB = maxlevelXA;
        //      }
        //      else
                {
                        maxlevelF = Homology (F, "F", X, "X", A, "A",
                                Y, "Y", B, "B", homXA, maxlevelXA,
                                homYB, maxlevelYB, homF, true);
                }
        }
        catch (...)
        {
                // if failed for two-layer map, then this is a bad error
                if (twolayer)
                        throw;

                // otherwise the map must be computed using the new method
                else
                {
                        sout << "*** Homology computation failed. "
                                "Using a careful two-layer method. ***\n";
                        error = true;
                }
        }
        
        // if the levels don't agree, the map must be computed again
        if (!error && ((maxlevelF != maxlevelXA) ||
                (maxlevelXA != maxlevelYB)))
        {
                if (twolayer)
                        throw "Wrong homology levels in the Conley index.";
                sout << "*** Homology levels don't agree. "
                        "Using a careful two-layer method. ***\n";
                error = true;
        }

        // show the time used for the computation of the index map
        sbug << "The map computed in " << compTime << ".\n";
        compTime = 0;

        // if an error occurred for the flat model, use the two-layer one
        if (error)
        {
                delete [] homXA;
                delete [] homYB;
                delete homF;
                X = emptyset;
                A = emptyset;
                Y = emptyset;
                B = emptyset;
                cubicalmap emptymap;
                F = emptymap;
                return this -> compute (true);
        }

        toplevel = dim; /* previously I had here "dim + 1" --- why? */
        for (int level = 0; level <= toplevel; ++ level)
        {
                std::vector<euclidom> coefficients;
                if (level <= maxlevelXA)
                {
                        for (int i = 0; i < homXA [level]. size (); ++ i)
                                coefficients. push_back
                                        (homXA [level]. coef (i));
                }
                coef. push_back (coefficients);
        }
        if (indmap)
                delete [] indmap;
        if (toplevel >= 0)
                indmap = new mmatrix<euclidom> [toplevel + 1];
        else
        {
                indmap = 0;
                nontrivialindex = 0;
        }
        for (int i = 0; i <= maxlevelXA; ++ i)
                indmap [i] = (*homF) [i];

        // DEBUG !!!
        if (false && homF)
        {
                debugfilename [3] = 'h';
                std::ofstream f (debugfilename);
                f << *homF;
        }

        delete [] homXA;
        delete [] homYB;
        delete homF;
        return 0;
} /* ConleyIndex::compute */

template <class IndexPair, class euclidom>
inline int ConleyIndex<IndexPair,euclidom>::reduce (bool shrink)
{
        // remove those generators whose image is zero or which do not
        // appear in any image, i.e., their preimage is zero
        const int twice_the_level = (toplevel + 1) << 1;
        for (int level = 0; level < twice_the_level; ++ level)
        {
                // retrieve the right matrix and determine: column or row?
                bool columns = (level > toplevel);
                chomp::homology::mmatrix<euclidom> &matr = indmap [columns ?
                        (level - toplevel - 1) : level];

                // get the lengths of columns/rows of the matrix
                int n = columns ? matr. getncols () : matr. getnrows ();
                if (n <= 0)
                        continue;
                int *lengths = new int [n];
                for (int i = 0; i < n; ++ i)
                        lengths [i] = columns ? matr. getcol (i). size () :
                                matr. getrow (i). size ();

                // remove from each column these elements which have
                // zero columns themselves
                for (int i = 0; i < n; ++ i)
                {
                        int length = lengths [i];
                        if (!length)
                                continue;
                        chomp::homology::chain<euclidom> col_row = columns ?
                                matr. getcol (i) : matr. getrow (i);
                        for (int j = 0; j < length; ++ j)
                        {
                                int row_col = col_row. num (j);
                                if (lengths [row_col] > 0)
                                        continue;
                                matr. add (columns ? row_col : i,
                                        columns ? i : row_col,
                                        -col_row. coef (j));
                        }
                }

                // forget the column lengths
                delete [] lengths;
        }

        // shrink the matrices and generators if requested to
        for (int level = toplevel; shrink && (level >= 0); -- level)
        {
                // retrieve the references to this level's data
                chomp::homology::mmatrix<euclidom> &matr = indmap [level];
                std::vector<euclidom> &coefs = coef [level];

                // if the matrix is empty
                if (!matr. getncols () || !coefs. size ())
                {
                        coefs. clear ();
                        chomp::homology::mmatrix<euclidom> emptymatr;
                        matr = emptymatr;
                //      if (level == toplevel)
                //              -- toplevel;
                        continue;
                }

                // select the numbers of generators/columns to keep
                std::vector<int> selected;
                int ncols = matr. getncols ();
                for (int i = 0; i < ncols; ++ i)
                        if (matr. getcol (i). size () >= 1)
                                selected. push_back (i);
                int length = selected. size ();

                // if the index at this level is trivial, skip the rest
                if (!length)
                {
                        coefs. clear ();
                        chomp::homology::mmatrix<euclidom> emptymatr;
                        matr = emptymatr;
                //      if (level == toplevel)
                //              -- toplevel;
                        continue;
                }

                // create a new matrix
                chomp::homology::mmatrix<euclidom> newmatr;
                newmatr. define (length, length);
                for (int i = 0; i < length; ++ i)
                {
                        for (int j = 0; j < length; ++ j)
                        {
                                newmatr. add (i, j, matr. get (selected [i],
                                        selected [j]));
                        }
                }
                matr = newmatr;

                // create a new sequence of coefficients
                std::vector<euclidom> newcoefs (length);
                for (int i = 0; i < length; ++ i)
                        newcoefs [i] = coefs [selected [i]];
                coefs = newcoefs;
        }

        // change the sign of negative matrix entries where possible
        // unless 1 == -1 in the Euclidean domain in question
        euclidom e1, minus1;
        e1 = 1;
        minus1 = -e1;
        if (e1 != minus1)
        {
                for (int level = 0; level <= toplevel; ++ level)
                {
                        // retrieve the reference to the matrix at this level
                        chomp::homology::mmatrix<euclidom> &matr =
                                indmap [level];

                        // find single -1's and change their sign
                        int n = matr. getncols ();
                        for (int i = 0; i < n; ++ i)
                        {
                                const chomp::homology::chain<euclidom> &column =
                                        matr. getcol (i);
                                if (column. size () != 1)
                                        continue;
                                if (column. coef (0) != minus1)
                                        continue;
                                if (column. num (0) <= i)
                                        continue;
                                matr. multiplyrow (i, minus1);
                                matr. multiplycol (i, minus1);
                        }
                }
        }

        return 0;
} /* ConleyIndex::reduce */

template <class IndexPair, class euclidom>
inline void ConleyIndex<IndexPair,euclidom>::truncate (int newlevel)
{
        if ((newlevel >= -1) && (newlevel <= toplevel))
        {
                toplevel = newlevel;
                coef. erase (coef. begin () + toplevel + 1, coef. end ());
        }
        return;
} /* ConleyIndex<IndexPair,euclidom>::truncate */

template <class IndexPair, class euclidom>
inline int ConleyIndex<IndexPair,euclidom>::eigenvalues (int level,
        std::vector<double> &Re, std::vector<double> &Im,
        bool nonzero, bool sorted) const
{
        // if the level is out of scope, return the error code
        if ((level < 0) || (level > toplevel))
                return -1;

        // retrieve the matrix and its size
        const chomp::homology::mmatrix<euclidom> &m = indmap [level];
        int fullsize = m. getncols ();
        if (fullsize <= 0)
                return 0;

        // retrieve the indices not related to torsion
        int size = fullsize;
        const std::vector<euclidom> &coefs = coef [level];
        std::vector<int> indices (size);
        int current = 0;
        for (int i = 0; i < fullsize; ++ i)
        {
                if (coefs [i]. delta () > 1)
                {
                        indices [i] = -1;
                        -- size;
                }
                else
                {
                        indices [i] = current;
                        ++ current;
                }
        }
        if (size <= 0)
                return 0;

        // prepare data for the LAPACK function to compute the eigenvalues
        char N = 'N';
        long int dimension = size;
        std::auto_ptr<double> APtr (new double [size * size]);
        double *A = APtr. get ();
        std::auto_ptr<double> rePtr (new double [size]);
        double *re = rePtr. get ();
        std::auto_ptr<double> imPtr (new double [size]);
        double *im = imPtr. get ();
        long int worksize = 3 * size + (2 * size + 128);
        std::auto_ptr<double> workPtr (new double [worksize]);
        double *work = workPtr. get ();
        long int result = 0;

        // translate the chain matrix to the matrix of doubles
        double *ptr = A;
        for (int i = 0; i < size; ++ i)
        {
                for (int j = 0; j < size; ++ j)
                {
                        *(ptr ++) = 0;
                }
        }
        for (int col = 0; col < fullsize; ++ col)
        {
                int icol = indices [col];
                if (icol < 0)
                        continue;
                const chomp::homology::chain<euclidom> c = m. getcol (col);
                int len = c. size ();
                for (int i = 0; i < len; ++ i)
                {
                        int irow = indices [c. num (i)];
                        if (irow < 0)
                                continue;
                        euclidom e = c. coef (i);
                        int n = e. delta ();
                        euclidom e1;
                        e1 = n;
                        if (e1 != e)
                                n = -n;
                        A [size * irow + icol] = n;
                }
        }

        // call the LAPACK function to compute the eigenvalues of the map
        // (note: the eigenvalues of A^T are the same as of the eigenvalues
        // of A, so I don't care that the matrices are reprsented
        // in a different way in Fortran)
        dgeev_ (&N, &N, &dimension, A, &dimension, re, im,
                0, &dimension, 0, &dimension, work, &worksize, &result);

        // throw an error message in case of failure
        if (result < 0)
                throw "Conley Index eigenvalues: LAPACK error.";
        else if (result > 0)
                throw "The QR algorithm failed in the Conley Index "
                        "Eigenvalues computation.";

        // gather the output data
        std::vector<ppComplex> values;
        const double epsilon = 0.001;
        const double threshold = 0.01;
        for (int i = 0; i < size; ++ i)
        {
                re [i] = round (re [i], epsilon, threshold);
                im [i] = round (im [i], epsilon, threshold);
                if (nonzero && !re [i] && !im [i])
                        continue;
                values. push_back (ppComplex (re [i], im [i]));
        }

        // make a note of any nontrivial eigenvalues
        if (values. size ())
                nontrivialindex = 1;

        // sort the eigenvalues if requested to
        if (sorted)
                std::sort (values. begin (), values. end ());

        // write the results to the provided data structures
        for (std::vector<ppComplex>::const_iterator it = values. begin ();
                it != values. end (); ++ it)
        {
                Re. push_back (it -> re);
                Im. push_back (it -> im);
        }

        return values. size ();
} /* ConleyIndex<IndexPair,euclidom>::eigenvalues */

template <class IndexPair, class euclidom>
inline bool ConleyIndex<IndexPair,euclidom>::trivial () const
{
        // if the answer has already been computed before, use it now
        if (nontrivialindex >= 0)
                return !nontrivialindex;

        for (int level = 0; level <= toplevel; ++ level)
        {
                std::vector<double> Re, Im;
                int n = eigenvalues (level, Re, Im, true, false);
                if (n)
                {
                        nontrivialindex = 1;
                        return false;
                }
        //      if (coef [level]. size () && (indmap [level]. getncols () > 0))
        //              return false;
        }
        nontrivialindex = 0;
        return true;
} /* ConleyIndex::trivial */

template <class IndexPair, class euclidom>
inline std::string ConleyIndex<IndexPair,euclidom>::EigenString (int n) const
{
        // make sure the level is correct
        if (!indmap || (n < 0) || (n > toplevel))
                return std::string ();

        // prepare the output string stream
        std::ostringstream out;

        // compute the nonzero eigenvalues and sort them
        std::vector<double> Re, Im;
        int count = eigenvalues (n, Re, Im);
        if (count > 0)
                nontrivialindex = 1;

        // output the eigenvalues
        out << "Eigenvalues " << n << ": (";
        if (count >= 0)
        {
                for (int i = 0; i < count; ++ i)
                {
                        if (i)
                                out << ", ";
                        if (Re [i])
                        {
                                out << Re [i];
                                if (Im [i] > 0)
                                        out << "+";
                        }
                        if (Im [i])
                                out << Im [i] << "i";
                }
        }
        else
                out << "ERROR";
        out << ").";

        return out. str ();
} /* ConleyIndex::EigenString */

template <class IndexPair, class euclidom>
inline void ConleyIndex<IndexPair,euclidom>::define (const int *betti,
        const chomp::homology::mmatrix<euclidom> *matr, int _toplevel)
{
        // reset the potentially remembered nontriviality
        nontrivialindex = -1;

        // set the top level of homology
        toplevel = _toplevel;

        // create the coefficients vector
        coef. clear ();
        euclidom e;
        e = 1;
        for (int i = 0; i <= toplevel; ++ i)
        {
                std::vector<euclidom> b;
                for (int j = betti [i]; j > 0; -- j)
                        b. push_back (e);
                coef. push_back (b);
        }

        // set the index maps
        if (indmap)
        {
                delete [] indmap;
                indmap = 0;
        }
        if (toplevel >= 0)
        {
                indmap = new chomp::homology::mmatrix<euclidom>
                        [toplevel + 1];
        }
        for (int level = 0; level <= toplevel; ++ level)
        {
                indmap [level] = matr [level];
        }
        return;
} /* ConleyIndex::define */

template <class IndexPair, class euclidom>
inline void ConleyIndex<IndexPair,euclidom>::define (const int *betti,
        int _toplevel)
{
        // create an array of the identity maps
        chomp::homology::mmatrix<euclidom> *matr = (_toplevel < 0) ? 0 :
                new chomp::homology::mmatrix<euclidom> [_toplevel + 1];
        for (int level = 0; level <= _toplevel; ++ level)
        {
                if (betti [level] > 0)
                        matr [level]. identity (betti [level]);
        }

        // define the index
        this -> define (betti, matr, _toplevel);

        // clean up and finalize
        if (matr)
                delete [] matr;
        return;
} /* ConleyIndex::define */

// --------------------------------------------------

template <class IndexPair, class euclidom>
std::ostream &operator << (std::ostream &out,
        const ConleyIndex<IndexPair,euclidom> &c)
{
        // write the top level
        out << c. toplevel << " ";
        if (c. toplevel < 0)
                return out;

        // write the coefficients
        for (int i = 0; i <= c. toplevel; ++ i)
        {
                int size = c. coef [i]. size ();
                out << size << " ";
                for (int j = 0; j < size; ++ j)
                        out << c. coef [i] [j] << " ";
        }

        // write the maps
        for (int i = 0; i <= c. toplevel; ++ i)
        {
                const chomp::homology::mmatrix<euclidom> &m = c. indmap [i];
                out << m. getnrows () << " " << m. getncols () << " ";
                for (int j = 0; j < m. getncols (); ++ j)
                {
                        const chomp::homology::chain<euclidom> &col =
                                m. getcol (j);
                        out << col. size () << " ";
                        for (int k = 0; k < col. size (); ++ k)
                        {
                                out << col. num (k) << " ";
                                out << col. coef (k) << " ";
                        }
                }
        }

        return out;
} /* operator << */

template <class IndexPair, class euclidom>
std::istream &operator >> (std::istream &in,
        ConleyIndex<IndexPair,euclidom> &c)
{
        // reset the data
        c. nontrivialindex = -1;
        if (c. indmap)
        {
                delete [] c. indmap;
                c. indmap = 0;
        }
        c. coef. clear ();

        // read the top level
        in >> c. toplevel;
        if (c. toplevel < 0)
                return in;

        // read the coefficients
        for (int i = 0; i <= c. toplevel; ++ i)
        {
                int size = -1;
                in >> size;
                if (size < 0)
                        throw "Cannot decode the Conley index from input.";
                std::vector<euclidom> coef;
                for (int j = 0; j < size; ++ j)
                {
                        euclidom e;
                        e = 0;
                        in >> e;
                        if (e. delta () <= 0)
                                throw "Wrong coefficient in the "
                                        "Conley index from input.";
                        coef. push_back (e);
                }
                c. coef. push_back (coef);
        }

        // read the maps
        if (c. toplevel >= 0)
                c. indmap = new chomp::homology::mmatrix<euclidom>
                        [c. toplevel + 1];
        for (int i = 0; i <= c. toplevel; ++ i)
        {
                chomp::homology::mmatrix<euclidom> &m = c. indmap [i];
                int nrows = 0, ncols = 0;
                in >> nrows >> ncols;
                if ((nrows < 0) || (ncols < 0))
                        throw "Wrong matrix size in the Conley index "
                                "from input.";
                if ((nrows > 0) && (ncols > 0))
                        m. define (nrows, ncols);
                for (int j = 0; j < ncols; ++ j)
                {
                        int length = -1;
                        in >> length;
                        if (length < 0)
                                throw "Negative column length in the "
                                        "Conley index from input.";
                        for (int k = 0; k < length; ++ k)
                        {
                                int num = -1;
                                in >> num;
                                if ((num < 0) || (num >= nrows))
                                        throw "Wrong row number in the "
                                                "Conley index from input.";
                                euclidom e;
                                e = 0;
                                in >> e;
                                if (e. delta () <= 0)
                                        throw "Wrong matrix entry in the "
                                                "Conley index from input.";
                                m. add (num, j, e);
                        }
                }
        }

        return in;
} /* operator >> */


#endif // _CMGRAPHS_CONINDEX_H_

---