partcrit.h

Go to the documentation of this file.

/// @addtogroup unifexp
/// @{

/////////////////////////////////////////////////////////////////////////////
///
/// @file partcrit.h
///
/// This file contains the definition of a special partition type in which
/// the size of the intervals close to a few images of critical points
/// is much smaller than in other areas.
///
/// @author Pawel Pilarczyk
///
/////////////////////////////////////////////////////////////////////////////

// Copyright (C) 2007 by Pawel Pilarczyk.
//
// This file is part of my research program 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 August 25, 2007. Last revision: August 29, 2007.

#ifndef _partcrit_h_
#define _partcrit_h_

#include <vector>
#include <string>
#include "chomp/system/textfile.h" // for debugging only
#include "parttype.h"


namespace unifexp {

// --------------------------------------------------
// --------------- critical partition ----------------
// --------------------------------------------------

/// A critical orbit based partition type.
/// The complement of the critical neighborhood is divided into intervals
/// of approximately the same length.
template <class numType>
class partCritical: public partType<numType>
{
public:
        /// Returns the name of the object.
        std::string name () const;

        /// Creates a partition based on the given map, the requested
        /// number of elements in the partition, and the width of the
        /// critical neighborhood.
        void create (const mapType<numType> &theMap, int partCount,
                const numType &delta);

private:
        /// A small helper functin that computes the distance
        /// between two numbers.
        static numType dist (const numType &x, const numType &y);

        /// Fills part of the partition table between the given entries
        /// in a uniform way based on the values at these ends.
        template <class vectType>
        static void fillUniform (vectType &tab, int first, int last);

}; /* class partCritical */

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

template <class numType>
std::string partCritical<numType>::name () const
{
        return std::string ("critical");
} /* partCritical::name */

template <class numType>
inline numType partCritical<numType>::dist (const numType &x,
        const numType &y)
{
        numType diff = x - y;
        return (diff < 0) ? -diff : diff;
} /* partCritical::dist */

template <class numType>
template <class vectType>
void partCritical<numType>::fillUniform (vectType &tab, int first, int last)
{
        const numType &numFirst = tab [first];
        const numType &numLast = tab [last];
        const numType numDiff = (numLast - numFirst) / (last - first);
        for (int i = first + 1; i < last; ++ i)
                tab [i] = numFirst + (i - first) * numDiff;
        return;
} /* partCritical::fillUniform */

template <class numType>
void partCritical<numType>::create (const mapType<numType> &theMap,
        int partCount, const numType &delta)
{
        using chomp::homology::sbug;

        // some parameters of this partition which can be tweaked
        const int coarseFraction = 10;
        const int nImages = 30;
        const int widthFraction = 100;
        const int iterCoef = 3;

        // make sure that the number of partition elements is large enough
        // to create a reasonable partition
        int nCrit = theMap. countCritical ();
        if (partCount < 2 * coarseFraction * (nCrit + 1))
                throw "Too small critical partition requested.";

        // prepare a coarse partition and fill in the two extreme bounds
        const int coarseCount = partCount / coarseFraction;
        std::vector<numType> coarse (coarseCount + 1);
        coarse [0] = theMap. leftBound ();
        coarse [coarseCount] = theMap. rightBound ();

        // compute the width of the entire area to partition, that is,
        // the sum of the widths of intervals in the domain of the map
        // after the critical neighborhood intervals have been removed
        numType width = coarse [coarseCount] - coarse [0] -
                2 * delta * nCrit;
        if (width <= 0)
                throw "Too wide critical neighborhood requested.";

        // compute several images of all the critical points
        std::vector<numType> images (nImages);
        std::vector<int> iterates (nImages);
        int maxIndex = 0;
        int curIndex = -nCrit;
        int iterate = 0;
        int iterateIndex = 0;
        while ((curIndex < maxIndex) && (maxIndex < nImages))
        {
                // update the iterate number
                if (curIndex == iterateIndex)
                {
                        ++ iterate;
                        iterateIndex = maxIndex;
                }

                // take one of the critical points or one of their images
                numType prevPoint = (curIndex < 0) ?
                        theMap. criticalPoint (-curIndex - 1) :
                        images [curIndex];
                ++ curIndex;

                // compute the image of this point
                numType image1, image2;
                theMap. image (prevPoint, prevPoint, image1, image2);
                numType image = (image1 == image2) ? image1 :
                        ((image1 + image2) / 2);

                // make sure that it is inside the domain bounds
                if ((image < theMap. leftBound ()) ||
                        (theMap. rightBound () < image))
                {
                        continue;
                }

                // make sure that it is not inside any critical neighborhood
                bool wrong = false;
                for (int j = 0; j < nCrit; ++ j)
                {
                        if (dist (image, theMap. criticalPoint (j)) < delta)
                        {
                                wrong = true;
                                break;
                        }
                }
                if (wrong)
                        continue;

                // add this image to the list
                images [maxIndex] = image;
                iterates [maxIndex] = iterate;
                ++ maxIndex;
        }

        // show the critical points and their iterates
        if (false && sbug. show)
        {
                sbug << "Critical points and their iterates:\n";
                for (int i = 0; i < maxIndex; ++ i)
                        sbug << std::setw (8) << iterates [i] << ": " <<
                                images [i] << "\n";
        }

        // create an initial coarse uniform partition
        std::vector<int> leftPoints;
        std::vector<int> rightPoints;
        int prev = 0;
        for (int i = 0; i < nCrit; ++ i)
        {
                const numType crit = theMap. criticalPoint (i);
                const numType left = crit - delta;
                const numType right = crit + delta;
                if ((left <= coarse [prev]) ||
                        (coarse [coarseCount] <= right))
                {
                        throw "Too large critical neighborhood requested.";
                }
                int n = prev + static_cast<int> (coarseCount / width *
                        (left - coarse [prev]));
                if (n <= prev)
                        n = prev + 1;
                if (coarseCount <= n + 1)
                        throw "Too few partition intervals requested.";
                coarse [n] = left;
                fillUniform (coarse, prev, n);
                leftPoints. push_back (prev);
                rightPoints. push_back (n);
                coarse [n + 1] = right;
                prev = n + 1;
        }
        fillUniform (coarse, prev, coarseCount);
        leftPoints. push_back (prev);
        rightPoints. push_back (coarseCount);

        // show the coarse partition
        if (false && sbug. show)
        {
                sbug << "Coarse partition (" << (coarseCount + 1) <<
                        " points):\n";
                for (int i = 0; i <= coarseCount; ++ i)
                        sbug << "\t" << coarse [i] << "\n";
        }

        // compute the weights for each interval in the coarse partition
        std::vector<numType> weights (coarseCount);
        int nPoints = leftPoints. size ();
        int cur = 0;
        numType sumWeights = 0;
        const numType widthScale = width / widthFraction;
        for (int n = 0; n < nPoints; ++ n)
        {
                int left = leftPoints [n];
                int right = rightPoints [n];
                for (int i = left; i < right; ++ i)
                {
                        numType midPoint = (coarse [i] + coarse [i + 1]) / 2;
                        numType weight = 0;
                        for (int curInd = 0; curInd < maxIndex; ++ curInd)
                        {
                                numType influence = exp (-dist (midPoint,
                                        images [curInd]) / widthScale);
                                if (iterates [curInd])
                                {
                                        influence /= iterCoef *
                                                iterates [curInd];
                                }
                                weight += influence;
                        }
                        sumWeights += weight;
                        weights [cur] = weight;
                        ++ cur;

                        // show the weight of the midpoint
                        if (false && sbug. show)
                        {
                                sbug << "Weight of " << midPoint << " = " <<
                                        weight << "\n";
                        }
                }
        }

        // allocate the full partition array
        this -> allocate (partCount);

        // create the full partition as a refinement of the coarse partition
        // by subdividing each coarse interval uniformly into a number of
        // intervals proportional to the weight associated with it
        int coarseCur = 0;
        int partCur = 0;
        int nRemaining = partCount - coarseCount;
        for (int n = 0; n < nPoints; ++ n)
        {
                int left = leftPoints [n];
                int right = rightPoints [n];
                for (int i = left; i < right; ++ i)
                {
                        // remember if this is the last interval to fill in
                        bool theLastInterval = (i == right - 1) &&
                                (n == nPoints - 1);
                        // compute the number of points to insert
                        int nInsert = static_cast<int>
                                ((partCount - coarseCount) *
                                weights [coarseCur] / sumWeights + 0.5);
                        if ((nRemaining < nInsert) || theLastInterval)
                                nInsert = nRemaining;
                        nRemaining -= nInsert;

                        // show how many points are inserted to this interval
                        if (false && sbug. show)
                        {
                                sbug << "Inserting " << nInsert <<
                                        " points into " <<
                                        i << " of [" << left << ", " <<
                                        right << "], " << nRemaining <<
                                        " remaining.\n";
                        }

                        // fill in the bounds of the partition segment
                        (*this) [partCur] = coarse [i];
                        (*this) [partCur + nInsert + 1] = coarse [i + 1];
                        fillUniform (*this, partCur, partCur + nInsert + 1);
                        partCur += nInsert + 1;
                        ++ coarseCur;
                }

                // add to the list of critical intervals
                if (n != nPoints - 1)
                {
                        if (false && sbug. show)
                        {
                                sbug << "*Critical: " <<
                                        (*this) [partCur] << "\n";
                        }

                        this -> addCritical (partCur);
                        ++ partCur;
                }
        }

        return;
} /* partCritical::create */


} // namespace unifexp

#endif // _partcrit_h_

/// @}