cubired1.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_CUBIRED1_H_
#define _CHOMP_HOMOLOGY_CUBIRED1_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 "chomp/homology/cubacycl.h"
 
#include <iostream>
#include <fstream>
#include <cstdlib>
 
namespace chomp {
namespace homology {
namespace reduction1 {
 
 
// --------------------------------------------------
// ----------------- add neighbors ------------------
// --------------------------------------------------
 
template <class tCube, class tCubeSet>
inline void addcubeneighbors (const tCube &q, int dim,
        const tCubeSet &cubset, bitfield *b, hashedset<tCube> &neighbors,
        hashedset<tCube> &queue, const hashedset<tCube> &notthese)
{
        // determine the neighbors of this cube (if not done yet)
        if (dim <= MaxBddDim)
                getneighbors (q, b, cubset, &neighbors, 0);
 
        // add the computed neighbors of this cube to the queue
        for (int_t j = 0; j < neighbors. size (); ++ j)
        {
                if (!notthese. check (neighbors [j]))
                        queue. add (neighbors [j]);
        }
 
        return;
} /* addcubeneighbors */
 
 
// --------------------------------------------------
// ------------------ reduce cubes ------------------
// --------------------------------------------------
 
template <class tCube, class tVerify>
int_t cubreducequiet (const hashedset<tCube> &maincset,
        hashedset<tCube> &cset, hashedset<tCube> &other,
        const hashedset<tCube> &keep, const tVerify &verify,
        bool quiet = true)
{
        // remember the initial number of cubes in both sets
        int_t prev = cset. size () + other. size ();
 
        // if the case is trivial, return the answer
        if (maincset. empty () && cset. empty ())
        {
                if (!other. empty ())
                {
                        hashedset<tCube> empty;
                        other = empty;
                }
                return prev;
        }
 
        // determine the space dimension
        int dim = (cset. empty ()) ? maincset [0]. dim () :
                cset [0]. dim ();
 
        // compute the maximal number of neighbors of a cube
        int_t maxneighbors = getmaxneighbors (dim);
 
        // prepare a bitfield for storing neighborhood configurations
        // and allocate it if necessary (this is a static variable)
        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);
        }
 
        // prepare a queue for cubes to check
        hashedset<tCube> queue;
 
        // prepare a counter for displaying the progress of computations
        int_t count = 0;
 
        // go through the other set, remove cubes which can be removed,
        // and add their neighbors (only in the other set) to the queue
        for (int_t i = 0; i < other. size (); ++ i)
        {
                // show the progress indicator
                ++ count;
                if (!quiet && !(count % 373))
                        scon << std::setw (10) << count <<
                                "\b\b\b\b\b\b\b\b\b\b";
 
                // if the cube should be kept, skip it
                if (keep. check (other [i]))
                        continue;
 
                // clear the neighborbits
                b. clearall (maxneighbors);
 
                // prepare a set for storing the neighbors of the cube
                hashedset<tCube> neighbors;
                hashedset<tCube> *none (0);
 
                // if the cube cannot be removed, skip it
                if (!acyclic_rel (other [i], dim,
                        makesetunion (maincset, cset), other, &b,
                        maxneighbors, none, &neighbors))
                {
                        continue;
                }
 
                // verify if the removal of the cube from A is OK
                if (!verify (other [i], other))
                        continue;
 
                // verify if the removal of the cube from X is OK
                if (!verify (other [i], makesetunion (maincset,
                        makesetunion (other, cset))))
                {
                        continue;
                }
 
                // remove this cube from the queue
                queue. remove (other [i]);
 
                // determine the neighbors of this cube (if not done yet)
                // and add the computed neighbors of this cube to the queue
                addcubeneighbors (other [i], dim, other, &b, neighbors,
                        queue, keep);
 
                // remove this cube from 'other'
                if (!quiet)
                        sseq << '0' << other [i] << '\n';
                other. removenum (i);
                -- i;
        }
 
        // show a temporary progress indicator and reset the counter
        if (!quiet)
                scon << '.';
        count = 0;
 
        // go through the main set (scan those cubes only that should
        // be considered for removal), remove cubes which can be removed,
        // and add their neighbors to the queue
        for (int_t i = 0; i < cset. size (); ++ i)
        {
                // show progress indicator
                if (!quiet && !(count % 373))
                        scon << std::setw (10) << count <<
                                "\b\b\b\b\b\b\b\b\b\b";
                ++ count;
 
                // if this cube should be kept, skip it
                if (keep. check (cset [i]))
                        continue;
 
                // clear the neighborbits
                b. clearall (maxneighbors);
 
                // prepare a set for storing the neighbors of the cube
                hashedset<tCube> neighbors;
 
                // if the neighborhood of the cube is not acyclic, skip it
                if (!acyclic (cset [i], dim, makesetunion (maincset,
                        makesetunion (cset, other)),
                        &b, maxneighbors, &neighbors))
                {
                        continue;
                }
 
                // verify if the removal of the cube from X is OK
                if (!verify (cset [i], makesetunion (maincset,
                        makesetunion (cset, other))))
                {
                        continue;
                }
 
                // remove this cube from the queue
                queue. remove (cset [i]);
 
                // determine the neighbors of this cube (if not done yet)
                // and add the computed neighbors of this cube to the queue;
                // neighbors in 'maincset' should be skipped, because
                // they are not considered for reduction anyway
                addcubeneighbors (cset [i], dim,
                        makesetunion (cset, other),
                        &b, neighbors, queue, keep);
 
                // remove this cube from 'cset'
                if (!quiet)
                        sseq << '0' << cset [i] << '\n';
                cset. removenum (i);
                -- i;
        }
 
        // update the temporary progress indicator and reset the counter
        if (!quiet)
                scon << "\b*";
        count = 0;
 
        // take cubes from the queue and remove them if possible
        while (!queue. empty ())
        {
                // update the progress indicator
                if (!quiet && !(count % 373))
                        scon << std::setw (10) << count <<
                                "\b\b\b\b\b\b\b\b\b\b";
                ++ count;
 
                // take a cube from the queue and remove it from the queue
                tCube c = queue [0];
                queue. removenum (0);
 
                // clean the neighborbits
                b. clearall (maxneighbors);
 
                // debug check: the cube must not be in 'keep' nor 'maincset'
                if (keep. check (c))
                        throw "A cube from 'keep' found in a reduction queue!";
                if (maincset. check (c))
                        throw "A cube from 'maincset' found in a reduction queue!";
 
                // if this cube belongs to the other set...
                if (other. check (c))
                {
                        // prepare a set for storing the neighbors
                        hashedset<tCube> neighbors;
 
                        if (!acyclic_rel (c, dim,
                                makesetunion (maincset, cset), other, &b,
                                maxneighbors, &neighbors, &neighbors))
                        {
                                continue;
                        }
 
                        // verify if the removal of the cube from A is OK
                        if (!verify (c, other))
                                continue;
 
                        // verify if the removal of the cube from X is OK
                        if (!verify (c, makesetunion (maincset,
                                makesetunion (other, cset))))
                        {
                                continue;
                        }
 
                        // determine the neighbors of this cube (if not yet)
                        // and add the computed neighbors to the queue
                        addcubeneighbors (c, dim,
                                makesetunion (cset, other),
                                &b, neighbors, queue, keep);
 
                        // remove the cube from the other set
                        if (!quiet)
                                sseq << '0' << c << '\n';
                        other. remove (c);
                }
 
                // otherwise, if this cube belongs to 'cset'...
                else
                {
                        // prepare a set for storing the neighbors
                        hashedset<tCube> neighbors;
 
                        // if the neighborhood of the cube is not acyclic
                        // then skip it
                        if (!acyclic (c, dim, makesetunion (maincset,
                                makesetunion (cset, other)),
                                &b, maxneighbors, &neighbors))
                        {
                                continue;
                        }
 
                        // verify if the removal of the cube from X is OK
                        if (!verify (c, makesetunion (maincset,
                                makesetunion (cset, other))))
                        {
                                continue;
                        }
 
                        // determine the neighbors of this cube (if not yet)
                        // and add the computed neighbors to the queue
                        addcubeneighbors (c, dim,
                                makesetunion (cset, other),
                                &b, neighbors, queue, keep);
 
                        // remove the cube from 'cset'
                        if (!quiet)
                                sseq << '0' << c << '\n';
                        cset. remove (c);
                }
        }
 
        // erase the temporary progress indicator
        if (!quiet)
                scon << "\b \b";
 
        // return the number of cubes removed from both sets
        return prev - cset. size () - other. size ();
} /* cubreducequiet */
 
 
// --------------------------------------------------
// ------------------ expand cubes ------------------
// --------------------------------------------------
 
template <class tCube, class tVerify>
int_t cubexpand (hashedset<tCube> &cset, hashedset<tCube> &other,
        const tVerify &verify, bool quiet = false)
{
        // if the case is trivial, return the answer
        if (cset. empty () || other. empty ())
                return 0;
 
        // remember the initial number of cubes in the main set
        int_t prev = cset. size ();
 
        // determine the space dimension
        int dim = cset [0]. dim ();
 
        // compute the maximal number of neighbors of a cube
        int_t maxneighbors = getmaxneighbors (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);
        }
 
        // prepare a queue for cubes to check
        hashedset<tCube> queue;
 
        // prepare a counter for displaying the progress of computations
        int_t count = 0;
 
        // go through the main set, move cubes to the other set
        // if possible, and add their neighbors to the queue
        for (int_t i = 0; i < cset. size (); ++ i)
        {
                // show progress indicator
                ++ count;
                if (!quiet && !(count % 373))
                        scon << std::setw (10) << count <<
                                "\b\b\b\b\b\b\b\b\b\b";
 
                // take a cube from 'cset'
                const tCube &c = cset [i];
 
                // clear the neighborbits
                b. clearall (maxneighbors);
 
                // if the neighborhood of the cube is not acyclic, skip it
                if (!acyclic (c, dim, other, &b, maxneighbors))
                        continue;
 
                // verify if moving this cube to 'other' is allowed
                if (!verify (c, cset, other))
                        continue;
 
                // remove this cube from the queue
                queue. remove (c);
 
                // add neighbors from 'cset' to the queue
                b. clearall (maxneighbors);
                getneighbors (c, &b, cset, &queue, 0);
 
                // move the cube to the other set
                if (!quiet)
                        sseq << "L\n" << '2' << c << '\n';
                other. add (c);
                cset. removenum (i);
                -- i;
        }
 
        // show a temporary progress indicator and reset the counter
        if (!quiet)
                scon << '*';
        count = 0;
 
        // take cubes from the queue and move them from 'cset' to 'other'
        while (!queue. empty ())
        {
                // show progress indicator
                if (!quiet && !(count % 373))
                        scon << std::setw (10) << count <<
                                "\b\b\b\b\b\b\b\b\b\b";
                ++ count;
 
                // take a cube from the queue
                const tCube c = queue [0];
                queue. removenum (0);
 
                // clear the neighborbits
                b. clearall (maxneighbors);
 
                // if the neighborhood of the cube is not acyclic, skip it
                if (!acyclic (c, dim, other, &b, maxneighbors))
                        continue;
 
                // verify if moving this cube to 'other' is allowed
                if (!verify (c, cset, other))
                        continue;
 
                // add the neighbors from 'cset' to the queue
                b. clearall (maxneighbors);
                getneighbors (c, &b, cset, &queue, 0);
 
                // move this cube from 'cset' to 'other'
                if (!quiet)
                        sseq << "L\n" << '2' << c << '\n';
                cset. remove (c);
                other. add (c);
        }
 
        // erase the temporary progress indicator
        if (!quiet)
                scon << "\b \b";
 
        // return the number of cubes removed from 'cset'
        return prev - cset. size ();
} /* cubexpand */
 
 
} // namespace reduction1
} // namespace homology
} // namespace chomp
 
#endif // _CHOMP_HOMOLOGY_CUBIRED1_H_