boundary.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// \file
///
/// Boundary computation at the level of chains of cells.
///
/////////////////////////////////////////////////////////////////////////////

// Copyright (C) 2009-2016 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 <http://www.gnu.org/licenses/>.

// Started on March 24, 2009. Last revision: February 5, 2016.


#ifndef _CHAINCON_BOUNDARY_H_
#define _CHAINCON_BOUNDARY_H_


// include some standard C++ header files
#include <vector>

// include relevant local header files
#include "chaincon/combchain.h"
#include "chaincon/chain.h"
#include "chaincon/combtensor.h"
#include "chaincon/tensor.h"
#include "chaincon/comblinmap.h"
#include "chaincon/linmap.h"
#include "chaincon/filtcomplex.h"

// include selected header files from the CHomP library
#include "chomp/system/config.h"
#include "chomp/struct/multitab.h"
//#include "chomp/system/textfile.h"


// --------------------------------------------------
// ----------------- helper classes -----------------
// --------------------------------------------------

/// A template of a simple class for unrestricted cell selection
/// in the boundary computation procedures.
class Unrestricted
{
public:
        /// Checks whether the given cell satisfies the restriction.
        /// In this class, always returns true.
        template <class CellT>
        bool check (const CellT &c) const;

}; /* class Unrestricted */

template <class CellT>
inline bool Unrestricted::check (const CellT &c) const
{
        return true;
} /* Unrestricted::check */

// --------------------------------------------------

/// A template of a simple class for negating the cell restriction
/// given by another object for the boundary computation procedures.
template <class CellRestrT>
class NegateRestr
{
public:
        /// Creates an object of the class based on another object
        /// whose decisions will be negated.
        NegateRestr (const CellRestrT &p);

        /// Checks whether the given cell satisfies the negation
        /// of the restriction given by the parent object.
        template <class CellT>
        bool check (const CellT &c) const;

private:
        /// The parent object whose decisions will be negated.
        const CellRestrT &parent;

}; /* class NegateRestr */

template <class CellRestrT>
inline NegateRestr<CellRestrT>::NegateRestr (const CellRestrT &p):
        parent (p)
{
        return;
} /* NegateRestr::NegateRestr */

template <class CellRestrT>
template <class CellT>
inline bool NegateRestr<CellRestrT>::check (const CellT &c) const
{
        return !parent. check (c);
} /* NegateRestr::check */


// --------------------------------------------------
// ------- boundary of a combinatorial chain --------
// --------------------------------------------------

/// Adds the boundary of a given cell to the provided chain
/// (takes boundary cells restricted by the given object,
/// preferably of the type hashedset<CellT>).
template <class CellT, class CellRestrT>
inline void computeBoundary (const CellT &c, tCombChain<CellT> &b,
        const CellRestrT &restr)
{
        // determine the length of the boundary of this cell
        int length = c. boundaryLength ();

        // process all the cells in its boundary
        for (int i = 0; i < length; ++ i)
        {
                // create the subsequent boundary cell
                CellT bc (c, i);

                // if the cell is not in the restriction then skip it
                if (!restr. check (bc))
                        continue;

                // add the boundary cell to the given chain
                b. add (bc);
        }
        return;
} /* computeBoundary */

/// Adds the boundary of a given chain to the provided chain
/// (takes boundary cells restricted by the given object,
/// preferably of the type hashedset<CellT>).
template <class CellT, class CellRestrT>
inline void computeBoundary (const tCombChain<CellT> &c,
        tCombChain<CellT> &b, const CellRestrT &restr)
{
        int_t size = c. size ();
        for (int_t i = 0; i < size; ++ i)
                computeBoundary (c. getCell (i), b, restr);
        return;
} /* computeBoundary */

/// Returns the boundary of a given chain (takes boundary cells restricted
/// by the given object, preferably of the type hashedset<CellT>).
template <class CellT, class CellRestrT>
inline tCombChain<CellT> boundary (const tCombChain<CellT> &c,
        const CellRestrT &restr)
{
        // prepare a chain that is going to contain the boundary of c
        tCombChain<CellT> b;

        // add the boundary of each cell in the chain to b
        computeBoundary (c, b, restr);

        // return the computed boundary chain
        return b;
} /* boundary */


// --------------------------------------------------
// -------------- boundary of a chain ---------------
// --------------------------------------------------

/// Adds the boundary of a given cell to the provided chain
/// with the given coefficient (takes boundary cells restricted
/// by the given object, preferably of the type hashedset<CellT>).
template <class CellT, class CoefT, class CellRestrT>
inline void computeBoundary (const CellT &c, tChain<CellT,CoefT> &b,
        const CoefT &coef, const CellRestrT &restr)
{
        // determine the length of the boundary of this cell
        int length = c. boundaryLength ();

        // process all the cells in its boundary
        for (int i = 0; i < length; ++ i)
        {
                // create the subsequent boundary cell
                CellT bc (c, i);

                // if the cell is not in the restriction then skip it
                if (!restr. check (bc))
                        continue;

                // determine the coefficient to use
                CoefT bcoef (c. boundaryCoef (i));
                bcoef *= coef;

                // add the boundary cell to the given chain
                b. add (bc, bcoef);
        }
        return;
} /* computeBoundary */

/// Adds the boundary of a given chain to the provided chain
/// (takes boundary cells restricted by the given object,
/// preferably of the type hashedset<CellT>).
template <class CellT, class CoefT, class CellRestrT>
inline void computeBoundary (const tChain<CellT,CoefT> &c,
        tChain<CellT,CoefT> &b, const CellRestrT &restr)
{
        int_t size = c. size ();
        for (int_t i = 0; i < size; ++ i)
                computeBoundary (c. getCell (i), b, c. getCoef (i), restr);
        return;
} /* computeBoundary */

/// Returns the boundary of a given chain (takes boundary cells restricted
/// by the given object, preferably of the type hashedset<CellT>).
template <class CellT, class CoefT, class CellRestrT>
inline tChain<CellT,CoefT> boundary (const tChain<CellT,CoefT> &c,
        const CellRestrT &restr)
{
        // prepare a chain that is going to contain the boundary of c
        tChain<CellT,CoefT> b;

        // add the boundary of each cell in the chain to b
        computeBoundary (c, b, restr);

        // return the computed boundary chain
        return b;
} /* boundary */

// Returns the boundary of a given chain (no restriction on boundary cells).
//template <class CellT, class CoefT>
//inline tChain<CellT,CoefT> boundary (const tChain<CellT,CoefT> &c)
//{
//      return boundary (c, Unrestricted<CellT> ());
//} /* boundary */


// --------------------------------------------------
// ------- boundary of a combinatorial tensor -------
// --------------------------------------------------

/// Adds the boundary of a given tensor to the provided tensor.
template <class Cell1T, class Cell2T, class CellRestrT>
inline void computeBoundary (const tCombTensor<Cell1T,Cell2T> &c,
        tCombTensor<Cell1T,Cell2T> &b, const CellRestrT &restr)
{
        int_t size = c. size ();
        for (int_t i = 0; i < size; ++ i)
        {
                tCombChain<Cell1T> left;
                computeBoundary (c. left (i), left, restr);
                tCombChain<Cell2T> right;
                computeBoundary (c. right (i), right, restr);
                b. add (left, right);
        }
        return;
} /* computeBoundary */

/// Returns the boundary of a given tensor.
template <class Cell1T, class Cell2T, class CellRestrT>
inline tCombTensor<Cell1T,Cell2T> boundary
        (const tCombTensor<Cell1T,Cell2T> &c, const CellRestrT &restr)
{
        // prepare a tensor that is going to contain the boundary of c
        tCombTensor<Cell1T,Cell2T> b;

        // add the boundary of each pair in the tensor c to the tensor b
        computeBoundary (c, b, restr);

        // return the computed boundary tensor
        return b;
} /* boundary */


// --------------------------------------------------
// -------------- boundary of a tensor --------------
// --------------------------------------------------

/// Adds the boundary of a given tensor to the provided tensor.
template <class Cell1T, class Cell2T, class CoefT, class CellRestrT>
inline void computeBoundary (const tTensor<Cell1T,Cell2T,CoefT> &c,
        tTensor<Cell1T,Cell2T,CoefT> &b, const CellRestrT &restr)
{
        int_t size = c. size ();
        for (int_t i = 0; i < size; ++ i)
        {
                tChain<Cell1T,CoefT> left;
                computeBoundary (c. left (i), left, restr);
                tChain<Cell2T,CoefT> right;
                computeBoundary (c. right (i), right, restr);
                b. add (left, right);
        }
        return;
} /* computeBoundary */

/// Returns the boundary of a given tensor.
template <class Cell1T, class Cell2T, class CoefT, class CellRestrT>
inline tTensor<Cell1T,Cell2T,CoefT> boundary
        (const tTensor<Cell1T,Cell2T,CoefT> &c, const CellRestrT &restr)
{
        // prepare a tensor that is going to contain the boundary of c
        tTensor<Cell1T,Cell2T,CoefT> b;

        // add the boundary of each pair in the tensor c to the tensor b
        computeBoundary (c, b, restr);

        // return the computed boundary tensor
        return b;
} /* boundary */


// --------------------------------------------------
// ----- boundary of a combinatorial linear map -----
// --------------------------------------------------

/// Computes the full boundary map for a cellular complex that contains
/// the cells in a provided set and all their boundary cells.
/// The map is assumed to be initially zero.
template <class CellArray, class CellT, class CellRestrT>
inline void computeBoundaryMap (const CellArray &cells,
        tCombLinMap<CellT,CellT> &f, const CellRestrT &restr)
{
        // prepare a set of cells whose boundaries yet have to be computed
        chomp::homology::hashedset<CellT> waiting;

        // process all the cells in the domain of the map
        // (note that the domain may increase during this process)
        int_t inputCellNumber = 0;
        const int_t inputCellCount = cells. size ();
        while (1)
        {
                // determine whether a cell from the array should be used
                bool useInput = (inputCellNumber < inputCellCount);

                // if the input array has been used and the set
                // of wating cells as well then exit the loop
                if (!useInput && waiting. empty ())
                        break;

                // take a cell for processing
                int_t waitingPosition = waiting. size () - 1;
                const CellT &c = useInput ? cells [inputCellNumber] :
                        waiting [waitingPosition];

                // determine the length of the boundary of this cell
                int length = c. boundaryLength ();

                // process all the cells in its boundary
                for (int b = 0; b < length; ++ b)
                {
                        // create the subsequent boundary cell
                        CellT bc (c, b);

                        // if the cell is not in the restriction then skip it
                        if (!restr. check (bc))
                                continue;

                        // add the boundary cell to the waiting list
                        if (bc. boundaryLength () &&
                                f. getImage (bc). empty ())
                        {
                                waiting. add (bc);
                        }

                        // add the boundary cell to the image of 'c'
                        f. add (c, bc);
                }

                // increase the counter of cells in the input array
                if (useInput)
                        ++ inputCellNumber;

                // or remove the waiting cell which was taken
                else
                        waiting. removenum (waitingPosition);
        }
        return;
} /* computeBoundaryMap */


// --------------------------------------------------
// ------------ boundary of a linear map ------------
// --------------------------------------------------

/// Computes the full boundary map for a cellular complex that contains
/// the cells in a provided set and all their boundary cells.
/// The map is assumed to be initially zero.
template <class CellArray, class CellT, class CoefT, class CellRestrT>
inline void computeBoundaryMap (const CellArray &cells,
        tLinMap<CellT,CellT,CoefT> &f, const CellRestrT &restr)
{
        // prepare a set of cells whose boundaries yet have to be computed
        chomp::homology::hashedset<CellT> waiting;

        // process all the cells in the domain of the map
        // (note that the domain may increase during this process)
        int_t inputCellNumber = 0;
        const int_t inputCellCount = cells. size ();
        while (1)
        {
                // determine whether a cell from the array should be used
                bool useInput = (inputCellNumber < inputCellCount);

                // if the input array has been used and the set
                // of wating cells as well then exit the loop
                if (!useInput && waiting. empty ())
                        break;

                // take a cell for processing
                int_t waitingPosition = waiting. size () - 1;
                const CellT &c = useInput ? cells [inputCellNumber] :
                        waiting [waitingPosition];

                // determine the length of the boundary of this cell
                int length = c. boundaryLength ();

                // process all the cells in its boundary
                for (int b = 0; b < length; ++ b)
                {
                        // create the subsequent boundary cell
                        CellT bc (c, b);

                        // if the cell is not in the restriction then skip it
                        if (!restr. check (bc))
                                continue;

                        // add the boundary cell to the waiting list
                        if (bc. boundaryLength () &&
                                f. getImage (bc). empty ())
                        {
                                waiting. add (bc);
                        }

                        // add the boundary cell to the image of 'c'
                        f. add (c, bc, CoefT (c. boundaryCoef (b)));
                }

                // increase the counter of cells in the input array
                if (useInput)
                        ++ inputCellNumber;

                // or remove the waiting cell which was taken
                else
                        waiting. removenum (waitingPosition);
        }
        return;
} /* computeBoundaryMap */


// --------------------------------------------------
// -------- boundaries in a filtered complex --------
// --------------------------------------------------

/// Adds boundaries to all the cells in the given filtered complex.
template <class CellT, class CellRestrT>
inline void addBoundaries (tFilteredComplex<CellT> &K,
        const CellRestrT &restr)
{
//      using chomp::homology::sbug;
        using chomp::homology::scon;

        // if the complex is empty then there is no work to do
        if (K. empty ())
                return;

        // determine the dimensions of all the cells in the input complex
        chomp::homology::multitable<std::vector<int_t> > cellsByDim (32);
        int topDim (-1);
        int_t Ksize = K. size ();
        for (int_t i = 0; i < Ksize; ++ i)
        {
                int dim (K [i]. dim ());
                if (dim >= 0)
                {
                        cellsByDim [dim]. push_back (i);
                        if (topDim < dim)
                                topDim = dim;
                }
        }

        // prepare another complex for the result and swap it with the input
        tFilteredComplex<CellT> L;
        L. swap (K);
        int_t Lsize (L. size ());
        int_t Lused (0);

        // prepare the current dimension of cells (starting with top dim)
        int curDim (topDim);

        // process all the cells in the filter (its size may be changing)
        int_t counter (0);
        for (int_t n = 0; (n < K. size ()) || (Lused < Lsize); ++ n)
        {
        //      sbug << "DEBUG: n = " << n << ", K. size () = " <<
        //              K. size () << ".\n";
                // add appropriate cells from L if the dimension
                // of the next cell to be processed allows for that
                if ((curDim >= 0) && ((n >= K. size ()) ||
                        (K [K. size () - n - 1]. dim () - 1 <= curDim)))
                {
                        size_t curDimSize (cellsByDim [curDim]. size ());
                //      sbug << "DEBUG: Adding " << curDimSize <<
                //              " cell(s) from L of dim " << curDim << ".\n";
                        for (size_t i = 0; i < curDimSize; ++ i)
                        {
                                K. add (L [cellsByDim [curDim] [i]]);
                                ++ counter;
                                if (!(counter % 5347))
                                        scon << std::setw (15) << counter <<
                                                "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b";
                        }
                        Lused += curDimSize;
                        -- curDim;
                        -- n;
                        continue;
                }

                // take the cell
                const CellT &c = K [K. size () - n - 1];
        //      sbug << "DEBUG: Computing bounary of " << c << ".\n";

                // process the boundary cells of the cell
                int len = c. boundaryLength ();
                for (int i = 0; i < len; ++ i)
                {
                        // create the subsequent boundary cell
                        CellT bc (c, i);

                        // if the cell is not in the restriction then skip it
                        if (!restr. check (bc))
                                continue;

                        // add the boundary cell to the filtered complex
                        int_t cn = K. add (bc);
                        ++ counter;
                        if (!(counter % 5347))
                                scon << std::setw (15) << counter <<
                                        "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b";

                        // make sure that the parent cell appears earlier
                        if (cn < n)
                        {
                                throw "A boundary cell appears too early "
                                        "in the filter.";
                        }
                }
        }

        return;
} /* addBoundaries */


#endif // _CHAINCON_BOUNDARY_H_