cubcell.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// \file
///
/// A cubical cell.
///
/////////////////////////////////////////////////////////////////////////////

// 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 6, 2016.


#ifndef _CHAINCON_CUBCELL_H_
#define _CHAINCON_CUBCELL_H_


// include relevant local header files
//#include "chaincon/emptycell.h"
//#include "chaincon/wrapping.h"

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

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


// --------------------------------------------------
// ------------------ cubical cell ------------------
// --------------------------------------------------

/// An elementary cubical cell with vertex coordinates of integer type.
/// Note that unless using the most standard built-in integers,
/// the vertex coordinate integer type must have hash key functions
/// (hashkey1 and hashkey2, see e.g. those in chomp/struct/hashsets.h).
/// The wrapping class must provide a static function "wrap" for wrapping
/// the coordinates, separately at each dimension.
/// The empty cell existence decision class must provide a static function
/// "exists" that returns true if the empty cell must be taken
/// into consideration or false otherwise.
template <class CoordT, class WrapT, class EmptyT>
class tCubCell
{
public:
        /// The type of coordinates.
        typedef CoordT CoordType;

        /// The type of the coord wrapping class.
        typedef WrapT WrapType;

        /// The type of the empty cell existence decision class.
        typedef EmptyT EmptyType;

        /// The default constructor of an empty cube.
        tCubCell ();

        /// The constructor of a full cube from an array of coordinates.
        template <class CoordArray>
        tCubCell (int dimension, const CoordArray &coordinates);

        /// The constructor of a full cube from two arrays of coordinates.
        template <class CoordArray>
        tCubCell (int dimension, const CoordArray &left,
                const CoordArray &right);

        /// The copy constructor.
        tCubCell (const tCubCell<CoordT,WrapT,EmptyT> &c);

        /// The constructor of the n-th boundary cube.
        tCubCell (const tCubCell<CoordT,WrapT,EmptyT> &c, int n);

        /// The constructor of the n-th component of the Alexander-Whitney
        /// diagonal, either the left one (if side = 0), or the right one
        /// (if side = 1).
        tCubCell (const tCubCell<CoordT,WrapT,EmptyT> &c, int n, int side);

        /// The assignment operator.
        tCubCell<CoordT,WrapT,EmptyT> &operator =
                (const tCubCell<CoordT,WrapT,EmptyT> &s);

        /// The destructor.
        ~tCubCell ();

        /// Returns the dimension of the embedding space.
        int spaceDim () const;

        /// Returns the dimension of the elementary cube.
        int dim () const;

        /// Returns the n-th left coordinate of the elementary cube.
        const CoordT left (int n) const;

        /// Returns the n-th right coordinate of the elementary cube.
        const CoordT right (int n) const;

        /// Returns the length of the boundary of the elementary cube.
        int boundaryLength () const;

        /// Returns the n-th coefficient in the boundary of the cube.
        int boundaryCoef (int n) const;

        /// Returns the number of terms in the Alexander-Whitneney diagonal.
        int diagonalLength () const;

        /// Returns hashing key no. 1, based on the internal data.
        int_t hashkey1 () const;

        /// Return shashing key no. 2, based on the internal data
        int_t hashkey2 () const;

        /// Returns true if and only if the cubes are the same.
        bool operator == (const tCubCell<CoordT,WrapT,EmptyT> &s) const;

        /// Swaps the data between two cubical cells.
        void swap (tCubCell<CoordT,WrapT,EmptyT> &s);

protected:
        /// The number that defines the left corner of the elementary cube.
        /// The dimension of the embedding space is stored in high bits
        /// of this number.
        /// As an exception, the empty cell has this number set to 0.
        int_t n1;

        /// The number that defines the right corner of the elementary cube.
        /// The dimension of the cubical cell is stored in the high bits
        /// of this number.
        /// As an exception, the empty cell has this number set to 0.
        int_t n2;

}; /* class tCubCell */

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

template <class CoordT, class WrapT, class EmptyT>
inline tCubCell<CoordT,WrapT,EmptyT>::tCubCell (): n1 (0), n2 (0)
{
        return;
} /* tCubCell::tCubCell */

template <class CoordT, class WrapT, class EmptyT>
template <class CoordArray>
inline tCubCell<CoordT,WrapT,EmptyT>::tCubCell (int dimension,
        const CoordArray &coordinates): n1 (0), n2 (0)
{
        // if the embedding space dimension is not positive
        // then create the empty cubical cell
        if (dimension <= 0)
                return;

        // if the requested dimension is too high then throw an error message
        if (dimension >= chomp::homology::MaxBasDim)
                throw "Too high cubical cell dimension requested.";

        // compute the number that represents the left corner of the cell
        CoordT c [chomp::homology::MaxBasDim];
        for (int i = 0; i < dimension; ++ i)
                c [i] = coordinates [i];
        WrapT::wrap (c, dimension);
        n1 = chomp::homology::tPointBase<CoordT>::number (c, dimension);
        n1 |= dimension << chomp::homology::NumBits;

        // compute the number that represents the right corner of the cell
        for (int i = 0; i < dimension; ++ i)
                ++ (c [i]);
        WrapT::wrap (c, dimension);
        n2 = chomp::homology::tPointBase<CoordT>::number (c, dimension);
        n2 |= dimension << chomp::homology::NumBits;

        return;
} /* tCubCell::tCubCell */

template <class CoordT, class WrapT, class EmptyT>
template <class CoordArray>
inline tCubCell<CoordT,WrapT,EmptyT>::tCubCell (int dimension,
        const CoordArray &left, const CoordArray &right): n1 (0), n2 (0)
{
        // if the embedding space dimension is not positive
        // then create the empty cubical cell
        if (dimension <= 0)
                return;

        // if the requested dimension is too high then throw an error message
        if (dimension >= chomp::homology::MaxBasDim)
                throw "Too high cubical cell dimension requested.";

        // compute the dimension of the cell
        int cellDimension = 0;
        for (int i = 0; i < dimension; ++ i)
        {
                if (left [i] != right [i])
                        ++ cellDimension;
        }

        // compute the number that represents the left corner of the cell
        CoordT c [chomp::homology::MaxBasDim];
        for (int i = 0; i < dimension; ++ i)
                c [i] = left [i];
        WrapT::wrap (c, dimension);
        n1 = chomp::homology::tPointBase<CoordT>::number (c, dimension);
        n1 |= dimension << chomp::homology::NumBits;

        // compute the number that represents the right corner of the cell
        for (int i = 0; i < dimension; ++ i)
                c [i] = right [i];
        WrapT::wrap (c, dimension);
        n2 = chomp::homology::tPointBase<CoordT>::number (c, dimension);
        n2 |= cellDimension << chomp::homology::NumBits;

        return;
} /* tCubCell::tCubCell */

template <class CoordT, class WrapT, class EmptyT>
inline tCubCell<CoordT,WrapT,EmptyT>::tCubCell
        (const tCubCell<CoordT,WrapT,EmptyT> &c):
        n1 (c. n1), n2 (c. n2)
{
        return;
} /* tCubCell::tCubCell */

template <class CoordT, class WrapT, class EmptyT>
inline tCubCell<CoordT,WrapT,EmptyT>::tCubCell
        (const tCubCell<CoordT,WrapT,EmptyT> &c, int n)
{
        using chomp::homology::sbug;

        // determine the dimension of the cell
        int cellDimension = c. dim ();

        // if the dimension of the cell is zero then the result should be empty
        if (!cellDimension)
        {
                n1 = 0;
                n2 = 0;
                return;
        }

        // determine the dimension of the embedding space
        int spaceDimension = c. spaceDim ();

        // prepare arrays for the coordinates of the corners of the cell
        CoordT newLeft [chomp::homology::MaxBasDim];
        CoordT newRight [chomp::homology::MaxBasDim];
        for (int i = 0; i < spaceDimension; ++ i)
        {
                newLeft [i] = c. left (i);
                newRight [i] = c. right (i);
        }

        // if this is the first set of cells then squeeze rightwards
        bool squeezeRight = (n < cellDimension);
        int nn ((n >= cellDimension) ? (n - cellDimension) : n);
        int counter = 0;
        bool squeezed = false;
        for (int i = 0; i < spaceDimension; ++ i)
        {
                if (newLeft [i] != newRight [i])
                {
                        if (nn == counter)
                        {
                                if (squeezeRight)
                                        newLeft [i] = newRight [i];
                                else
                                        newRight [i] = newLeft [i];
                                squeezed = true;
                                break;
                        }
                        ++ counter;
                }
        }
        if (!squeezed)
        {
                sbug << "DEBUG: c = " << c << ".\n";
                sbug << "DEBUG: cellDimension = " << cellDimension << ".\n";
                sbug << "DEBUG: spaceDimension = " << spaceDimension << ".\n";
                sbug << "DEBUG: squeezeRight = " <<
                    (squeezeRight ? "true" : "false") << ".\n";
                sbug << "DEBUG: n = " << n << ".\n";
                sbug << "DEBUG: nn = " << nn << ".\n";
                throw "Wrong boundary element requested.";
        }
        n1 = chomp::homology::tPointBase<CoordT>::number
                (newLeft, spaceDimension);
        n1 |= spaceDimension << chomp::homology::NumBits;
        n2 = chomp::homology::tPointBase<CoordT>::number
                (newRight, spaceDimension);
        n2 |= (cellDimension - 1) << chomp::homology::NumBits;
        return;
} /* tCubCell::tCubCell */

/// The version of the AW diagonal. For testing only.
static int diagVersion = 0;

template <class CoordT, class WrapT, class EmptyT>
inline tCubCell<CoordT,WrapT,EmptyT>::tCubCell
        (const tCubCell<CoordT,WrapT,EmptyT> &c, int n, int side)
{
        int cellDimension = c. dim ();
        if (cellDimension != 2)
                throw "AW diagonal supported for cubes only in dim 2.";
        CoordT newLeft [chomp::homology::MaxBasDim];
        CoordT newRight [chomp::homology::MaxBasDim];
        int spaceDimension = c. spaceDim ();
        int indices [2];
        int index = 0;
        for (int i = 0; i < spaceDimension; ++ i)
        {
                newLeft [i] = c. left (i);
                newRight [i] = c. right (i);
                if (newLeft [i] != newRight [i])
                        indices [index ++] = i;
        }

        if (diagVersion == 1)
        {

        if (side == 0)
        {
                switch (n)
                {
                case 0:
                        newRight [indices [0]] = newLeft [indices [0]];
                        newLeft [indices [1]] = newRight [indices [1]];
                        break;
                case 1:
                        break;
                case 2:
                case 3:
                        newLeft [indices [0]] = newRight [indices [0]];
                        break;
                case 4:
                        newLeft [indices [1]] = newRight [indices [1]];
                        break;
                case 5:
                        newRight [indices [0]] = newLeft [indices [0]];
                        break;
                }
        }
        else
        {
                switch (n)
                {
                case 0:
                        break;
                case 1:
                        newRight [indices [0]] = newLeft [indices [0]];
                        newRight [indices [1]] = newLeft [indices [1]];
                        break;
                case 2:
                        newLeft [indices [1]] = newRight [indices [1]];
                        break;
                default:
                        newRight [indices [0]] = newLeft [indices [0]];
                        break;
                }
        }

        }

        else if (diagVersion == 2)
        {

        if (side == 0)
        {
                switch (n)
                {
                case 0:
                        newRight [indices [1]] = newLeft [indices [1]];
                        break;
                case 1:
                        newLeft [indices [0]] = newRight [indices [0]];
                        break;
                case 2:
                        newRight [indices [0]] = newLeft [indices [0]];
                        break;
                }
        }
        else
        {
                switch (n)
                {
                case 0:
                case 1:
                        newLeft [indices [0]] = newRight [indices [0]];
                        break;
                case 2:
                        newLeft [indices [1]] = newRight [indices [1]];
                        break;
                }
        }

        }

        else if (diagVersion == 3)
        {

        if (side == 0)
        {
                switch (n)
                {
                case 0:
                case 1:
                        newRight [indices [1]] = newLeft [indices [1]];
                        break;
                case 2:
                        newLeft [indices [0]] = newRight [indices [0]];
                        break;
                }
        }
        else
        {
                switch (n)
                {
                case 0:
                        newLeft [indices [0]] = newRight [indices [0]];
                        break;
                case 1:
                case 2:
                        newLeft [indices [1]] = newRight [indices [1]];
                        break;
                }
        }

        }

        else
                throw "Undefined version of A-W diagonal requested.\n";

        int newDimension = 0;
        for (int i = 0; i < spaceDimension; ++ i)
        {
                if (newLeft [i] != newRight [i])
                        ++ newDimension;
        }
        n1 = chomp::homology::tPointBase<CoordT>::number
                (newLeft, spaceDimension);
        n1 |= spaceDimension << chomp::homology::NumBits;
        n2 = chomp::homology::tPointBase<CoordT>::number
                (newRight, spaceDimension);
        n2 |= newDimension << chomp::homology::NumBits;
        return;
} /* tCubCell::tCubCell */

template <class CoordT, class WrapT, class EmptyT>
inline tCubCell<CoordT,WrapT,EmptyT> &
        tCubCell<CoordT,WrapT,EmptyT>::operator =
        (const tCubCell<CoordT,WrapT,EmptyT> &c)
{
        n1 = c. n1;
        n2 = c. n2;
        return *this;
} /* tCubCell::operator = */

template <class CoordT, class WrapT, class EmptyT>
inline tCubCell<CoordT,WrapT,EmptyT>::~tCubCell ()
{
        return;
} /* tCubCell::~tCubCell */

template <class CoordT, class WrapT, class EmptyT>
inline int tCubCell<CoordT,WrapT,EmptyT>::spaceDim () const
{
        return (n1 >> chomp::homology::NumBits);
} /* tCubCell::spaceDim */

template <class CoordT, class WrapT, class EmptyT>
inline int tCubCell<CoordT,WrapT,EmptyT>::dim () const
{
        if ((n1 == 0) && (n2 == 0))
                return -1;
        return (n2 >> chomp::homology::NumBits);
} /* tCubCell::dim */

template <class CoordT, class WrapT, class EmptyT>
inline const CoordT tCubCell<CoordT,WrapT,EmptyT>::left (int n) const
{
        return chomp::homology::tPointBase<CoordT>::coord
                (n1 & chomp::homology::NumMask, spaceDim ()) [n];
} /* tCubCell::left */

template <class CoordT, class WrapT, class EmptyT>
inline const CoordT tCubCell<CoordT,WrapT,EmptyT>::right (int n) const
{
        return chomp::homology::tPointBase<CoordT>::coord
                (n2 & chomp::homology::NumMask, spaceDim ()) [n];
} /* tCubCell::right */

template <class CoordT, class WrapT, class EmptyT>
inline int tCubCell<CoordT,WrapT,EmptyT>::boundaryLength () const
{
        int dimension = dim ();
        if (EmptyT::exists ())
                return (dimension > 0) ? (dimension << 1) : (dimension + 1);
        else
                return dimension << 1;
} /* tCubCell::boundaryLength */

template <class CoordT, class WrapT, class EmptyT>
inline int tCubCell<CoordT,WrapT,EmptyT>::boundaryCoef (int n) const
{
        if (!n)
                return 1;
        int dimension = dim ();
        if (n >= dimension)
                n -= dimension - 1;
        return (n & 1) ? -1 : 1;
} /* tCubCell::boundaryCoef */

template <class CoordT, class WrapT, class EmptyT>
inline int tCubCell<CoordT,WrapT,EmptyT>::diagonalLength () const
{
        if (dim () != 2)
                throw "A-W diagonal supported only in dim 2.";
        return 3; //6;
} /* tCubCell::diagonalLength */

template <class CoordT, class WrapT, class EmptyT>
inline int_t tCubCell<CoordT,WrapT,EmptyT>::hashkey1 () const
{
        return ((n1 ^ 0x55555555u) << 21) ^ ((n1 ^ 0xAAAAAAAAu) << 9) ^
                ((n2 ^ 0xAAAAAAAAu) >> 1);
} /* tCubCell::hashkey1 */

template <class CoordT, class WrapT, class EmptyT>
inline int_t tCubCell<CoordT,WrapT,EmptyT>::hashkey2 () const
{
        return ((n2 ^ 0xAAAAAAAAu) << 22) ^ ((n2 ^ 0x55555555u) << 10) ^
                (n1 ^ 0x55555555u);
} /* tCubCell::hashkey2 */

template <class CoordT, class WrapT, class EmptyT>
inline bool tCubCell<CoordT,WrapT,EmptyT>::operator ==
        (const tCubCell<CoordT,WrapT,EmptyT> &s) const
{
        return (n1 == s. n1) && (n2 == s. n2);
} /* tCubCell::operator == */

template <class CoordT, class WrapT, class EmptyT>
inline void tCubCell<CoordT,WrapT,EmptyT>::swap
        (tCubCell<CoordT,WrapT,EmptyT> &s)
{
        std::swap (n1, s. n1);
        std::swap (n2, s. n2);
        return;
} /* tCubCell::swap */

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

/// Generates a hashing key no. 1 for a cubical cell,
/// composed of hashing keys for the coordinates.
/// This key is to be used in a hashed set.
template <class CoordT, class WrapT, class EmptyT>
inline int_t hashkey1 (const tCubCell<CoordT,WrapT,EmptyT> &c)
{
        return c. hashkey1 ();
/*
//      using chomp::homology::hashkey1;
        using ::hashkey1;
        int d = c. spaceDim ();
        if (d <= 0)
                return 0;
        else if (d == 1)
                return hashkey1 (c. left (0)) << 2;
        else if (d == 2)
        {
                return (((hashkey1 (c. left (0)) & 0x65555555u) << 10) ^
                        ((hashkey1 (c. right (1)) & 0xAAAAAAAAu)));
        }
        else if (d == 3)
        {
                return (((hashkey1 (c. left (0)) & 0x65555555u) << 10) ^
                        ((hashkey1 (c. right (1)) & 0xAAAAAAAAu) << 2) ^
                        ((hashkey1 (c. left (2)) & 0xA5294252u) >> 1));
        }
        else
        {
                return (((hashkey1 (c. left ((d >> 2) + 1)) & 0x65555555u) << 16) ^
                        ((hashkey1 (c. right (d >> 3)) & 0x52942529u) << 10) ^
                        ((hashkey1 (c. right (d >> 1)) & 0xAAAAAAAAu) << 4) ^
                        ((hashkey1 (c. left (d - 1)) & 0xA5294252u) >> 1));
        }
*/
} /* hashkey1 */

/// Generates a hashing key no. 2 for a cubical cell,
/// composed of hashing keys for the coordinates.
/// This key is to be used in a hashed set.
template <class CoordT, class WrapT, class EmptyT>
inline int_t hashkey2 (const tCubCell<CoordT,WrapT,EmptyT> &c)
{
        return c. hashkey2 ();
/*
//      using chomp::homology::hashkey2;
        using ::hashkey2;
        int d = c. spaceDim ();
        if (d <= 0)
                return 0;
        else if (d == 1)
                return hashkey2 (c. right (0)) << 4;
        else if (d == 2)
        {
                return (((hashkey2 (c. right (0)) & 0xAAAAAAAAu) >> 1) ^
                        ((hashkey2 (c. left (1)) & 0x55555555u) << 7));
        }
        else if (d == 3)
        {
                return (((hashkey2 (c. right (2)) & 0x55555555u) << 9) ^
                        ((hashkey2 (c. right (0)) & 0xAAAAAAAAu) << 3) ^
                        ((hashkey2 (c. left (1)) & 0xA5294A52u) >> 1));
        }
        else
        {
                return (((hashkey2 (c. right (d - 2)) & 0x55555555u) << 11) ^
                        ((hashkey1 (c. left ((d >> 3) + 1)) & 0x52942529u) << 19) ^
                        ((hashkey2 (c. right (d >> 2)) & 0xAAAAAAAAu) << 5) ^
                        ((hashkey2 (c. left ((d >> 1) + 1)) & 0xA5294A52u) >> 1));
        }
*/
} /* hashkey2 */

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

/// Writes an elementary cube in the text mode to an output stream
/// as the cartesian product of elementary intervals.
template <class CoordT, class WrapT, class EmptyT>
std::ostream &operator << (std::ostream &out,
        const tCubCell<CoordT,WrapT,EmptyT> &c)
{
        int dim = c. spaceDim ();
        out << '[';
        for (int i = 0; i < dim; ++ i)
        {
                if (i)
                        out << "]x[";
                CoordT left = c. left (i);
                CoordT right = c. right (i);
                out << left;
                if (left != right)
                        out << ',' << right;
        }
        out << ']';
        return out;
} /* operator << */

/// Reads an elementary cube from the input stream in the text format.
/// Skips all the comments and spaces in the input stream
/// before checking for an opening parenthesis or bracket.
/// Then reads the cartesian product that defines the cube.
/// If no elementary cube is found in the input then the cube
/// is not modified.
template <class CoordT, class WrapT, class EmptyT>
std::istream &operator >> (std::istream &in,
        tCubCell<CoordT,WrapT,EmptyT> &c)
{
        using chomp::homology::ignorecomments;

        // ignore any comments, spaces, tabs and new-line characters
        ignorecomments (in);

        // check for an opening parenthesis of bracket
        if ((in. peek () != '[') && (in. peek () != '('))
                return in;

        // read the opening parenthesis or bracket if found
        int ch = in. get ();
        ignorecomments (in);

        // consider the case of the empty cubical cell or the empty cube
        if ((in. peek () == ']') || (in. peek () == ')'))
        {
                in. get ();
                const CoordT *null = 0;
                c = tCubCell<CoordT,WrapT,EmptyT> (-1, null, null);
                return in;
        }

        int spaceDim = 0;

        // read a full cube with coordinates in parentheses
        if (ch == '(')
        {
                CoordT coordinates [chomp::homology::MaxBasDim];
                while (1)
                {
                        // read a consecutive coordinate of the cube
                        CoordT &coordinate = coordinates [spaceDim];
                        coordinate = 0;
                        in >> coordinate;
                        ignorecomments (in);
                        ++ spaceDim;

                        // read the separator if any
                        if (in. peek () == ',')
                        {
                                ch = in. get ();
                                ignorecomments (in);
                        }

                        // if this is a closing parenthesis then finish
                        if (in. peek () == ')')
                        {
                                c = tCubCell<CoordT,WrapT,EmptyT> (spaceDim,
                                        coordinates);
                                return in;
                        }
                }
        }

        // prepare arrays for the left and right endpoints of the cell
        CoordT leftArray [chomp::homology::MaxBasDim];
        CoordT rightArray [chomp::homology::MaxBasDim];

        // read the cartesian product of intervals
        while (1)
        {
                CoordT &left = leftArray [spaceDim];
                CoordT &right = rightArray [spaceDim];
                left = 0;
                right = 0;
                in >> left;
                ignorecomments (in);
                if (in. peek () == ']')
                        right = left;
                else if (in. get () == ',')
                {
                        in >> right;
                        ignorecomments (in);
                }
                else
                        throw "Can't read a cubical cell.";
                in. get ();
                ignorecomments (in);
                ++ spaceDim;
                if (in. peek () != 'x')
                        break;
                in. get ();
                ignorecomments (in);
                if (in. get () != '[')
                        throw "An interval expected for a cubical cell.";
                if (spaceDim >= chomp::homology::MaxBasDim)
                        throw "Too high dimension of a cubical cell.";
        }

        // define the cubical cell
        c = tCubCell<CoordT,WrapT,EmptyT> (spaceDim, leftArray, rightArray);

        return in;
} /* operator >> */


#endif // _CHAINCON_CUBCELL_H_