cubmaps.h

Go to the documentation of this file.
00001 
00002 
00003 
00015 
00016 // Copyright (C) 1997-2010 by Pawel Pilarczyk.
00017 //
00018 // This file is part of the Homology Library.  This library is free software;
00019 // you can redistribute it and/or modify it under the terms of the GNU
00020 // General Public License as published by the Free Software Foundation;
00021 // either version 2 of the License, or (at your option) any later version.
00022 //
00023 // This library is distributed in the hope that it will be useful,
00024 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00025 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00026 // GNU General Public License for more details.
00027 //
00028 // You should have received a copy of the GNU General Public License along
00029 // with this software; see the file "license.txt".  If not, write to the
00030 // Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
00031 // MA 02111-1307, USA.
00032 
00033 // Started in January 2002. Last revision: October 3, 2008.
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 // ------------------- generators -------------------
00058 // --------------------------------------------------
00059 
00062 // Note: The main loop(s) of this function are copied from the function
00063 // "createprojection".
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 } /* writegenerators */
00138 
00139 
00140 // --------------------------------------------------
00141 // ------------------- map images -------------------
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                                 // compensate for space wrapping
00171                                 if (right [j] < left [j])
00172                                         right [j] = left [j] + 1;
00173                                 // compute the bounds
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 } /* createimages */
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 } /* createimages */
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 } /* createimages */
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 } /* createimages */
00228 
00229 
00230 // --------------------------------------------------
00231 // ------------------ projections -------------------
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         // go through the list of all the dimensions which are of concern
00245         for (int d = 0; d <= Ycompl. dim (); ++ d)
00246         {
00247                 if ((!level || level [d]) && (Fcompl. dim () >= d))
00248                 {
00249                         // take sets of cells of this dimension
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                         // go through the list of cells in Fcompl of dim. d
00258                         for (int_t i = 0; i < Fset. size (); ++ i)
00259                         {
00260                                 // get this cell and its coordinates
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                                 // check if this cell has no width in the
00268                                 // directions that are discarded
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                                 // create the projected cell
00289                                 if (!(tCell::PointBase::check
00290                                         (left + offset, outdim)))
00291                                         continue;
00292                                 if (!(tCell::PointBase::check
00293                                         (right + offset, outdim)))
00294                                         continue;
00295                         //      tCell projected (left + offset,
00296                         //              right + offset, outdim);
00297                                 tCell projected (Fcell, offset, outdim);
00298 
00299                                 // find its number in Y
00300                                 int_t nr = Yset. getnumber (projected);
00301 
00302                                 // if not found, discard it
00303                                 if (nr < 0)
00304                                         continue;
00305 
00306                                 // add the pair to the projection map
00307                                 euclidom e;
00308                                 e = 1;
00309                                 cmap. add (d, nr, i, e);
00310                         }
00311                 }
00312         }
00313         return;
00314 } /* createprojection */
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         // go through the list of all the dimensions which are of concern
00330         for (int d = 0; d <= c. dim (); ++ d)
00331         {
00332                 if ((!level || level [d]) && (c. dim () >= d))
00333                 {
00334                         // take sets of cells of this dimension
00335                         const hashedset<tCell> &cset = c [d];
00336                         if (cset. empty ())
00337                                 continue;
00338 
00339                         // go through the list of cells in c of dimension d
00340                         for (int_t i = 0; i < cset. size (); ++ i)
00341                         {
00342                                 // get this cell and its coordinates
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                                 // check if this cell has no width in the
00350                                 // directions that are discarded
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                                 // add the projected cell to the complex
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 } /* project */
00397 
00398 
00399 // --------------------------------------------------
00400 // ------------------- read cubes -------------------
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 } /* readdomain */
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 } /* readimage */
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         // the general format: [x,y,z] -> {[a,b,c], [d,e,f]}
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                 //      in. get (); // '}'
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 } /* readimage */
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; // (ignored)
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 } /* readselective */
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 } /* readselective */
00660 
00661 
00662 } // namespace homology
00663 } // namespace chomp
00664 
00665 #endif // _CHOMP_CUBES_CUBMAPS_H_
00666 
00668