29 #ifndef _CHAINCON_ATMODEL3_H_ 30 #define _CHAINCON_ATMODEL3_H_ 38 #include "chomp/system/config.h" 39 #include "chomp/system/textfile.h" 40 #include "chomp/system/timeused.h" 41 #include "chomp/struct/hashsets.h" 42 #include "chomp/struct/multitab.h" 43 #include "chomp/struct/integer.h" 44 #include "chomp/homology/chains.h" 55 #ifndef _product_is_identity_ 56 #define _product_is_identity_ 60 (
const chomp::homology::mmatrix<chomp::homology::integer> &A,
61 const chomp::homology::mmatrix<chomp::homology::integer> &G)
63 if (A. getncols () != G. getnrows ())
65 if (G. getncols () != A. getnrows ())
68 using chomp::homology::integer;
69 using chomp::homology::chain;
73 int_t nCols = A. getncols ();
74 for (int_t nCol = 0; nCol < nCols; ++ nCol)
76 const chain<integer> &col (A. getcol (nCol));
78 int_t size = col. size ();
79 for (int_t pos = 0; pos < size; ++ pos)
81 img. add (G. getcol (col. num (pos)),
84 if ((img. size () != 1) ||
85 (img. num (0) != nCol) ||
86 (img. coef (0) != one))
93 #endif // _product_is_identity_ 101 template <
class CellT,
class CellArray1,
class CellArray2,
class CellRestrT>
106 const CellRestrT &restr)
108 using chomp::homology::chaincomplex;
109 using chomp::homology::chain;
110 using chomp::homology::mmatrix;
111 using chomp::homology::timeused;
112 using chomp::homology::sout;
114 using chomp::homology::sbug;
115 using chomp::homology::integer;
116 using chomp::homology::hashedset;
117 using chomp::homology::multitable;
128 sbug <<
"[AT Model Algorithm 3.]\n";
129 sout <<
"Constructing the chain complex of K...\n";
134 integer::initialize (2);
141 const int_t KSize = K. size ();
142 const int dim = K [KSize - 1]. dim ();
143 const int trace_generators = 1;
144 const int trace_bases = 1;
145 chaincomplex<integer> cx (dim, trace_generators, trace_bases);
148 multitable<hashedset<CellT> > cells;
149 for (int_t i = 0; i < KSize; ++ i)
152 const CellT &c = K [i];
158 throw "Reduced homology not supported " 159 "in combinatorial Algorithm 3.\n" 160 "Please, use Algorithm 4 with Z_2 instead.";
164 int_t cNum = cells [d]. size ();
168 int_t bLen = c. boundaryLength ();
169 for (int_t j = 0; j < bLen; ++ j)
178 int_t bNum = cells [d - 1]. getnumber (bCell);
185 cx. add (d, bNum, cNum, one);
190 for (
int d = 0; d <= dim; ++ d)
192 cx. def_gen (d, cells [d]. size ());
196 sout <<
"The chain complex created in " << constrTime <<
".\n";
199 sout <<
"Computing the SNF of the boundary operator...\n";
205 cx. simple_form (level, quiet);
208 sout <<
"The SNF computed in " << snfTime <<
".\n";
215 sout <<
"Change of basis check: ";
216 for (
int q = 0; q <= dim; ++ q)
218 const mmatrix<integer> &G (cx. gethomgen (q));
219 const mmatrix<integer> &A (cx. getchgbasis (q));
223 sout <<
"Failed at level " << q <<
".\n";
224 throw "Wrong change of basis detected.";
227 sout <<
"Completed successfully in " << verifTime <<
".\n";
231 sout <<
"Composing an algebraic topological model from the SNF...\n";
240 multitable <mmatrix<integer> > psi;
243 const mmatrix<integer> zeroMatrix;
244 for (
int q = 0; q <= dim; ++ q)
247 int_t gen_number = 0;
251 const mmatrix<integer> &Dq = cx. getboundary (q);
253 const mmatrix<integer> &Dq1 = (q < dim) ?
254 cx. getboundary (q + 1) : zeroMatrix;
256 const mmatrix<integer> &Gq = cx. gethomgen (q);
258 const mmatrix<integer> &Aq = cx. getchgbasis (q);
263 psi [q - 1]. define (cells [q]. size (),
264 cells [q - 1]. size ());
268 int ncells = cells [q]. size ();
269 for (
int i = 0; i < ncells; ++ i)
272 const chain<integer> &column = Dq. getcol (i);
275 if (!column. empty ())
277 int b = column. num (0);
279 throw "Negative psi index.";
280 psi [q - 1]. add (i, b, one);
284 else if ((q == dim) || (Dq1. getrow (i). empty ()))
291 const chain<integer> &genChain =
295 const CellT &repr = cells [q] [gen_number];
300 int_t genSize = genChain. size ();
301 for (int_t j = 0; j < genSize; ++ j)
303 int_t n = genChain. num (j);
304 const CellT &genCell = cells [q] [n];
305 incl. add (repr, genCell);
309 const chain<integer> &piRow =
311 int_t piRowSize = piRow. size ();
312 for (int_t j = 0; j < piRowSize; ++ j)
314 int_t num = piRow. num (j);
315 const CellT &cell = cells [q] [num];
316 pi. add (cell, repr);
324 sout <<
"Transferring the computed maps to the original basis...\n";
325 for (
int q = 1; q <= dim; ++ q)
328 const mmatrix<integer> &M_A = cx. getchgbasis (q - 1);
329 const mmatrix<integer> &M_psi = psi [q - 1];
330 const mmatrix<integer> &M_A1 = cx. gethomgen (q);
333 int_t sizeA = M_A. getncols ();
334 int_t sizePsi = M_psi. getncols ();
335 int_t sizeA1 = M_A1. getncols ();
338 for (int_t i = 0; i < sizeA; ++ i)
341 const chain<integer> &column = M_A. getcol (i);
345 int_t sizeC = column. size ();
346 for (int_t j = 0; j < sizeC; ++ j)
348 int_t n = column. num (j);
351 psiC. add (M_psi. getcol (n), one);
355 chain<integer> phiColumn;
356 int_t sizePsiC = psiC. size ();
357 for (
int j = 0; j < sizePsiC; ++ j)
359 int_t n = psiC. num (j);
362 phiColumn. add (M_A1. getcol (n), one);
366 int_t sizePhiColumn = phiColumn. size ();
367 for (
int j = 0; j < sizePhiColumn; ++ j)
369 int_t n = phiColumn. num (j);
370 phi. add (cells [q - 1] [i],
377 sout <<
"AT model composed from the SNF in " << atModTime <<
".\n";
380 sout <<
"AT model computed in " << compTime <<
".\n";
386 #endif // _CHAINCON_ATMODEL3_H_ void algTopModel3(const CellArray1 &K, CellArray2 &H, tCombLinMap< CellT, CellT > &pi, tCombLinMap< CellT, CellT > &incl, tCombLinMap< CellT, CellT > &phi, const CellRestrT &restr)
Computes an algebraic topological model for the given filtered finite cell complex "K"...
A combinatorial chain, that is, a chain with Z_2 coefficients.
A combinatorial linear map (for coefficients in Z_2).
A combinatorial linear map.
bool product_is_identity(const chomp::homology::mmatrix< chomp::homology::integer > &A, const chomp::homology::mmatrix< chomp::homology::integer > &G)
Verifies if the product of the given two matrices is an identity matrix.