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-2014 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 3 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,
// please, see <https://www.gnu.org/licenses/>.
 
// Started on February 11, 2008. Last revision: March 29, 2022.
 
 
#ifndef _CMGRAPHS_WORKER_H_
#define _CMGRAPHS_WORKER_H_
 
 
// include some standard C++ header files
#include <string>
#include <sstream>
#include <cstdio>
#include <vector>
#include <string>
 
// 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 "typedyns.h"
#include "utils.h"
#include "compmdec.h"
#include "plotmdec.h"
#include "eigenval.h"
#include "dataconv.h"
#include "datavector.h"
#include "matchdec.h"
#include "mapopt.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);
 
        /// Destructor part for the worker specialization.
        ~Worker ();
 
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;
 
        /// Waiting time measurement between consecutive pieces.
        chomp::homology::timeused waitingTime;
 
        /// Is the worker waiting for anything?
        bool waiting;
 
}; /* class Worker */
 
// --------------------------------------------------
 
inline Worker::Worker (bool _local): local (_local), waiting (true)
{
        waitingTime = 0;
        return;
} /* Worker::Worker */
 
// --------------------------------------------------
 
inline Worker::~Worker ()
{
        using namespace chomp::homology;
 
//DEBUG if (waiting && !local && (static_cast<double> (waitingTime) > 1.0))
        if (waiting && !local)
        {
                sbug << "It took " << waitingTime <<
                        " to decide to quit.\n";
        }
        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;
 
        // display waiting time if considerable
        if (static_cast<double> (waitingTime) > 1.0)
        {
                sbug << "It took " << waitingTime <<
                        " to receive data.\n";
        }
        waiting = false;
 
        // prepare a stopwatch to measure the time of processing
        timeused processingTime;
        processingTime = 0;
 
        // 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 for ODE integration optimization data
        std::string mapOptPrefix;
        data >> mapOptPrefix;
 
        // read the file name prefix to which Morse decompositions are saved
        std::string sharePrefix;
        data >> sharePrefix;
 
        // read the file name prefix to which pictures of the phase space
        // are saved
        std::string phaseSpacePrefix;
        data >> phaseSpacePrefix;
 
        // read the file name prefix to which the list of cubes are saved
        std::string cubesPrefix;
        data >> cubesPrefix;
 
        // read the file name prefix for cached Morse decompositions
        std::string morseDecPrefix;
        data >> morseDecPrefix;
 
        // read the file name prefix for Conley-Morse graphs
        std::string graphsPrefix;
        data >> graphsPrefix;
 
        // read the file name prefix for post-processing information
        std::string procPrefix;
        data >> procPrefix;
 
        // should the full phase space be plotted in a PNG file?
        bool fullPhaseSpace (false);
        data >> fullPhaseSpace;
 
        // receive map optimization data
        std::vector<std::vector<std::string> > mapOptHints;
        data >> mapOptHints;
 
        // 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 an object for storing the ranges of coordinates
        // in the Morse sets
        CoordMinMax coordMinMax;
 
        // 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;
        int countWritten = 0;
        int 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, the phase space PNG,
                // and the sets of cubes to be saved to files
                std::string pictureFileName = phaseSpacePrefix +
                        coord2str (curCoord, paramDim) + ".png";
                std::string morseDecFileName = morseDecPrefix. empty () ?
                        "" : (morseDecPrefix +
                        coord2str (curCoord, paramDim) + ".bz2");
                std::string graphFileName = graphsPrefix. empty () ?
                        "" : (graphsPrefix +
                        coord2str (curCoord, paramDim) + ".txt");
                std::string cacheFileName ((intBoundaryPoint &&
                        !sharePrefix. empty ()) ? (sharePrefix +
                        coord2str (curCoord, paramDim) + ".ind") : "");
                std::string cubesFilePrefix = cubesPrefix. empty () ? "" :
                        (cubesPrefix + coord2str (curCoord, paramDim) +
                        "c");
                std::string procFilePrefix = procPrefix. empty () ? "" :
                        (procPrefix + coord2str (curCoord, paramDim));
                std::string mapOptFileName = mapOptPrefix. empty () ? "" :
                        (mapOptPrefix + coord2str (curCoord, paramDim));
 
                // 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 ())) &&
                        (morseDecPrefix. empty () ||
                        fileExists (morseDecFileName. c_str ())) &&
                        (graphsPrefix. empty () ||
                        fileExists (graphFileName. c_str ())))
                {
                        bool exists = false;
                        if (cacheFileName. empty ())
                                sout << "No Conley index cache to compute. ";
                        else if (fileExists (cacheFileName. c_str ()))
                        {
                                sout << "Conley index cache already exists. ";
                                exists = true;
                        }
                        if (phaseSpacePrefix. empty ())
                                sout << "No phase space picture to save. ";
                        else if (fileExists (pictureFileName. c_str ()))
                        {
                                sout << "Phase space picture already exists. ";
                                exists = true;
                        }
                        if (morseDecPrefix. empty ())
                                sout << "No Morse decomp cache to save. ";
                        else if (fileExists (morseDecFileName. c_str ()))
                        {
                                sout << "Morse decomp cache already exists. ";
                                exists = true;
                        }
                        else if (fileExists (graphFileName. c_str ()))
                        {
                                sout << "C-M graph txt file already exists. ";
                                exists = true;
                        }
#ifndef CONFIG_ONLYPREPROCESSING
                        if (exists)
                                sout << "Skipping the box, "
                                        "because at least one file exists.\n";
                        else
                                sout << "Skipping the box, "
                                        "because the result would be lost.\n";
                        break;
#else
                        if (exists)
                        {
                                sout << "Skipping the box.\n";
                                break;
                        }
                        sout << "Processing the box anyway.\n";
#endif
                }
 
                // 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);
 
                // apply map optimization data
                if (!mapOptHints [curNumber]. empty ())
                {
                        sbug << "Using " <<
                                mapOptHints [curNumber]. size () <<
                                " map optimization hints.\n";
                        for (std::vector<std::string>::const_iterator it =
                                mapOptHints [curNumber]. begin ();
                                it != mapOptHints [curNumber]. end (); ++ it)
                        {
                                theMap [curNumber] -> useOptInfo (*it);
                        }
                }
 
                // 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, finalDepth, *(theMap [curNumber]));
                theCubMap0 [curNumber] -> cache = true;
                theCubMap1 [curNumber] = new theCubMapType (offset, width,
                        1 << (finalDepth + 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, pictureFileName, cubesFilePrefix,
                        morseDecFileName, graphFileName,
                        procFilePrefix, mapOptFileName,
                        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
                int_t nMorseSets = morseDec [curNumber] -> count ();
                for (int_t nSet = 0; nSet < nMorseSets; ++ nSet)
                        coordMinMax ((*(morseDec [curNumber])) [nSet]);
 
                // save a picture of the computed Morse decomposition
                if (
#if (!defined (PLOT_X_COORD) && !defined (PLOT_Y_COORD))
                        (spaceDim == 2) &&
#endif
                        !phaseSpacePrefix. empty () &&
                        !fileExists (pictureFileName. c_str ()))
                {
                        sbug << "Saving '" << pictureFileName << "'...\n";
#ifdef PLOT_X_COORD
                        const int xCoord (PLOT_X_COORD);
#else
                        const int xCoord (0);
#endif
#ifdef PLOT_Y_COORD
                        const int yCoord (PLOT_Y_COORD);
#else
                        const int yCoord (1);
#endif
                        plotMorseDecompositionPNG (*(morseDec [curNumber]),
                                pictureFileName. c_str (), xCoord, yCoord,
                                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);
 
        // 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 (coordMinMax. coordLeftMin, spaceDim);
        data << spcCube (coordMinMax. 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 (static_cast<int_t> (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;
 
        // send the map optimization data
        std::vector<std::string> mapOptFinal;
        for (int i = 0; i < countBoxes; ++ i)
        {
                std::string line = theMap [i] -> getOptInfo ();
                if (!line. empty ())
                        mapOptFinal. push_back (line);
        }
        data << mapOptFinal;
 
        // delete the allocated Morse decompositions
        for (int i = 0; i < countBoxes; ++ i)
        {
                sbug << "Deleting the map for " << paramBoxes [i] << "...\n";
                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";
 
        // reset waiting timer
        waitingTime. reset ();
        waiting = true;
 
        // return the result
        return mwOk;
} /* Worker::Process */
 
 
#endif // _CMGRAPHS_WORKER_H_