chaincon/cubcell.h

Go to the documentation of this file.

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

// Copyright (C) 2009-2011 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: March 4, 2011.


#ifndef _CHAINCON_CUBCELL_H_
#define _CHAINCON_CUBCELL_H_


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

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

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


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

/// An elementary cubical cell with vertices of integer type.
template <class CoordT>
class tCubCell
{
public:
        /// The type of coordinates.
        typedef CoordT CoordType;

        /// 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> &c);

        /// The constructor of the n-th boundary cube.
        tCubCell (const tCubCell<CoordT> &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> &c, int n, int side);

        /// The assignment operator.
        tCubCell<CoordT> &operator = (const tCubCell<CoordT> &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;

        /// The equality operator.
        bool operator == (const tCubCell<CoordT> &s) const;

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

private:
        /// 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.
        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.
        int_t n2;

}; /* class tCubCell */

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

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

template <class CoordT>
template <class CoordArray>
inline tCubCell<CoordT>::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;

        // 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];
        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]);
        n2 = chomp::homology::tPointBase<CoordT>::number (c, dimension);
        n2 |= dimension << chomp::homology::NumBits;

        return;
} /* tCubCell::tCubCell */

template <class CoordT>
template <class CoordArray>
inline tCubCell<CoordT>::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;

        // 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];
        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];
        n2 = chomp::homology::tPointBase<CoordT>::number (c, dimension);
        n2 |= cellDimension << chomp::homology::NumBits;

        return;
} /* tCubCell::tCubCell */

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

template <class CoordT>
inline tCubCell<CoordT>::tCubCell (const tCubCell<CoordT> &c, int n)
{
        // determine the dimension of the embedding space and of the cell
        int spaceDimension = c. spaceDim ();
        int cellDimension = c. dim ();

        // 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);
        if (n >= cellDimension)
                n -= cellDimension;
        int counter = 0;
        bool squeezed = false;
        for (int i = 0; i < spaceDimension; ++ i)
        {
                if (newLeft [i] != newRight [i])
                {
                        if (n == counter)
                        {
                                if (squeezeRight)
                                        newLeft [i] = newRight [i];
                                else
                                        newRight [i] = newLeft [i];
                                squeezed = true;
                                break;
                        }
                        ++ counter;
                }
        }
        if (!squeezed)
                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>
inline tCubCell<CoordT>::tCubCell (const tCubCell<CoordT> &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>
inline tCubCell<CoordT> &tCubCell<CoordT>::operator =
        (const tCubCell<CoordT> &c)
{
        n1 = c. n1;
        n2 = c. n2;
        return *this;
} /* tCubCell::operator = */

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

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

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

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

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

template <class CoordT>
inline int tCubCell<CoordT>::boundaryLength () const
{
        int dimension = dim ();
#ifdef NO_EMPTY_CELL
        return dimension << 1;
#else
        return (dimension > 0) ? (dimension << 1) : 1;
#endif
} /* tCubCell::boundaryLength */

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

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

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

template <class CoordT>
inline void tCubCell<CoordT>::swap (tCubCell<CoordT> &s)
{
        swap (n1, s. n1);
        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>
inline int_t hashkey1 (const tCubCell<CoordT> &c)
{
        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)) & 0x655555u) << 10) ^
                        ((hashkey1 (c. right (1)) & 0xAAAAAAu));
        }
        else
        {
                return ((hashkey1 (c. left (0)) & 0x655555u) << 10) ^
                        ((hashkey1 (c. right (d >> 1)) & 0xAAAAAAu) << 4) ^
                        ((hashkey1 (c. left (d - 1)) & 0xAAAAAAu) >> 6);
        }
} /* 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>
inline int_t hashkey2 (const tCubCell<CoordT> &c)
{
        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)) & 0xAAAAAAu) >> 5) ^
                        ((hashkey2 (c. left (1)) & 0x555555u) << 7);
        }
        else
        {
                return ((hashkey2 (c. right (d - 1)) & 0x555555u) << 9) ^
                        ((hashkey2 (c. right (0)) & 0xAAAAAAu) << 1) ^
                        ((hashkey2 (c. left (d >> 1)) & 0xAAAAAAu) >> 5);
        }
} /* hashkey2 */

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

/// Writes an elementary cube in the text mode to an output stream
/// as the cartesian product of elementary intervals.
template <class CoordT>
std::ostream &operator << (std::ostream &out, const tCubCell<CoordT> &c)
{
        int dim = c. spaceDim ();
        out << '[';
        for (int i = 0; i < dim; ++ i)
        {
                if (i)
                        out << "]x[";
                int left = c. left (i);
                int 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 characters in the input stream until
/// an opening parenthesis is found.
/// 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>
std::istream &operator >> (std::istream &in, tCubCell<CoordT> &c)
{
        // ignore any comments, spaces, tabs and new-line characters
        chomp::homology::ignorecomments (in);

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

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

        // 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> (-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;
                        ++ spaceDim;

                        // read the next character
                        ch = in. get ();

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

                        // if this is an invalid separator then complain
                        if ((ch != ',') && (ch != ' '))
                                throw "Can't read a full cube.";
                }
        }

        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;
                if (in. peek () == ']')
                        right = left;
                else if (in. get () == ',')
                        in >> right;
                else
                        throw "Can't read a cubical cell.";
                in. get ();
                ++ spaceDim;
                if (in. peek () != 'x')
                        break;
                in. get ();
                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> (spaceDim, leftArray, rightArray);

        return in;
} /* operator >> */


#endif // _CHAINCON_CUBCELL_H_