simplex.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// \file
///
/// A simplex class with arbitrary vertices.
///
/////////////////////////////////////////////////////////////////////////////

// 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: August 4, 2015.


#ifndef _CHAINCON_SIMPLEX_H_
#define _CHAINCON_SIMPLEX_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 relevant local header files
#include "chaincon/emptycell.h"


// --------------------------------------------------
// --------------------- simplex --------------------
// --------------------------------------------------

/// A simplex with vertices of arbitrary type.
/// Note that the vertex type must have hash key functions
/// (defined in chomp/struct/hashsets.h for most built-in integers,
/// defined in chaincon/stringhash.h for std::string, or your own).
/// 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.
/// Important: The vertices are assumed to be given in the ascending order;
/// otherwise the operator == may give false negatives!
template <class VertexT, class EmptyT>
class tSimplex
{
public:
        /// The type of the vertex.
        typedef VertexT VertexType;

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

        /// The default constructor of an empty simplex.
        tSimplex ();

        /// The constructor of a simplex from an array of vertices.
        template <class VertexArray>
        tSimplex (int dimension, const VertexArray &vertices);

        /// The copy constructor.
        tSimplex (const tSimplex<VertexT,EmptyT> &s);

        /// The constructor of the n-th face of a simplex (if n >= 0)
        /// or the k-th degeneracy operator (if n = -1 - k < 0).
        tSimplex (const tSimplex<VertexT,EmptyT> &s, int n);

        /// The assignment operator.
        tSimplex<VertexT,EmptyT> &operator =
                (const tSimplex<VertexT,EmptyT> &s);

        /// The destructor.
        ~tSimplex ();

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

        /// Returns the n-th vertex of the simplex.
        const VertexT &operator [] (int n) const;

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

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

        /// Checks if the simplex is degenerate, that is, if it has two
        /// identical vertices one after another. If this is the case
        /// then returns the position of the second vertex.
        /// Returns 0 if the simplex is not degenerate.
        int degenerate () const;

        /// The equality operator.
        /// Verifies if the vertices are the same and in the same order.
        /// Note that this may result in false negatives if the vertices
        /// are stored in various orders.
        bool operator == (const tSimplex<VertexT,EmptyT> &s) const;

        /// Swaps the data between two simplices.
        void swap (tSimplex<VertexT,EmptyT> &s);

private:
        /// The dimension of the simplex.
        int dimension;

        /// An array of vertices of the simplex.
        VertexT *vertices;
}; /* class tSimplex */

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

template <class VertexT, class EmptyT>
inline tSimplex<VertexT,EmptyT>::tSimplex ():
        dimension (-1), vertices (0)
{
        return;
} /* tSimplex::tSimplex */

template <class VertexT, class EmptyT>
template <class VertexArray>
inline tSimplex<VertexT,EmptyT>::tSimplex (int dimension,
        const VertexArray &vertices)
{
        // make sure that the vertices are in the ascending order
        for (int i = 0; i < dimension; ++ i)
        {
                if (vertices [i + 1] < vertices [i])
                        throw "Wrong order of vertices of a simplex.";
        }

        // allocate memory for the simplex and store its dimension
        this -> dimension = (dimension >= 0) ? dimension : -1;
        this -> vertices = (dimension >= 0) ? new VertexT [dimension + 1] : 0;

        // copy the values of the vertices to the inernal array
        for (int i = 0; i <= dimension; ++ i)
                this -> vertices [i] = vertices [i];
        return;
} /* tSimplex::tSimplex */

template <class VertexT, class EmptyT>
inline tSimplex<VertexT,EmptyT>::tSimplex (const tSimplex<VertexT,EmptyT> &s)
{
        dimension = s. dimension;
        vertices = (dimension >= 0) ? new VertexT [dimension + 1] : 0;
        for (int i = 0; i <= dimension; ++ i)
                vertices [i] = s. vertices [i];
        return;
} /* tSimplex::tSimplex */

template <class VertexT, class EmptyT>
inline tSimplex<VertexT,EmptyT>::tSimplex (const tSimplex<VertexT,EmptyT> &s,
        int n)
{
        if (n >= 0)
        {
                dimension = (s. dimension >= 0) ? (s. dimension - 1) : -1;
                vertices = (dimension >= 0) ? new VertexT [dimension + 1] : 0;
                for (int i = 0; i < n; ++ i)
                        vertices [i] = s. vertices [i];
                for (int i = n; i <= dimension; ++ i)
                        vertices [i] = s. vertices [i + 1];
        }
        else
        {
                dimension = (s. dimension >= 0) ? (s. dimension + 1) : -1;
                vertices = (dimension >= 0) ? new VertexT [dimension + 1] : 0;
                for (int i = 0; i < -n; ++ i)
                        vertices [i] = s. vertices [i];
                for (int i = -n; i <= dimension; ++ i)
                        vertices [i] = s. vertices [i - 1];
        }
        return;
} /* tSimplex::tSimplex */

template <class VertexT, class EmptyT>
inline tSimplex<VertexT,EmptyT> &tSimplex<VertexT,EmptyT>::operator =
        (const tSimplex<VertexT,EmptyT> &s)
{
        if (dimension != s. dimension)
        {
                if (vertices)
                        delete [] vertices;
                dimension = s. dimension;
                vertices = (dimension >= 0) ?
                        new VertexT [dimension + 1] : 0;
        }
        for (int i = 0; i <= dimension; ++ i)
                this -> vertices [i] = s. vertices [i];
        return *this;
} /* tSimplex::operator = */

template <class VertexT, class EmptyT>
inline tSimplex<VertexT,EmptyT>::~tSimplex ()
{
        if (vertices)
                delete [] vertices;
        return;
} /* tSimplex::~tSimplex */

template <class VertexT, class EmptyT>
inline int tSimplex<VertexT,EmptyT>::dim () const
{
        return dimension;
} /* tSimplex::dim */

template <class VertexT, class EmptyT>
inline const VertexT &tSimplex<VertexT,EmptyT>::operator [] (int n) const
{
        return vertices [n];
} /* tSimplex::operator [] */

template <class VertexT, class EmptyT>
inline int tSimplex<VertexT,EmptyT>::boundaryLength () const
{
        if (EmptyT::exists ())
                return dimension + 1;
        else
                return (dimension > 0) ? (dimension + 1) : 0;
} /* tSimplex::boundaryLength */

template <class VertexT, class EmptyT>
inline int tSimplex<VertexT,EmptyT>::boundaryCoef (int n) const
{
        return (n & 1) ? -1 : 1;
} /* tSimplex::boundaryCoef */

template <class VertexT, class EmptyT>
inline int tSimplex<VertexT,EmptyT>::degenerate () const
{
        for (int i = 1; i <= dimension; ++ i)
        {
                if (vertices [i - 1] == vertices [i])
                        return i;
        }
        return 0;
} /* tSimplex::degenerate */

template <class VertexT, class EmptyT>
inline bool tSimplex<VertexT,EmptyT>::operator ==
        (const tSimplex<VertexT,EmptyT> &s) const
{
        if (dimension != s. dimension)
                return false;
        for (int i = 0; i <= dimension; ++ i)
        {
                if (!(vertices [i] == s. vertices [i]))
                        return false;
        }
        return true;
} /* tSimplex::operator == */

template <class VertexT, class EmptyT>
inline void tSimplex<VertexT,EmptyT>::swap (tSimplex<VertexT,EmptyT> &s)
{
        std::swap (dimension, s. dimension);
        std::swap (vertices, s. vertices);
        return;
} /* tSimplex::swap */

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

/// Generates a hashing key no. 1 for a simplex,
/// composed of hashing keys for the vertices.
/// This key is to be used in a hashed set.
template <class VertexT, class EmptyT>
inline int_t hashkey1 (const tSimplex<VertexT,EmptyT> &s)
{
//      using chomp::homology::hashkey1;
        using ::hashkey1;
        int d = s. dim ();
        if (d < 0)
                return 0;
        else if (d == 0)
                return hashkey1 (s [0]) << 2;
        else if (d == 1)
        {
                return ((hashkey1 (s [0]) & 0x655555u) << 11) ^
                        ((hashkey1 (s [1]) & 0xAA00AAu) >> 1);
        }
        else
        {
                return ((hashkey1 (s [0]) & 0x655555u) << 15) ^
                        ((hashkey1 (s [d >> 1]) & 0xAA00AAu) << 1) ^
                        ((hashkey1 (s [d]) & 0xAAAAAAu) >> 7);
        }
} /* hashkey1 */

/// Generates a hashing key no. 2 for a simplex,
/// composed of hashing keys for the vertices.
/// This key is to be used in a hashed set.
template <class VertexT, class EmptyT>
inline int_t hashkey2 (const tSimplex<VertexT,EmptyT> &s)
{
//      using chomp::homology::hashkey2;
        using ::hashkey2;
        int d = s. dim ();
        if (d < 0)
                return 0;
        else if (d == 0)
                return hashkey2 (s [0]) << 4;
        else if (d == 1)
        {
                return ((hashkey2 (s [0]) & 0xAAAAAAu) >> 2) ^
                        ((hashkey2 (s [1]) & 0x555555u) << 8);
        }
        else
        {
                return ((hashkey2 (s [d]) & 0x555555u) << 7) ^
                        ((hashkey2 (s [0]) & 0xAA00AAu) << 5) ^
                        ((hashkey2 (s [d >> 1]) & 0xAAAAu) >> 5);
        }
} /* hashkey2 */

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

/// Writes a simplex in the text mode to an output stream.
/// All the vertices are written using their operator <<,
/// and are enclosed in parentheses in a comma-separated list.
template <class VertexT, class EmptyT>
std::ostream &operator << (std::ostream &out, const tSimplex<VertexT,EmptyT> &s)
{
        int dim = s. dim ();
        out << "(";
        for (int i = 0; i <= dim; ++ i)
        {
                if (i)
                        out << ',';
                out << s [i];
        }
        out << ")";
        return out;
} /* operator << */

/// Reads a simplex 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 vertices
/// of the simplex using their operator >> as a comma-separated list.
/// The closing parenthesis indicates the end of the simplex definition.
/// If the closing parenthesis is not found then an exception is thrown.
/// If no simplex is found in the input then the simplex is not modified.
template <class VertexT, class EmptyT>
std::istream &operator >> (std::istream &in, tSimplex<VertexT,EmptyT> &s)
{
        // ignore any comments, spaces, tabs and new-line characters
        chomp::homology::ignorecomments (in);

        // read the opening parenthesis
        int ch = in. get ();
        if (ch != '(')
                return in;

        // read a comma-separated list of vertices of the simplex
        std::vector<VertexT> vertices;
        ch = in. peek ();
        while ((ch != EOF) && (ch != ')'))
        {
                VertexT vertex;
                in >> vertex;
                vertices. push_back (vertex);
                while ((in. peek () == ' ') || (in. peek () == ','))
                        ch = in. get ();
                if (in. peek () == ')')
                        ch = in. get ();
        }
        if (ch != ')')
                throw "Simplex reading error: ')' expected.\n";

        // make sure that the vertices are in the ascending order
        for (unsigned i = 1; i < vertices. size (); ++ i)
        {
                if (vertices [i - 1] >= vertices [i])
                        throw "Simplex reading error: wrong vertices.\n";
        }

        // define the simplex
        s = tSimplex<VertexT,EmptyT> (vertices. size () - 1, vertices);

        return in;
} /* operator >> */


#endif // _CHAINCON_SIMPLEX_H_