indpair.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file indpair.h
///
/// The Conley index pair computation.
/// This file contains the definition of routines for the computation
/// of the Conley index pair using a combinatorial cubical multivalued map.
///
/// @author Pawel Pilarczyk
///
/////////////////////////////////////////////////////////////////////////////
 
// Copyright (C) 1997-2014 by Pawel Pilarczyk.
//
// This file is part of my research software package. This is free software:
// you can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// This software is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this software; see the file "license.txt". If not,
// please, see <https://www.gnu.org/licenses/>.
 
// Started on July 19, 2006. Last revision: March 19, 2013.
 
 
#ifndef _CMGRAPHS_INDPAIR_H_
#define _CMGRAPHS_INDPAIR_H_
 
 
// include some standard C++ header files
#include <iostream>
#include <algorithm>
#include <new>
#include <memory>
 
// include selected header files from the CHomP library
#include "chomp/system/textfile.h"
#include "chomp/system/timeused.h"
#include "chomp/struct/hashsets.h"
#include "chomp/cubes/cube.h"
 
// include local header files
#include "config.h"
#include "typedefs.h"
#include "spacewrap.h"
 
 
// --------------------------------------------------
// ------------------- 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:
        /// The type of a cube.
        typedef chomp::homology::Cube CubeType;
 
        /// The type of a set of cubes.
        typedef chomp::homology::SetOfCubes CubSetType;
 
        /// Returns the dimension of the underlying space.
        int dim () const {return 1;}
 
        /// Returns the number of cubes in the isolating neighborhood.
        int_t countInv () const {return 0;}
 
        /// Returns the number of cubes in the exit set.
        int_t countExit () const {return 0;}
 
        /// Returns the coordinates of the n-th cube
        /// (first Inv, then Exit, followed by outside).
        template <class IntType>
        int getcube (int_t n, IntType *coord) const {return 0;}
 
        /// Returns the number of cubes in the image of the n-th cube.
        int_t imgcount (int_t n) const {return 0;}
 
        /// Returns the number of the n-th cube in the image.
        int_t getimgcube (int_t n, int_t 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 type of a cube.
        typedef cubetype CubeType;
 
        /// The type of a set of cubes.
        typedef cubsettype CubSetType;
 
        /// The only allowed constructor.
        IndexPair (const mapcomp &_theCubMap = 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_t 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_t countExit () const;
 
        /// Returns the exit set as a set of cubes.
        const cubsettype &getExit () const;
 
        /// Retrieves the coordinates of the n-th cube.
        template <class IntType>
        int getcube (int_t n, IntType *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_t imgcount (int_t n) const;
 
        /// Returns the number of the i-th cube
        /// in the image of the n-th cube.
        int_t getimgcube (int_t n, int_t i) const;
 
        /// Computes the map on the entire set S, adds the necessary cubes
        /// to the exit set, and also computes the map on the exit set,
        /// so that the entire 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 from the 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 &theCubMap;
 
        /// Adds the image of the map computed on one cube.
        /// If the cube is in S, then the part of its image
        /// that sticks out of S is added to the exit set.
        int addimage (int_t n, const cubsettype &img);
 
}; /* class IndexPair */
 
// --------------------------------------------------
 
template <class mapcomp, class cubetype, class cubsettype>
inline IndexPair<mapcomp,cubetype,cubsettype>::IndexPair
        (const mapcomp &_theCubMap): theCubMap (_theCubMap)
{
        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_t 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_t 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>
template <class IntType>
inline int IndexPair<mapcomp,cubetype,cubsettype>::getcube
        (int_t n, IntType *coord) const
{
        int_t countS = S. size ();
        if (n < countS)
        {
                S [n]. coord (coord);
                return 0;
        }
        int_t 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_t IndexPair<mapcomp,cubetype,cubsettype>::imgcount (int_t n) const
{
        int_t 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_t IndexPair<mapcomp,cubetype,cubsettype>::getimgcube
        (int_t n, int_t i) const
{
        int_t countS = S. size ();
        const cubetype &q = (n < countS) ? (F (S [n])) [i] :
                (F (exitSet [n - countS])) [i];
        int_t num = S. getnumber (q);
        if (num >= 0)
                return num;
        num = exitSet. getnumber (q);
        if (num >= 0)
                return (countS + num);
        int_t 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>::addimage (int_t n,
        const cubsettype &img)
{
        // add the cubes in the image to the exit set or to the outside set
        int_t countS = S. size ();
        int_t imgSize = img. size ();
        for (int_t i = 0; i < imgSize; ++ i)
        {
                const cubetype &qImg (img [i]);
                if (n < countS)
                {
                        if (!S. check (qImg))
                                exitSet. add (qImg);
                }
                else // if (n >= countS)
                {
                        if (S. check (qImg))
                                throw "Wrong index pair detected.";
                        else if (!exitSet. check (qImg))
                                outside. add (qImg);
                }
        }
 
        // 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::addimage */
 
template <class mapcomp, class cubetype, class cubsettype>
inline int IndexPair<mapcomp,cubetype,cubsettype>::compute ()
{
        // set space wrapping if necessary and make sure that the wrapping
        // will be unset upon exiting this procedure;
        // if the servants are launched in separate processes then this is
        // not necessary, I only put it here for testing the local servant
        localSpaceWrapping wrapping (theCubMap. intWidth ());
 
        // determine the size of the isolating neighborhood
        int_t countS = S. size ();
 
        chomp::homology::sbug << "Computing the index pair... ";
 
        // compute the images of all the cubes
        int_t countSX = countS;
        for (int_t n = 0; n < countSX; ++ n)
        {
                // determine the cube whose image is to be computed
                const cubetype &q ((n < countS) ? S [n] :
                        exitSet [n - countS]);
 
                // compute the image of this cube
                cubsettype &img = F [q];
                theCubMap (q, &img, 0, 0, 0, !ignoreIsolationForConleyIndex);
 
                // add the image of the cube to respective sets if necessary
                addimage (n, img);
 
                // update the end-of-the-loop bound to include the exit set
                if (n == countS - 1)
                {
                        chomp::homology::sbug << "* ";
                        countSX = countS + exitSet. size ();
                }
        }
 
        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 */
 
 
#endif // _CMGRAPHS_INDPAIR_H_