matcharr.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
/// @file matcharr.h
///
/// Array for matching Morse decompositions.
/// This file defines a class for storing a multi-dimensional array
/// whose entries can be efficiently joined into chains which define
/// separate classes of matching entries.
///
/// @author Pawel Pilarczyk
/////////////////////////////////////////////////////////////////////////////

// Copyright (C) 1997-2008 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 2 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, write to the
// Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
// MA 02111-1307, USA.

// Started on November 16, 2006. Last revision: February 19, 2008.


#ifndef _CMGRAPHS_MATCHARR_H_
#define _CMGRAPHS_MATCHARR_H_


// include some standard C++ header files
#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <vector>
#include <map>
#include <algorithm>
#include <cstring>

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


// --------------------------------------------------
// ------------------- MATCHARRAY -------------------
// --------------------------------------------------

/// A multi-dimensional array with links forward/backward between
/// the entries. Note that the type of integers
/// provided as template's argument must be large enough
/// to allow indexing of all the entries in a linear array.
template <class CoordType, class IntType>
class MatchArray
{
public:
        /// The only allowed constructor in which the size of the array
        /// and its dimension must be given.
        /// Optionally the coordinates of the minimal corner may be provided.
        MatchArray (int _dim, const CoordType *_sizes,
                const CoordType *_corner = 0);

        /// The destructor.
        ~MatchArray ();

        /// Joins the chains corresponding to the two parameter n-tuples.
        void join (const CoordType *x, const CoordType *y);

        /// Checks if the two parameters have already been joined, that is,
        /// if they belong to the same chain.
        /// Returns "true" if yes, "false" if not.
        bool joined (const CoordType *x, const CoordType *y) const;

        /// Saves the equivalence classes larger than 1 element
        /// to the provided vector of vectors (the classes are appended).
        /// The type of coordinates of the cube type must match "CoordType".
        template <class CubeType>
        void saveClasses (std::vector<std::vector<CubeType> > &matchClasses)
                const;

        /// Saves the sets to separate files. Sorts the sets by their size.
        void saveSets (const char *prefix) const;

private:
        /// The copy constructor is not allowed.
        MatchArray (const MatchArray<CoordType,IntType> &a);

        /// The assignment operator is not allowed.
        MatchArray &operator = (const MatchArray<CoordType,IntType> &a);

        /// Backward links for chains.
        IntType *link1;

        /// Forward links for chains.
        IntType *link2;

        /// The unique number for each chain. This is stored here
        /// to speed up the verification if two boxes are in the same chain.
        IntType *num;

        /// The size of the multi-dimensional array in each direction.
        CoordType *sizes;

        /// The coordinates of the minimal corner of the array if any.
        CoordType *corner;

        /// The dimension of the array.
        int dim;

        /// The length of the allocated memory arrays.
        int length;

        /// Translates a flat offset array into a multidimensional one.
        void link2coords (IntType link, CoordType *x) const;

        /// Computes the memory offset that corresponds to the given
        /// multidimensional position in the array.
        IntType coords2link (const CoordType *x) const;

        /// A small class whose objects store an IntType object
        /// and an int value. The objects are compared by the value.
        struct pair2
        {
                pair2 (IntType _x, int _value): x (_x), value (_value) {}
                IntType x;
                int value;
                bool operator < (const pair2 &other) const
                {
                        return (this -> value < other. value);
                }
        }; /* struct pair2 */

        /// Retrieves a set of first elements of all the chains
        /// which have at least two elements.
        /// The list is sorted according to the lengths of the chains
        /// in the descending order.
        void getChains (std::vector <MatchArray<CoordType,IntType>::pair2>
                &chains) const;

        /// Saves one set to a file.
        void saveSet (const char *prefix, int counter, int pos,
                CoordType *workspace, int ndigits) const;

        /// Joins two sets, given the offsets in the flat arrays.
        void join (IntType xpos, IntType ypos);

}; /* class MatchArray */

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

template <class CoordType, class IntType>
inline MatchArray<CoordType,IntType>::MatchArray (int _dim,
        const CoordType *_sizes, const CoordType *_corner):
        link1 (0), link2 (0), num (0), sizes (0), corner (0),
        dim (_dim), length (0)
{
        if (!_sizes)
                throw "Undefined array sizes.";
        if (_dim <= 0)
                throw "Non-positive dimension.";
        sizes = new CoordType [dim];
        length = 1;
        for (int i = 0; i < dim; ++ i)
        {
                sizes [i] = _sizes [i];
                long prevLength = length;
                length *= sizes [i];
                if (length < prevLength)
                        throw "Too large matching array requested.";
        }
        if (_corner)
        {
                corner = new CoordType [dim];
                for (int i = 0; i < dim; ++ i)
                        corner [i] = _corner [i];
        }
        link1 = new IntType [length];
        link2 = new IntType [length];
        num = new IntType [length];
        for (int offset = 0; offset < length; ++ offset)
        {
                num [offset] = offset;
                link1 [offset] = -1;
                link2 [offset] = -1;
        }

        if (false) // DEBUG
        {
                CoordType *c = new CoordType [dim];
                for (int i = 0; i < dim; ++ i)
                        c [i] = corner ? corner [i] : 0;
                for (int offset = 0; offset < length; ++ offset)
                {
                        // verify the correctness of the entry
                        if (coords2link (c) != offset)
                                chomp::homology::sbug << "Mismatch at "
                                        "array offset " << offset <<
                                        ": " << coords2link (c) << ".\n";
                        // increase the multi-counter
                        int i = 0;
                        while (1)
                        {
                                if (++ c [i] < sizes [i] +
                                        (corner ? corner [i] : 0))
                                        break;
                                c [i] = corner ? corner [i] : 0;
                                ++ i;
                                if (i >= dim)
                                        break;
                        }
                }
                delete [] c;
        //      chomp::homology::sbug << "Note: Array offsets verified "
        //              "successfully.\n";
        }

        return;
} /* MatchArray::MatchArray */

template <class CoordType, class IntType>
inline MatchArray<CoordType,IntType>::~MatchArray ()
{
        if (corner)
                delete [] corner;
        delete [] sizes;
        delete [] num;
        delete [] link2;
        delete [] link1;
        return;
} /* MatchArray::~MatchArray */

template <class CoordType, class IntType>
inline MatchArray<CoordType,IntType> &MatchArray<CoordType,IntType>::
	operator = (const MatchArray &)
{
        throw "Array assignment operator not allowed.";
        return *this;
} /* MatchArray::operator = */

template <class CoordType, class IntType>
inline MatchArray<CoordType,IntType>::MatchArray
        (const MatchArray<CoordType,IntType> &)
{
        throw "Array copy constructor not allowed.";
        return;
} /* MatchArray::MatchArray */

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

template <class CoordType, class IntType>
inline void MatchArray<CoordType,IntType>::link2coords (IntType link,
        CoordType *x) const
{
        int n = link;
        for (int i = 0; i < dim - 1; ++ i)
        {
                x [i] = n % sizes [i] + (corner ? corner [i] : 0);
                n /= sizes [i];
        }
        x [dim - 1] = n % sizes [dim - 1] + (corner ? corner [dim - 1] : 0);
        return;
} /* MatchArray::link2coords */

template <class CoordType, class IntType>
inline IntType MatchArray<CoordType,IntType>::coords2link
        (const CoordType *x) const
{
        IntType n = x [dim - 1] - (corner ? corner [dim - 1] : 0);
        for (int i = dim - 2; i >= 0; -- i)
        {
                n *= sizes [i];
                n += x [i] - (corner ? corner [i] : 0);
        }
        return n;
} /* MatchArray::coords2link */

template <class CoordType, class IntType>
inline void MatchArray<CoordType,IntType>::join (IntType xpos, IntType ypos)
{
        if (num [xpos] == num [ypos])
                return;
        IntType newnum = num [xpos];
        while (link1 [xpos] >= 0)
                xpos = link1 [xpos];
        while (link2 [ypos] >= 0)
                ypos = link2 [ypos];
        link1 [xpos] = ypos;
        link2 [ypos] = xpos;
        while (1)
        {
                num [ypos] = newnum;
                if (link1 [ypos] < 0)
                        break;
                ypos = link1 [ypos];
        }
        return;
} /* MatchArray::join */

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

template <class CoordType, class IntType>
inline void MatchArray<CoordType,IntType>::saveSet (const char *prefix,
        int counter, int pos, CoordType *workspace, int ndigits) const
{
        // create the file name and open the file
        std::ostringstream fname;
        fname << prefix;
        int ten = 1;
        for (int i = 1; i < ndigits; ++ i)
                ten *= 10;
        while (counter < ten)
        {
                fname << '0';
                ten /= 10;
        }
        fname << counter << ".cub";
        std::ofstream out (fname. str (). c_str ());

        // go to the beginning of the list
        while (link1 [pos] >= 0)
                pos = link1 [pos];

        // write all the elements of the list to the file
        while (pos >= 0)
        {
                link2coords (pos, workspace);
                for (int i = 0; i < dim; ++ i)
                        out << (i ? "," : "(") << workspace [i];
                out << ")\n";
        //      link1 [x] [y] = -1;
        //      IntType link = link2 [pos];
        //      link2 [pos] = -1;
        //      pos = link;
                pos = link2 [pos];
        }
        return;
} /* MatchArray::saveSet */

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

//template <class CoordType, class IntType>
//inline bool operator <
//      (const typename MatchArray<CoordType,IntType>::pair2 &x,
//      const typename MatchArray<CoordType,IntType>::pair2 &y)
//{
//      return (x. value < y. value);
//} /* operator < */

template <class CoordType, class IntType>
inline void MatchArray<CoordType,IntType>::join (const CoordType *x,
        const CoordType *y)
{
        join (coords2link (x), coords2link (y));
        return;
} /* MatchArray::join */

template <class CoordType, class IntType>
inline bool MatchArray<CoordType,IntType>::joined (const CoordType *x,
        const CoordType *y) const
{
        int xpos = coords2link (x);
        if (xpos < 0)
                return false;
        int ypos = coords2link (y);
        if (ypos < 0)
                return false;
        return (num [xpos] == num [ypos]);
} /* MatchArray::joined */

template <class CoordType, class IntType>
inline void MatchArray<CoordType,IntType>::getChains
        (std::vector<MatchArray<CoordType,IntType>::pair2> &chains) const
{
        // create a list of chain beginnings
        chomp::homology::sbug << "Analyzing the sets of matched boxes... ";
        for (IntType cur = 0; cur < length; ++ cur)
        {
                if ((link1 [cur] < 0) && (link2 [cur] >= 0))
                {
                        int setsize = 0;
                        IntType pos = cur;
                        while (pos >= 0)
                        {
                                ++ setsize;
                                pos = link2 [pos];
                        }
                        chains. push_back (pair2 (cur, setsize));
                }
        }
        chomp::homology::sbug << chains. size () << " nontrivial classes.\n";

        // sort the list of chain beginnings
        std::sort (chains. begin (), chains. end ());

        return;
} /* MatchArray::getChains */

template <class CoordType, class IntType>
inline void MatchArray<CoordType,IntType>::saveSets (const char *prefix)
        const
{
        // create a list of chain beginnings
        std::vector<pair2> sets;
        this -> getChains (sets);

        // compute the number of digits to be used in the file names
        int ndigits = 1;
        unsigned ten = 10;
        while (sets. size () >= ten)
        {
                ++ ndigits;
                ten *= 10;
        }

        // write the sets to files
        int counter = 0;
        chomp::homology::sbug << "Saving the nontrivial "
                "equivalence classes... ";
        CoordType *workspace = new CoordType [dim];
        for (int i = sets. size () - 1; i >= 0; -- i)
        {
                saveSet (prefix, ++ counter, sets [i]. x, workspace,
                        ndigits);
        }
        delete [] workspace;
        chomp::homology::sout << "Done.\n";
        
        return;
} /* MatchArray::saveSets */

template <class CoordType, class IntType>
template <class CubeType>
inline void MatchArray<CoordType,IntType>::saveClasses
        (std::vector<std::vector<CubeType> > &matchClasses) const
{
        // create a list of chain beginnings
        std::vector<pair2> sets;
        getChains (sets);

        // create vectors of cubes and add them to the given one
        CoordType *coord = new CoordType [dim];
        std::vector<CubeType> emptyVector;
        for (int i = sets. size () - 1; i >= 0; -- i)
        {
                // go to the beginning of the list
                IntType pos = sets [i]. x;
                while (link1 [pos] >= 0)
                        pos = link1 [pos];

                // create a vector with all the elements of the list
                matchClasses. push_back (emptyVector);
                std::vector<CubeType> &vect =
                        matchClasses [matchClasses. size () - 1];
                vect. resize (sets [i]. value);
                int j = 0;
                while (pos >= 0)
                {
                        link2coords (pos, coord);
                        vect [j ++] = CubeType (coord, dim);
                        pos = link2 [pos];
                }
        }
        delete [] coord;
        return;
} /* MatchArray::saveClasses */


#endif // _CMGRAPHS_MATCHARR_H

---