worker.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file worker.h
///
/// The worker class for the Conley-Morse graphs computation program.
/// This file contains the definition of the worker class which processes
/// a single chunk of the parameter range.
///
/// @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 February 11, 2008. Last revision: April 16, 2008.


#ifndef _CMGRAPHS_WORKER_H_
#define _CMGRAPHS_WORKER_H_


// include some standard C++ header files
#include <string>
#include <sstream>
#include <cstdio>

// include selected header files from the CHomP library
#include "chomp/system/config.h"
#include "chomp/system/textfile.h"
#include "chomp/system/timeused.h"
#include "chomp/cubes/cube.h"
#include "chomp/multiwork/mw.h"
#include "chomp/struct/digraph.h"

// include local header files
#include "config.h"
#include "typedefs.h"
#include "utils.h"
#include "compmdec.h"
#include "plotmdec.h"
#include "eigenval.h"
#include "dataconv.h"
#include "matchdec.h"


// --------------------------------------------------
// ---------------- the worker class -----------------
// --------------------------------------------------

/// The worker class that processes single chunks of data
/// which contain the definition of a selected rectangular area
/// in the set of parameters to process.
class Worker: public chomp::multiwork::mwWorker
{
public:
        /// The default constructor
        Worker (bool _local = false);

private:
        /// A function for the initialization of a worker.
        int Initialize (chomp::multiwork::mwData &data);

        /// A function for processing a piece of data by a worker.
        int Process (chomp::multiwork::mwData &data);

        /// The copy constructor is not allowed.
        Worker (const Worker &) {return;}

        /// The assignment operator is not allowed.
        Worker &operator = (const Worker &) {return *this;}

private:
        /// Is the work being done locally? That is, in the same process
        /// as the coordinator? If so, then the global variables are not
        /// reset, such as pointbase, because this might affect the
        /// coordinator's data.
        bool local;

}; /* class Worker */

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

inline Worker::Worker (bool _local): local (_local)
{
        return;
} /* Worker::Worker */

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

inline int Worker::Initialize (chomp::multiwork::mwData &data)
{
        // ignore the received data and say it was Ok
        return chomp::multiwork::mwOk;
} /* Worker::Initialize */

inline int Worker::Process (chomp::multiwork::mwData &data)
{
        using namespace chomp::homology;
        using namespace chomp::multiwork;

        // prepare a stopwatch to measure the time of processing
        timeused processingTime;

        // reset the variables for computing max image diameter and volume
        theCubMapType::maxImgDiam = 0;
        theCubMapType::maxImgVol = 0;

        // decode the number of the data pack to process
        int dataNumber = 0;
        data >> dataNumber;

        // the parameter space dimension
        int dim = 0;
        data >> dim;
        if (dim != paramDim)
        {
                sout << "! Wrong parameter space dimension (" << dim <<
                        "). Rejecting the data.\n";
                data. Reset ();
                return mwReject;
        }

        // the ranges of the parameter boxes
        parCoord parLeft [paramDim];
        parCoord parRight [paramDim];
        for (int i = 0; i < paramDim; ++ i)
        {
                int n = 0;
                data >> n;
                parLeft [i] = static_cast<parCoord> (n);
                data >> n;
                parRight [i] = static_cast<parCoord> (n);
        }

        // the size of Morse sets to skip the Conley index computation above
        int skipIndices = 0;
        data >> skipIndices;

        // read the file name prefix to which Morse decompositions are saved
        std::string sharePrefix;
        data >> sharePrefix;

        // read the file name prefix to which Morse decompositions are saved
        std::string phaseSpacePrefix;
        data >> phaseSpacePrefix;

        // should the full phase space be plotted in a PNG file?
        bool fullPhaseSpace (false);
        data >> fullPhaseSpace;

        // the ending control number
        int ctrl = 0;
        data >> ctrl;
        if (ctrl != controlNumber)
        {
                data. Reset ();
                sout << "! Data incomplete. Rejecting it.\n";
                return mwReject;
        }

        // compute the number of parameter boxes to process
        int volume = 1;
        for (int i = 0; i < paramDim; ++ i)
        {
                int prevVolume = volume;
                volume *= parRight [i] - parLeft [i];
                if (volume < prevVolume)
                        throw "Too many boxes for processing by a worker.";
        }

        // show a message on which boxes are going to be processed
        if (volume > 1)
        {
                sout << "Processing " << dataNumber << ": " << volume <<
                        " boxes in ";
                for (int i = 0; i < paramDim; ++ i)
                {
                        sout << (i ? "x" : "") << "[" << parLeft [i] <<
                                "," << parRight [i] << "]";
                }
                sout << "...\n";
        }
        // show a message about processing a single parameter box
        else
        {
                sout << "Processing " << dataNumber << ": the box " <<
                        parCube (parLeft, paramDim) << "...\n";
        }


        // prepare the offset and the width of the phase space region
        double offset [spaceDim];
        for (int i = 0; i < spaceDim; ++ i)
                offset [i] = spaceOffset [i];
        double width [spaceDim];
        for (int i = 0; i < spaceDim; ++ i)
                width [i] = spaceWidth [i];

        // prepare arrays for storing the ranges of coordinates
        // in the Morse sets
        spcCoord coordLeftMin [spaceDim];
        spcCoord coordRightMax [spaceDim];
        for (int i = 0; i < spaceDim; ++ i)
        {
                coordLeftMin [i] = 1 << finalDepth;
                coordRightMax [i] = 0;
        }

        // prepare the data for storing results of computations
        // for individual boxes
        parCubes paramBoxes;
        std::vector<theMapType *> theMap (volume);
        std::vector<theCubMapType *> theCubMap0 (volume);
        std::vector<theCubMapType *> theCubMap1 (volume);
        std::vector<theMorseDecompositionType *> morseDec (volume);
        std::vector<bool> morseDecComputed (volume, false);
        std::vector<double> timeUsed (volume, 0);
        std::vector<std::vector<int> > wrongIndices (volume);
        std::vector<std::vector<int> > skippedIndices (volume);
        std::vector<std::vector<int> > attractors (volume);

        // run one set of computations
        int countRead = 0, countWritten = 0, countComputed = 0;
        parRect rect (parLeft, parRight, paramDim);
        const parCoord *curCoord;
        while ((curCoord = rect. get ()) != 0)
        {
                // start measuring the time of processing this box
                chomp::homology::timeused stopwatch;

                // determine if this is an internal boundary point
                bool intBoundaryPoint = internalBoundaryPoint
                        (curCoord, parLeft, parRight, paramSubdiv, paramDim);

                // force reading/writing cached Conley indices
                intBoundaryPoint = true;

                // prepare filenames for cached indices and phase space PNG
                std::string pictureFileName = phaseSpacePrefix +
                        coord2str (curCoord, paramDim) + ".png";
                std::string cacheFileName ((intBoundaryPoint &&
                        !sharePrefix. empty ()) ? (sharePrefix +
                        coord2str (curCoord, paramDim) + ".ind") : "");

                // skip this parameter box if the corresponding index cache
                // file exists and the task was to generate this cache only
                if ((volume == 1) &&
                        (cacheFileName. empty () ||
                        fileExists (cacheFileName. c_str ())) &&
                        (phaseSpacePrefix. empty () ||
                        fileExists (pictureFileName. c_str ())))
                {
                        if (cacheFileName. empty ())
                                sout << "No Conley index cache to compute.";
                        else if (fileExists (cacheFileName. c_str ()))
                                sout << "Conley index cache already exists.";
                        sout << " Skipping the box.\n";
                        break;
                }

                // add the analyzed box to the set of parameter boxes
                int curNumber = paramBoxes. size ();
                parCube paramBox (curCoord, paramDim);
                paramBoxes. add (paramBox);
                if (paramBoxes. size () == curNumber)
                        throw "A repeated box found in the iterations.";

                // prepare the actual parameters for the map
                double leftMapParam [paramCount];
                double rightMapParam [paramCount];
                computeParam (curCoord, leftMapParam, rightMapParam);
                sbug << "Box " << paramBox << ":";
                for (int i = 0; i < paramCount; ++ i)
                {
                        sbug << " [" << leftMapParam [i] << "," <<
                                rightMapParam [i] << "]";
                }
                sbug << ".\n";

                // create the map object and set its parameters
                theMap [curNumber] = new theMapType ();
                theMap [curNumber] -> setParam (leftMapParam, rightMapParam,
                        paramCount);

                // prepare an object of the cubical map
                // at the highest subdivision level to be used
                // in the definition of the Morse decomposition
                theCubMap0 [curNumber] = new theCubMapType (offset, width,
                        1 << finalDepth, *(theMap [curNumber]));
                theCubMap0 [curNumber] -> cache = true;
                theCubMap1 [curNumber] = new theCubMapType (offset, width,
                        1 << (finalDepth + 1), *(theMap [curNumber]));
                theCubMap1 [curNumber] -> cache = true;

                // compute the Morse decomposition
                sbug << "Computing Morse decomposition for " <<
                        paramBox << "...\n";
                bool computed = false;
                morseDec [curNumber] = computeMorseDecomposition
                        (*(theMap [curNumber]),
                        *(theCubMap0 [curNumber]), *(theCubMap1 [curNumber]),
                        offset, width, skipIndices, cacheFileName,
                        computed, wrongIndices [curNumber],
                        skippedIndices [curNumber], attractors [curNumber],
                        false, leftMapParam, rightMapParam, paramCount);
                morseDecComputed [curNumber] = computed;

                // update a counter of computed Morse decompositions
                if (computed)
                {
                        ++ countComputed;
                        if (intBoundaryPoint)
                                ++ countWritten;
                }
                else
                        ++ countRead;

                // update the min and max coordinates of Morse sets
                coordMinMax (*(morseDec [curNumber]),
                        coordLeftMin, coordRightMax);

                // save a picture of the computed Morse decomposition
                if ((spaceDim == 2) && !phaseSpacePrefix. empty () &&
                        !fileExists (pictureFileName. c_str ()))
                {
                        sbug << "Saving '" << pictureFileName << "'...\n";
                        plotMorseDecompositionPNG (*(morseDec [curNumber]),
                                pictureFileName. c_str (),
                                fullPhaseSpace ? (1 << finalDepth) : 0);
                }

                // remember the time used for processing this cube
                timeUsed [curNumber] = stopwatch;

                // indicate the percentage of work completed
                sbug << timeUsed [curNumber] << " sec, " <<
                        ((curNumber + 1) * 100 / volume) << "% done.\n";
                sbug << "--------------------------------\n";
        }

        // show statistics of the computations
        if (volume > 1)
        {
                sout << countRead << " sets of Conley indices read, " <<
                        countComputed << " computed (" <<
                        countWritten << " written).\n";
        }

        // process the computed Morse decompositions
        // to determine the continuation relations between them
        std::vector<std::vector<parCube> > matchClasses;
        if (volume > 1)
        {
                sbug << "Matching " << volume <<
                        " Morse decompositions...\n";
                std::vector<std::vector<int> > noWrongIndices (volume);
                matchMorseDecompositions (paramBoxes, morseDec,
                //      theMap, offset, width, finalDepth, theCubMap0,
                        theCubMap1,
                        ignoreIsolationForContinuation ? noWrongIndices :
                        wrongIndices, parLeft, parRight, matchClasses);

                // show a summary of what has just been computed
                int nClasses = matchClasses. size ();
                if (nClasses == 1)
                        sbug << "There is only one nontrivial class";
                else
                        sbug << "There are " <<  nClasses << " nontrivial classes";
                sbug << " of equivalent Morse decompositions:\n";
                if (nClasses == 0)
                        sbug << 0;
                for (int n = 0; n < nClasses; ++ n)
                        sbug << (n ? "+" : "") << matchClasses [n]. size ();
                sbug << " elements.\n";
        }

        // show computed max image diameter and volume
        sbug << "Max img diam = " << theCubMapType::maxImgDiam <<
                ", max img vol = " << theCubMapType::maxImgVol << ".\n";

        // send back the data containing the result of the processing
        data. Reset ();
        data << dataNumber;

        // send the processing time
        data << static_cast<double> (processingTime);
        processingTime = 0;

        // send the lists of numbers of Morse sets with wrong indices
        data << wrongIndices;

        // send the lists of numbers of Morse sets whose indices were skipped
        data << skippedIndices;

        // send the lists of attractors
        data << attractors;

        // send the cubes with the coordinate ranges
        data << spcCube (coordLeftMin, spaceDim);
        data << spcCube (coordRightMax, spaceDim);

        // send the number of boxes in the data chunk
        int countBoxes = paramBoxes. size ();
        data << countBoxes;

        // send data for each box
        for (int curNumber = 0; curNumber < countBoxes; ++ curNumber)
        {
                // send the coordinates of the box
                data << paramBoxes [curNumber];

                // send the time used for processing this box
                data << timeUsed [curNumber];

                // send the information on whether this Morse decomposition
                // was computed or if it was known before
                data << morseDecComputed [curNumber];

                // send the transitive reduction of the Morse graph
                diGraph<> connGraph;
                morseDec [curNumber] -> makegraph (connGraph);
                diGraph<> morseGraph;
                transitiveReduction (connGraph, morseGraph);
                data << morseGraph;

                // send the sizes of each Morse set
                int nSets = morseDec [curNumber] -> count ();
                for (int n = 0; n < nSets; ++ n)
                        data << (*(morseDec [curNumber])) [n]. size ();

                // send the Conley indices of the Morse sets
                // and the nonzero eigenvalues of each index
                if (nSets != morseGraph. countVertices ())
                        throw "Wrong number of Morse sets encountered.";
                for (int n = 0; n < nSets; ++ n)
                {
                        const theConleyIndexType &index =
                                morseDec [curNumber] -> index (n);
                        data << index;
                        IndexEigenValues eigenValues (index);
                        data << eigenValues;
                }
        }

        // send data regarding the continuation relation between the boxes
        data << matchClasses;

        // send computed max image diameter and volume
        data << theCubMapType::maxImgDiam;
        data << theCubMapType::maxImgVol;

        // delete the allocated Morse decompositions
        for (int i = 0; i < countBoxes; ++ i)
        {
                delete morseDec [i];
                delete theCubMap1 [i];
                delete theCubMap0 [i];
                delete theMap [i];
        }

        // reset the environment of computations
        // unless the work is being done locally
        if (!local)
        {
                sbug << "(Resetting the pointset data.)\n";
                parCube::PointBase::reset ();
                spcCube::PointBase::reset ();
                Cube::PointBase::reset ();
        }

        // add a data chunk separator
        sout << "================================\n";

        // return the result
        return mwOk;
} /* Worker::Process */


#endif // _CMGRAPHS_WORKER_H_

---