invpart.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file invpart.h
///
/// Computation of the invariant part of a set of cubes.
/// This file contains the definition of a function
/// for the computation of the invariant part of a set of cubes
/// with respect to a given multivalued map.
/// It also contains the procedure for verifying whether the invariant part
/// is nonempty after subdividing the set of cubes a few times.
///
/// @author Pawel Pilarczyk
///
/////////////////////////////////////////////////////////////////////////////

// Copyright (C) 1997-2008 by Pawel Pilarczyk.
//
// This file is part of my research software package.  This is free software;
// you can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation;
// either version 2 of the License, or (at your option) any later version.
//
// This software is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this software; see the file "license.txt".  If not, write to the
// Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
// MA 02111-1307, USA.

// Started on February 12, 2007. Last revision: March 6, 2008.


#ifndef _CMGRAPHS_INVPART_H_
#define _CMGRAPHS_INVPART_H_


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

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


// --------------------------------------------------
// ----------------- invariant part -----------------
// --------------------------------------------------

/// Computes X := Inv (X) using the algorithm by Bill Kalies and Hyunju Ban.
/// If the graph 'gInv' is given, then the resulting graph is created,
/// otherwise only the set X is modified.
template <class typeCubes>
inline void invariantPart (typeCubes &X, const chomp::homology::diGraph<> &g,
        chomp::homology::diGraph<> *gInv = 0)
{
        // do nothing if the set X is empty
        if (X. empty ())
                return;

        // compute the chain recurrent set of the graph
        // and remember the transposed graph
        chomp::homology::multitable<int> compVertices, compEnds;
        chomp::homology::diGraph<> gT;
        int count = chomp::homology::SCC (g, compVertices, compEnds,
                static_cast<chomp::homology::diGraph<> *> (0), &gT);

        // retrieve vertices that can be reached from the chain recurrent set
        int nVertices = X. size ();
        char *forward = new char [nVertices];
        for (int i = 0; i < nVertices; ++ i)
                forward [i] = 0;
        for (int cur = 0; cur < count; ++ cur)
                g. DFScolor (forward, '\1',
                        compVertices [cur ? compEnds [cur - 1] : 0]);

        // retrieve vertices that can reach the chain recurrent set
        char *backward = new char [nVertices];
        for (int i = 0; i < nVertices; ++ i)
                backward [i] = 0;
        for (int cur = 0; cur < count; ++ cur)
                gT. DFScolor (backward, '\1',
                        compVertices [cur ? compEnds [cur - 1] : 0]);

        // compute the new set of cubes and exit if gInv is null
        typeCubes invX;
        if (!gInv)
        {
                for (int i = 0; i < nVertices; ++ i)
                        if (forward [i] && backward [i])
                                invX. add (X [i]);
        }
        else
        {
                // join both tables and compute the new set of cubes
                for (int i = 0; i < nVertices; ++ i)
                {
                        forward [i] &= backward [i];
                        if (forward [i])
                                invX. add (X [i]);
                }

                // restrict the graph to the colored vertices
                g. subgraph (*gInv, forward);
        }

        // swap the set of cubes and the corresponding graph
        X. swap (invX);

        // clean the memory and return
        delete [] backward;
        delete [] forward;
        return;
} /* invariantPart */


/// Computes X := Inv (X) using the algorithm by Bill Kalies and Hyunju Ban.
/// Uses the combinatorial cubical multivalued map provided.
/// Returns 0 on success or -1 if failed. Shows messages to 'sbug'.
template <class typeCubes, class theCubMapType>
inline int invariantPart (typeCubes &X, const theCubMapType &theCubMap)
{
        using chomp::homology::sbug;

        // compute the graph representation of the cubical map
        sbug << "F ";
        chomp::homology::diGraph<> g;
        try
        {
                // compute the graph that represents the map restricted to X
                computeMapGraph (X, g, theCubMap);

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

        }
        catch (const char *msg)
        {
                // return -1 in case of failure
                sbug << "Failed: " << msg << "\n";
                return -1;
        }

        // compute the invariant part of X in case of success
        sbug << "Inv... ";
        invariantPart (X, g);
        sbug << X. size () << " left.\n";

        return 0;
} /* invariantPart */


/// Verifies whether the invariant part of the given set is empty
/// by applying a limited number of refinements until the set is empty
/// or the number of cubes becomes too large.
/// Returns "true" if it is empty, "false" if this could not be proved.
inline bool checkEmptyInv (const spcCubes &X, int subdivDepth,
        const double *offset, const double *width, const theMapType &theMap,
        const theCubMapType &theCubMap0, const theCubMapType &theCubMap1)
{
        using chomp::homology::sbug;
        using chomp::homology::diGraph;

        if (X. empty () || (X. size () > maxRefineSize0))
                return false;
        spcCubes Z;
        int subdivisions = 0;
        while (++ subdivisions <= refineDepth)
        {
                // subdivide X or Z to Y
                spcCubes Y;
                const spcCubes &source = Z. empty () ? X : Z;
                sbug << "* Subdiv " << source. size () << " -> ";
                subdivideCubes (source, Y);
                sbug << Y. size () << ". ";
                spcCubes emptySet;
                Z = emptySet;

                // compute the invariant part of Y
                theCubMapType subdivCubMap (offset, width,
                        1 << (subdivDepth + subdivisions), theMap);
                subdivCubMap. cache = false;
                const theCubMapType &theCubMap =
                        (subdivisions > 1) ? subdivCubMap :
                        (subdivisions ? theCubMap1 : theCubMap0);
                int result = invariantPart (Y, theCubMap);

                // if the computation failed then return the negative result
                if (result < 0)
                        return false;

                // return the positive result if the invariant part is empty
                if (Y. empty ())
                {
                        sbug << "* Empty after " << subdivisions <<
                                " subdivisions).\n";
                        return true;
                }

                // remember the current set of cubes as Z
                Z. swap (Y);

                // exit the loop if there are too many cubes
                if (Z. size () > maxRefineSize1)
                        break;
        }
        sbug << "* Nonempty after " << subdivisions << " subdivisions (" <<
                Z. size () << " cubes left).\n";
        return false;
} /* checkEmptyInv */


#endif // _CMGRAPHS_INVPART_H_

---