The Original CHomP Software
gcomplex.h
Go to the documentation of this file.
1
3
16
17// Copyright (C) 1997-2020 by Pawel Pilarczyk.
18//
19// This file is part of my research software package. This is free software:
20// you can redistribute it and/or modify it under the terms of the GNU
21// General Public License as published by the Free Software Foundation,
22// either version 3 of the License, or (at your option) any later version.
23//
24// This software is distributed in the hope that it will be useful,
25// but WITHOUT ANY WARRANTY; without even the implied warranty of
26// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27// GNU General Public License for more details.
28//
29// You should have received a copy of the GNU General Public License
30// along with this software; see the file "license.txt". If not,
31// please, see <https://www.gnu.org/licenses/>.
32
33// Started in January 2002. Last revision: January 23, 2010.
34
35
36#ifndef _CHOMP_HOMOLOGY_GCOMPLEX_H_
37#define _CHOMP_HOMOLOGY_GCOMPLEX_H_
38
39#include "chomp/system/config.h"
40#include "chomp/system/textfile.h"
41#include "chomp/struct/integer.h"
42#include "chomp/struct/hashsets.h"
43#include "chomp/homology/chains.h"
44#include "chomp/homology/gcomplex.h"
45
46#include <iostream>
47#include <fstream>
48#include <iomanip>
49#include <cstdlib>
50
51namespace chomp {
52namespace homology {
53
54
55// class templates defined within this header file (in this order):
56template <class cell, class euclidom>
57class gcomplex;
58template <class cell, class euclidom, class element>
59class mvcellmap;
60
61
62// --------------------------------------------------
63// --------------- geometric complex ----------------
64// --------------------------------------------------
65
66template <class cell>
67int boundarycoef (const cell &, int i);
68
69template <class cell>
70int boundarylength (const cell &);
71
72template <class cell>
73cell boundarycell (const cell &, int i);
74
83template <class cell, class euclidom>
84class gcomplex
85{
86public:
88 gcomplex ();
89
91 gcomplex (const gcomplex<cell,euclidom> &c);
92
94 ~gcomplex ();
95
97 gcomplex &operator = (const gcomplex<cell,euclidom> &c);
98
101 int dim () const;
102
104 int_t size () const;
105
107 bool empty () const;
108
110 const hashedset<cell> &operator [] (int d) const;
111
113 const hashedset<cell> &getcob (const cell &c) const;
114
117 const hashedset<cell> &getcob (const cell &c, int d) const;
118
120 gcomplex<cell,euclidom> &add (const cell &c);
121
124 gcomplex<cell,euclidom> &add (const hashedset<cell> &c, int d);
125
127 gcomplex<cell,euclidom> &add (const hashedset<cell> &c);
128
130 gcomplex<cell,euclidom> &add (const gcomplex<cell,euclidom> &c);
131
133 gcomplex<cell,euclidom> &remove (const cell &c);
134
137 gcomplex<cell,euclidom> &remove (const hashedset<cell> &c, int d);
138
140 gcomplex<cell,euclidom> &remove (const hashedset<cell> &c);
141
143 gcomplex<cell,euclidom> &remove (const gcomplex<cell,euclidom> &c);
144
146 gcomplex<cell,euclidom> &removeall (int d);
147
150 bool check (const cell &c) const;
151
155 int_t addboundaries (int d, bool addcob = false);
156
160 int_t addboundaries (bool addcob = false);
161
168 int_t addboundaries (int d, gcomplex<cell,euclidom> &notthese,
169 bool keepused = false);
170
177 int_t addboundaries (gcomplex<cell,euclidom> &notthese,
178 bool keepused = false);
179
181 int_t addboundaries (int d, chaincomplex<euclidom> &c);
182
184 int_t addboundaries (chaincomplex<euclidom> &c);
185
187 int_t addboundaries (int d, chaincomplex<euclidom> &c,
188 gcomplex<cell,euclidom> &notthese,
189 bool keepused = false);
190
192 int_t addboundaries (chaincomplex<euclidom> &c,
193 gcomplex<cell,euclidom> &notthese,
194 bool keepused = false);
195
198 int_t addboundaries (int d, chaincomplex<euclidom> *c,
199 gcomplex<cell,euclidom> *notthese,
200 bool dontadd, bool keepused, bool addcob);
201
204 int_t addboundaries (chaincomplex<euclidom> *c,
205 gcomplex<cell,euclidom> *notthese,
206 bool dontadd, bool keepused, bool addcob);
207
213 int_t collapse (int d,
214 gcomplex<cell,euclidom> &other,
215 const gcomplex<cell,euclidom> &keep,
216 bool addbd, bool addcob, bool disjoint,
217 bool quiet = false);
218
226 int_t collapse (gcomplex<cell,euclidom> &other,
227 gcomplex<cell,euclidom> &keep,
228 bool addbd, bool addcob, bool disjoint,
229 const int *level = NULL, bool quiet = false);
230
231private:
233 hashedset<cell> **tab;
234
236 mvmap<cell,cell> **cob;
237
239 int n;
240
242 void increasedimension (int d);
243
245 void decreasedimension ();
246
247}; /* class gcomplex */
248
249// --------------------------------------------------
250
251template <class cell, class euclidom>
252inline gcomplex<cell,euclidom>::gcomplex (): tab (NULL), cob (NULL), n (0)
253{
254 return;
255} /* gcomplex<cell,euclidom>::gcomplex */
256
257template <class cell, class euclidom>
258gcomplex<cell,euclidom>::gcomplex (const gcomplex<cell,euclidom> &c)
259{
260 n = c. n;
261 if (n > 0)
262 {
263 tab = new hashedset<cell> *[n];
264 cob = new mvmap<cell,cell> *[n];
265 if (!tab || !cob)
266 throw "Cannot copy a geometric complex.";
267 }
268 else
269 {
270 tab = NULL;
271 cob = NULL;
272 }
273 for (int i = 0; i < n; ++ i)
274 {
275 tab [i] = new hashedset<cell> (*(c. tab [i]));
276 cob [i] = new mvmap<cell,cell> (*(c. cob [i]));
277 if (!tab [i] || !cob [i])
278 throw "Cannot copy part of a geometric complex.";
279 }
280 return;
281} /* gcomplex<cell,euclidom>::gcomplex */
282
283template <class cell, class euclidom>
284gcomplex<cell,euclidom> &gcomplex<cell,euclidom>::operator =
285 (const gcomplex<cell,euclidom> &c)
286{
287 // if the sizes of the tables are not the same, re-allocate the table
288 // and use the copying constructor to create an identical table
289 if (n != c. n)
290 {
291 if (tab)
292 {
293 for (int i = 0; i < n; ++ i)
294 delete tab [i];
295 delete [] tab;
296 }
297 if (cob)
298 {
299 for (int i = 0; i < n; ++ i)
300 delete cob [i];
301 delete [] cob;
302 }
303 n = c. n;
304 if (n > 0)
305 {
306 tab = new hashedset<cell> *[n];
307 cob = new mvmap<cell,cell> *[n];
308 if (!tab || !cob)
309 throw "Cannot copy a geometric complex.";
310 }
311 else
312 {
313 tab = NULL;
314 cob = NULL;
315 }
316 for (int i = 0; i < n; ++ i)
317 {
318 tab [i] = new hashedset<cell> (*(c. tab [i]));
319 cob [i] = new mvmap<cell,cell> (*(c. cob [i]));
320 if (!tab [i] || !cob [i])
321 throw "Cannot copy part of a geom. complex.";
322 }
323 }
324 // otherwise copy the source table to this complex
325 else
326 {
327 for (int i = 0; i < n; ++ i)
328 {
329 *(tab [i]) = *(c. tab [i]);
330 *(cob [i]) = *(c. cob [i]);
331 }
332 }
333 return *this;
334} /* gcomplex<cell,euclidom>::operator = */
335
336template <class cell, class euclidom>
337gcomplex<cell,euclidom>::~gcomplex ()
338{
339 for (int i = 0; i < n; ++ i)
340 {
341 if (tab [i])
342 delete tab [i];
343 if (cob [i])
344 delete cob [i];
345 }
346 if (tab)
347 delete [] tab;
348 if (cob)
349 delete [] cob;
350 return;
351} /* gcomplex<cell,euclidom>::~gcomplex */
352
353template <class cell, class euclidom>
354inline int gcomplex<cell,euclidom>::dim () const
355{
356 return n - 1;
357} /* gcomplex<cell,euclidom>::dim */
358
359template <class cell, class euclidom>
360int_t gcomplex<cell,euclidom>::size () const
361{
362 int_t count = 0;
363 for (int i = 0; i < n; ++ i)
364 count += tab [i] -> size ();
365 return count;
366} /* gcomplex<cell,euclidom>::size */
367
368template <class cell, class euclidom>
369bool gcomplex<cell,euclidom>::empty () const
370{
371 if (!n)
372 return true;
373 for (int i = 0; i < n; ++ i)
374 if (!(*(tab [i])). empty ())
375 return false;
376 return true;
377} /* gcomplex<cell,euclidom>::empty */
378
379template <class cell, class euclidom>
380inline const hashedset<cell> &gcomplex<cell,euclidom>::operator [] (int d)
381 const
382{
383 if ((d < 0) || (d >= n))
384 throw "Dimension out of range for retrieving cells.";
385 return *(tab [d]);
386} /* gcomplex<cell,euclidom>::operator [] */
387
388template <class cell, class euclidom>
389inline const hashedset<cell> &gcomplex<cell,euclidom>::getcob (const cell &c)
390 const
391{
392 return getcob (c, c. dim ());
393} /* gcomplex<cell,euclidom>::getcob */
394
395template <class cell, class euclidom>
396inline const hashedset<cell> &gcomplex<cell,euclidom>::getcob (const cell &c,
397 int d) const
398{
399 if ((d < 0) || (d >= n))
400 throw "Dimension out of range for coboundary.";
401 return (*(cob [d])) (c);
402} /* gcomplex<cell,euclidom>::getcob */
403
404template <class cell, class euclidom>
405inline gcomplex<cell,euclidom> &gcomplex<cell,euclidom>::add (const cell &c)
406{
407 // get the dimension of the cell
408 int d = c. dim ();
409 if (d < 0)
410 throw "Negative dimension of a cell to be added.";
411
412 // if the dimension of the cell is beyond the allocated table,
413 // then increase the table
414 if (n <= d)
415 increasedimension (d);
416
417 // add the cell to the suitable set of cells
418 tab [d] -> add (c);
419
420 return *this;
421} /* gcomplex<cell,euclidom>::add */
422
423template <class cell, class euclidom>
424inline gcomplex<cell,euclidom> &gcomplex<cell,euclidom>::add
425 (const hashedset<cell> &c, int d)
426{
427 // increase the dimension of the complex if necessary
428 if (d > n)
429 increasedimension (d);
430
431 // add the set of cells to the suitable set
432 tab [d] -> add (c);
433
434 return *this;
435} /* gcomplex<cell,euclidom>::add */
436
437template <class cell, class euclidom>
438gcomplex<cell,euclidom> &gcomplex<cell,euclidom>::add
439 (const hashedset<cell> &c)
440{
441 int_t size = c. size ();
442 for (int_t i = 0; i < size; ++ i)
443 add (c [i]);
444
445 return *this;
446} /* gcomplex<cell,euclidom>::add */
447
448template <class cell, class euclidom>
449gcomplex<cell,euclidom> &gcomplex<cell,euclidom>::add
450 (const gcomplex<cell,euclidom> &c)
451{
452 // increase the dimension of the complex if necessary
453 if (c. n > n)
454 increasedimension (c. dim ());
455
456 // add the sets of cells
457 for (int i = 0; i < c. n; ++ i)
458 {
459 tab [i] -> add (*(c. tab [i]));
460 /* const hashedset<cell> &cset = *(c. tab [i]);
461 for (int j = 0; j < cset. size (); ++ j)
462 {
463 const cell &thecell = cset [j];
464 tab [i] -> add (thecell);
465 if (cob && c. cob)
466 (*(cob [i])) [thecell] =
467 ((*(c. cob [i])) (thecell));
468 }
469 */
470 }
471
472 return *this;
473} /* gcomplex<cell,euclidom>::add */
474
475template <class cell, class euclidom>
476inline gcomplex<cell,euclidom> &gcomplex<cell,euclidom>::remove
477 (const cell &c)
478{
479 // get the dimension of the cell
480 int d = c. dim ();
481 if (d < 0)
482 throw "Negative dimension of a cell to be removed.";
483
484 // if the dimension of the cell is beyond the allocated table,
485 // then ignore it
486 if (n <= d)
487 return *this;
488
489 // remove the cell from the suitable set of cells
490 tab [d] -> remove (c);
491
492 // decrease the dimension of the complex if no cells remain
493 if ((d == n - 1) && tab [d] -> empty ())
494 decreasedimension ();
495
496 return *this;
497} /* gcomplex<cell,euclidom>::remove */
498
499template <class cell, class euclidom>
500gcomplex<cell,euclidom> &gcomplex<cell,euclidom>::remove
501 (const gcomplex<cell,euclidom> &c)
502{
503 // figure out the number of tables to work on
504 int m = c. n;
505 if (m > n)
506 m = n;
507
508 // remove the sets of cells
509 for (int i = 0; i < m; ++ i)
510 tab [i] -> remove (*(c. tab [i]));
511
512 // decrease the dimension of the complex if no highdim cells remain
513 decreasedimension ();
514
515 return *this;
516} /* gcomplex<cell,euclidom>::remove */
517
518template <class cell, class euclidom>
519inline gcomplex<cell,euclidom> &gcomplex<cell,euclidom>::remove
520 (const hashedset<cell> &c, int d)
521{
522 // remove the set of cells to the suitable set
523 tab [d] -> remove (c);
524
525 // decrease the dimension of the complex if no highdim cells remain
526 decreasedimension ();
527
528 return *this;
529} /* gcomplex<cell,euclidom>::remove */
530
531template <class cell, class euclidom>
532gcomplex<cell,euclidom> &gcomplex<cell,euclidom>::remove
533 (const hashedset<cell> &c)
534{
535 int_t size = c. size ();
536 for (int_t i = 0; i < size; ++ i)
537 remove (c [i]);
538
539 return *this;
540} /* gcomplex<cell,euclidom>::remove */
541
542template <class cell, class euclidom>
543gcomplex<cell,euclidom> &gcomplex<cell,euclidom>::removeall (int d)
544{
545 if (d == n - 1)
546 {
547 -- n;
548 delete tab [n];
549 delete cob [n];
550 decreasedimension ();
551 }
552 else if ((d >= 0) && (d < n))
553 {
554 delete tab [d];
555 tab [d] = new hashedset<cell>;
556 }
557
558 return *this;
559} /* gcomplex<cell,euclidom>::removeall */
560
561template <class cell, class euclidom>
562bool gcomplex<cell,euclidom>::check (const cell &c) const
563{
564 // get the dimension of the cell
565 int d = c. dim ();
566 if (d < 0)
567 throw "Negative dimension of a cell to be checked.";
568
569 // if the dimension of the cell is beyond the allocated table,
570 // then it is not there
571 if (n <= d)
572 return false;
573
574 // check for the existence of the cell in the suitable set of cells
575 return tab [d] -> check (c);
576} /* gcomplex<cell,euclidom>::check */
577
578template <class cell, class euclidom>
579int_t gcomplex<cell,euclidom>::addboundaries (int d,
580 chaincomplex<euclidom> *c, gcomplex<cell,euclidom> *notthese,
581 bool dontadd, bool keepused, bool addcob)
582{
583 // if the dimension is inappropriate, do nothing
584 if ((d <= 0) || (d >= n))
585 return 0;
586
587 // first add boundaries to the other cell complex
588 if (notthese && !dontadd)
589 notthese -> addboundaries (d);
590
591 int_t prevsize = tab [d - 1] -> size ();
592 hashedset<cell> &cset = *(tab [d]);
593 for (int_t i = 0; i < cset. size (); ++ i)
594 {
595 const cell &thecell = cset [i];
596 int len = boundarylength (thecell);
597 for (int j = 0; j < len; ++ j)
598 {
599 // take the j-th boundary cell
600 cell bcell = boundarycell (thecell, j);
601
602 // add it to the cell complex unless it is unwanted
603 if (!notthese || !notthese -> check (bcell))
604 {
605 tab [d - 1] -> add (bcell);
606 if (c)
607 {
608 int icoef = boundarycoef (thecell, j);
609 euclidom coef;
610 if (icoef < 0)
611 {
612 coef = -icoef;
613 coef = -coef;
614 }
615 else
616 coef = icoef;
617 c -> add (d, tab [d - 1] ->
618 getnumber (bcell), i, coef);
619 }
620 if (addcob)
621 (*(cob [d - 1])) [bcell].
622 add (thecell);
623 }
624 }
625 }
626
627 // clean the used level in 'notthese'
628 if (notthese && !keepused)
629 notthese -> removeall (d);
630
631 return tab [d - 1] -> size () - prevsize;
632} /* gcomplex<cell,euclidom>::addboundaries */
633
634template <class cell, class euclidom>
635int_t gcomplex<cell,euclidom>::addboundaries (chaincomplex<euclidom> *c,
636 gcomplex<cell,euclidom> *notthese, bool dontadd, bool keepused,
637 bool addcob)
638{
639 int_t countadded = 0;
640 for (int d = n - 1; d > 0; -- d)
641 countadded += addboundaries (d, c, notthese,
642 dontadd, keepused, addcob);
643 return countadded;
644} /* gcomplex<cell,euclidom>::addboundaries */
645
646template <class cell, class euclidom>
647inline int_t gcomplex<cell,euclidom>::addboundaries (int d, bool addcob)
648{
649 return addboundaries (d, NULL, NULL, false, false, addcob);
650} /* gcomplex<cell,euclidom>::addboundaries */
651
652template <class cell, class euclidom>
653inline int_t gcomplex<cell,euclidom>::addboundaries (bool addcob)
654{
655 return addboundaries (NULL, NULL, false, false, addcob);
656} /* gcomplex<cell,euclidom>::addboundaries */
657
658template <class cell, class euclidom>
659inline int_t gcomplex<cell,euclidom>::addboundaries (int d,
660 gcomplex<cell,euclidom> &notthese, bool keepused)
661{
662 return addboundaries (d, NULL, &notthese, false, keepused, false);
663} /* gcomplex<cell,euclidom>::addboundaries */
664
665template <class cell, class euclidom>
666int_t gcomplex<cell,euclidom>::addboundaries
667 (gcomplex<cell,euclidom> &notthese, bool keepused)
668{
669 return addboundaries (NULL, &notthese, false, keepused, false);
670} /* gcomplex<cell,euclidom>::addboundaries */
671
672template <class cell, class euclidom>
673inline int_t gcomplex<cell,euclidom>::addboundaries (int d,
674 chaincomplex<euclidom> &c)
675{
676 return addboundaries (d, &c, NULL, false, false, false);
677} /* gcomplex<cell,euclidom>::addboundaries */
678
679template <class cell, class euclidom>
680inline int_t gcomplex<cell,euclidom>::addboundaries (chaincomplex<euclidom> &c)
681{
682 return addboundaries (&c, NULL, false, false, false);
683} /* gcomplex<cell,euclidom>::addboundaries */
684
685template <class cell, class euclidom>
686inline int_t gcomplex<cell,euclidom>::addboundaries (int d,
687 chaincomplex<euclidom> &c, gcomplex<cell,euclidom> &notthese,
688 bool keepused)
689{
690 return addboundaries (d, &c, &notthese, false, keepused, false);
691} /* gcomplex<cell,euclidom>::addboundaries */
692
693template <class cell, class euclidom>
694inline int_t gcomplex<cell,euclidom>::addboundaries (chaincomplex<euclidom> &c,
695 gcomplex<cell,euclidom> &notthese, bool keepused)
696{
697 return addboundaries (&c, &notthese, false, keepused, false);
698} /* gcomplex<cell,euclidom>::addboundaries */
699
700// --------------------------------------------------
701
704template <class element>
705int_t findelem (const multitable<element> &tab, const element &e, int_t len)
706{
707 if (len <= 0)
708 return -1;
709
710 // prepare a random search starting point and step increase
711 int_t i = static_cast<int_t> (std::rand ()) % len;
712 if (i < 0)
713 i = static_cast<int_t> (1) - i;
714 int_t step = static_cast<int_t> (std::rand () + 17) % len;
715 if (step < 0)
716 step = static_cast<int_t> (1) - step;
717 if (step < 17)
718 step = 17;
719 if (step > len)
720 step = len >> 1;
721 if (step < 1)
722 step = 1;
723
724 // jump randomly in the table to find some element if possible
725 int_t jumping = len >> 1;
726 while (jumping --)
727 {
728 // if ((i < 0) || (i >= len))
729 // throw "Wrong random number.";
730 if (tab (i) == e)
731 return i;
732 // if ((i + 1 < len) && (tab (i + 1) == e))
733 // return i + 1;
734 if (jumping)
735 {
736 i += step;
737 if (i >= len)
738 i -= len;
739 }
740 }
741
742 // if not found, try finding the element thoroughly
743 for (int_t i = 0; i < len; ++ i)
744 {
745 if (tab (i) == e)
746 return i;
747 }
748 return -1;
749} /* findelem */
750
751template <class cell, class euclidom>
752int_t gcomplex<cell,euclidom>::collapse (int d,
753 gcomplex<cell,euclidom> &other, const gcomplex<cell,euclidom> &keep,
754 bool addbd, bool addcob, bool disjoint, bool quiet)
755{
756 if ((d <= 0) || (d > dim ()))
757 return 0;
758
759 // prepare references to the sets of cells of interest
760 hashedset<cell> &cset = *(tab [d]);
761 hashedset<cell> &bset = *(tab [d - 1]);
762 hashedset<cell> empty;
763 hashedset<cell> &cother = (other. dim () >= d) ?
764 *(other. tab [d]) : empty;
765 hashedset<cell> &bother = (other. dim () >= d - 1) ?
766 *(other. tab [d - 1]) : empty;
767 const hashedset<cell> &ckeep = (keep. dim () >= d) ?
768 *(keep. tab [d]) : empty;
769 const hashedset<cell> &bkeep = (keep. dim () >= d - 1) ?
770 *(keep. tab [d - 1]) : empty;
771
772 // show the currently processed dimension
773 if (!quiet)
774 {
775 if (d > 9)
776 scon << "#\b";
777 else
778 scon << d << '\b';
779 }
780
781 // go through all the cells from A and generate their boundaries
782 if (addbd)
783 {
784 if (!quiet)
785 {
786 sseq << "\"Adding boundaries of cells in A...\"\n";
787 sseq << "D 0\n";
788 }
789 for (int_t i = 0; i < cother. size (); ++ i)
790 {
791 // select the given cell in 'cother'
792 const cell &thecell = cother [i];
793
794 // detect the length of the boundary of this cell
795 int blen = boundarylength (thecell);
796
797 // add to 'bother' all the cells in this boundary
798 for (int j = 0; j < blen; ++ j)
799 bother. add (boundarycell (thecell, j));
800
801 // write the boundary cells to the sequential file
802 if (!quiet && (sseq. show || sseq. log))
803 {
804 for (int j = 0; j < blen; ++ j)
805 sseq << '2' << boundarycell
806 (thecell, j) << '\n';
807 }
808 }
809 }
810 if (!quiet)
811 sseq << "D 100\n";
812
813 // if the complexes should be disjoint, remove 'bother' from 'bset'
814 if (disjoint)
815 bset. remove (bother);
816
817 // prepare tables for coboundaries of cells in X and their lengths
818 multitable<int> coblen;
819 multitable<hashedset <cell> > coboundary;
820 if (!bset. empty ())
821 coblen. fill (0, bset. size ());
822
823 // go through the list of all the cells of dimension 'd' in the set,
824 // add their boundaries to 'bset' and create the coboundary links;
825 // note: these cells may appear in the 'other' complex
826 if (!quiet)
827 sseq << "\"Adding boundaries of cells of dimension " <<
828 d << "\"\nD 0\n";
829 bool maximal = (d == dim ());
830 for (int_t i = 0; i < cset. size (); ++ i)
831 {
832 // select the i-th d-dimensional cell
833 const cell &thecell = cset [i];
834
835 // check if this cell belongs to 'cother'
836 bool cbelongs = cother. check (thecell);
837
838 // if this cell should be removed, do it
839 if (cbelongs && (disjoint || maximal))
840 {
841 // if (!quiet && (sseq. show || sseq. log))
842 // sseq << '0' << thecell << '\n';
843 cother. remove (thecell);
844 cset. removenum (i --);
845 continue;
846 }
847
848 // detect the length of the boundary of this cell
849 int blen = boundarylength (thecell);
850
851 // should this cell be kept secure from the collapses?
852 bool keepit = ckeep. check (thecell);
853
854 // go through all the cells in this boundary
855 for (int j = 0; j < blen; ++ j)
856 {
857 // take the j-th boundary cell
858 cell bcell = boundarycell (thecell, j);
859
860 // check if this cell belongs to the other complex
861 bool bbelongs = bother. check (bcell);
862
863 // if it is in 'bother', then skip it if disjoint
864 if (bbelongs && disjoint)
865 continue;
866
867 // add this cell to 'bset' or get its number
868 int_t prev = bset. size ();
869 int_t number = addbd ? bset. add (bcell) :
870 bset. getnumber (bcell);
871
872 // if the cell is not there, skip it
873 if (number < 0)
874 continue;
875
876 // if this is the first occurrence of the cell
877 if ((prev < bset. size ()) || (coblen [number] == 0))
878 {
879 // write it to the sequence if necessary
880 if (!quiet && !bbelongs)
881 sseq << '1' << bcell << '\n';
882
883 // if this cell should be kept, mark it
884 if (keepit || bkeep. check (bset [number]) ||
885 bother. check (bset [number]))
886 coblen [number] = 13;
887 else
888 coblen [number] = 1;
889 coboundary [number]. add (thecell);
890 }
891
892 // otherwise add the corresponding coboundary link
893 else
894 {
895 ++ (coblen [number]);
896 coboundary [number]. add (thecell);
897 }
898 }
899 }
900 if (!quiet)
901 sseq << "D 100\n";
902
903 // show a dot
904 if (!quiet)
905 scon << "*\b";
906
907 // prepare tables for cells to be removed
908 hashedset<cell> cremove, bremove;
909 int_t nremove = 0;
910
911 // go through all the free faces
912 if (!quiet)
913 sseq << "\"Collapsing free faces...\"\n";
914 while (1)
915 {
916 // find a free face
917 int_t ncell = findelem (coblen, 1, bset. size ());
918
919 // if not found then finish
920 if (ncell < 0)
921 break;
922
923 // collapse this cell with its parent cell
924 const cell &parent = coboundary [ncell] [0];
925 int blen = boundarylength (parent);
926 coblen [ncell] = 0;
927
928 // remove the parent cell from coboundaries
929 for (int j = 0; j < blen; ++ j)
930 {
931 cell thecell = boundarycell (parent, j);
932 int_t number = bset. getnumber (thecell);
933 if ((number >= 0) && (coblen [number] > 0))
934 {
935 coboundary [number]. remove (parent);
936 -- (coblen [number]);
937 }
938 }
939
940 // write these cells to the sequence file
941 if (!quiet)
942 {
943 sseq << '0' << bset [ncell] << '\n';
944 sseq << '0' << parent << '\n';
945 }
946
947 // mark these cells for removal
948 cremove. add (parent);
949 bremove. add (bset [ncell]);
950 ++ nremove;
951
952 // clear the coboundary, because it is no longer used
953 coboundary [ncell] = empty;
954 }
955
956 // add the computed coboundaries if required
957 if (addcob)
958 {
959 for (int_t i = 0; i < bset. size (); ++ i)
960 if (!bremove. check (bset [i]))
961 {
962 coboundary [i]. remove (cremove);
963 (*(cob [d - 1])) [bset [i]]. add
964 (coboundary [i]);
965 }
966 }
967
968 // remove cells that are scheduled for removal from 'cset'
969 if (nremove == cset. size ())
970 {
971 removeall (d);
972 other. removeall (d);
973 }
974 else
975 {
976 for (int_t i = 0; i < nremove; ++ i)
977 cset. remove (cremove [i]);
978 }
979
980 // remove from the set of boundary cells these scheduled for removal
981 for (int_t i = 0; i < nremove; ++ i)
982 bset. remove (bremove [i]);
983
984 // update the dimension of the cell complex - is this necessary?
985 decreasedimension ();
986
987 // show a dot
988 if (!quiet)
989 scon << '.';
990
991 return nremove;
992} /* gcomplex<cell,euclidom>::collapse */
993
994template <class cell, class euclidom>
995int_t gcomplex<cell,euclidom>::collapse (gcomplex<cell,euclidom> &other,
996 gcomplex<cell,euclidom> &keep, bool addbd, bool addcob,
997 bool disjoint, const int *level, bool quiet)
998{
999 // determine the lower bound for the adding boundaries levels
1000 int dmin = 0;
1001 if (level)
1002 {
1003 while ((dmin < dim ()) && (!level [dmin]))
1004 ++ dmin;
1005 if (dmin && level [dmin])
1006 -- dmin;
1007 }
1008
1009 // add boundaries to high-dimensional cells in 'keep'
1010 while (keep. dim () > dim ())
1011 {
1012 keep. addboundaries (keep. dim ());
1013 keep. removeall (keep. dim ());
1014 }
1015
1016 // add boundaries to high-dimensional cells in 'other'
1017 if (other. dim () > dim ())
1018 {
1019 other. addboundaries (other. dim ());
1020 other. removeall (other. dim ());
1021 }
1022
1023 // add boundaries and collapse in all the dimensions of interest
1024 int_t counter = 0;
1025 for (int d = dim (); d > dmin; -- d)
1026 {
1027 // add boundaries to the other cell complexes
1028 keep. addboundaries (d);
1029
1030 // add boundaries and collapse this level
1031 counter += collapse (d, other, keep, addbd, addcob, disjoint,
1032 quiet);
1033
1034 // remove unnecessary cells from 'other'
1035 if (disjoint)
1036 other. removeall (d);
1037
1038 // remove unnecessary cells from 'keep'
1039 keep. removeall (d);
1040 }
1041
1042 // forget all the other cells of minimal dimension which are not used
1043 if (!disjoint && (dim () >= dmin) && (other. dim () >= dmin))
1044 {
1045 hashedset<cell> &cset = *(tab [dmin]);
1046 hashedset<cell> &cother = *(other. tab [dmin]);
1047 for (int_t i = 0; i < cother. size (); ++ i)
1048 {
1049 if (!(cset. check (cother [i])))
1050 cother. removenum (i --);
1051 }
1052 }
1053
1054 // remove unused cells which were supposed to be kept
1055 if (!keep. empty ())
1056 {
1057 gcomplex<cell,euclidom> empty;
1058 keep = empty;
1059 }
1060
1061 if (!quiet)
1062 scon << ' ';
1063
1064 return counter;
1065} /* gcomplex<cell,euclidom>::collapse */
1066
1067template <class cell, class euclidom>
1068void gcomplex<cell,euclidom>::increasedimension (int d)
1069{
1070 // create a new table for sets of cells
1071 hashedset<cell> **newtab = new hashedset<cell> *[d + 1];
1072 mvmap<cell,cell> **newcob = new mvmap<cell,cell> *[d + 1];
1073
1074 // copy anything that was in the old table
1075 for (int i = 0; i < n; ++ i)
1076 {
1077 newtab [i] = tab [i];
1078 newcob [i] = cob [i];
1079 }
1080
1081 // delete the old table if it was allocated
1082 if (tab)
1083 delete [] tab;
1084 if (cob)
1085 delete [] cob;
1086
1087 // initialize the remaining portion of the new table
1088 tab = newtab;
1089 cob = newcob;
1090 for (int i = n; i < d + 1; ++ i)
1091 {
1092 tab [i] = new hashedset<cell>;
1093 cob [i] = new mvmap<cell,cell>;
1094 }
1095
1096 // set the new dimension
1097 n = d + 1;
1098
1099 return;
1100} /* gcomplex<cell,euclidom>::increasedimension */
1101
1102template <class cell, class euclidom>
1103void gcomplex<cell,euclidom>::decreasedimension ()
1104{
1105 while (n && tab [n - 1] -> empty ())
1106 {
1107 -- n;
1108 delete tab [n];
1109 delete cob [n];
1110 }
1111 if (!n)
1112 {
1113 delete [] tab;
1114 tab = NULL;
1115 delete [] cob;
1116 cob = NULL;
1117 }
1118 return;
1119} /* gcomplex<cell,euclidom>::decreasedimension */
1120
1121// --------------------------------------------------
1122
1127template <class cell, class euclidom>
1128chaincomplex<euclidom> &createchaincomplex (chaincomplex<euclidom> &c,
1129 const gcomplex<cell,euclidom> &g, bool quiet = false)
1130{
1131 if (g. dim () < 0)
1132 return c;
1133 c. def_gen (0, g [0]. size ());
1134 for (int d = 1; d <= g. dim (); ++ d)
1135 {
1136 c. def_gen (d, g [d]. size ());
1137 for (int_t i = 0; i < g [d]. size (); ++ i)
1138 {
1139 int len = boundarylength (g [d] [i]);
1140 for (int j = 0; j < len; ++ j)
1141 {
1142 // take the j-th boundary cell
1143 cell thecell = boundarycell (g [d] [i], j);
1144
1145 // add it to the chain complex
1146 if (g. check (thecell))
1147 {
1148 int icoef = boundarycoef
1149 (g [d] [i], j);
1150 euclidom coef;
1151 if (icoef < 0)
1152 {
1153 coef = -icoef;
1154 coef = -coef;
1155 }
1156 else
1157 coef = icoef;
1158 c. add (d, g [d - 1]. getnumber
1159 (thecell), i, coef);
1160 }
1161 }
1162 }
1163 if (!quiet)
1164 {
1165 if (d < g. dim ())
1166 scon << '.';
1167 else
1168 scon << ". ";
1169 }
1170 }
1171 return c;
1172} /* createchaincomplex */
1173
1179template <class cell, class euclidom>
1180std::ostream &writechaincomplex (std::ostream &out,
1181 const gcomplex<cell,euclidom> &g, bool symbolicnames = false,
1182 bool quiet = false)
1183{
1184 if (g. dim () < 0)
1185 return out;
1186 out << "chaincomplex\n\n";
1187 out << "maxdimension " << g. dim () << "\n\n";
1188 out << "dimension 0: " << g [0]. size () << "\n\n";
1189 for (int d = 1; d <= g. dim (); ++ d)
1190 {
1191 out << "dimension " << d << ": " << g [d]. size () << "\n";
1192 for (int_t i = 0; i < g [d]. size (); ++ i)
1193 {
1194 bool cellwritten = false;
1195 int len = boundarylength (g [d] [i]);
1196 for (int j = 0; j < len; ++ j)
1197 {
1198 // take the j-th boundary cell
1199 cell thecell = boundarycell (g [d] [i], j);
1200
1201 // add it to the chain complex
1202 if (g. check (thecell))
1203 {
1204 int icoef = boundarycoef
1205 (g [d] [i], j);
1206 euclidom coef;
1207 if (icoef < 0)
1208 {
1209 coef = -icoef;
1210 coef = -coef;
1211 }
1212 else
1213 coef = icoef;
1214 if (!cellwritten)
1215 {
1216 out << "\t# ";
1217 if (symbolicnames)
1218 out << g [d] [i];
1219 else
1220 out << (i + 1);
1221 out << " = ";
1222 if (-coef == 1)
1223 out << "- ";
1224 }
1225 else if (coef == 1)
1226 out << " + ";
1227 else if (-coef == 1)
1228 out << " - ";
1229 else
1230 out << " + " << coef <<
1231 " * ";
1232 if (symbolicnames)
1233 out << thecell;
1234 else
1235 out << (1 + g [d - 1].
1236 getnumber (thecell));
1237 cellwritten = true;
1238 }
1239 }
1240 if (cellwritten)
1241 out << '\n';
1242 }
1243 if (!quiet)
1244 {
1245 if (d < g. dim ())
1246 scon << '.';
1247 else
1248 scon << ". ";
1249 }
1250 out << '\n';
1251 }
1252 return out;
1253} /* writechaincomplex */
1254
1259template <class cell, class euclidom>
1260chaincomplex<euclidom> &createchaincomplex (chaincomplex<euclidom> &c,
1261 const gcomplex<cell,euclidom> &g, const gcomplex<cell,euclidom> &rel,
1262 bool quiet = false)
1263{
1264 // if the geometric complex is empty, don't modify the chain complex
1265 if (g. dim () < 0)
1266 return c;
1267
1268 // prepare a table of numbers of ignored cells in the geom. complex
1269 multitable<int_t> *numbers0 = new multitable<int_t>;
1270 int_t count0 = g [0]. size ();
1271 if (rel. dim () >= 0)
1272 {
1273 count0 -= rel [0]. size ();
1274 int_t j = 0;
1275 for (int_t i = 0; i < rel [0]. size (); ++ i)
1276 {
1277 (*numbers0) [j] = g [0]. getnumber (rel [0] [i]);
1278 if ((*numbers0) [j] < 0)
1279 ++ count0;
1280 else
1281 ++ j;
1282 }
1283 }
1284
1285 // set the number of generators of the given dimension
1286 c. def_gen (0, count0);
1287
1288 // create boundary matrices
1289 for (int d = 1; d <= g. dim (); ++ d)
1290 {
1291 // prepare a table of numbers to be ignored
1292 multitable<int_t> *numbers1 = new multitable<int_t>;
1293 int_t count1 = g [d]. size ();
1294 if (rel. dim () >= d)
1295 {
1296 count1 -= rel [d]. size ();
1297 int_t j = 0;
1298 for (int_t i = 0; i < rel [d]. size (); ++ i)
1299 {
1300 (*numbers1) [j] =
1301 g [d]. getnumber (rel [d] [i]);
1302 if ((*numbers1) [j] < 0)
1303 ++ count1;
1304 else
1305 ++ j;
1306 }
1307 }
1308
1309 // set the number of generators of this dimension
1310 c. def_gen (d, count1);
1311
1312 // create the boundary connections with coefficients
1313 for (int_t i = 0; i < g [d]. size (); ++ i)
1314 {
1315 // if this cell is in the other complex, ignore it
1316 if ((rel. dim () >= d) &&
1317 (rel [d]. check (g [d] [i])))
1318 continue;
1319
1320 // determine the number of this cell
1321 int_t n1 = i;
1322 if (n1 >= count1)
1323 n1 = (*numbers1) [n1 - count1];
1324
1325 // get the length of the cell
1326 int len = boundarylength (g [d] [i]);
1327
1328 // retrieve boundary cells and make boundary formula
1329 for (int j = 0; j < len; ++ j)
1330 {
1331 // take the j-th boundary cell
1332 cell thecell = boundarycell (g [d] [i], j);
1333
1334 // add it to the chain complex if relevant
1335 if (g [d - 1]. check (thecell) &&
1336 ((rel. dim () < d - 1) ||
1337 (!rel [d - 1]. check (thecell))))
1338 {
1339 // determine the number of the cell
1340 int_t n0 = g [d - 1].
1341 getnumber (thecell);
1342
1343 // if out of range, translate it
1344 if (n0 >= count0)
1345 n0 = (*numbers0)
1346 [n0 - count0];
1347
1348 // determine the right coefficient
1349 int icoef = boundarycoef
1350 (g [d] [i], j);
1351 euclidom coef;
1352 if (icoef < 0)
1353 {
1354 coef = -icoef;
1355 coef = -coef;
1356 }
1357 else
1358 coef = icoef;
1359
1360 // add this link to the boundary
1361 c. add (d, n0, n1, coef);
1362 }
1363 }
1364 }
1365
1366 // forget unnecessary tables and prepare for the next loop
1367 delete numbers0;
1368 numbers0 = numbers1;
1369 count0 = count1;
1370
1371 // show progress indicator
1372 if (!quiet)
1373 {
1374 if (d < g. dim ())
1375 scon << '.';
1376 else
1377 scon << ". ";
1378 }
1379 }
1380
1381 // release used memory
1382 delete numbers0;
1383
1384 // finish
1385 return c;
1386} /* createchaincomplex */
1387
1389template <class cell, class euclidom>
1390std::ostream &writegenerators (std::ostream &out, const chain<euclidom> *hom,
1391 const chaincomplex<euclidom> &c,
1392 const gcomplex<cell,euclidom> &g, const int *level = NULL)
1393{
1394 bool firstlist = true;
1395 for (int d = 0; d <= c. dim (); ++ d)
1396 {
1397 if ((!level || level [d]) && !hom [d]. empty ())
1398 {
1399 if (firstlist)
1400 firstlist = false;
1401 else
1402 out << '\n';
1403 if (hom [d]. size () == 1)
1404 out << "The generator of H_" << d <<
1405 " follows:" << '\n';
1406 else
1407 out << "The " << hom [d]. size () <<
1408 " generators of H_" << d <<
1409 " follow:" << '\n';
1410 const hashedset<cell> &cset = g [d];
1411 for (int_t i = 0; i < hom [d]. size (); ++ i)
1412 {
1413 if (hom [d]. size () > 1)
1414 out << "generator " << (i + 1) <<
1415 '\n';
1416 const chain<euclidom> &lst =
1417 c. gethomgen (d, hom [d]. num (i));
1418 for (int_t j = 0; j < lst. size (); ++ j)
1419 out << lst. coef (j) << " * " <<
1420 cset [lst. num (j)] << '\n';
1421 }
1422 }
1423 }
1424 return out;
1425} /* writegenerators */
1426
1428template <class cell, class euclidom>
1429gcomplex<cell,euclidom> &creategraph (const mvmap<cell,cell> &m,
1430 gcomplex<cell,euclidom> &graph)
1431{
1432 for (int_t i = 0; i < m. size (); ++ i)
1433 {
1434 const cell &e = m. get (i);
1435 const hashedset<cell> &f = m (i);
1436 for (int_t j = 0; j < f. size (); ++ j)
1437 graph. add (e * f [j]);
1438 }
1439 return graph;
1440} /* creategraph */
1441
1446template <class cell, class euclidom>
1447bool acyclic (gcomplex<cell,euclidom> &c)
1448{
1449 // collapse the geometric complex and check if you can get one elem.
1450 gcomplex<cell,euclidom> empty;
1451 c. collapse (empty, empty, true, false, false, NULL, true);
1452 if (c. size () == 1)
1453 return true;
1454
1455 // create the corresponding chain complex and compute its homology
1456 chaincomplex<euclidom> cx (c. dim ());
1457 createchaincomplex (cx, c, true);
1458 chain<euclidom> *hom;
1459 cx. simple_form (static_cast<int *> (0), true);
1460 cx. simple_homology (hom, 0);
1461
1462 // if there is anything above the zero level, the set is not acyclic
1463 for (int i = 1; i <= cx. dim (); ++ i)
1464 {
1465 if (!hom [i]. empty ())
1466 {
1467 delete [] hom;
1468 return false;
1469 }
1470 }
1471
1472 // if there is more than one connected component, it is not acyclic
1473 if (hom [0]. size () != 1)
1474 {
1475 delete [] hom;
1476 return false;
1477 }
1478
1479 // if the coefficient is not 1, the set is not acyclic
1480 if (hom [0]. getcoefficient (0) == 1)
1481 {
1482 delete [] hom;
1483 return true;
1484 }
1485 else
1486 {
1487 delete [] hom;
1488 return false;
1489 }
1490} /* acyclic */
1491
1492// --------------------------------------------------
1493
1495template <class cell, class euclidom>
1496std::ostream &operator << (std::ostream &out,
1497 const gcomplex<cell,euclidom> &c)
1498{
1499 out << "; Geometric complex of dimension " << c. dim () <<
1500 " (" << c. size () << " cells total)." << '\n';
1501 for (int i = 0; i <= c. dim (); ++ i)
1502 out << "; Cells of dimension " << i << ':' << '\n' << c [i];
1503 return out;
1504} /* operator << */
1505
1507template <class cell, class euclidom>
1508std::istream &operator >> (std::istream &in, gcomplex<cell,euclidom> &c)
1509{
1510 ignorecomments (in);
1511 while (closingparenthesis (in. peek ()) != EOF)
1512 {
1513 cell x;
1514 in >> x;
1515 c. add (x);
1516 ignorecomments (in);
1517 }
1518 return in;
1519} /* operator >> */
1520
1521
1522// --------------------------------------------------
1523// ------------------ mv cell map -------------------
1524// --------------------------------------------------
1525
1528template <class cell, class euclidom, class element>
1529class mvcellmap
1530{
1531public:
1535 mvcellmap (gcomplex<cell,euclidom> *_g = 0);
1536
1540 mvcellmap (gcomplex<cell,euclidom> &_g);
1541
1543 mvcellmap (const mvcellmap<cell,euclidom,element> &m);
1544
1546 mvcellmap &operator = (const mvcellmap<cell,euclidom,element> &m);
1547
1549 ~mvcellmap ();
1550
1552 int dim () const;
1553
1555 const hashedset<cell> &get (int d) const;
1556
1558 const gcomplex<cell,euclidom> &getdomain () const;
1559
1561 const hashedset<element> &operator () (const cell &c) const;
1562
1565 void add (int d, const cell &c,
1566 const hashedset<element> &set);
1567
1569 void add (const cell &c, const hashedset<element> &set);
1570
1573 void add (int d, int_t n, const hashedset<element> &set);
1574
1577 void add (int d, const cell &c, const element &e);
1578
1580 void add (const cell &c, const element &e);
1581
1584 void add (int d, int_t n, const element &e);
1585
1586private:
1588 gcomplex<cell,euclidom> *g;
1589
1591 multitable <hashedset<element> > *images;
1592
1594 int dimension;
1595
1596}; /* class mvcellmap */
1597
1598// --------------------------------------------------
1599
1600template <class cell, class euclidom, class element>
1601inline mvcellmap<cell,euclidom,element>::mvcellmap
1602 (gcomplex<cell,euclidom> *_g)
1603{
1604 g = _g;
1605 if (!g || (g -> dim () < 0))
1606 {
1607 images = NULL;
1608 dimension = -1;
1609 return;
1610 }
1611 dimension = g -> dim ();
1612 images = new multitable <hashedset<element> > [dimension + 1];
1613 if (!images)
1614 throw "Cannot create a multivalued cellular map.";
1615 return;
1616} /* mvcellmap<cell,euclidom,element>::mvcellmap */
1617
1618template <class cell, class euclidom, class element>
1619inline mvcellmap<cell,euclidom,element>::mvcellmap
1620 (gcomplex<cell,euclidom> &_g)
1621{
1622 g = &_g;
1623 if (!g || (g -> dim () < 0))
1624 {
1625 images = NULL;
1626 dimension = -1;
1627 return;
1628 }
1629 dimension = g -> dim ();
1630 images = new multitable <hashedset<element> > [dimension + 1];
1631 if (!images)
1632 throw "Cannot create a multivalued cellular map.";
1633 return;
1634} /* mvcellmap<cell,euclidom,element>::mvcellmap */
1635
1636template <class cell, class euclidom, class element>
1637mvcellmap<cell,euclidom,element>::mvcellmap
1638 (const mvcellmap<cell,euclidom,element> &m)
1639{
1640 g = m. g;
1641 if (!g || (g -> dim () < 0))
1642 {
1643 images = NULL;
1644 dimension = -1;
1645 return;
1646 }
1647 dimension = g -> dim ();
1648 images = new multitable <hashedset<element> > [dimension + 1];
1649 if (!images)
1650 throw "Unable to construct a copy of a multivalued "
1651 "cellular map.";
1652 for (int i = 0; i < dimension + 1; ++ i)
1653 images [i] = m. images [i];
1654 return;
1655} /* mvcellmap<cell,euclidom,element>::mvcellmap */
1656
1657template <class cell, class euclidom, class element>
1658mvcellmap<cell,euclidom,element> &mvcellmap<cell,euclidom,element>::operator =
1659 (const mvcellmap<cell,euclidom,element> &m)
1660{
1661 if (images)
1662 delete [] images;
1663 g = m. g;
1664 dimension = m. dimension;
1665 if (!g || (dimension < 0))
1666 {
1667 images = NULL;
1668 return *this;
1669 }
1670 images = new multitable <hashedset<element> > [dimension + 1];
1671 if (!images)
1672 throw "Cannot copy a multivalued cellular map.";
1673 for (int i = 0; i < dimension + 1; ++ i)
1674 images [i] = m. images [i];
1675 return *this;
1676} /* mvcellmap<cell,euclidom,element>::operator = */
1677
1678template <class cell, class euclidom, class element>
1679inline mvcellmap<cell,euclidom,element>::~mvcellmap ()
1680{
1681 if (images)
1682 delete [] images;
1683 return;
1684} /* mvcellmap<cell,euclidom,element>::~mvcellmap */
1685
1686template <class cell, class euclidom, class element>
1687int mvcellmap<cell,euclidom,element>::dim () const
1688{
1689 return dimension;
1690} /* mvcellmap<cell,euclidom,element>::dim */
1691
1692template <class cell, class euclidom, class element>
1693const hashedset<cell> &mvcellmap<cell,euclidom,element>::get (int d) const
1694{
1695 if ((d < 0) || (d > dimension))
1696 throw "Wrong dimension to retrieve from mvcellmap.";
1697 return (*g) [d];
1698} /* mvcellmap<cell,euclidom,element>::get */
1699
1700template <class cell, class euclidom, class element>
1701const gcomplex<cell,euclidom> &mvcellmap<cell,euclidom,element>::getdomain ()
1702 const
1703{
1704 return *g;
1705} /* mvcellmap<cell,euclidom,element>::getdomain */
1706
1707template <class cell, class euclidom, class element>
1708const hashedset<element> &mvcellmap<cell,euclidom,element>::operator ()
1709 (const cell &c) const
1710{
1711 int d = c. dim ();
1712 if (d > dimension)
1713 throw "Trying to get the image of a cell of dim too high.";
1714 int_t n = (*g) [d]. getnumber (c);
1715 if (n < 0)
1716 throw "Trying to get the image of a cell not in the domain.";
1717 return images [d] [n];
1718} /* mvcellmap<cell,euclidom,element>::operator () */
1719
1720template <class cell, class euclidom, class element>
1721inline void mvcellmap<cell,euclidom,element>::add (int d, const cell &c,
1722 const hashedset<element> &set)
1723{
1724 if ((d < 0) || (d > dimension))
1725 throw "Wrong dimension for adding a cell to a map.";
1726 if (!g -> check (c))
1727 // throw "The cell does not belong to the domain of the map.";
1728 g -> add (c);
1729 images [d] [(*g) [d]. getnumber (c)]. add (set);
1730 return;
1731} /* mvcellmap<cell,euclidom,element>::add */
1732
1733template <class cell, class euclidom, class element>
1734inline void mvcellmap<cell,euclidom,element>::add (const cell &c,
1735 const hashedset<element> &set)
1736{
1737 add (c. dim (), c, set);
1738 return;
1739} /* mvcellmap<cell,euclidom,element>::add */
1740
1741template <class cell, class euclidom, class element>
1742inline void mvcellmap<cell,euclidom,element>::add (int d, int_t n,
1743 const hashedset<element> &set)
1744{
1745 if ((d < 0) || (d > dimension))
1746 throw "Wrong dimension for adding an element to the image.";
1747 images [d] [n]. add (set);
1748 return;
1749} /* mvcellmap<cell,euclidom,element>::add */
1750
1751template <class cell, class euclidom, class element>
1752inline void mvcellmap<cell,euclidom,element>::add (int d, const cell &c,
1753 const element &e)
1754{
1755 if ((d < 0) || (d > dimension))
1756 throw "Wrong dimension for adding a cell to a map.";
1757 if (!g -> check (c))
1758 // throw "The cell does not belong to the domain of the map.";
1759 g -> add (c);
1760 images [d] [(*g) [d]. getnumber (c)]. add (e);
1761 return;
1762} /* mvcellmap<cell,euclidom,element>::add */
1763
1764template <class cell, class euclidom, class element>
1765inline void mvcellmap<cell,euclidom,element>::add (const cell &c,
1766 const element &e)
1767{
1768 add (c. dim (), c, e);
1769 return;
1770} /* mvcellmap<cell,euclidom,element>::add */
1771
1772template <class cell, class euclidom, class element>
1773inline void mvcellmap<cell,euclidom,element>::add (int d, int_t n,
1774 const element &set)
1775{
1776 if ((d < 0) || (d > dimension))
1777 throw "Wrong dimension for adding an element to the image.";
1778 images [d] [n]. add (set);
1779 return;
1780} /* mvcellmap<cell,euclidom,element>::add */
1781
1782// --------------------------------------------------
1783
1787template <class cell, class euclidom, class element>
1788void creategraph (const mvcellmap<cell,euclidom,element> &m,
1789 gcomplex<cell,euclidom> &c, bool addbd)
1790{
1791 // go through all the dimensions of the domain cell complex
1792 for (int d1 = 0; d1 <= m. dim (); ++ d1)
1793 {
1794 // go through all the cells of this dimension
1795 const hashedset<cell> &cset = m. get (d1);
1796 for (int_t i = 0; i < cset. size (); ++ i)
1797 {
1798 // create a cell complex of the image of this cell
1799 const cell &thecell = cset [i];
1800 const hashedset<element> &set = m (thecell);
1801 gcomplex<cell,euclidom> img;
1802 for (int_t j = 0; j < set. size (); ++ j)
1803 img. add (cell (set [j]));
1804 if (addbd)
1805 img. addboundaries ();
1806
1807 // go through all the dimensions of this cell complex
1808 for (int d2 = 0; d2 <= img. dim (); ++ d2)
1809 {
1810 // add all the products to the graph
1811 const hashedset<cell> &cs = img [d2];
1812 for (int_t j = 0; j < cs. size (); ++ j)
1813 c. add (thecell * cs [j]);
1814 }
1815 }
1816 if (d1 < m. dim ())
1817 scon << '.';
1818 else
1819 scon << ' ';
1820 }
1821 return;
1822} /* creategraph */
1823
1827template <class cell, class euclidom, class element>
1828void creategraph (const mvcellmap<cell,euclidom,element> &m,
1829 const gcomplex<cell,euclidom> &rel,
1830 gcomplex<cell,euclidom> &c, bool addbd)
1831{
1832 // go through all the dimensions of the domain cell complex
1833 for (int d1 = 0; d1 <= m. dim (); ++ d1)
1834 {
1835 // go through all the cells of this dimension
1836 const hashedset<cell> &cset = m. get (d1);
1837 for (int_t i = 0; i < cset. size (); ++ i)
1838 {
1839 // create a cell complex of the image of this cell
1840 const cell &thecell = cset [i];
1841 const hashedset<element> &set = m (thecell);
1842 gcomplex<cell,euclidom> img;
1843 for (int_t j = 0; j < set. size (); ++ j)
1844 img. add (cell (set [j]));
1845 if (addbd)
1846 img. addboundaries ();
1847
1848 // go through all the dimensions of this cell complex
1849 for (int d2 = 0; d2 <= img. dim (); ++ d2)
1850 {
1851 // add all the products to the graph
1852 const hashedset<cell> &cs = img [d2];
1853 for (int_t j = 0; j < cs. size (); ++ j)
1854 {
1855 cell bigcell = thecell * cs [j];
1856 if (!rel. check (bigcell))
1857 c. add (bigcell);
1858 }
1859 }
1860 }
1861 if (d1 < m. dim ())
1862 scon << '.';
1863 else
1864 scon << ' ';
1865 }
1866 return;
1867} /* creategraph */
1868
1873template <class cell, class euclidom, class element>
1874void createcellmap (const mvcellmap<cell,euclidom,element> &m,
1875 mvcellmap<cell,euclidom,cell> &cm)
1876{
1877 // go through all the dimensions of the domain cell complex
1878 int spacedim = -1;
1879 const gcomplex<cell,euclidom> &dom = m. getdomain ();
1880 for (int d = m. dim (); d >= 0; -- d)
1881 {
1882 // go through all the cells of this dimension
1883 const hashedset<cell> &cset = m. get (d);
1884 for (int_t i = 0; i < cset. size (); ++ i)
1885 {
1886 // extract the cell of interest and its image
1887 const cell &thecell = cset [i];
1888 const hashedset<element> &set = m (thecell);
1889
1890 // the image must not be empty!
1891 if (set. empty ())
1892 throw "An empty image of a cell found.";
1893
1894 // correct the dimension of the space if necessary
1895 if (spacedim < 0)
1896 spacedim = set [0]. dim ();
1897
1898 // if this is maximal dimension, extract one point
1899 if (d == spacedim)
1900 {
1901 // cm. add (d, thecell, cell (set [0]. coord (),
1902 // set [0]. coord (), spacedim));
1903 cm. add (d, thecell, cell (set [0], 0));
1904 continue;
1905 }
1906
1907 // create a cell complex of the image of this cell
1908 gcomplex<cell,euclidom> img;
1909 for (int_t j = 0; j < set. size (); ++ j)
1910 img. add (cell (set [j]));
1911
1912 // create a cell complex to keep while collapsing
1913 gcomplex<cell,euclidom> keep;
1914 const hashedset<cell> &cob =
1915 dom. getcob (thecell, d);
1916 for (int_t j = 0; j < cob. size (); ++ j)
1917 keep. add (cm (cob [j]));
1918
1919 // reduce 'img' towards 'keep'
1920 gcomplex<cell,euclidom> empty;
1921 img. collapse (empty, keep, 1, 0, 0, NULL, 1);
1922
1923 // set this to be the image of this cell
1924 // *** THIS MUST BE IMPROVED ***
1925 for (int j = 0; j <= img. dim (); ++ j)
1926 cm. add (d, thecell, img [j]);
1927
1928 // show progress indicator
1929 if (i && !(i % 373))
1930 scon << std::setw (12) << i <<
1931 "\b\b\b\b\b\b\b\b\b\b\b\b";
1932 }
1933 if (cset. size () > 373)
1934 scon << " \b\b\b\b\b\b\b\b\b\b\b\b";
1935 scon << '.';
1936 }
1937 scon << ' ';
1938 return;
1939} /* createcellmap */
1940
1947template <class cell, class euclidom, class element>
1948bool createcellmap (const mvcellmap<cell,euclidom,element> &m,
1949 const mvcellmap<cell,euclidom,element> &rel,
1950 mvcellmap<cell,euclidom,cell> &cm, bool verifyacyclicity)
1951{
1952 // prepare the default result
1953 bool result = true;
1954
1955 // go through all the dimensions of the domain cell complex
1956 const gcomplex<cell,euclidom> &dom = m. getdomain ();
1957 const gcomplex<cell,euclidom> &reldom = rel. getdomain ();
1958
1959 int maxdim = m. dim ();
1960 for (int d = maxdim; d >= 0; -- d)
1961 {
1962 // go through all the cells of this dimension
1963 const hashedset<cell> &cset = m. get (d);
1964 for (int_t i = 0; i < cset. size (); ++ i)
1965 {
1966 // extract the cell of interest and its image
1967 const cell &thecell = cset [i];
1968 const hashedset<element> &set = m (thecell);
1969
1970 // the image must not be empty!
1971 if (set. empty ())
1972 throw "An empty image of a cell found.";
1973
1974 // if this is maximal dimension, extract one point
1975 if (d == maxdim)
1976 {
1977 if (reldom. check (thecell))
1978 continue;
1979 // throw "This cell shouldn't be here.";
1980 // cm. add (d, thecell, cell (set [0]. coord (),
1981 // set [0]. coord (), set [0]. dim ()));
1982 cm. add (d, thecell, cell (set [0], 0));
1983 continue;
1984 }
1985
1986 // create a cell complex of the image of this cell
1987 gcomplex<cell,euclidom> img;
1988 for (int_t j = 0; j < set. size (); ++ j)
1989 img. add (cell (set [j]));
1990
1991 // create a cell complex to keep while collapsing
1992 gcomplex<cell,euclidom> keep;
1993 const hashedset<cell> &cob =
1994 dom. getcob (thecell, d);
1995
1996 for (int_t j = 0; j < cob. size (); ++ j)
1997 keep. add (cm (cob [j]));
1998
1999 // create a cell complex from 'rel'
2000 gcomplex<cell,euclidom> other;
2001 if (reldom. check (thecell))
2002 {
2003 const hashedset<element> &relset =
2004 rel (thecell);
2005 for (int_t j = 0; j < relset. size (); ++ j)
2006 {
2007 cell c = cell (relset [j]);
2008 keep. add (c);
2009 other. add (c);
2010 // img. remove (c);
2011 }
2012 }
2013
2014 // reduce 'img' towards 'keep'
2015 gcomplex<cell,euclidom> empty;
2016 img. collapse (empty, keep, 1, 0, 0, NULL, 1);
2017
2018 // verify the acyclicity of this image if requested
2019 if (verifyacyclicity)
2020 {
2021 gcomplex<cell,euclidom> imgdupl (img);
2022 if (!acyclic (imgdupl))
2023 {
2024 result = false;
2025 verifyacyclicity = false;
2026 }
2027 }
2028
2029 // remove the other cells and their faces from img
2030 other. addboundaries ();
2031 img. remove (other);
2032
2033 // verify the acyclicity of the other complex
2034 if (verifyacyclicity && other. size ())
2035 {
2036 // note: acycl. check destroys 'other'
2037 if (!acyclic (other))
2038 {
2039 result = false;
2040 verifyacyclicity = false;
2041 }
2042 }
2043
2044 // set this to be the image of this cell
2045 // *** THIS MUST BE IMPROVED ***
2046 for (int j = 0; j <= img. dim (); ++ j)
2047 cm. add (d, thecell, img [j]);
2048
2049 // show progress indicator
2050 if (i && !(i % 373))
2051 scon << std::setw (12) << i <<
2052 "\b\b\b\b\b\b\b\b\b\b\b\b";
2053 }
2054 if (cset. size () > 373)
2055 scon << " \b\b\b\b\b\b\b\b\b\b\b\b";
2056 scon << '.';
2057 }
2058 scon << ' ';
2059 return result;
2060} /* createcellmap */
2061
2062// --------------------------------------------------
2063
2065template <class cell, class euclidom, class element>
2066std::ostream &operator << (std::ostream &out,
2067 const mvcellmap<cell,euclidom,element> &m)
2068{
2069 for (int d = 0; d <= m. dim (); ++ d)
2070 {
2071 out << "; Dimension " << d << ':' << '\n';
2072 const hashedset<cell> &cset = m. get (d);
2073 for (int_t i = 0; i < cset. size (); ++ i)
2074 {
2075 out << cset [i] << " -> ";
2076 write (out, m (cset [i]), SMALL_SIZE) << '\n';
2077 }
2078 }
2079 return out;
2080} /* operator << */
2081
2082
2083} // namespace homology
2084} // namespace chomp
2085
2086#endif // _CHOMP_HOMOLOGY_GCOMPLEX_H_
2087
2089
chains.h
This file contains classes and functions related to algebraic chain complexes and chain maps,...
chomp::homology::chain
This class defines objects which represent chains as finite sequences of elements identified by integ...
Definition: chains.h:94
chomp::homology::chaincomplex
This is an implementation of a chain complex over an arbitrary ring.
Definition: chains.h:2662
chomp::homology::gcomplex
The class that defines a geometric complex - a set of cells (cubes, simplices, etc).
Definition: gcomplex.h:85
chomp::homology::gcomplex::increasedimension
void increasedimension(int d)
Increases the dimension of the complex to take this cell.
Definition: gcomplex.h:1068
chomp::homology::gcomplex::n
int n
The number of tables.
Definition: gcomplex.h:239
chomp::homology::gcomplex::size
int_t size() const
Returns the number of cells in the complex.
Definition: gcomplex.h:360
chomp::homology::gcomplex::remove
gcomplex< cell, euclidom > & remove(const cell &c)
Remove a cell from the geometric complex.
Definition: gcomplex.h:477
chomp::homology::gcomplex::check
bool check(const cell &c) const
Check whether the given cell is in the complex.
Definition: gcomplex.h:562
chomp::homology::gcomplex::collapse
int_t collapse(int d, gcomplex< cell, euclidom > &other, const gcomplex< cell, euclidom > &keep, bool addbd, bool addcob, bool disjoint, bool quiet=false)
Adds boundaries of all the cells of the given dimension and then performs free face collapses.
Definition: gcomplex.h:752
chomp::homology::gcomplex::operator[]
const hashedset< cell > & operator[](int d) const
Returns the set of cells of the given dimension.
Definition: gcomplex.h:380
chomp::homology::gcomplex::empty
bool empty() const
Returns 'true' iff the cell complex is empty.
Definition: gcomplex.h:369
chomp::homology::gcomplex::getcob
const hashedset< cell > & getcob(const cell &c) const
Returns the coboundary of the given cell.
Definition: gcomplex.h:389
chomp::homology::gcomplex::cob
mvmap< cell, cell > ** cob
The tables with coboundaries of cells of these dimensions.
Definition: gcomplex.h:236
chomp::homology::gcomplex::tab
hashedset< cell > ** tab
The tables with cells of dimension 0, 1, 2, etc.
Definition: gcomplex.h:233
chomp::homology::gcomplex::gcomplex
gcomplex()
The default constructor.
Definition: gcomplex.h:252
chomp::homology::gcomplex::dim
int dim() const
Returns the dimension of the complex, that is, the highest dimension of any cell contained in the com...
Definition: gcomplex.h:354
chomp::homology::gcomplex::add
gcomplex< cell, euclidom > & add(const cell &c)
Add a cell to the geometric complex.
Definition: gcomplex.h:405
chomp::homology::gcomplex::decreasedimension
void decreasedimension()
Frees empty sets if there are no cells of high dimensions.
Definition: gcomplex.h:1103
chomp::homology::gcomplex::operator=
gcomplex & operator=(const gcomplex< cell, euclidom > &c)
The assignment operator.
Definition: gcomplex.h:285
chomp::homology::gcomplex::~gcomplex
~gcomplex()
The destructor.
Definition: gcomplex.h:337
chomp::homology::gcomplex::addboundaries
int_t addboundaries(int d, bool addcob=false)
Adds boundaries of all the cells of given dimension to the geometric complex.
Definition: gcomplex.h:647
chomp::homology::gcomplex::removeall
gcomplex< cell, euclidom > & removeall(int d)
Remove all the cells of the given dimension.
Definition: gcomplex.h:543
chomp::homology::multitable
A container for elements placed in a table (like a vector) that is actually built of many smaller tab...
Definition: multitab.h:65
chomp::homology::mvcellmap
This class represents a multivalued map whose domain is a geometric complex.
Definition: gcomplex.h:1530
chomp::homology::mvcellmap::operator()
const hashedset< element > & operator()(const cell &c) const
Returns the image of a given cell.
Definition: gcomplex.h:1709
chomp::homology::mvcellmap::mvcellmap
mvcellmap(gcomplex< cell, euclidom > *_g=0)
The constructor of a map with its domain.
Definition: gcomplex.h:1602
chomp::homology::mvcellmap::add
void add(int d, const cell &c, const hashedset< element > &set)
Adds a set to the image of a given cell, provided the dimension of the cell is known.
Definition: gcomplex.h:1721
chomp::homology::mvcellmap::get
const hashedset< cell > & get(int d) const
Returns the given level of the geometric complex.
Definition: gcomplex.h:1693
chomp::homology::mvcellmap::dim
int dim() const
Returns the dimension of the domain of the map.
Definition: gcomplex.h:1687
chomp::homology::mvcellmap::g
gcomplex< cell, euclidom > * g
A pointer to the domain of the map.
Definition: gcomplex.h:1588
chomp::homology::mvcellmap::getdomain
const gcomplex< cell, euclidom > & getdomain() const
Returns a reference of the domain cell complex of the map.
Definition: gcomplex.h:1701
chomp::homology::mvcellmap::~mvcellmap
~mvcellmap()
The destructor.
Definition: gcomplex.h:1679
chomp::homology::mvcellmap::dimension
int dimension
The dimension of the domain of the map.
Definition: gcomplex.h:1594
chomp::homology::mvcellmap::images
multitable< hashedset< element > > * images
The array of images of the elements of each dimension.
Definition: gcomplex.h:1591
chomp::homology::mvcellmap::operator=
mvcellmap & operator=(const mvcellmap< cell, euclidom, element > &m)
The assignment operator.
Definition: gcomplex.h:1659
config.h
This file contains some precompiler definitions which indicate the operating system and/or compiler u...
gcomplex.h
This file contains a definition of a general geometric complex which represents a polyhedron.
int_t
int int_t
Index type for indexing arrays, counting cubes, etc.
Definition: config.h:115
hashsets.h
This file contains the definition of the container "hashedset" which can be used to represent a set o...
SMALL_SIZE
#define SMALL_SIZE
This constant passed to the function 'write' makes the hashed set be displayed in a way that is appro...
Definition: hashsets.h:813
integer.h
This file defines a class "integer" which represents the ring of integers or the field of integers mo...
chomp::homology::sseq
outputstream sseq
An auxiliary stream which captures sequences of processed data.
chomp::homology::writegenerators
std::ostream & writegenerators(std::ostream &out, const chain< euclidom > *hom, const chaincomplex< euclidom > &c, const gcomplex< tCell, euclidom > &g, const int *level, int xdim, int format=0)
Writes projected homology generators of a cubical complex to a file.
Definition: cubmaps.h:64
chomp::homology::acyclic
bool acyclic(int dim, BitField &b)
Checks whether this cube's nieghbors form an acyclic set.
Definition: cubacycl.h:70
chomp::homology::closingparenthesis
int closingparenthesis(int ch)
Returns the matching closing parenthesis for the given opening one or EOF if none.
Definition: textfile.h:411
chomp::homology::operator<<
std::ostream & operator<<(std::ostream &out, const bincube< Dim, twoPower > &b)
Definition: bincube.h:907
chomp::homology::write
int write(std::ostream &out, const coordtype *c, int dim, char parenthesis=40, char closing=0)
Definition: pointset.h:2148
chomp::homology::addboundaries
void addboundaries(gcomplex< cell, euclidom > &Xcompl, gcomplex< cell, euclidom > &Acompl, int minlevel, bool bothsets, const char *Xname, const char *Aname)
Adds boundaries to the geometric complex X or to both X and A.
Definition: homtools.h:835
chomp::homology::boundarylength
int boundarylength(const tCellBase< coordtype > &q)
Returns the length of the boundary of a cell.
Definition: cellbase.h:501
chomp::homology::collapse
void collapse(gcomplex< cell, euclidom > &Xcompl, gcomplex< cell, euclidom > &Acompl, gcomplex< cell, euclidom > &Xkeepcompl, const char *Xname, const char *Aname, bool addbd=true, bool addcob=false, bool disjoint=true, const int *level=NULL)
Collapses a pair of geometric complexes.
Definition: homtools.h:742
chomp::homology::boundarycell
tCellBase< coordtype > boundarycell(const tCellBase< coordtype > &q, int i, bool onlyexisting)
Computes the i-th boundary element of a cell.
Definition: cellbase.h:485
chomp::homology::createcellmap
void createcellmap(const mvcellmap< cell, euclidom, element > &m, mvcellmap< cell, euclidom, cell > &cm)
Creates the graph of the map as a cell complex while reducing as possible.
Definition: gcomplex.h:1874
chomp::homology::findelem
int_t findelem(const multitable< element > &tab, const element &e, int_t len)
Finds the given element in the table of given length.
Definition: gcomplex.h:705
chomp::homology::operator>>
std::istream & operator>>(std::istream &in, bincube< Dim, twoPower > &b)
Definition: bincube.h:914
chomp::homology::writechaincomplex
std::ostream & writechaincomplex(std::ostream &out, const gcomplex< cell, euclidom > &g, bool symbolicnames=false, bool quiet=false)
Writes out a chain complex of the geometric cell complex.
Definition: gcomplex.h:1180
chomp::homology::decreasedimension
void decreasedimension(gcomplex< cell, euclidom > &Acompl, int dim, const char *name)
Decreases the dimension of the geometric complex by adding boundary cells to all the cells on higher ...
Definition: homtools.h:813
chomp::homology::scon
outputstream scon
The console output stream to which one should put all the junk that spoils the log file,...
chomp::homology::createchaincomplex
chaincomplex< euclidom > & createchaincomplex(chaincomplex< euclidom > &c, const gcomplex< cell, euclidom > &g, bool quiet=false)
Creates an algebraic chain complex based on the data from the given geometric cell complex.
Definition: gcomplex.h:1128
chomp::homology::ignorecomments
void ignorecomments(std::istream &in)
Ignores white characters (spaces, tabulators, CRs and LFs), as well as comments from the input text f...
Definition: textfile.h:385
chomp::homology::creategraph
gcomplex< cell, euclidom > & creategraph(const mvmap< cell, cell > &m, gcomplex< cell, euclidom > &graph)
Add a graph of a multivalued cell map to the cell complex.
Definition: gcomplex.h:1429
chomp::homology::boundarycoef
int boundarycoef(const tCellBase< coordtype > &q, int i)
Returns the i-th coefficient in the boundary of a cell.
Definition: cellbase.h:508
chomp
This namespace contains the entire CHomP library interface.
Definition: bitmaps.h:51
textfile.h
This file contains some useful functions related to the text input/output procedures.
Generated by doxygen 1.9.4