graphs.h

Go to the documentation of this file.

/// @addtogroup unifexp
/// @{

/////////////////////////////////////////////////////////////////////////////
///
/// @file graphs.h
///
/// This file contains the definition of the main procedures for the
/// computations which use graph algorithms.
///
/// @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 March 8, 2007. Last revision: August 22, 2007.

#ifndef _graphs_h_
#define _graphs_h_

#include "rounding.h" // (this header must be included first - BOOST)

#include "chomp/system/config.h"
#include "chomp/system/textfile.h"
#include "chomp/struct/digraph.h"

#include "maptype.h"
#include "parttype.h"


namespace unifexp {

// --------------------------------------------------
// ------------------ create graph ------------------
// --------------------------------------------------

/// Computes the images of intervals in the partition
/// and fills in the corresponding weighted graph.
/// Fills out the table of images of critical neighborhoods.
template <class wType, class numType>
inline void createGraph (chomp::homology::diGraph<wType> &g,
        const mapType<numType> &theMap, const partType<numType> &thePart,
        int partCount, std::vector<int> &critImages)
{
        using namespace chomp::homology;

        // compute the images of the intervals
        for (int i = 0; i < partCount; ++ i)
        {
                // compute the image of the interval [i,i+1]
                const numType &x1 = thePart [i];
                const numType &x2 = thePart [i + 1];
                if (false && sbug. show)
                {
                        sbug << "Img of [" << x1 << "," << x2 << "]... ";
                }
                numType y1, y2;
                theMap. image (x1, x2, y1, y2);

                // cover this image with a series of intervals
                int i1 = thePart. find (y1);
                if (i1 < 0)
                        i1 = 0;
                if ((i1 > 0) && (thePart [i1] == y1))
                        -- i1;
                int i2 = thePart. find (y2);
                if (i2 < partCount)
                        ++ i2;
                if (false && sbug. show)
                {
                        sbug << "(" << y1 << "," << y2 << "): " <<
                                i1 << "-" << i2 << " (";
                        sbug << thePart [i1] << "," <<
                                thePart [i2] << ").\n";
                }

                // add the vertex that corresponds to this interval
                g. addVertex ();

                // if this is part of a critical neighborhood,
                // then add the image intervals to the sources list
                if (thePart. isCritical (i))
                {
                        for (int j = i1; j < i2; ++ j)
                                critImages. push_back (j);
                        continue;
                }

                // add the edges to the graph
                for (int j = i1; j < i2; ++ j)
                {
                        g. addEdge (j, theMap. minLogDerivative
                                (x1, x2, thePart [j], thePart [j + 1]));
                }
        }

        return;
} /* createGraph */


// --------------------------------------------------
// ------------ find lambda(s) and log C ------------
// --------------------------------------------------

/// Computes lambda, C, and lambda0 with the given partition
/// and for the given map.
/// Does not compute those values whose pointers are set to 0.
template <class numType>
inline void findLambdaC (const numType &delta,
        const mapType<numType> &theMap, partType<numType> &thePart,
        int partCount, numType *lambda, numType *logC, numType *lambda0,
        int rigorous, int sparseGraph)
{
        using namespace chomp::homology;

        // create a partition
        thePart. create (theMap, partCount, delta);
        if (false && sbug. show)
        {
                const partType<numType> &thePartConst = thePart;
                numType minWidth = -1, maxWidth = 0;
                bool firstTime = true;
                for (int i = 0; i <= partCount; ++ i)
                {
                        sbug << (i ? "\t" : "Partition:\n\t") <<
                                thePartConst [i] << "\n";
                        if (!i || thePartConst. isCritical (i - 1))
                                continue;
                        const numType width = thePartConst [i] -
                                thePartConst [i - 1];
                        if (firstTime)
                        {
                                minWidth = maxWidth = width;
                                firstTime = false;
                        }
                        else if (width < minWidth)
                                minWidth = width;
                        else if (maxWidth < width)
                                maxWidth = width;
                }
                sbug << "Partition interval widths: from " << minWidth <<
                        " to " << maxWidth << ".\n";
        }

        // create a graph that represents the map on this partition
        sbug << "Creating the graph...\n";
        diGraph<numType> g;
        std::vector<int> critImages;
        createGraph (g, theMap, thePart, partCount, critImages);

        // show the information about the graph that has been created
        sbug << "The graph has " << g. countVertices () << " vertices, " <<
                g. countEdges () << " edges.\n";
        if (false && sbug. show)
        {
                sbug << "The graph:\n";
                g. show (sbug, true);
                sbug << "\n";
        }

        // create a substitute lambda if necessary
        numType lambdaTemp;
        if (!lambda && (lambda0 || logC))
                lambda = &lambdaTemp;

        // prepare a transposed graph if necessary
        diGraph<numType> gT;

        // create an object for rounding in graph algorithms
        tRounding<numType> rounding;

        // compute lambda if necessary
        if (lambda)
        {
                sbug << "Computing the exponent lambda... ";
                if (rigorous)
                {
                        *lambda = g. minMeanCycleWeight (rounding,
                                lambda0 ? &gT : 0);
                }
                else
                        *lambda = g. minMeanCycleWeight (lambda0 ? &gT : 0);
                sbug << *lambda << ".\n";
        }

        // compute lambda0 if requested to
        if (lambda0)
        {
                // find the minimum mean path weight starting at the image
                // of the critical neighborhood
                sbug << "Computing path weights... ";
                numType minStart = 0;
                if (rigorous)
                {
                        minStart = g. minMeanPathWeight (rounding,
                                critImages, critImages. size ());
                }
                else
                {
                        minStart = g. minMeanPathWeight (critImages,
                                critImages. size ());
                }
                sbug << minStart << " and... ";

                // find the minimum mean path weight ending at the critical
                // neighborhood
                std::vector<int> critIntervals;
                int countCrit = theMap. countCritical ();
                for (int i = 0; i < countCrit; ++ i)
                        critIntervals. push_back (thePart. getCritical (i));
                numType minEnd = 0;
                if (rigorous)
                {
                        minEnd = gT. minMeanPathWeight (rounding,
                                critIntervals, critIntervals. size ());
                }
                else
                {
                        minEnd = gT. minMeanPathWeight (critIntervals,
                                critIntervals. size ());
                }
                sbug << minEnd << ".\n";
                numType minPath = (minStart < minEnd) ? minStart : minEnd;
                *lambda0 = (minPath < *lambda) ? minPath : *lambda;
        }

        // compute C if requested to
        if (logC)
        {
                // decrease the weights of the edges by lambda if necessary
                numType zero (0);
                if (!(*lambda == zero))
                {
                        // decrease lambda a tiny bit to compensate for
                        // rounding errors which otherwise cause the creation
                        // of negative loops in the graph
                        if (*lambda < 0)
                                *lambda *= 1.00000000001;
                        else
                                *lambda *= 0.99999999999;
                        sout << "Decreasing lambda to " << *lambda <<
                                " to avoid spurious negative loops.\n";

                        // decrease the weights by the value of lambda
                        int nEdges = g. countEdges ();
                        for (int edge = 0; edge < nEdges; ++ edge)
                        {
                                g. setWeight (edge, rounding. sub_down
                                        (g. getWeight (edge), *lambda));
                        }
                }

                // use Floyd-Warshall or Johnson's algorithm
                // to compute the minimal path weight in the graph
                sbug << "Computing log C... ";
                if (rigorous)
                {
                        *logC = g. minPathWeight (rounding,
                                false, sparseGraph);
                }
                else
                {
                        *logC = g. minPathWeight (false, sparseGraph);
                }
                sbug << *logC << "\n";
        }

        return;
} /* findLambdaC */


// --------------------------------------------------
// ------------------- find delta -------------------
// --------------------------------------------------

/// Finds a possibly small value of delta using the bisection method.
/// The initial values of the arguments of this function are ignored.
/// The actual computations are done for log(1/delta), and the accuracy
/// is  applied to these rescaled numbers.
/// The bisection method is used until the requested accuracy has been
/// reached, and the two returned values contain delta that is too small,
/// and delta that is still good.
template <class numType>
inline void findDeltaBisection (numType lambdaMin, numType resolution,
        const mapType<numType> &theMap, partType<numType> &thePart,
        int partCount, numType &deltaBad, numType &deltaGood,
        bool considerPaths, int rigorous, int sparseGraph)
{
        using namespace chomp::homology;

        sout << "Finding min delta for which " <<
                (considerPaths ? "lambda0" : "lambda") << " > " <<
                lambdaMin << "...\n";

        // compute a reasonable guess for the initial delta
        double width = theMap. rightBound () - theMap. leftBound ();
        int nCrit = theMap. countCritical ();
        for (int i = 0; i <= nCrit; ++ i)
        {
                double left = (i == 0) ? theMap. leftBound () :
                        theMap. criticalPoint (i - 1);
                double right = (i == nCrit) ? theMap. rightBound () :
                        theMap. criticalPoint (i);
                if (right - left < width)
                        width = right - left;
        }

        // WARNING: delta stands for log (1/delta)
        double deltaMin = log (1 / (width / 8));
        double deltaMax = log (1 / (width / 16));
        sbug << "Initial guess for delta: [" << exp (-deltaMax) << ", " <<
                exp (-deltaMin) << "].\n";
        sbug << "Initial guess for log(1/delta): [" << deltaMin << ", " <<
                deltaMax << "].\n";

        // decrease the lower limit for delta until it is good,
        // then increase the upper limit for delta until it is good,
        // and finally use bisection to narrow down the delta scope
        // until the given accuracy is accomplished
        bool deltaMinOK = false;
        bool deltaMaxOK = false;
        numType accuracy = 0;
        numType delta = deltaMin;
        while (1)
        {
                // if the two bounds for delta are already known...
                if (deltaMinOK && deltaMaxOK)
                {
                        // determine the accuracy for computing delta
                        if (!accuracy)
                        {
                                accuracy = deltaMin * resolution;
                                sbug << "Computing log(1/delta) "
                                        "with the accuracy of " <<
                                        accuracy << "...\n";
                        }
                        if (deltaMax - deltaMin < accuracy)
                                break;
                }

                // indicate the current value of delta under consideration
                scon << "[" << deltaMin << "," << deltaMax << "]       \r";
                sbug << "(Current estimate for log(1/delta): " << deltaMin <<
                        " - " << deltaMax << ".)\n";
                sbug << "Checking log(1/delta) = " << delta <<
                        " (delta = " << exp (-delta) << ")...\n";

                // compute the appropriate lambda for this delta
                numType lambda = 0;
                numType *none = 0;
                findLambdaC (exp (-delta), theMap, thePart, partCount,
                        considerPaths ? none : &lambda, none,
                        considerPaths ? &lambda : none,
                        rigorous, sparseGraph);

                // prepare a new delta to verify
                if (!deltaMinOK)
                {
                        if (lambdaMin < lambda)
                        {
                                sbug << "Good min delta. :)\n";
                                deltaMinOK = true;
                                if (deltaMaxOK)
                                        delta = (deltaMin + deltaMax) / 2;
                                else
                                        delta = deltaMax;
                        }
                        else
                        {
                                sbug << "Bad min delta. %(\n";
                                deltaMaxOK = true;
                                // increase the lower limit by twice the
                                // current difference between min and max
                                deltaMin = deltaMin * 3 - deltaMax * 2;
                                deltaMax = delta;
                                delta = deltaMin;
                        }
                }
                else if (!deltaMaxOK)
                {
                        if (lambdaMin < lambda)
                        {
                                sbug << "Good max delta. %)\n";
                                // increase the upper limit by twice the
                                // current difference between min and max
                                deltaMax = deltaMax * 3 - deltaMin * 2;
                                deltaMin = delta;
                                delta = deltaMax;
                        }
                        else
                        {
                                sbug << "Bad max delta. :(\n";
                                deltaMaxOK = true;
                                deltaMax = delta;
                                delta = (deltaMin + deltaMax) / 2;
                        }
                }
                else // if (deltaMinOK && deltaMaxOK)
                {
                        if (lambdaMin < lambda)
                        {
                                sbug << "Good delta. :)\n";
                                deltaMin = delta;
                                delta = (deltaMin + deltaMax) / 2;
                        }
                        else
                        {
                                sbug << "Bad delta. :(\n";
                                deltaMax = delta;
                                delta = (deltaMin + deltaMax) / 2;
                        }
                }
        }

        // show the answer for the log
        sout << "Log(1/delta) is between " << deltaMin << " and " <<
                deltaMax << ".\n";

        // make the right numbers
        deltaBad = exp (-deltaMax);
        deltaGood = exp (-deltaMin);

        // show the actual answer
        sout << "Delta is between " << deltaBad << " and " << deltaGood <<
                ".\n";

        return;
} /* findDeltaBisection */


} // namespace unifexp

#endif // _graphs_h_

/// @}