partderiv.h

Go to the documentation of this file.

/// @addtogroup unifexp
/// @{

/////////////////////////////////////////////////////////////////////////////
///
/// @file partderiv.h
///
/// This file contains the definition of a special partition type in which
/// the size of the intervals on which the derivative of f is large
/// 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 November 14, 2007. Last revision: November 14, 2007.

#ifndef _partderiv_h_
#define _partderiv_h_

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


namespace unifexp {

// --------------------------------------------------
// -------------- derivative partition --------------
// --------------------------------------------------

/// A derivative-based partition type.
/// The complement of the critical neighborhood is divided into intervals
/// of approximately the same length.
template <class numType>
class partDerivative: 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:
        /// 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 partDerivative */

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

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

template <class numType>
template <class vectType>
void partDerivative<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;
} /* partDerivative::fillUniform */

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

        // the fraction of all partition points to use for a coarse partition
        const int coarseFraction = 10;

        // 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 derivative-based 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.";

        // 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;
        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 y1 = 0, y2 = 0;
                        theMap. image (midPoint, midPoint, y1, y2);
                        numType deriv = std::exp (theMap. minLogDerivative
                                (midPoint, midPoint, y1, y2));
                        numType weight = std::exp (deriv);
                        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;
} /* partDerivative::create */


} // namespace unifexp

#endif // _partderiv_h_

/// @}