cubmaps.h

Go to the documentation of this file.

// Copyright (C) 1997-2010 by Pawel Pilarczyk.
//
// This file is part of the Homology Library.  This library is free software;
// you can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation;
// either version 2 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this software; see the file "license.txt".  If not, write to the
// Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
// MA 02111-1307, USA.

// Started in January 2002. Last revision: October 3, 2008.


#ifndef _CHOMP_CUBES_CUBMAPS_H_
#define _CHOMP_CUBES_CUBMAPS_H_

#include "chomp/system/config.h"
#include "chomp/system/textfile.h"
#include "chomp/cubes/pointset.h"
#include "chomp/struct/integer.h"
#include "chomp/struct/hashsets.h"
#include "chomp/homology/gcomplex.h"
#include "chomp/cubes/pointbas.h"
#include "chomp/cubes/cube.h"
#include "chomp/cubes/cell.h"

#include <iostream>
#include <fstream>
#include <cstdlib>

namespace chomp {
namespace homology {

// --------------------------------------------------
// ------------------- generators -------------------
// --------------------------------------------------

// Note: The main loop(s) of this function are copied from the function
// "createprojection".
template <class euclidom, class tCell>
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)
{
        typedef typename tCell::CoordType coordType;
        bool firstlist = true;
        for (int d = 0; d <= c. dim (); ++ d)
        {
                if ((level && !level [d]) || hom [d]. empty ())
                        continue;
                if (firstlist)
                        firstlist = false;
                else if (format == 0)
                        out << '\n';
                if (format == 0)
                {
                        if (hom [d]. size () == 1)
                        {
                                out << "The generator of H_" << d <<
                                        " follows:" << '\n';
                        }
                        else
                        {
                                out << "The " << hom [d]. size () <<
                                        " generators of H_" << d <<
                                        " follow:" << '\n';
                        }
                }
                const hashedset<tCell> &cset = g [d];
                for (int_t i = 0; i < hom [d]. size (); ++ i)
                {
                        if (format == 0)
                        {
                                if (hom [d]. size () > 1)
                                        out << "generator " <<
                                                (i + 1) << '\n';
                        }
                        const chain<euclidom> &lst =
                                c. gethomgen (d, hom [d]. num (i));
                        for (int_t j = 0; j < lst. size (); ++ j)
                        {
                                coordType left [tCell::MaxDim];
                                cset [lst. num (j)]. leftcoord (left);
                                coordType right [tCell::MaxDim];
                                cset [lst. num (j)]. rightcoord (right);
                                int projdim = 0;
                                for (int k = 0; k < xdim; ++ k)
                                {
                                        if (left [k] != right [k])
                                                ++ projdim;
                                }
                                if (projdim != d)
                                        continue;
                                if (format == 0)
                                {
                                        out << lst. coef (j) << " * " <<
                                                tCell (left, right, xdim) <<
                                                '\n';
                                }
                                else
                                {
                                        for (int k = 0; k < xdim; ++ k)
                                                out << left [k] << " ";
                                        for (int k = 0; k < xdim; ++ k)
                                                out << right [k] << " ";
                                        out << lst. coef (j) << " " <<
                                                (i + 1) << " " << d << '\n';
                                }
                        }
                }
        }
        return out;
} /* writegenerators */


// --------------------------------------------------
// ------------------- map images -------------------
// --------------------------------------------------

template <class euclidom, class tCell, class tCube>
int_t createimages (mvcellmap<tCell,euclidom,tCube> &m,
        const mvmap<tCube,tCube> &f1, const mvmap<tCube,tCube> &f2,
        const hashedset<tCube> &dom1, const hashedset<tCube> &dom2)
{
        typedef typename tCell::CoordType coordType;
        int_t countimages = 0;
        coordType leftbound [tCell::MaxDim];
        coordType rightbound [tCell::MaxDim];
        coordType left [tCell::MaxDim];
        coordType right [tCell::MaxDim];
        for (int d = 0; d <= m. dim (); ++ d)
        {
                const hashedset<tCell> &cset = m. get (d);
                if (cset. empty ())
                        continue;
                const int spacedim = cset [0]. spacedim ();
                for (int_t i = 0; i < cset. size (); ++ i)
                {
                        cset [i]. leftcoord (left);
                        cset [i]. rightcoord (right);
                        for (int_t j = 0; j < spacedim; ++ j)
                        {
                                // compensate for space wrapping
                                if (right [j] < left [j])
                                        right [j] = left [j] + 1;
                                // compute the bounds
                                leftbound [j] = static_cast<coordType>
                                        (left [j] - (left [j] == right [j]));
                                rightbound [j] = static_cast<coordType>
                                        (right [j] +
                                        (left [j] == right [j]));
                        }
                        tRectangle<coordType> r (leftbound, rightbound,
                                spacedim);
                        const coordType *c;
                        while ((c = r. get ()) != NULL)
                        {
                                if (!tCube::PointBase::check (c, spacedim))
                                        continue;
                                tCube q (c, spacedim);
                                if (dom1. check (q))
                                        m. add (d, i, f1 (q));
                                if (dom2. check (q))
                                        m. add (d, i, f2 (q));
                        }
                        countimages += m (cset [i]). size ();
                }
        }
        return countimages;
} /* createimages */

template <class euclidom, class tCell, class tCube>
inline int_t createimages (mvcellmap<tCell,euclidom,tCube> &m,
        const mvmap<tCube,tCube> &f, const hashedset<tCube> &dom)
{
        mvmap<tCube,tCube> emptymap;
        hashedset<tCube> emptyset;
        return createimages (m, f, emptymap, dom, emptyset);
} /* createimages */

template <class euclidom, class tCell, class tCube>
int_t createimages (mvcellmap<tCell,euclidom,tCube> &m,
        const mvmap<tCube,tCube> &f1, const mvmap<tCube,tCube> &f2)
{
        const hashedset<tCube> &dom1 = f1. getdomain ();
        const hashedset<tCube> &dom2 = f2. getdomain ();
        return createimages (m, f1, f2, dom1, dom2);
} /* createimages */

template <class euclidom, class tCell, class tCube>
inline int_t createimages (mvcellmap<tCell,euclidom,tCube> &m,
        const mvmap<tCube,tCube> &f)
{
        mvmap<tCube,tCube> emptymap;
        return createimages (m, f, emptymap);
} /* createimages */


// --------------------------------------------------
// ------------------ projections -------------------
// --------------------------------------------------

template <class euclidom, class tCell>
void createprojection (const gcomplex<tCell,euclidom> &Fcompl,
        const gcomplex<tCell,euclidom> &Ycompl,
        chainmap<euclidom> &cmap,
        int offset, int outdim, int discarddim, int *level = NULL)
{
        typedef typename tCell::CoordType coordType;
        // go through the list of all the dimensions which are of concern
        for (int d = 0; d <= Ycompl. dim (); ++ d)
        {
                if ((!level || level [d]) && (Fcompl. dim () >= d))
                {
                        // take sets of cells of this dimension
                        const hashedset<tCell> &Fset = Fcompl [d];
                        if (Fset. empty ())
                                continue;
                        const hashedset<tCell> &Yset = Ycompl [d];
                        if (Yset. empty ())
                                continue;

                        // go through the list of cells in Fcompl of dim. d
                        for (int_t i = 0; i < Fset. size (); ++ i)
                        {
                                // get this cell and its coordinates
                                const tCell &Fcell = Fset [i];
                                coordType left [tCell::MaxDim];
                                Fcell. leftcoord (left);
                                coordType right [tCell::MaxDim];
                                Fcell. rightcoord (right);

                                // check if this cell has no width in the
                                // directions that are discarded
                                register int j;
                                for (j = 0; j < offset; ++ j)
                                        if (left [j] != right [j])
                                        {
                                                j = offset + 33;
                                                break;
                                        }
                                if (j > offset)
                                        continue;
                                for (j = 0; j < discarddim; ++ j)
                                        if (left [offset + outdim + j] !=
                                                right [offset + outdim + j])
                                        {
                                                j = discarddim + 33;
                                                break;
                                        }
                                if (j > discarddim)
                                        continue;

                                // create the projected cell
                                if (!(tCell::PointBase::check
                                        (left + offset, outdim)))
                                        continue;
                                if (!(tCell::PointBase::check
                                        (right + offset, outdim)))
                                        continue;
                        //      tCell projected (left + offset,
                        //              right + offset, outdim);
                                tCell projected (Fcell, offset, outdim);

                                // find its number in Y
                                int_t nr = Yset. getnumber (projected);

                                // if not found, discard it
                                if (nr < 0)
                                        continue;

                                // add the pair to the projection map
                                euclidom e;
                                e = 1;
                                cmap. add (d, nr, i, e);
                        }
                }
        }
        return;
} /* createprojection */

template <class euclidom, class tCell>
void project (const gcomplex<tCell,euclidom> &c,
        gcomplex<tCell,euclidom> &img,
        const gcomplex<tCell,euclidom> &only,
        int offset, int outdim, int discarddim, const int *level,
        bool watchforimages)
{
        typedef typename tCell::CoordType coordType;

        // go through the list of all the dimensions which are of concern
        for (int d = 0; d <= c. dim (); ++ d)
        {
                if ((!level || level [d]) && (c. dim () >= d))
                {
                        // take sets of cells of this dimension
                        const hashedset<tCell> &cset = c [d];
                        if (cset. empty ())
                                continue;

                        // go through the list of cells in c of dimension d
                        for (int_t i = 0; i < cset. size (); ++ i)
                        {
                                // get this cell and its coordinates
                                const tCell &thecell = cset [i];
                                coordType left [tCell::MaxDim];
                                thecell. leftcoord (left);
                                coordType right [tCell::MaxDim];
                                thecell. rightcoord (right);

                                // check if this cell has no width in the
                                // directions that are discarded
                                int j;
                                for (j = 0; j < offset; ++ j)
                                        if (left [j] != right [j])
                                        {
                                                j = offset + 33;
                                                break;
                                        }
                                if (j > offset)
                                        continue;
                                for (j = 0; j < discarddim; ++ j)
                                        if (left [offset + outdim + j] !=
                                                right [offset + outdim + j])
                                        {
                                                j = discarddim + 33;
                                                break;
                                        }
                                if (j > discarddim)
                                        continue;

                                // add the projected cell to the complex
                                if (!(tCell::PointBase::check
                                        (left + offset, outdim))
                                        || !(tCell::PointBase::check
                                        (right + offset, outdim)))
                                {
                                        if (watchforimages)
                                                throw "Inclusion undefined!";
                                        else
                                                continue;
                                }
                                tCell projected (left + offset,
                                        right + offset, outdim);
                                if ((only. dim () >= 0) &&
                                        !only. check (projected))
                                {
                                        if (watchforimages)
                                                throw "Inclusion undefined.";
                                        else
                                                continue;
                                }
                                img. add (projected);
                        }
                }
        }
        return;
} /* project */


// --------------------------------------------------
// ------------------- read cubes -------------------
// --------------------------------------------------

template <class tCube>
std::istream &readdomain (std::istream &in, hashedset<tCube> &dom)
{
        ignorecomments (in);
        if ((in. peek () != 'S') && (in. peek () != 's'))
        {
                mvmap<tCube,tCube> Fdummy;
                return readdomain (in, dom, Fdummy);
        }

        while (in. peek () != EOF)
        {
                if (in. peek () == '[')
                {
                        in. get ();
                        tCube q;
                        in >> q;
                        if (!in)
                                throw "Can't read the file.";
                        dom. add (q);
                }
                ignoreline (in);
                ignorecomments (in);
        }
        
        return in;
} /* readdomain */

template <class tCube>
std::istream &readimage (std::istream &in, hashedset<tCube> &img)
{
        ignorecomments (in);
        if ((in. peek () != 'S') && (in. peek () != 's'))
        {
                mvmap<tCube,tCube> Fdummy;
                return readimage (in, img, Fdummy);
        }

        typename tCube::CoordType left [tCube::MaxDim];
        typename tCube::CoordType right [tCube::MaxDim];
        while (in. peek () != EOF)
        {
                if (in. peek () == '[')
                {
                        in. get ();
                        tCube dummy;
                        in >> dummy;
                        in >> dummy;
                        ignorecomments (in);
                        in. get ();
                        ignorecomments (in);
                        in. get ();
                        tCube q1, q2;
                        in >> q1 >> q2;
                        if (!in)
                                throw "Can't read the file.";
                        ignorecomments (in);
                        in. get ();
                        ignorecomments (in);
                        int dim = q1. dim ();
                        q1. coord (left);
                        q2. coord (right);
                        tRectangle<typename tCube::CoordType> r
                                (left, right, dim);
                        const typename tCube::CoordType *c;
                        while ((c = r. get ()) != NULL)
                                img. add (tCube (c, dim));
                }
                else
                        ignoreline (in);
                ignorecomments (in);
        }

        return in;
} /* readimage */

template <class tCube>
std::istream &readimage (std::istream &in, const hashedset<tCube> &dom,
        hashedset<tCube> &img)
{
        ignorecomments (in);
        // the general format: [x,y,z] -> {[a,b,c], [d,e,f]}
        if ((in. peek () != 'S') && (in. peek () != 's'))
        {
                while (in. peek () != EOF)
                {
                        if ((closingparenthesis (in. peek ()) == EOF) &&
                                !std::isdigit (in. peek ()))
                        {
                                ignoreline (in);
                                ignorecomments (in);
                                continue;
                        }
                        tCube q;
                        in >> q;
                        bool ignore = !dom. check (q);
                        ignorecomments (in);
                        while (in. peek () == '-')
                                in. get ();
                        in. get (); // '>'
                        ignorecomments (in);
                        int opening = in. get ();
                        int closing = closingparenthesis (opening);
                        if (closing == EOF)
                                throw "An opening brace '{' expected.";
                        while ((in. peek () != EOF) &&
                                (in. peek () != closing))
                        {
                                if (ignore)
                                {
                                        if (in. get () == opening)
                                        {
                                                while ((in. peek () != EOF) &&
                                                        (in. peek () != closing))
                                                {
                                                        in. get ();
                                                }
                                                in. get (); // '}'
                                        }
                                }
                                else
                                {
                                        in >> q;
                                        img. add (q);
                                        ignorecomments (in);
                                        if (in. peek () == ',')
                                        {
                                                in. get ();
                                                ignorecomments (in);
                                        }
                                }
                        }
                //      in. get (); // '}'
                        ignoreline (in);
                        ignorecomments (in);
                }
        }

        typename tCube::CoordType left [tCube::MaxDim];
        typename tCube::CoordType right [tCube::MaxDim];
        while (in. peek () != EOF)
        {
                if (in. peek () == '[')
                {
                        in. get ();
                        tCube domcube1, domcube2;
                        in >> domcube1;
                        if (!dom. check (domcube1))
                        {
                                ignoreline (in);
                                ignorecomments (in);
                                continue;
                        }
                        in >> domcube2;
                        ignorecomments (in);
                        in. get ();
                        ignorecomments (in);
                        in. get ();
                        tCube q1, q2;
                        in >> q1 >> q2;
                        if (!in)
                                throw "Can't read the file.";
                        ignorecomments (in);
                        in. get ();
                        ignorecomments (in);
                        if (dom. check (domcube1))
                        {
                                int dim = q1. dim ();
                                q1. coord (left);
                                q2. coord (right);
                                tRectangle<typename tCube::CoordType> r
                                        (left, right, dim);
                                const typename tCube::CoordType *c;
                                while ((c = r. get ()) != NULL)
                                        img. add (tCube (c, dim));
                        }
                }
                else
                        ignoreline (in);
                ignorecomments (in);
        }

        return in;
} /* readimage */

template <class tCube>
std::istream &readselective (std::istream &in, const hashedset<tCube> &dom1,
        const hashedset<tCube> &dom2, mvmap<tCube,tCube> &m)
{
        if (dom1. empty () && dom2. empty ())
        {
                sout << "Warning: The domain of the map is empty.\n";
                return in;
        }

        ignorecomments (in);
        if ((in. peek () != 'S') && (in. peek () != 's'))
                return readselective (in, m, dom1, dom2);

        typename tCube::CoordType left [tCube::MaxDim];
        typename tCube::CoordType right [tCube::MaxDim];
        while (in. peek () != EOF)
        {
                if (in. peek () == '[')
                {
                        in. get ();
                        tCube domcube;
                        in >> domcube;
                        int dim = domcube. dim ();
                        if (dom1. check (domcube) || dom2. check (domcube))
                        {
                                tCube q1, q2;
                                in >> q1; // (ignored)
                                ignorecomments (in);
                                in. get ();
                                ignorecomments (in);
                                in. get ();
                                in >> q1 >> q2;
                                if (!in)
                                        throw "Can't read the file.";
                                hashedset<tCube> &img = m [domcube];
                                q1. coord (left);
                                q2. coord (right);
                                tRectangle<typename tCube::CoordType> r
                                        (left, right, dim);
                                const typename tCube::CoordType *c;
                                while ((c = r. get ()) != NULL)
                                        img. add (tCube (c, dim));
                        }
                }
                ignoreline (in);
                ignorecomments (in);
        }

        return in;
} /* readselective */

template <class tCube>
inline std::istream &readselective (std::istream &in,
        const hashedset<tCube> &dom, mvmap<tCube,tCube> &m)
{
        hashedset<tCube> empty;
        return readselective (in, dom, empty, m);
} /* readselective */


} // namespace homology
} // namespace chomp

#endif // _CHOMP_CUBES_CUBMAPS_H_