matchdec.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file matchdec.h
///
/// Matching Morse decompositions.
/// This file contains the definitions of a procedure for matching
/// Morse decompositions. A relation between the sets of Morse sets
/// for two Morse decompositions is determined by checking which Morse sets
/// intersect their counterparts. If this relation is a bijection then
/// the edges in Morse graphs is also verified.
///
/// @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 16, 2008. Last revision: March 19, 2008.


#ifndef _CMGRAPHS_MATCHDEC_H_
#define _CMGRAPHS_MATCHDEC_H_


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

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

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


// --------------------------------------------------
// --------- intersections of cubical sets ----------
// --------------------------------------------------

/// Verifies whether two sets of cubes intersect each other.
/// Returns "true" if they share at least one common cube
/// or "false" otherwise.
template <class typeCubes>
inline bool intersect (const typeCubes &X, const typeCubes &Y)
{
        // choose the smaller set and the larger one
        bool XY = (X. size () < Y. size ());
        const typeCubes &smaller = XY ? X : Y;
        const typeCubes &larger = XY? Y : X;

        // check each cube in the smaller set if it is in the larger one
        int smallerSize = smaller. size ();
        for (int i = 0; i < smallerSize; ++ i)
        {
                if (larger. check (smaller [i]))
                        return true;
        }

        return false;
} /* intersect */

/// Verifies whether two sets of cubes intersect each other.
/// The first set is with respect to a twice finer grid that the other one.
/// Returns "true" if they share at least one common cube
/// or "false" otherwise.
template <class typeCubes>
inline bool intersectSub (const typeCubes &fine, const typeCubes &coarse)
{
        typedef typename typeCubes::value_type typeCube;
        typedef typename typeCube::CoordType typeCoord;

        // determine the number of cubes in the first set
        int size = fine. size ();
        if (!size)
                return false;

        // determine the dimension of the cubes and allocate working memory
        int dim = fine [0]. dim ();
        std::auto_ptr<typeCoord> coordPtr (new typeCoord [dim]);
        typeCoord *coord = coordPtr. get ();

        // check each cube in the finer set if it is contained
        // in some cube in the coarser set
        for (int i = 0; i < size; ++ i)
        {
                fine [i]. coord (coord);
                for (int j = 0; j < dim; ++ j)
                        coord [j] /= 2;
                if (coarse. check (typeCube (coord, dim)))
                        return true;
        }

        return false;
} /* intersectSub */


// --------------------------------------------------
// -------------- matching Morse sets ---------------
// --------------------------------------------------

/// Verifies if the two given Morse decompositions are equivalent
/// by matching Morse sets which intersect each other.
/// If requested, matches also the ordering in the Morse decomposition.
inline bool matchMorseSets (const theMorseDecompositionType &morseDecM,
        const std::vector<int> &wrongIndM, const theCubMapType &theCubMap1M,
        std::vector<spcCubes> &refinedM,
        const theMorseDecompositionType &morseDecN,
        const std::vector<int> &wrongIndN, const theCubMapType &theCubMap1N,
        std::vector<spcCubes> &refinedN,
        bool matchOrdering)
{
        using chomp::homology::sbug;

        // check the numbers of Morse sets in both Morse decompositions
        int nSets = morseDecM. count ();
        if (nSets != morseDecN. count ())
                return false;
        if (!nSets)
                return true;

        // try to determine which Morse sets in one Morse decomposition
        // intersect which Morse sets in the other Morse decomposition
        std::vector<std::vector<int> > MtoN (nSets);
        std::vector<std::vector<int> > NtoM (nSets);
        for (int m = 0; m < nSets; ++ m)
        {
                for (int n = 0; n < nSets; ++ n)
                {
                        if (!intersect (morseDecM [m], morseDecN [n]))
                                continue;
                        MtoN [m]. push_back (n);
                        NtoM [n]. push_back (m);
                }
        }

        // check the bijection of intersections
        bool bijection = true;
        for (int i = 0; i < nSets; ++ i)
        {
                if (MtoN [i]. empty () || NtoM [i]. empty ())
                        return false;
                if ((MtoN [i]. size () > 1) || (NtoM [i]. size () > 1))
                        bijection = false;
        }

        // in case of no bijection, try refining some of the Morse sets
        if (!bijection)
        {
                // show a message on what is going to be done
                sbug << "Note: No bijection in matching. Refining...\n";

                // increase the arrays of refined Morse sets
                if (refinedM. empty ())
                        refinedM. resize (nSets);
                if (refinedN. empty ())
                        refinedN. resize (nSets);

                // prepare a list of Morse sets which take part
                // in the ambiguous intersections
                std::vector<bool> refineM (nSets, false);
                std::vector<bool> refineN (nSets, false);
                for (int n = 0; n < nSets + nSets; ++ n)
                {
                        // prepare variables for M -> N or N -> M
                        int index = (n < nSets) ? n : (n - nSets);
                        std::vector<std::vector<int> > &AtoB =
                                (n < nSets) ? MtoN : NtoM;
                        if (AtoB [index]. size () <= 1)
                                continue;
                        std::vector<bool> &refineA =
                                (n < nSets) ? refineM : refineN;
                        std::vector<bool> &refineB =
                                (n < nSets) ? refineN : refineM;

                        // mark the numbers of sets which should be refined
                        refineA [index] = true;
                        for (unsigned i = 0; i < AtoB [index]. size (); ++ i)
                                refineB [AtoB [index] [i]] = true;
                }

                // refine all the Morse sets for which ambiguity was found
                for (int n = 0; n < nSets + nSets; ++ n)
                {
                        // prepare variables for M -> N or N -> M
                        int index = (n < nSets) ? n : (n - nSets);
                        std::vector<bool> &refineA =
                                (n < nSets) ? refineM : refineN;
                        if (!refineA [index])
                                continue;
                        std::vector<spcCubes> &refinedA =
                                (n < nSets) ? refinedM : refinedN;
                        if (!refinedA [index]. empty ())
                                continue;
                        const theMorseDecompositionType &morseDecA =
                                (n < nSets) ? morseDecM : morseDecN;
                        if (morseDecA [index]. size () > maxRefineSize0)
                                continue;
                        const theCubMapType &theCubMap1 =
                                (n < nSets) ? theCubMap1M : theCubMap1N;
                        sbug << "- match subdiv " <<
                                morseDecA [index]. size () << " -> ";
                        subdivideCubes (morseDecA [index], refinedA [index]);
                        sbug << refinedA [index]. size () << ". ";
                        int result = invariantPart (refinedA [index],
                                theCubMap1);
                        if (result < 0)
                        {
                                spcCubes emptySet;
                                refinedA [index] = emptySet;
                        }
                }

                // compare the refined Morse sets to avoid the ambiguity
                for (int n = 0; n < nSets + nSets; ++ n)
                {
                        // prepare variables for M -> N or N -> M
                        int indA = (n < nSets) ? n : (n - nSets);
                        std::vector<std::vector<int> > &AtoB =
                                (n < nSets) ? MtoN : NtoM;
                        if (AtoB. size () <= 1)
                                continue;
                        std::vector<std::vector<int> > &BtoA =
                                (n < nSets) ? NtoM : MtoN;
                        const theMorseDecompositionType &morseDecA =
                                (n < nSets) ? morseDecM : morseDecN;
                        const theMorseDecompositionType &morseDecB =
                                (n < nSets) ? morseDecN : morseDecM;
                        std::vector<spcCubes> &refinedA =
                                (n < nSets) ? refinedM : refinedN;
                        std::vector<spcCubes> &refinedB =
                                (n < nSets) ? refinedN : refinedM;

                        // modify the relation A -> B
                        std::vector<int>::iterator begin =
                                AtoB [indA]. begin ();
                        for (unsigned i = 0; i < AtoB [indA]. size (); ++ i)
                        {
                                // determine the index in B
                                int indB = AtoB [indA] [i];

                                // check which Morse sets are refined
                                bool refA = !refinedA [indA]. empty ();
                                bool refB = !refinedB [indB]. empty ();

                                // verify the intersection of the sets
                                if (refA && refB && intersect
                                        (refinedA [indA], refinedB [indB]))
                                {
                                        continue;
                                }
                                else if (refA && !refB && intersectSub
                                        (refinedA [indA], morseDecB [indB]))
                                {
                                        continue;
                                }
                                else if (!refA && refB && intersectSub
                                        (refinedB [indB], morseDecA [indA]))
                                {
                                        continue;
                                }
                                else if (!refA && !refB)
                                        continue;

                                // remove the link A <-> B
                                BtoA [indB]. erase (std::find
                                        (BtoA [indB]. begin (),
                                        BtoA [indB]. end (), indA));
                                AtoB [indA]. erase (begin + i);
                                -- i;
                        }
                }

                // check the bijection of intersections for the second time
                for (int i = 0; i < nSets; ++ i)
                {
                        if (MtoN [i]. empty () || NtoM [i]. empty ())
                        {
                                sbug << "WARNING: The intersection became "
                                        "empty!!!\n";
                                return false;
                        }
                        if ((MtoN [i]. size () > 1) ||
                                (NtoM [i]. size () > 1))
                        {
                                sbug << "Note: No bijection in matching, "
                                        "even after refinements.\n";
                                return false;
                        }
                }
        }

        // check if the Morse sets with wrong indices match each other
        if (wrongIndM. size () != wrongIndN. size ())
                return false;
        int wrongIndSize = wrongIndM. size ();
        for (int i = 0; i < wrongIndSize; ++ i)
        {
                int number = MtoN [wrongIndM [i]] [0];
                bool found = false;
                for (int j = 0; j < wrongIndSize; ++ j)
                {
                        if (wrongIndN [j] != number)
                                continue;
                        found = true;
                        break;
                }
                if (!found)
                        return false;
        }

        // compare the graphs of the Morse decompositions
        if (matchOrdering)
        {
                // create the arrays with connections for both Morse decomps
                chomp::homology::flatMatrix<char> orderM (nSets);
                chomp::homology::flatMatrix<char> orderN (nSets);
                for (int i = 0; i < nSets; ++ i)
                {
                        for (int j = 0; j < nSets; ++ j)
                        {
                                orderM [i] [j] = morseDecM.
                                        connected (i, j) ? 1 : 0;
                                orderN [i] [j] = morseDecN. connected
                                        (MtoN [i] [0], MtoN [j] [0]) ? 1 : 0;
                        }
                }

                // compute the transitive closure of the ordering
                chomp::homology::transitiveClosure (orderM, nSets);
                chomp::homology::transitiveClosure (orderN, nSets);

                // compare the arrays if they are the same indeed
                const char *memoryM = orderM. memory ();
                const char *memoryN = orderN. memory ();
                int memSize = nSets * nSets;
                for (int i = 0; i < memSize; ++ i)
                {
                        if (*(memoryM ++) != *(memoryN ++))
                        {
                                sbug << "Note: Different Morse ordering.\n";
                                return false;
                        }
                }
        //      sbug << "Note: Morse order matched.\n";
        }

        return true;
} /* matchMorseSets */


// --------------------------------------------------
// ---------- matching Morse decompositions ---------
// --------------------------------------------------

/// Matches Morse decompositions corresponding to parameter cubes
/// in the entire range. The set of cubes is used to determine
/// the numbers in the indexing of Morse decompositions.
inline void matchMorseDecompositions (const parCubes &paramBoxes,
        const std::vector<theMorseDecompositionType *> &morseDec,
//      const std::vector<theMapType *> &theMap,
//      const double *offset, const double *width, int subdivDepth,
//      const std::vector<theCubMapType *> &theCubMap0,
        const std::vector<theCubMapType *> &theCubMap1,
        const std::vector<std::vector<int> > &wrongIndices,
        const parCoord *parLeft, const parCoord *parRight,
        std::vector<std::vector<parCube> > &matchClasses)
{
        using chomp::homology::sbug;

        // create an array for matching the set of Morse decompositions
        parCoord sizes [paramDim];
        for (int i = 0; i < paramDim; ++ i)
                sizes [i] = parRight [i] - parLeft [i];
        MatchArray<parCoord,int> matchArray (paramDim, sizes, parLeft);

        // determine the total number of boxes to process
        int nBoxes = paramBoxes. size ();

        // prepare an array of refined Morse sets
        std::vector<std::vector<spcCubes> > refined (nBoxes);

        // process all the neighbors of each Morse set while avoiding
        // repetitive verification for the same pairs
        for (int boxNumber = 0; boxNumber < nBoxes; ++ boxNumber)
        {
                // retrieve the coordinates of the box
                parCoord coord [paramDim];
                paramBoxes [boxNumber]. coord (coord);

                // prepare the coordinates of the corners of the neighborhood
                parCoord coordLeft [paramDim];
                parCoord coordRight [paramDim];
                for (int i = 0; i < paramDim; ++ i)
                {
                        coordLeft [i] = coord [i] - 1;
                        if (coordLeft [i] < parLeft [i])
                                coordLeft [i] = parLeft [i];
                        coordRight [i] = coord [i] + 2;
                        if (coordRight [i] > parRight [i])
                                coordRight [i] = parRight [i];
                }

                // iterate the neighbors of the cube (including itself)
                parRect neighbors (coordLeft, coordRight, paramDim);
                const parCoord *neighborCoord;
                while ((neighborCoord = neighbors. get ()) != 0)
                {
                        // make sure that the neighbor precedes
                        // the current box in the lexicographical order
                        int compare = 0;
                        for (int i = 0; i < paramDim; ++ i)
                        {
                                if (neighborCoord [i] < coord [i])
                                {
                                        compare = -1;
                                        break;
                                }
                                else if (neighborCoord [i] > coord [i])
                                {
                                        compare = 1;
                                        break;
                                }
                        }
                        if (compare >= 0)
                                continue;

                        // if the boxes have already been joined
                        // then skip the verification of continuation
                        if (matchArray. joined (coord, neighborCoord))
                                continue;

                        // check matching between the corresponding
                        // Morse decompositions
                        int neighborNumber = paramBoxes. getnumber
                                (parCube (neighborCoord, paramDim));
                        if (neighborNumber < 0)
                                throw "Matching a non-existent neighbor.";
                        sbug << "Matching " << parCube (coord, paramDim) <<
                                " & " << parCube (neighborCoord, paramDim) <<
                                "...\n";
                        bool matchResult = matchMorseSets
                                (*(morseDec [boxNumber]),
                                wrongIndices [boxNumber],
                                *(theCubMap1 [boxNumber]),
                                refined [boxNumber],
                                *(morseDec [neighborNumber]),
                                wrongIndices [neighborNumber],
                                *(theCubMap1 [neighborNumber]),
                                refined [neighborNumber],
                                compareMorseOrdering);
                        if (matchResult == true)
                        {
                                sbug << "Success. :)\n";
                                matchArray. join (coord, neighborCoord);
                        }
                        else
                        {
                                sbug << "Failure. :(\n";
                        }
                }
        }

        // retrieve classes larger than 1 from the array
        matchArray. saveClasses (matchClasses);

        return;
} /* matchMorseDecompositions */


#endif // _CMGRAPHS_MATCHDEC_H_

---