35#ifndef _CMGRAPHS_MATCHDEC_H_
36#define _CMGRAPHS_MATCHDEC_H_
45#include "chomp/system/config.h"
46#include "chomp/system/textfile.h"
47#include "chomp/struct/digraph.h"
48#include "chomp/struct/flatmatr.h"
49#include "chomp/struct/autoarray.h"
65template <
class typeCubes>
66inline bool intersect (
const typeCubes &X,
const typeCubes &Y)
69 bool XY = (X. size () < Y. size ());
70 const typeCubes &smaller = XY ? X : Y;
71 const typeCubes &larger = XY? Y : X;
74 int smallerSize = smaller. size ();
75 for (
int i = 0; i < smallerSize; ++ i)
77 if (larger. check (smaller [i]))
88template <
class typeCubes>
89inline bool intersectSub (
const typeCubes &fine,
const typeCubes &coarse)
91 using chomp::homology::auto_array;
93 typedef typename typeCubes::value_type typeCube;
94 typedef typename typeCube::CoordType typeCoord;
97 int size = fine. size ();
102 int dim = fine [0]. dim ();
103 auto_array<typeCoord> coordPtr (
new typeCoord [dim]);
104 typeCoord *coord = coordPtr. get ();
108 for (
int i = 0; i < size; ++ i)
110 fine [i]. coord (coord);
111 for (
int j = 0; j < dim; ++ j)
113 if (coarse. check (typeCube (coord, dim)))
129 const std::vector<int> &wrongIndM,
const theCubMapType &theCubMap1M,
130 std::vector<spcCubes> &refinedM,
132 const std::vector<int> &wrongIndN,
const theCubMapType &theCubMap1N,
133 std::vector<spcCubes> &refinedN,
136 using chomp::homology::sbug;
139 int nSets = morseDecM. count ();
140 if (nSets != morseDecN. count ())
147 std::vector<std::vector<int> > MtoN (nSets);
148 std::vector<std::vector<int> > NtoM (nSets);
149 for (
int m = 0; m < nSets; ++ m)
151 for (
int n = 0; n < nSets; ++ n)
153 if (!
intersect (morseDecM [m], morseDecN [n]))
155 MtoN [m]. push_back (n);
156 NtoM [n]. push_back (m);
161 bool bijection =
true;
162 for (
int i = 0; i < nSets; ++ i)
164 if (MtoN [i]. empty () || NtoM [i]. empty ())
166 if ((MtoN [i]. size () > 1) || (NtoM [i]. size () > 1))
174 sbug <<
"Note: No bijection in matching. Refining...\n";
177 if (refinedM. empty ())
178 refinedM. resize (nSets);
179 if (refinedN. empty ())
180 refinedN. resize (nSets);
184 std::vector<bool> refineM (nSets,
false);
185 std::vector<bool> refineN (nSets,
false);
186 for (
int n = 0; n < nSets + nSets; ++ n)
189 int index = (n < nSets) ? n : (n - nSets);
190 std::vector<std::vector<int> > &AtoB =
191 (n < nSets) ? MtoN : NtoM;
192 if (AtoB [index]. size () <= 1)
194 std::vector<bool> &refineA =
195 (n < nSets) ? refineM : refineN;
196 std::vector<bool> &refineB =
197 (n < nSets) ? refineN : refineM;
200 refineA [index] =
true;
201 for (
unsigned i = 0; i < AtoB [index]. size (); ++ i)
202 refineB [AtoB [index] [i]] =
true;
206 for (
int n = 0; n < nSets + nSets; ++ n)
209 int index = (n < nSets) ? n : (n - nSets);
210 std::vector<bool> &refineA =
211 (n < nSets) ? refineM : refineN;
212 if (!refineA [index])
214 std::vector<spcCubes> &refinedA =
215 (n < nSets) ? refinedM : refinedN;
216 if (!refinedA [index]. empty ())
219 (n < nSets) ? morseDecM : morseDecN;
223 (n < nSets) ? theCubMap1M : theCubMap1N;
224 sbug <<
"- match subdiv " <<
225 morseDecA [index]. size () <<
" -> ";
227 sbug << refinedA [index]. size () <<
". ";
233 refinedA [index] = emptySet;
238 for (
int n = 0; n < nSets + nSets; ++ n)
241 int indA = (n < nSets) ? n : (n - nSets);
242 std::vector<std::vector<int> > &AtoB =
243 (n < nSets) ? MtoN : NtoM;
244 if (AtoB. size () <= 1)
246 std::vector<std::vector<int> > &BtoA =
247 (n < nSets) ? NtoM : MtoN;
249 (n < nSets) ? morseDecM : morseDecN;
251 (n < nSets) ? morseDecN : morseDecM;
252 std::vector<spcCubes> &refinedA =
253 (n < nSets) ? refinedM : refinedN;
254 std::vector<spcCubes> &refinedB =
255 (n < nSets) ? refinedN : refinedM;
258 std::vector<int>::iterator begin =
259 AtoB [indA]. begin ();
260 for (
unsigned i = 0; i < AtoB [indA]. size (); ++ i)
263 int indB = AtoB [indA] [i];
266 bool refA = !refinedA [indA]. empty ();
267 bool refB = !refinedB [indB]. empty ();
271 (refinedA [indA], refinedB [indB]))
276 (refinedA [indA], morseDecB [indB]))
281 (refinedB [indB], morseDecA [indA]))
285 else if (!refA && !refB)
289 BtoA [indB]. erase (std::find
290 (BtoA [indB]. begin (),
291 BtoA [indB]. end (), indA));
292 AtoB [indA]. erase (begin + i);
298 for (
int i = 0; i < nSets; ++ i)
300 if (MtoN [i]. empty () || NtoM [i]. empty ())
302 sbug <<
"WARNING: The intersection became "
306 if ((MtoN [i]. size () > 1) ||
307 (NtoM [i]. size () > 1))
309 sbug <<
"Note: No bijection in matching, "
310 "even after refinements.\n";
317 if (wrongIndM. size () != wrongIndN. size ())
319 int wrongIndSize = wrongIndM. size ();
320 for (
int i = 0; i < wrongIndSize; ++ i)
322 int number = MtoN [wrongIndM [i]] [0];
324 for (
int j = 0; j < wrongIndSize; ++ j)
326 if (wrongIndN [j] != number)
339 chomp::homology::flatMatrix<char> orderM (nSets);
340 chomp::homology::flatMatrix<char> orderN (nSets);
341 for (
int i = 0; i < nSets; ++ i)
343 for (
int j = 0; j < nSets; ++ j)
345 orderM [i] [j] = morseDecM.
346 connected (i, j) ? 1 : 0;
347 orderN [i] [j] = morseDecN. connected
348 (MtoN [i] [0], MtoN [j] [0]) ? 1 : 0;
353 chomp::homology::transitiveClosure (orderM, nSets);
354 chomp::homology::transitiveClosure (orderN, nSets);
357 const char *memoryM = orderM. memory ();
358 const char *memoryN = orderN. memory ();
359 int memSize = nSets * nSets;
360 for (
int i = 0; i < memSize; ++ i)
362 if (*(memoryM ++) != *(memoryN ++))
364 sbug <<
"Note: Different Morse ordering.\n";
383 const std::vector<theMorseDecompositionType *> &morseDec,
387 const std::vector<theCubMapType *> &theCubMap1,
388 const std::vector<std::vector<int> > &wrongIndices,
390 std::vector<std::vector<parCube> > &matchClasses)
392 using chomp::homology::sbug;
397 sizes [i] = parRight [i] - parLeft [i];
401 int nBoxes = paramBoxes. size ();
404 std::vector<std::vector<spcCubes> > refined (nBoxes);
408 for (
int boxNumber = 0; boxNumber < nBoxes; ++ boxNumber)
412 paramBoxes [boxNumber]. coord (coord);
419 coordLeft [i] = coord [i] - 1;
420 if (coordLeft [i] < parLeft [i])
421 coordLeft [i] = parLeft [i];
422 coordRight [i] = coord [i] + 2;
423 if (coordRight [i] > parRight [i])
424 coordRight [i] = parRight [i];
430 while ((neighborCoord = neighbors. get ()) != 0)
437 if (neighborCoord [i] < coord [i])
442 else if (neighborCoord [i] > coord [i])
453 if (matchArray. joined (coord, neighborCoord))
458 int neighborNumber = paramBoxes. getnumber
460 if (neighborNumber < 0)
461 throw "Matching a non-existent neighbor.";
466 (*(morseDec [boxNumber]),
467 wrongIndices [boxNumber],
468 *(theCubMap1 [boxNumber]),
470 *(morseDec [neighborNumber]),
471 wrongIndices [neighborNumber],
472 *(theCubMap1 [neighborNumber]),
473 refined [neighborNumber],
475 if (matchResult ==
true)
477 sbug <<
"Success. :)\n";
478 matchArray. join (coord, neighborCoord);
482 sbug <<
"Failure. :(\n";
488 matchArray. saveClasses (matchClasses);
A generic map computation routine that computes a rigorous cubical multivalued map based on a functio...
The Morse decoposition class.
Choice of configuration settings.
void invariantPart(typeCubes &X, const chomp::homology::diGraph<> &g, chomp::homology::diGraph<> *gInv=0)
Computes X := Inv (X) using the algorithm by Bill Kalies and Hyunju Ban.
Array for matching Morse decompositions.
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)
Verifies if the two given Morse decompositions are equivalent by matching Morse sets which intersect ...
bool intersect(const typeCubes &X, const typeCubes &Y)
Verifies whether two sets of cubes intersect each other.
bool intersectSub(const typeCubes &fine, const typeCubes &coarse)
Verifies whether two sets of cubes intersect each other.
void matchMorseDecompositions(const parCubes ¶mBoxes, const std::vector< theMorseDecompositionType * > &morseDec, const std::vector< theCubMapType * > &theCubMap1, const std::vector< std::vector< int > > &wrongIndices, const parCoord *parLeft, const parCoord *parRight, std::vector< std::vector< parCube > > &matchClasses)
Matches Morse decompositions corresponding to parameter cubes in the entire range.
const int paramDim
The dimension of the parameter space to iterate.
const int maxRefineSize0
The maximal allowed size of a set of cubes in the phase space which can be refined at the initial sub...
const bool compareMorseOrdering
Should the ordering between the Morse sets be taken into consideration while determining whether two ...
Customizable data types for the Conley-Morse graphs computation program.
Data types for the dynamical systems data structures.
chomp::homology::tCubeFix< paramDim, parCoord > parCube
The type of a cube in the set of parameters.
short int parCoord
The type of coordinates of cubes in the set of parameters.
chomp::homology::tRectangle< parCoord > parRect
The type of a rectangle used to iterate sets of cubes of parameters.
chomp::homology::hashedset< parCube > parCubes
The type of a set of cubes in the set of parameters.
chomp::homology::hashedset< spcCube > spcCubes
The type of a set of cubes in the phase space.
void subdivideCubes(const typeCubes &X, typeCubes &result)
Subdivides one set of cubes to another set of cubes with respect to twice finer grid.