covboxes.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file covboxes.h
///
/// A box cover type.
/// This file contains the definition of a class that represents
/// a cover of a bounded rectangular region in R^n by means of open boxes.
/// The class includes a method for computing a cover of any open box
/// or a point by means of the boxes in this class.
///
/// @author Pawel Pilarczyk
///
/////////////////////////////////////////////////////////////////////////////

// Copyright (C) 1997-2010 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 March 13, 2009. Last revision: June 18, 2009.


#ifndef _FINRESDYN_COVBOXES_H_
#define _FINRESDYN_COVBOXES_H_


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

// include some header files from the CHomP library
#include "chomp/struct/hashsets.h"
#include "chomp/cubes/pointset.h"
#include "chomp/cubes/cube.h"

// include local header files
#include "rounding.h"
#include "streams.h"


// --------------------------------------------------
// --------------- dummy box numbers ----------------
// --------------------------------------------------

/// A dummy class for adding box numbers.
class DummyBoxNumbers
{
public:
        /// Forgets the box number.
        void add (int) {};

}; /* class DummyBoxNumbers */


// --------------------------------------------------
// ------------------- box cover --------------------
// --------------------------------------------------

/// A cover of a subset in R^n consisting of open boxes.
template <class NumType>
class tCoverBoxes
{
public:
        /// The type of numbers in the Henon map.
        typedef NumType NumberType;

        /// The constructor of an open cover object.
        tCoverBoxes (const int spaceDim, const NumType *leftCorner,
                const NumType *boxWidths, const int *partCounts,
                const NumType *overlap = 0);

        /// The destructor of an open cover object.
        ~tCoverBoxes ();

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

        /// Returns the number of boxes in the cover.
        int size () const;

        /// Returns the coordinates of the given box.
        /// Throws an exception if the box with this number is undefined.
        template <class ResultType>
        void get (const int n,
                ResultType *leftBounds, ResultType *rightBounds) const;

        /// Returns the coordinates of the left corner of the range box.
        template <class ResultType>
        void getLeftCorner (ResultType *leftCorner) const;

        /// Returns the width of the range box in each direction.
        template <class ResultType>
        void getBoxWidths (ResultType *boxWidths) const;

        /// Covers the given box and adds numbers of all the boxes
        /// in the cover to a list using the provided class object
        /// with its method "add". Adds boxes necessary to cover
        /// the rectangle to the cover if they are not there yet.
        /// Throws an exception if the predefined bounds of the covered
        /// rectangular box have been exceeded.
        /// If requested, comptues a rigorous upper bound for the square
        /// of the diameter of the cover.
        template <class InputType, class BoxNumbers>
        void cover (const InputType *leftBounds,
                const InputType *rightBounds, BoxNumbers &boxNumbers,
                NumType *diameter2 = 0);

        /// Covers a given point with open boxes.
        template <class InputType>
        void cover (const InputType *x);

        /// Forgets the cover elements.
        void forget ();

private:
        /// The dimension of the space.
        int spaceDim;

        /// The lower left corner of the covered rectangle.
        NumType *leftCorner;

        /// The widths of the covered rectangle in all the directions.
        NumType *boxWidths;

        /// The upper right corner of the covered rectangle.
        NumType *rightCorner;

        /// The numbers of parts into which the box is subdivided
        /// in each direction.
        int *partCounts;

        /// The minimal overlap width.
        NumType *overlap;

        /// The type of integer coordinates for boxes that represent
        /// cover elements.
        typedef int CoordType;

        /// A temporary buffer for the coordinates of two cubes.
        CoordType *buffer;

        /// The type of a cube that represents a cover element.
        typedef chomp::homology::tCubeBase<CoordType> CubeType;

        /// The set of boxes in the cover.
        chomp::homology::hashedset<CubeType> boxes;

        /// Computes the real left coordinate in the given direction
        /// of a box represented by integer coordinates.
        NumType leftBound (int coord, int dir) const;

        /// Computes the real right coordinate in the given direction
        /// of a box represented by integer coordinates.
        NumType rightBound (int coord, int dir) const;

        /// Computes the left coordinate in the given direction
        /// of the leftmost box to cover the given real coordinate.
        CoordType leftCoord (const NumType &bound, int dir) const;

        /// Computes the right coordinate in the given direction
        /// of the rightmost box to cover the given real coordinate.
        CoordType rightCoord (const NumType &bound, int dir) const;

        /// The copy constructor should not be used.
        tCoverBoxes (const tCoverBoxes<NumType> &);

        /// The assignment operator should not be used.
        tCoverBoxes<NumType> &operator = (const tCoverBoxes<NumType> &);

}; /* class tCoverBoxes */

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

template <class NumType>
inline tCoverBoxes<NumType>::tCoverBoxes (const int spaceDim,
        const NumType *leftCorner, const NumType *boxWidths,
        const int *partCounts, const NumType *overlap)
{
        // make sure the space dimension is within the correct range
        if ((spaceDim <= 0) || (spaceDim > 100000))
                throw "Wrong space dimension for a box open cover.";
        this -> spaceDim = spaceDim;

        // verify the data and throw an exception if necessary
        for (int i = 0; i < spaceDim; ++ i)
        {
                if (boxWidths [i] <= 0)
                        throw "Trying to use a non-positive box width.";
                if (partCounts [i] <= 0)
                        throw "Trying to use a non-positive no. of parts.";
                if (overlap && (overlap [i] <= 0))
                        throw "Trying to use a too small overlap.";
                if (overlap && (2.00001 * overlap [i] >=
                        boxWidths [i] / partCounts [i]))
                {
                        throw "Trying to use too large overlap.";
                }
        }

        // copy the lower left corner of the covered area
        this -> leftCorner = new NumType [spaceDim];
        for (int i = 0; i < spaceDim; ++ i)
                this -> leftCorner [i] = leftCorner [i];

        // copy the widths of the covered area in each direction
        this -> boxWidths = new NumType [spaceDim];
        for (int i = 0; i < spaceDim; ++ i)
                this -> boxWidths [i] = boxWidths [i];

        // compute an upper bound for the region
        this -> rightCorner = new NumType [spaceDim];
        for (int i = 0; i < spaceDim; ++ i)
        {
                this -> rightCorner [i] = tRounding<NumberType>::add_up
                        (this -> leftCorner [i], this -> boxWidths [i]);
        }

        // copy the number of parts in each direction
        this -> partCounts = new int [spaceDim];
        for (int i = 0; i < spaceDim; ++ i)
                this -> partCounts [i] = partCounts [i];

        // copy or create the overlap size in each direction
        this -> overlap = new NumType [spaceDim];
        if (!overlap)
        {
                const NumType min_number =
                        tRounding<NumType>::min_number ();
                for (int i = 0; i < spaceDim; ++ i)
                        this -> overlap [i] = min_number;
        }
        else
        {
                for (int i = 0; i < spaceDim; ++ i)
                        this -> overlap [i] = overlap [i];
        }

        // allocate a buffer for temporary point coordinates
        this -> buffer = new CoordType [2 * spaceDim];

        return;
} /* tCoverBoxes::tCoverBoxes */

template <class NumType>
inline tCoverBoxes<NumType>::tCoverBoxes (const tCoverBoxes<NumType> &)
{
        throw "Trying to use the copy constructor for an open cover.";
        return;
} /* tCoverBoxes::tCoverBoxes */

template <class NumType>
inline tCoverBoxes<NumType> &tCoverBoxes<NumType>::operator =
        (const tCoverBoxes<NumType> &)
{
        throw "Trying to use the assignment operator for an open cover.";
        return *this;
} /* tCoverBoxes::operator = */

template <class NumType>
inline tCoverBoxes<NumType>::~tCoverBoxes ()
{
        delete [] leftCorner;
        delete [] boxWidths;
        delete [] rightCorner;
        delete [] partCounts;
        delete [] overlap;
        delete [] buffer;
        return;
} /* tCoverBoxes::~tCoverBoxes */

template <class NumType>
inline int tCoverBoxes<NumType>::dim () const
{
        return spaceDim;
} /* tCoverBoxes::dim */

template <class NumType>
inline int tCoverBoxes<NumType>::size () const
{
        return boxes. size ();
} /* tCoverBoxes::size */

template <class NumType>
inline void tCoverBoxes<NumType>::forget ()
{
        chomp::homology::hashedset<CubeType> empty;
        boxes = empty;
        return;
} /* tCoverBoxes::forget */

template <class NumType>
inline NumType tCoverBoxes<NumType>::leftBound (int coord, int dir) const
{
        if (coord == 0)
                return leftCorner [dir];
        typedef tRounding<NumType> rnd;
        NumType left1 = rnd::mul_down (boxWidths [dir], coord);
        NumType left2 = rnd::div_down (left1, partCounts [dir]);
        return rnd::add_down (leftCorner [dir], left2);
} /* tCoverBoxes::leftBound */

template <class NumType>
inline NumType tCoverBoxes<NumType>::rightBound (int coord, int dir) const
{
        typedef tRounding<NumType> rnd;
        if (coord == partCounts [dir])
                return rnd::add_up (leftCorner [dir], boxWidths [dir]);
        NumType right1 = rnd::mul_up (boxWidths [dir], coord + 1);
        NumType right2 = rnd::div_up (right1, partCounts [dir]);
        NumType right3 = rnd::add_up (right2, overlap [dir]);
        return rnd::add_up (leftCorner [dir], right3);
} /* tCoverBoxes::rightBound */

template <class NumType>
inline typename tCoverBoxes<NumType>::CoordType
        tCoverBoxes<NumType>::leftCoord
        (const NumType &bound, int dir) const
{
        // make sure that the left end is within the region
        if (bound < leftCorner [dir])
                throw "Trying to cover a box that is too large (left).";

        // return the left corner of the region if this is the case
        if (bound == leftCorner [dir])
                return 0;

        // compute a candidate for the left corner coordinate
        NumType offset = bound - leftCorner [dir];
        int coord = static_cast<CoordType>
                (offset * partCounts [dir] / boxWidths [dir]);
        if (coord < 0)
                return 0;

        // verify if the bound could be improved
        if ((coord + 1 < partCounts [dir]) &&
                (leftBound (coord + 1, dir) <= bound))
        {
                sbug << "* Left bound inc improve.\n";
                ++ coord;
        }

        // if the real lower bound that corresponds to this box
        // is too large then decrease the coordinate of the box
        if ((coord > 0) && (leftBound (coord, dir) > bound))
        {
                sbug << "* Left bound dec cover.\n";
                -- coord;
        }

        // if the box to the left overlaps this bound then take that box
        if ((coord > 0) && (rightBound (coord - 1, dir) > bound))
        {
                sbug << "* Left bound dec overlap.\n";
                -- coord;
        }

        return coord;
} /* tCoverBoxes::leftCoord */

template <class NumType>
inline typename tCoverBoxes<NumType>::CoordType
        tCoverBoxes<NumType>::rightCoord
        (const NumType &bound, int dir) const
{
        // make sure that the right end is within the region
        if (bound > rightCorner [dir])
                throw "Trying to cover a box that is too large (right).";

        // return the right corner of the region if this is the case
        if (bound == rightCorner [dir])
                return partCounts [dir];

        // compute a candidate for the right corner coordinate
        NumType offset = bound - leftCorner [dir] - overlap [dir];
        int coord = (offset > 0) ? (static_cast<CoordType>
                (offset * partCounts [dir] / boxWidths [dir]) + 1) : 1;
        if (coord >= partCounts [dir])
                return partCounts [dir];

        // if the bound could be improved then do it
        if ((coord > 1) && (rightBound (coord - 2, dir) >= bound))
        {
                sbug << "* Right bound inc improve.\n";
                -- coord;
        }

        // if the actual right bound of this box is below the bound
        // then take the next box to the right
        if ((coord < partCounts [dir]) &&
                (rightBound (coord - 1, dir) < bound))
        {
                sbug << "* Right bound inc cover.\n";
                ++ coord;
        }

        // if the box to the right overlaps this bound then take that box
        if ((coord < partCounts [dir]) &&
                (leftBound (coord, dir) < bound))
        {
                sbug << "* Right bound inc overlap.\n";
                ++ coord;
        }

        return coord;
} /* tCoverBoxes::rightCoord */

template <class NumType>
template <class ResultType>
inline void tCoverBoxes<NumType>::get (const int n,
                ResultType *leftBounds, ResultType *rightBounds) const
{
        // extract the rounding class to shorten the notation
        typedef tRounding<ResultType> rnd;

        // make sure that the box number is within the correct range
        if ((n < 0) || (n >= boxes. size ()))
                throw "Wrong cover element requested.";

        // extract the integer coordinates of the requested box
        CoordType *boxCoord = buffer;
        boxes [n]. coord (boxCoord);

        // compute the real coordinates of the box
        for (int i = 0; i < spaceDim; ++ i)
        {
                leftBounds [i] = leftBound (boxCoord [i], i);
                rightBounds [i] = rightBound (boxCoord [i], i);
        }
        return;
} /* tCoverBoxes::get */

template <class NumType>
template <class ResultType>
inline void tCoverBoxes<NumType>::getLeftCorner (ResultType *leftCorner)
        const
{
        for (int i = 0; i < spaceDim; ++ i)
                leftCorner [i] = this -> leftCorner [i];
        return;
} /* tCoverBoxes::getLeftCorner */

template <class NumType>
template <class ResultType>
inline void tCoverBoxes<NumType>::getBoxWidths (ResultType *boxWidths)
        const
{
        for (int i = 0; i < spaceDim; ++ i)
                boxWidths [i] = this -> boxWidths [i];
        return;
} /* tCoverBoxes::getBoxWidths */

template <class NumType>
template <class InputType, class BoxNumbers>
inline void tCoverBoxes<NumType>::cover (const InputType *leftBounds,
        const InputType *rightBounds, BoxNumbers &boxNumbers,
        NumType *diameter2)
{
        // prepare arrays of box coordinates
        CoordType *left = buffer;
        CoordType *right = buffer + spaceDim;

        // compute the lower left corner of the integer box coordinates
        for (int i = 0; i < spaceDim; ++ i)
        {
                left [i] = leftCoord (leftBounds [i], i);
                right [i] = rightCoord (rightBounds [i], i);
        }

        // compute an upper bound for the diameter of the covering set
        if (diameter2)
        {
                NumType squares = 0;
                for (int i = 0; i < spaceDim; ++ i)
                {
                        NumType diff = tRounding<NumType>::sub_up
                                (rightBound (right [i], i),
                                leftBound (left [i], i));
                        NumType diff2 = tRounding<NumType>::mul_up
                                (diff, diff);
                        squares = tRounding<NumType>::add_up
                                (squares, diff2);
                }
                //using std::sqrt;
                *diameter2 = squares;
        }

        // add all the boxes in the computed range to the set
        chomp::homology::tRectangle<CoordType> r (left, right, spaceDim);
        const CoordType *c;
        while ((c = r. get ()) != 0)
        {
                CubeType q (c, spaceDim);
                int n = boxes. getnumber (q);
                if (n >= 0)
                        boxNumbers. add (n);
                else
                {
                        boxNumbers. add (boxes. size ());
                        boxes. add (q);
                }
        }

        return;
} /* tCoverBoxes::cover */

template <class NumType>
template <class InputType>
inline void tCoverBoxes<NumType>::cover (const InputType *x)
{
        // extract the rounding class to shorten the notation
        typedef tRounding<NumType> rnd;

        // prepare a possibly small open interval that contains the point
        NumType *xLeft = new NumType [spaceDim];
        NumType *xRight = new NumType [spaceDim];
        NumType min_number = rnd::min_number ();
        for (int i = 0; i < spaceDim; ++ i)
        {
                xLeft [i] = rnd::add_down (x [i], min_number);
                xRight [i] = rnd::add_up (x [i], min_number);
        }

        // cover this open interval
        DummyBoxNumbers d;
        this -> cover (xLeft, xRight, d);

        // clean up the memory and return
        delete [] xRight;
        delete [] xLeft;
        return;
} /* tCoverBoxes::cover */

template <class NumType>
inline std::ostream &operator << (std::ostream &out,
        const tCoverBoxes<NumType> &cover)
{
        int dim = cover. dim ();
        NumType *buffer = new NumType [dim << 1];
        NumType *leftBound = buffer;
        NumType *rightBound = buffer + dim;
        for (int n = 0; n < cover. size (); ++ n)
        {
                cover. get (n, leftBound, rightBound);
                for (int i = 0; i < dim; ++ i)
                {
                        out << (i ? " x [" : "[") << leftBound [i] <<
                                "," << rightBound [i] << "]";
                }
                out << "\n";
        }
        delete [] buffer;
        return out;
} /* operator << */


#endif // _FINRESDYN_COVBOXES_H_