cubacycl.h

Go to the documentation of this file.

 
 
 
// Copyright (C) 1997-2020 by Pawel Pilarczyk.
//
// This file is part of my research software package. This 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 3 of the License, or (at your option) any later version.
//
// This software 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,
// please, see <https://www.gnu.org/licenses/>.
 
// Started in January 2002. Last revision: February 2, 2018.
 
 
#ifndef _CHOMP_HOMOLOGY_CUBACYCL_H_
#define _CHOMP_HOMOLOGY_CUBACYCL_H_
 
#include "chomp/system/config.h"
#include "chomp/system/textfile.h"
#include "chomp/cubes/pointset.h"
#include "chomp/homology/chains.h"
#include "chomp/struct/bitfield.h"
#include "chomp/struct/integer.h"
#include "chomp/struct/hashsets.h"
#include "chomp/struct/setunion.h"
#include "chomp/struct/autoarray.h"
#include "chomp/homology/gcomplex.h"
#include "chomp/cubes/pointbas.h"
#include "chomp/cubes/cube.h"
#include "chomp/cubes/cell.h"
#include "chomp/cubes/cubmaps.h"
#include "chomp/cubes/neighbor.h"
#include "chomp/homology/bddacycl.h"
#include "chomp/homology/tabulate.h"
#include "chomp/homology/known.h"
 
#include <iostream>
#include <fstream>
#include <cstdlib>
 
namespace chomp {
namespace homology {
 
// --------------------------------------------------
// ---------- acyclicity of neighborhoods -----------
// --------------------------------------------------
 
inline bool acyclic (int dim, BitField &b)
{
        // look for the answer in the tabulated data
        const char *table = tabulated [dim];
        if (table)
        {
                return Tabulated::get (table,
                        bits2int (b, getmaxneighbors (dim)));
        }
 
        // look for the answer among the known combinations
        SetOfBitFields *known = knownbits [dim];
        int answerKnown = known ? known -> check (b) : -1;
 
        // if the answer is known then return it
        if (answerKnown >= 0)
                return (!!answerKnown);
 
        // create a model cube for building the neighborhood
        typedef tCubeBase<short> tCube;
        if (dim > tCube::MaxDim)
                throw "Cube dimension too high for acyclicity verification.";
        short c [tCube::MaxDim];
        for (int i = 0; i < dim; ++ i)
                c [i] = 0;
        tCube qzero (c, dim);
 
        // find out whether the set of neighbors is acyclic or not
        typedef tCube::CellType tCell;
        gcomplex<tCell,integer> neighbors;
        addneighbors (qzero, b, neighbors, true);
        bool answerBool (acyclic (neighbors));
 
        // add the answer to the list of known ones
        if (known)
                known -> add (b, answerBool);
 
        return answerBool;
} /* acyclic */
 
template <class tCube>
inline bool acyclic (const hashedset<tCube> &cset)
{
        // the empty set is not acyclic
        if (cset. empty ())
                return false;
 
        // a singleton is acyclic
        if (cset. size () == 1)
                return true;
 
        // reduce the set and see if there is only one cube remaining
        hashedset<tCube> emptycubes;
        hashedset<tCube> theset (cset);
        cubreducequiet (emptycubes, theset, emptycubes); // !!!
        if (theset. size () == 1)
                return true;
 
        // create a cubical complex from this set of cubes
        gcomplex<typename tCube::CellType,integer> ccompl;
        int_t size = cset. size ();
        for (int_t i = 0; i < size; ++ i)
                ccompl. add (typename tCube::CellType (cset [i]));
 
        // check whether this geometric complex is acyclic or not
        return acyclic (ccompl);
} /* acyclic */
 
template <class tCube, class tCubeSet1, class tCubeSet2>
inline bool acyclic (const tCube &q, int dim, const tCubeSet1 &cset,
        BitField *b, int_t maxneighbors, tCubeSet2 *neighbors)
{
        // use a binary decision diagram code if possible
        if (dim <= MaxBddDim)
                return bddacyclic (q, dim, cset, *b);
 
        // get the neighbors of the cube
        int_t n = getneighbors (q, b, cset, neighbors, 0);
 
        // if the answer is obvious from the number of neighbors, return it
        if ((n == 0) || (n == maxneighbors))
                return false;
        if (n == 1)
                return true;
 
        // return the information on whether this set of neighbors is acyclic
        return acyclic (dim, *b);
} /* acyclic */
 
template <class tCube, class tCubeSet>
inline bool acyclic (const tCube &q, int dim, const tCubeSet &cset,
        BitField *b, int_t maxneighbors)
{
        hashedset<tCube> *neighbors = 0;
        return acyclic (q, dim, cset, b, maxneighbors, neighbors);
} /* acyclic */
 
template <class tCube, class tSet1, class tSet2>
bool acyclic_rel (const tCube &q, int dim, const tSet1 &cset,
        const tSet2 &other, BitField *b, int_t maxneighbors,
        hashedset<tCube> *neighbors_main, hashedset<tCube> *neighbors_other)
{
        // use a binary decision diagram code if possible
        if (dim <= MaxBddDim)
        {
                if (!getneighbors (q, b, cset, 1))
                        return true;
                if (!bddacyclic (q, dim, other, *b))
                        return false;
                return bddacyclic (q, dim, makesetunion (cset, other), *b);
        }
 
        // get the neighbors from the other set
        int_t nother = getneighbors (q, b, other, neighbors_other, 0);
 
        // verify if this set of neighbors is acyclic
        bool otheracyclic = (nother < maxneighbors) &&
                ((nother == 1) || (nother && acyclic (dim, *b)));
 
        // add the neighbors from 'cset'
        int_t ncset = getneighbors (q, b, cset, neighbors_main, 0);
 
        // if there are no cubes in 'cset', then this cube is ok
        if (!ncset)
                return true;
 
        // if there are neighbors in 'cset' and the neighbors
        // from 'other' are not acyclic, the cube canot be removed
        if (!otheracyclic)
                return false;
 
        // if there are neighbors from 'cset' then check if the neighbors
        // from both 'cset' and 'other' taken together form an acyclic set
        if ((ncset + nother > 1) && ((ncset + nother == maxneighbors) ||
                !acyclic (dim, *b)))
        {
                return false;
        }
        return true;
} /* acyclic_rel */
 
 
// --------------------------------------------------
// -------------- acyclicity for maps ---------------
// --------------------------------------------------
 
template <class tCube, class tCell, class tSet>
int_t computeimage (hashedset<tCube> &img, const tCell &face,
        const mvmap<tCube,tCube> &map, const tSet &cset,
        const tCube &ignore)
{
        // get the coordinates of the cubical cell
        typename tCell::CoordType left [tCell::MaxDim];
        face. leftcoord (left);
        typename tCell::CoordType right [tCell::MaxDim];
        face. rightcoord (right);
 
        // find the images of all the cubes containing this cell
        int spacedim = face. spacedim ();
        typename tCell::CoordType leftb [tCell::MaxDim];
        typename tCell::CoordType rightb [tCell::MaxDim];
        for (int j = 0; j < spacedim; ++ j)
        {
                typename tCell::CoordType diff =
                        (left [j] == right [j]) ? 1 : 0;
                leftb [j] = (left [j] - diff);
                rightb [j] = (right [j] + diff);
        }
        tRectangle<typename tCell::CoordType> r (leftb, rightb, spacedim);
        const typename tCell::CoordType *c;
        int_t countimages = 0;
        while ((c = r. get ()) != NULL)
        {
                if (!tCube::PointBase::check (c, spacedim))
                        continue;
                tCube q (c, spacedim);
                if (q == ignore)
                        continue;
                if (!cset. check (q))
                        continue;
                img. add (map (q));
                ++ countimages;
        }
 
        return countimages;
} /* computeimage */
 
template <class tCube, class tSet>
bool remainsacyclic (const mvmap<tCube,tCube> &map, const tCube &q,
        const tSet &cset)
{
        // if the map is trivial then the answer is obvious
        if (map. getdomain (). empty ())
                return true;
 
        // compute the maximal number of neighbors of a cube
        int_t maxneighbors = getmaxneighbors (q. dim ());
 
        // prepare a bitfield and allocate it if necessary
        static BitField b;
        static int_t _maxneighbors = 0;
        static auto_array<unsigned char> buffer;
        if (maxneighbors != _maxneighbors)
        {
                _maxneighbors = maxneighbors;
                buffer. reset (new unsigned char [(maxneighbors + 7) >> 3]);
                b. define (buffer. get (), maxneighbors);
        }
 
        // clear the neighborbits
        b. clearall (maxneighbors);
 
        // get the bitfield representing the set of the neighbors of the cube
        getneighbors (q, &b, cset, 0);
 
        // create all the faces of the cube
        gcomplex<typename tCube::CellType,integer> faces;
        addneighbors (q, b, faces);
        faces. addboundaries ();
 
        // compute the new images of all the faces
        // and determine if they are acyclic
        int startdim = faces. dim ();
        for (int d = startdim; d >= 0; -- d)
        {
                for (int_t i = 0; i < faces [d]. size (); ++ i)
                {
                        // compute the image of the face in the set
                        hashedset<tCube> img;
                        int_t n = computeimage (img, faces [d] [i], map,
                                cset, q);
 
                        // if this is the image of only one cube, it is Ok
                        if (n == 1)
                                continue;
 
                        // verify whether the large image (with 'q')
                        // can be reduced towards the small one (without 'q')
                        hashedset<tCube> imgsurplus = map (q);
                        imgsurplus. remove (img);
                        cubreducequiet (img, imgsurplus);
                        if (!imgsurplus. empty ())
                                return false;
                }
        }
 
        return true;
} /* remainsacyclic */
 
template <class tCube>
class MapRemainsAcyclic
{
public:
        MapRemainsAcyclic (const mvmap<tCube,tCube> &the_map);
 
        template <class tSet>
        bool operator () (const tCube &q, const tSet &cset) const;
 
protected:
        const mvmap<tCube,tCube> &map;
 
}; /* class MapRemainsAcyclic */
 
template <class tCube>
inline MapRemainsAcyclic<tCube>::MapRemainsAcyclic (const mvmap<tCube,tCube> &the_map):
        map (the_map)
{
        // determine if the acyclicity of the map should be considered
//      bool checkacyclic = !map. getdomain (). empty ();
 
        return;
} /* MapRemainsAcyclic<tCube>::MapRemainsAcyclic */
 
template <class tCube>
template <class tSet>
inline bool MapRemainsAcyclic<tCube>::operator () (const tCube &q,
        const tSet &cset) const
{
        return remainsacyclic (map, q, cset);
} /* MapRemainsAcyclic<tCube>::operator () */
 
// --------------------------------------------------
 
template <class tCube>
class MapCanExpand
{
public:
        MapCanExpand (const mvmap<tCube,tCube> &_map,
                hashedset<tCube> &_imgsrc, hashedset<tCube> &_img,
                bool _indexmap, bool _checkacyclic, bool _quiet);
 
        template <class tSet1, class tSet2>
        bool operator () (const tCube &c, const tSet1 &cset,
                const tSet2 &other) const;
 
protected:
        const mvmap<tCube,tCube> &map;
 
        hashedset<tCube> &imgsrc;
 
        hashedset<tCube> &img;
 
        bool indexmap;
 
        bool checkacyclic;
 
        bool quiet;
 
}; /* class MapCanExpand */
 
template <class tCube>
inline MapCanExpand<tCube>::MapCanExpand (const mvmap<tCube,tCube> &_map,
        hashedset<tCube> &_imgsrc, hashedset<tCube> &_img,
        bool _indexmap, bool _checkacyclic, bool _quiet):
        map (_map), imgsrc (_imgsrc), img (_img),
        indexmap (_indexmap), checkacyclic (_checkacyclic),
        quiet (_quiet)
{
        return;
} /* MapCanExpand<tCube>::MapCanExpand */
 
template <class tCube>
template <class tSet1, class tSet2>
inline bool MapCanExpand<tCube>::operator () (const tCube &c,
        const tSet1 &cset, const tSet2 &other) const
{
        // take the image and reduce what is outside 'img'
        const hashedset<tCube> &cubimage = map (c);
        hashedset<tCube> ximage = cubimage;
        if (indexmap)
                ximage. add (c);
        ximage. remove (img);
        cubreducequiet (img, ximage);
 
        // if the image could not be reduced, skip the cube
        if (!ximage. empty ())
                return false;
 
        // make sure that the acyclicity of the map is not spoiled
        if (checkacyclic && !remainsacyclic (map, c, other))
                return false;
 
        // add the image of this cube to the image of the map
        if (!quiet && !cubimage. empty ())
        {
                sseq << "R\n";
                for (int_t j = 0; j < cubimage. size (); ++ j)
                        sseq << '2' << cubimage [j] << '\n';
                if (indexmap)
                        sseq << '2' << c << '\n';
        }
        img. add (cubimage);
        if (indexmap)
                img. add (c);
 
        // remove from 'imgsrc' all the cubes added to 'img'
        imgsrc. remove (cubimage);
        if (indexmap)
                imgsrc. remove (c);
 
        return true;
} /* MapCanExpand<tCube>::operator () */
 
// --------------------------------------------------
 
template <class tCube>
class MapCanExpandDummy
{
public:
        MapCanExpandDummy ();
 
        template <class tSet1, class tSet2>
        bool operator () (const tCube &q, const tSet1 &cset,
                const tSet2 &other) const;
 
}; /* class MapCanExpandDummy */
 
template <class tCube>
inline MapCanExpandDummy<tCube>::MapCanExpandDummy ()
{
        return;
} /* MapCanExpandDummy<tCube>::MapCanExpandDummy */
 
template <class tCube>
template <class tSet1, class tSet2>
inline bool MapCanExpandDummy<tCube>::operator () (const tCube &,
        const tSet1 &, const tSet2 &) const
{
        return true;
} /* MapCanExpandDummy<tCube>::operator () */
 
 
} // namespace homology
} // namespace chomp
 
#endif // _CHOMP_HOMOLOGY_CUBACYCL_H_