00001
00002
00003
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #ifndef _CHOMP_CUBES_CUBMAPS_H_
00037 #define _CHOMP_CUBES_CUBMAPS_H_
00038
00039 #include "chomp/system/config.h"
00040 #include "chomp/system/textfile.h"
00041 #include "chomp/cubes/pointset.h"
00042 #include "chomp/struct/integer.h"
00043 #include "chomp/struct/hashsets.h"
00044 #include "chomp/homology/gcomplex.h"
00045 #include "chomp/cubes/pointbas.h"
00046 #include "chomp/cubes/cube.h"
00047 #include "chomp/cubes/cell.h"
00048
00049 #include <iostream>
00050 #include <fstream>
00051 #include <cstdlib>
00052
00053 namespace chomp {
00054 namespace homology {
00055
00056
00057
00058
00059
00062
00063
00064 template <class euclidom, class tCell>
00065 std::ostream &writegenerators (std::ostream &out, const chain<euclidom> *hom,
00066 const chaincomplex<euclidom> &c, const gcomplex<tCell,euclidom> &g,
00067 const int *level, int xdim, int format = 0)
00068 {
00069 typedef typename tCell::CoordType coordType;
00070 bool firstlist = true;
00071 for (int d = 0; d <= c. dim (); ++ d)
00072 {
00073 if ((level && !level [d]) || hom [d]. empty ())
00074 continue;
00075 if (firstlist)
00076 firstlist = false;
00077 else if (format == 0)
00078 out << '\n';
00079 if (format == 0)
00080 {
00081 if (hom [d]. size () == 1)
00082 {
00083 out << "The generator of H_" << d <<
00084 " follows:" << '\n';
00085 }
00086 else
00087 {
00088 out << "The " << hom [d]. size () <<
00089 " generators of H_" << d <<
00090 " follow:" << '\n';
00091 }
00092 }
00093 const hashedset<tCell> &cset = g [d];
00094 for (int_t i = 0; i < hom [d]. size (); ++ i)
00095 {
00096 if (format == 0)
00097 {
00098 if (hom [d]. size () > 1)
00099 out << "generator " <<
00100 (i + 1) << '\n';
00101 }
00102 const chain<euclidom> &lst =
00103 c. gethomgen (d, hom [d]. num (i));
00104 for (int_t j = 0; j < lst. size (); ++ j)
00105 {
00106 coordType left [tCell::MaxDim];
00107 cset [lst. num (j)]. leftcoord (left);
00108 coordType right [tCell::MaxDim];
00109 cset [lst. num (j)]. rightcoord (right);
00110 int projdim = 0;
00111 for (int k = 0; k < xdim; ++ k)
00112 {
00113 if (left [k] != right [k])
00114 ++ projdim;
00115 }
00116 if (projdim != d)
00117 continue;
00118 if (format == 0)
00119 {
00120 out << lst. coef (j) << " * " <<
00121 tCell (left, right, xdim) <<
00122 '\n';
00123 }
00124 else
00125 {
00126 for (int k = 0; k < xdim; ++ k)
00127 out << left [k] << " ";
00128 for (int k = 0; k < xdim; ++ k)
00129 out << right [k] << " ";
00130 out << lst. coef (j) << " " <<
00131 (i + 1) << " " << d << '\n';
00132 }
00133 }
00134 }
00135 }
00136 return out;
00137 }
00138
00139
00140
00141
00142
00143
00147 template <class euclidom, class tCell, class tCube>
00148 int_t createimages (mvcellmap<tCell,euclidom,tCube> &m,
00149 const mvmap<tCube,tCube> &f1, const mvmap<tCube,tCube> &f2,
00150 const hashedset<tCube> &dom1, const hashedset<tCube> &dom2)
00151 {
00152 typedef typename tCell::CoordType coordType;
00153 int_t countimages = 0;
00154 coordType leftbound [tCell::MaxDim];
00155 coordType rightbound [tCell::MaxDim];
00156 coordType left [tCell::MaxDim];
00157 coordType right [tCell::MaxDim];
00158 for (int d = 0; d <= m. dim (); ++ d)
00159 {
00160 const hashedset<tCell> &cset = m. get (d);
00161 if (cset. empty ())
00162 continue;
00163 const int spacedim = cset [0]. spacedim ();
00164 for (int_t i = 0; i < cset. size (); ++ i)
00165 {
00166 cset [i]. leftcoord (left);
00167 cset [i]. rightcoord (right);
00168 for (int_t j = 0; j < spacedim; ++ j)
00169 {
00170
00171 if (right [j] < left [j])
00172 right [j] = left [j] + 1;
00173
00174 leftbound [j] = static_cast<coordType>
00175 (left [j] - (left [j] == right [j]));
00176 rightbound [j] = static_cast<coordType>
00177 (right [j] +
00178 (left [j] == right [j]));
00179 }
00180 tRectangle<coordType> r (leftbound, rightbound,
00181 spacedim);
00182 const coordType *c;
00183 while ((c = r. get ()) != NULL)
00184 {
00185 if (!tCube::PointBase::check (c, spacedim))
00186 continue;
00187 tCube q (c, spacedim);
00188 if (dom1. check (q))
00189 m. add (d, i, f1 (q));
00190 if (dom2. check (q))
00191 m. add (d, i, f2 (q));
00192 }
00193 countimages += m (cset [i]). size ();
00194 }
00195 }
00196 return countimages;
00197 }
00198
00200 template <class euclidom, class tCell, class tCube>
00201 inline int_t createimages (mvcellmap<tCell,euclidom,tCube> &m,
00202 const mvmap<tCube,tCube> &f, const hashedset<tCube> &dom)
00203 {
00204 mvmap<tCube,tCube> emptymap;
00205 hashedset<tCube> emptyset;
00206 return createimages (m, f, emptymap, dom, emptyset);
00207 }
00208
00211 template <class euclidom, class tCell, class tCube>
00212 int_t createimages (mvcellmap<tCell,euclidom,tCube> &m,
00213 const mvmap<tCube,tCube> &f1, const mvmap<tCube,tCube> &f2)
00214 {
00215 const hashedset<tCube> &dom1 = f1. getdomain ();
00216 const hashedset<tCube> &dom2 = f2. getdomain ();
00217 return createimages (m, f1, f2, dom1, dom2);
00218 }
00219
00221 template <class euclidom, class tCell, class tCube>
00222 inline int_t createimages (mvcellmap<tCell,euclidom,tCube> &m,
00223 const mvmap<tCube,tCube> &f)
00224 {
00225 mvmap<tCube,tCube> emptymap;
00226 return createimages (m, f, emptymap);
00227 }
00228
00229
00230
00231
00232
00233
00237 template <class euclidom, class tCell>
00238 void createprojection (const gcomplex<tCell,euclidom> &Fcompl,
00239 const gcomplex<tCell,euclidom> &Ycompl,
00240 chainmap<euclidom> &cmap,
00241 int offset, int outdim, int discarddim, int *level = NULL)
00242 {
00243 typedef typename tCell::CoordType coordType;
00244
00245 for (int d = 0; d <= Ycompl. dim (); ++ d)
00246 {
00247 if ((!level || level [d]) && (Fcompl. dim () >= d))
00248 {
00249
00250 const hashedset<tCell> &Fset = Fcompl [d];
00251 if (Fset. empty ())
00252 continue;
00253 const hashedset<tCell> &Yset = Ycompl [d];
00254 if (Yset. empty ())
00255 continue;
00256
00257
00258 for (int_t i = 0; i < Fset. size (); ++ i)
00259 {
00260
00261 const tCell &Fcell = Fset [i];
00262 coordType left [tCell::MaxDim];
00263 Fcell. leftcoord (left);
00264 coordType right [tCell::MaxDim];
00265 Fcell. rightcoord (right);
00266
00267
00268
00269 register int j;
00270 for (j = 0; j < offset; ++ j)
00271 if (left [j] != right [j])
00272 {
00273 j = offset + 33;
00274 break;
00275 }
00276 if (j > offset)
00277 continue;
00278 for (j = 0; j < discarddim; ++ j)
00279 if (left [offset + outdim + j] !=
00280 right [offset + outdim + j])
00281 {
00282 j = discarddim + 33;
00283 break;
00284 }
00285 if (j > discarddim)
00286 continue;
00287
00288
00289 if (!(tCell::PointBase::check
00290 (left + offset, outdim)))
00291 continue;
00292 if (!(tCell::PointBase::check
00293 (right + offset, outdim)))
00294 continue;
00295
00296
00297 tCell projected (Fcell, offset, outdim);
00298
00299
00300 int_t nr = Yset. getnumber (projected);
00301
00302
00303 if (nr < 0)
00304 continue;
00305
00306
00307 euclidom e;
00308 e = 1;
00309 cmap. add (d, nr, i, e);
00310 }
00311 }
00312 }
00313 return;
00314 }
00315
00320 template <class euclidom, class tCell>
00321 void project (const gcomplex<tCell,euclidom> &c,
00322 gcomplex<tCell,euclidom> &img,
00323 const gcomplex<tCell,euclidom> &only,
00324 int offset, int outdim, int discarddim, const int *level,
00325 bool watchforimages)
00326 {
00327 typedef typename tCell::CoordType coordType;
00328
00329
00330 for (int d = 0; d <= c. dim (); ++ d)
00331 {
00332 if ((!level || level [d]) && (c. dim () >= d))
00333 {
00334
00335 const hashedset<tCell> &cset = c [d];
00336 if (cset. empty ())
00337 continue;
00338
00339
00340 for (int_t i = 0; i < cset. size (); ++ i)
00341 {
00342
00343 const tCell &thecell = cset [i];
00344 coordType left [tCell::MaxDim];
00345 thecell. leftcoord (left);
00346 coordType right [tCell::MaxDim];
00347 thecell. rightcoord (right);
00348
00349
00350
00351 int j;
00352 for (j = 0; j < offset; ++ j)
00353 if (left [j] != right [j])
00354 {
00355 j = offset + 33;
00356 break;
00357 }
00358 if (j > offset)
00359 continue;
00360 for (j = 0; j < discarddim; ++ j)
00361 if (left [offset + outdim + j] !=
00362 right [offset + outdim + j])
00363 {
00364 j = discarddim + 33;
00365 break;
00366 }
00367 if (j > discarddim)
00368 continue;
00369
00370
00371 if (!(tCell::PointBase::check
00372 (left + offset, outdim))
00373 || !(tCell::PointBase::check
00374 (right + offset, outdim)))
00375 {
00376 if (watchforimages)
00377 throw "Inclusion undefined!";
00378 else
00379 continue;
00380 }
00381 tCell projected (left + offset,
00382 right + offset, outdim);
00383 if ((only. dim () >= 0) &&
00384 !only. check (projected))
00385 {
00386 if (watchforimages)
00387 throw "Inclusion undefined.";
00388 else
00389 continue;
00390 }
00391 img. add (projected);
00392 }
00393 }
00394 }
00395 return;
00396 }
00397
00398
00399
00400
00401
00402
00406 template <class tCube>
00407 std::istream &readdomain (std::istream &in, hashedset<tCube> &dom)
00408 {
00409 ignorecomments (in);
00410 if ((in. peek () != 'S') && (in. peek () != 's'))
00411 {
00412 mvmap<tCube,tCube> Fdummy;
00413 return readdomain (in, dom, Fdummy);
00414 }
00415
00416 while (in. peek () != EOF)
00417 {
00418 if (in. peek () == '[')
00419 {
00420 in. get ();
00421 tCube q;
00422 in >> q;
00423 if (!in)
00424 throw "Can't read the file.";
00425 dom. add (q);
00426 }
00427 ignoreline (in);
00428 ignorecomments (in);
00429 }
00430
00431 return in;
00432 }
00433
00437 template <class tCube>
00438 std::istream &readimage (std::istream &in, hashedset<tCube> &img)
00439 {
00440 ignorecomments (in);
00441 if ((in. peek () != 'S') && (in. peek () != 's'))
00442 {
00443 mvmap<tCube,tCube> Fdummy;
00444 return readimage (in, img, Fdummy);
00445 }
00446
00447 typename tCube::CoordType left [tCube::MaxDim];
00448 typename tCube::CoordType right [tCube::MaxDim];
00449 while (in. peek () != EOF)
00450 {
00451 if (in. peek () == '[')
00452 {
00453 in. get ();
00454 tCube dummy;
00455 in >> dummy;
00456 in >> dummy;
00457 ignorecomments (in);
00458 in. get ();
00459 ignorecomments (in);
00460 in. get ();
00461 tCube q1, q2;
00462 in >> q1 >> q2;
00463 if (!in)
00464 throw "Can't read the file.";
00465 ignorecomments (in);
00466 in. get ();
00467 ignorecomments (in);
00468 int dim = q1. dim ();
00469 q1. coord (left);
00470 q2. coord (right);
00471 tRectangle<typename tCube::CoordType> r
00472 (left, right, dim);
00473 const typename tCube::CoordType *c;
00474 while ((c = r. get ()) != NULL)
00475 img. add (tCube (c, dim));
00476 }
00477 else
00478 ignoreline (in);
00479 ignorecomments (in);
00480 }
00481
00482 return in;
00483 }
00484
00488 template <class tCube>
00489 std::istream &readimage (std::istream &in, const hashedset<tCube> &dom,
00490 hashedset<tCube> &img)
00491 {
00492 ignorecomments (in);
00493
00494 if ((in. peek () != 'S') && (in. peek () != 's'))
00495 {
00496 while (in. peek () != EOF)
00497 {
00498 if ((closingparenthesis (in. peek ()) == EOF) &&
00499 !std::isdigit (in. peek ()))
00500 {
00501 ignoreline (in);
00502 ignorecomments (in);
00503 continue;
00504 }
00505 tCube q;
00506 in >> q;
00507 bool ignore = !dom. check (q);
00508 ignorecomments (in);
00509 while (in. peek () == '-')
00510 in. get ();
00511 in. get ();
00512 ignorecomments (in);
00513 int opening = in. get ();
00514 int closing = closingparenthesis (opening);
00515 if (closing == EOF)
00516 throw "An opening brace '{' expected.";
00517 while ((in. peek () != EOF) &&
00518 (in. peek () != closing))
00519 {
00520 if (ignore)
00521 {
00522 if (in. get () == opening)
00523 {
00524 while ((in. peek () != EOF) &&
00525 (in. peek () != closing))
00526 {
00527 in. get ();
00528 }
00529 in. get ();
00530 }
00531 }
00532 else
00533 {
00534 in >> q;
00535 img. add (q);
00536 ignorecomments (in);
00537 if (in. peek () == ',')
00538 {
00539 in. get ();
00540 ignorecomments (in);
00541 }
00542 }
00543 }
00544
00545 ignoreline (in);
00546 ignorecomments (in);
00547 }
00548 }
00549
00550 typename tCube::CoordType left [tCube::MaxDim];
00551 typename tCube::CoordType right [tCube::MaxDim];
00552 while (in. peek () != EOF)
00553 {
00554 if (in. peek () == '[')
00555 {
00556 in. get ();
00557 tCube domcube1, domcube2;
00558 in >> domcube1;
00559 if (!dom. check (domcube1))
00560 {
00561 ignoreline (in);
00562 ignorecomments (in);
00563 continue;
00564 }
00565 in >> domcube2;
00566 ignorecomments (in);
00567 in. get ();
00568 ignorecomments (in);
00569 in. get ();
00570 tCube q1, q2;
00571 in >> q1 >> q2;
00572 if (!in)
00573 throw "Can't read the file.";
00574 ignorecomments (in);
00575 in. get ();
00576 ignorecomments (in);
00577 if (dom. check (domcube1))
00578 {
00579 int dim = q1. dim ();
00580 q1. coord (left);
00581 q2. coord (right);
00582 tRectangle<typename tCube::CoordType> r
00583 (left, right, dim);
00584 const typename tCube::CoordType *c;
00585 while ((c = r. get ()) != NULL)
00586 img. add (tCube (c, dim));
00587 }
00588 }
00589 else
00590 ignoreline (in);
00591 ignorecomments (in);
00592 }
00593
00594 return in;
00595 }
00596
00598 template <class tCube>
00599 std::istream &readselective (std::istream &in, const hashedset<tCube> &dom1,
00600 const hashedset<tCube> &dom2, mvmap<tCube,tCube> &m)
00601 {
00602 if (dom1. empty () && dom2. empty ())
00603 {
00604 sout << "Warning: The domain of the map is empty.\n";
00605 return in;
00606 }
00607
00608 ignorecomments (in);
00609 if ((in. peek () != 'S') && (in. peek () != 's'))
00610 return readselective (in, m, dom1, dom2);
00611
00612 typename tCube::CoordType left [tCube::MaxDim];
00613 typename tCube::CoordType right [tCube::MaxDim];
00614 while (in. peek () != EOF)
00615 {
00616 if (in. peek () == '[')
00617 {
00618 in. get ();
00619 tCube domcube;
00620 in >> domcube;
00621 int dim = domcube. dim ();
00622 if (dom1. check (domcube) || dom2. check (domcube))
00623 {
00624 tCube q1, q2;
00625 in >> q1;
00626 ignorecomments (in);
00627 in. get ();
00628 ignorecomments (in);
00629 in. get ();
00630 in >> q1 >> q2;
00631 if (!in)
00632 throw "Can't read the file.";
00633 hashedset<tCube> &img = m [domcube];
00634 q1. coord (left);
00635 q2. coord (right);
00636 tRectangle<typename tCube::CoordType> r
00637 (left, right, dim);
00638 const typename tCube::CoordType *c;
00639 while ((c = r. get ()) != NULL)
00640 img. add (tCube (c, dim));
00641 }
00642 }
00643 ignoreline (in);
00644 ignorecomments (in);
00645 }
00646
00647 return in;
00648 }
00649
00653 template <class tCube>
00654 inline std::istream &readselective (std::istream &in,
00655 const hashedset<tCube> &dom, mvmap<tCube,tCube> &m)
00656 {
00657 hashedset<tCube> empty;
00658 return readselective (in, dom, empty, m);
00659 }
00660
00661
00662 }
00663 }
00664
00665 #endif // _CHOMP_CUBES_CUBMAPS_H_
00666
00668