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-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 February 12, 2007. Last revision: October 9, 2022.
 
 
#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"
#include "typedyns.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_t> compVertices (8192);
        chomp::homology::multitable<int_t> compEnds (8192);
        chomp::homology::diGraph<> gT;
        int_t 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_t nVertices = X. size ();
        char *forward = new char [nVertices];
        for (int_t i = 0; i < nVertices; ++ i)
                forward [i] = 0;
        for (int_t cur = 0; cur < count; ++ cur)
        {
                g. DFScolor (forward, '\1',
                        compVertices [cur ? compEnds [cur - 1] :
                        static_cast<int_t> (0)]);
        }
 
        // retrieve vertices that can reach the chain recurrent set
        char *backward = new char [nVertices];
        for (int_t i = 0; i < nVertices; ++ i)
                backward [i] = 0;
        for (int_t cur = 0; cur < count; ++ cur)
        {
                gT. DFScolor (backward, '\1',
                        compVertices [cur ? compEnds [cur - 1] :
                        static_cast<int_t> (0)]);
        }
 
        // compute the new set of cubes and exit if gInv is null
        typeCubes invX;
        if (!gInv)
        {
                for (int_t 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_t 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,
        bool cropping = true)
{
        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, cropping);
 
                // show the average size of the image of a cube
                int_t avgImgSize = g. countEdges () /
                        (X. size () ? X. size () : static_cast<int_t> (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;
                ++ subdivisions;
 
                // compute the invariant part of Y
                theCubMapType subdivCubMap (offset, width,
                        1 << (subdivDepth + subdivisions),
                        subdivDepth + subdivisions, theMap);
                subdivCubMap. cache = false;
                if (subdivisions > 1)
                {
                        subdivCubMap. setAllowedImgSize
                                (maxImageDiameter << 3, maxImageVolume << 3);
                }
                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 == 1) ? " subdivision.\n" :
                                " 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 == 1) ? " subdivision (" :
                " subdivisions (") << Z. size () << " cubes left).\n";
        return false;
} /* checkEmptyInv */
 
/// Computes the strongly connected path components
/// of the graph representation of the multivalued map on the given set.
/// Returns the number of components and fills in the array provided
/// or returns -1 in case of failure.
/// Shows messages to 'sbug'.
template <class CubSetType, class CubSetArrayType, class CubMapType>
inline int_t findSCCs (const CubSetType &X, const CubMapType &theCubMap,
        CubSetArrayType &components)
{
        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_t avgImgSize = g. countEdges () /
                        (X. size () ? X. size () : static_cast<int_t> (1));
                sbug << "[avg " << avgImgSize << "] ";
 
        }
        catch (const char *msg)
        {
                // return -1 in case of failure
                sbug << "Failed: " << msg << "\n";
                return -1;
        }
 
        // compute the components of the chain recurrent set of the graph
        sbug << "SCC... ";
        chomp::homology::multitable<int_t> compVertices (8192);
        chomp::homology::multitable<int_t> compEnds (8192);
        chomp::homology::diGraph<> *sccAddr = 0;
        chomp::homology::diGraph<> *transpAddr = 0;
        int_t nSets = chomp::homology::SCC (g, compVertices, compEnds,
                sccAddr, transpAddr);
        sbug << nSets << ".\n";
 
        // put the corresponding cubes into the array of sets of cubes
        int_t cur = 0;
        for (int_t n = 0; n < nSets; ++ n)
        {
                int_t end = compEnds [n];
                while (cur < end)
                {
                        components [n]. add (X [compVertices [cur]]);
                        ++ cur;
                }
        }
 
        return nSets;
} /* findSCCs */
 
 
 
#endif // _CMGRAPHS_INVPART_H_