cmsingle.cpp

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file cmsingle.cpp
///
/// Computing a single Conley-Morse graph.
/// This program constitutes an interface to computing a Conley-Morse graph
/// for a single parameter box the same way it would be computed by the
/// program "cmgraphs". Its purpose is to extract specific data
/// for a single parameter box after having located it with the aid
/// of the results obtained by the program "cmgraphs", in particular,
/// by having analyzed the continuation diagram.
/// Please, use the program "cmgraphs" to run computations for a range
/// of parameters and to do continuation matching
/// between the Morse decompositions.
///
/// @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 26, 2008. Last revision: July 15, 2014.
 
 
// include some standard C++ header files
#include <exception>
#include <sstream>
#include <fstream>
#include <algorithm>
#include <new>
 
// 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/system/arg.h"
 
// include local header files
#include "config.h"
#include "confinfo.h"
#include "typedefs.h"
#include "typedyns.h"
#include "utils.h"
#include "compmdec.h"
#include "plotmdec.h"
#include "eigenval.h"
#include "dotgraph.h"
#include "mapopt.h"
 
 
// --------------------------------------------------
// -------------------- OVERTURE --------------------
// --------------------------------------------------
 
/// The title of the program which is displayed every time
/// the program is launcued.
const char *title = "\
This is a program which computes a single Conley-Morse graph.\n\
Ver. 0.06, August 22, 2014. Copyright (C) 2005-2014 by Pawel Pilarczyk.\n\
This is free software. No warranty. Consult 'license.txt' for details.";
 
/// The help information about the program which is displayed
/// if the program is launched without command-line arguments
/// or with incorrect arguments.
const char *helpinfo = "\
This is the front-end program for the computation of a Conley-Morse graph\n\
for a single parameter box. It is a supplementary program to 'cmgraphs',\n\
and is capable of saving the results of all its computations to files.\n\
The parameter box is identified by the integer coordinates of its minimal\n\
vertex and corresponds to a box in the subdivision of the parameter area,\n\
for example, if the entire parameter range [-1,1]x[-1,1] is divided into\n\
20x20 boxes then '-b 0,19' corresponds to [-1.0,-0.9]x[0.9,1.0].\n\
Command line arguments (the argument '-b' is mandatory):\n\
-b x1,x2,...,xN - the integer coordinates of a parameter box to process,\n\
-g file.txt - file name to save the text code for the Conley-Morse graph,\n\
-s file.ind - file name to save/load cached Conley index information,\n\
-p file.png - file name to save a PNG picture of the Morse decomposition,\n\
-q prefix - prefix for file names to save Morse sets to as lists of cubes,\n\
-o prefix - file name prefix for post-processing Morse sets information,\n\
-t filename - file name with optimization information (read & append),\n\
-u prefix - file name prefix for ODE integration optimization information,\n\
-z file.bz2 - name of a temporary file with Morse decomposition cache,\n\
--conn - compute connecting orbits between Morse sets (memory-consuming),\n\
--full - use the full range of the phase space to plot Morse sets,\n\
--log filename - file name to save all the text output of the program to,\n\
--quiet - suppress any output written to the terminal,\n\
--debug - display additional information useful for debugging,\n\
--help - make the program show this brief help information and exit.\n\
For more information ask the author at http://www.pawelpilarczyk.com/.";
 
 
// --------------------------------------------------
// --------------- Single Computation ---------------
// --------------------------------------------------
 
int runSingleComp (const char *boxCoordinates,
        const char *mapOptFileName, const char *mapOptFileLocal,
        const char *graphName, const char *shareName,
        const char *phaseSpaceName, bool fullPhaseSpace,
        const char *cubesFilePrefix, const char *morseDecName,
        const char *procFilePrefix,
        bool connOrbits, int skipIndices)
{
        using chomp::homology::sout;
        using chomp::homology::sbug;
 
        // show the global configuration information
        std::ostringstream info;
        showConfigInfo (info, 0, 0);
        sout << "\n" << info. str () << "\n";
 
        // reset the variables for computing max image diameter and volume
        theCubMapType::maxImgDiam = 0;
        theCubMapType::maxImgVol = 0;
 
        // decode the parameter box
        parCube paramBox;
        std::string boxCoordStr ("(" + std::string (boxCoordinates) + ")");
        std::istringstream boxCoordStream (boxCoordStr);
        boxCoordStream >> paramBox;
        sout << "Computing Morse decomposition for " << paramBox << ".\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];
 
        // set up an array of coordinates of the box
        parCoord curCoord [paramDim];
        paramBox. coord (curCoord);
 
        // prepare the actual parameters for the map
        double leftMapParam [paramCount];
        double rightMapParam [paramCount];
        computeParam (curCoord, leftMapParam, rightMapParam);
        sout << "Parameter values:";
        for (int i = 0; i < paramCount; ++ i)
        {
                sout << " [" << leftMapParam [i] << "," <<
                        rightMapParam [i] << "]";
        }
        sout << "\n";
 
        // create the map object and set its parameters
        theMapType theMap;
        theMap. setParam (leftMapParam, rightMapParam, paramCount);
 
        // read map optimization data and apply it to the map
        MapOptimization mapOptimization;
        size_t prevSize (mapOptimization. size ());
        if (mapOptFileName && *mapOptFileName)
        {
                std::ifstream in (mapOptFileName);
                if (in)
                {
                        sout << "Reading map optimization data... ";
                        mapOptimization. loadData (in);
                        sout << (mapOptimization. size () - prevSize) <<
                                " items added.\n";
 
                        sout << "Optimizing the map... ";
                        std::vector<std::string> data;
                        mapOptimization. prepareData (data,
                                leftMapParam, rightMapParam, paramCount);
                        sout << data. size () << " hints... ";
                        for (std::vector<std::string>::const_iterator it =
                                data. begin (); it != data. end (); ++ it)
                        {
                                theMap. useOptInfo (*it);
                        }
                        sout << "Done.\n";
                }
        }
 
        // prepare an object of the cubical map
        // at the highest subdivision level to be used
        // in the definition of the Morse decomposition
        theCubMapType theCubMap0 (offset, width, 1 << finalDepth,
                finalDepth, theMap);
        theCubMap0. cache = false;
        theCubMapType theCubMap1 (offset, width, 1 << (finalDepth + 1),
                finalDepth + 1, theMap);
        theCubMap1. cache = false;
 
        // prepare a boolean value to store the information on whether
        // the Conley indices in the Morse decomposition were computed
        // or read from cache
        bool morseDecComputed (false);
 
        // prepare an array in which wrong Conley indices will be recorded
        std::vector<int> wrongIndices;
 
        // prepare an array in which skipped Conley indices will be recorded
        std::vector<int> skippedIndices;
 
        // prepare an array in which attractors will be recorded
        std::vector<int> attractors;
 
        // compute the Morse decomposition
        theMorseDecompositionType *morseDec = computeMorseDecomposition
                (theMap, theCubMap0, theCubMap1, offset, width, skipIndices,
                shareName ? shareName : "",
                phaseSpaceName ? phaseSpaceName : "",
                cubesFilePrefix ? cubesFilePrefix : "",
                morseDecName ? morseDecName : "",
                graphName ? graphName : "",
                procFilePrefix ? procFilePrefix : "",
                mapOptFileLocal ? mapOptFileLocal : "",
                morseDecComputed, wrongIndices,
                skippedIndices, attractors, connOrbits,
                leftMapParam, rightMapParam, paramCount);
 
        // compute the min and max coordinates of Morse sets
        if (!sbug. show)
        {
                CoordMinMax coordMinMax;
                int_t nMorseSets = morseDec -> count ();
                for (int_t nSet = 0; nSet < nMorseSets; ++ nSet)
                        coordMinMax ((*morseDec) [nSet]);
                if (coordMinMax. defined ())
                {
                        sout << "The computed Morse sets are in " <<
                                coordMinMax << ".\n";
                        sout << "Real coordinates: ";
                        showRealCoords (sout, coordMinMax) << ".\n";
                }
        }
 
        // show computed max image diameter and volume
        sout << "Max img diam = " << theCubMapType::maxImgDiam <<
                ", max img vol = " << theCubMapType::maxImgVol << ".\n";
 
        // save a picture of the computed Morse decomposition
        if (
#if (!defined (PLOT_X_COORD) && !defined (PLOT_Y_COORD))
                (spaceDim == 2) &&
#endif
                phaseSpaceName && *phaseSpaceName)
        {
                sout << "Saving a picture of Morse sets to '" <<
                        phaseSpaceName << "'...\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, phaseSpaceName,
                        xCoord, yCoord,
                        fullPhaseSpace ? (1 << finalDepth) : 0);
        }
 
        // save the computed Conley-Morse graph for the "dot" program
        // NOTE: This procedure is contained in "computeMorseDecomposition".
        if (false && (graphName && *graphName))
        {
                // create the Morse graph from the Morse decomposition
                chomp::homology::diGraph<> connGraph;
                morseDec -> makegraph (connGraph);
                chomp::homology::diGraph<> morseGraph;
                transitiveReduction (connGraph, morseGraph);
 
                // gather the sizes of each Morse set into a vector
                int nSets = morseGraph. countVertices ();
                std::vector<int_t> sizes (nSets);
                for (int n = 0; n < nSets; ++ n)
                        sizes [n] = (*morseDec) [n]. size ();
 
                // gather the Conley indices and the eigenvalues
                std::vector<theConleyIndexType> indices (nSets);
                std::vector<IndexEigenValues> eigenValues (nSets);
                for (int n = 0; n < nSets; ++ n)
                {
                        indices [n] = morseDec -> index (n);
                        eigenValues [n] = IndexEigenValues (indices [n]);
                }
 
                sout << "Saving the Conley-Morse graph to '" <<
                        graphName << "'...\n";
                std::ofstream graphFile (graphName);
                writeDotGraph (graphFile, morseGraph, sizes,
                        indices, eigenValues,
                        wrongIndices, skippedIndices, attractors);
                graphFile << "\n";
                graphFile. close ();
        }
 
        // save the computed Morse sets to separate text files
        // NOTE: This procedure is contained in "computeMorseDecomposition".
//      if (cubesFilePrefix && *cubesFilePrefix)
//      {
//              sout << "Saving the Morse decomposition to '" <<
//                      cubesFilePrefix << "*'...\n";
//              morseDec -> savetofiles (cubesFilePrefix);
//      }
 
        // delete the Morse deocmposition
        delete morseDec;
 
        // save new map optimization data if any
        if (mapOptFileName && *mapOptFileName)
        {
                std::string line = theMap. getOptInfo ();
                if (!line. empty ())
                {
                        mapOptimization. addData (line, true);
                        std::ofstream mapOptFile;
                        mapOptFile. open (mapOptFileName,
                                std::ios::out | std::ios::app);
                        mapOptimization. saveData (mapOptFile, true);
                }
        }
 
        // release pointset data to make fewer valgrind complaints
        sbug << "(Releasing the pointset data.)\n";
        parCube::PointBase::reset ();
        spcCube::PointBase::reset ();
        chomp::homology::Cube::PointBase::reset ();
 
        return 0;
} /* runSingleComp */
 
 
// --------------------------------------------------
// ---------------------- MAIN ----------------------
// --------------------------------------------------
 
/// The main procedure of the program.
/// Returns: 0 = Ok, -1 = Error, 1 = Help displayed, 2 = Wrong arguments.
int main (int argc, char *argv [])
{
        using namespace chomp::homology;
 
        // prepare user-configurable data
        int skipIndices = 0;
        char *boxCoordinates = 0;
        char *mapOptFileName = 0;
        char *mapOptFileLocal = 0;
        char *graphName = 0;
        char *shareName = 0;
        char *phaseSpaceName = 0;
        char *cubesFilePrefix = 0;
        char *morseDecName = 0;
        char *procFilePrefix = 0;
        bool fullPhaseSpace = false;
        bool connOrbits = false;
        const int maxComments = 16;
        int nComments = 0;
        char *comments [maxComments];
 
        // interprete the command-line arguments
        arguments a;
        arg (a, "b", boxCoordinates);
        arg (a, "t", mapOptFileName);
        arg (a, "u", mapOptFileLocal);
        arg (a, "g", graphName);
        arg (a, "s", shareName);
        arg (a, "p", phaseSpaceName);
        arg (a, "z", morseDecName);
        arg (a, "c", cubesFilePrefix);//deprecated
        arg (a, "q", cubesFilePrefix);
        arg (a, "o", procFilePrefix);
        arg (a, "-comment", comments, nComments, maxComments);
        arg (a, "-skip-indices", skipIndices, 1);
        argswitch (a, "-full", fullPhaseSpace, true);
        argswitch (a, "-conn", connOrbits, true);
        arghelp (a);
 
        argstreamprepare (a);
        int argresult = a. analyze (argc, argv);
        argstreamset ();
 
        // if no box coordinates defined then show the help message
        if (!boxCoordinates)
                argresult = 1;
 
        // show the program's main title
        if (argresult >= 0)
                sout << title << '\n';
 
        // if something was incorrect, show an additional message and exit
        if (argresult < 0)
        {
                sout << "Call with '--help' for help.\n";
                return 2;
        }
 
        // if help requested, show help information
        if (argresult > 0)
        {
                sout << helpinfo << '\n';
                return 1;
        }
 
        // try running the main function and catch an error message if thrown
        try
        {
                // set an appropriate program time message
                program_time = "Aborted after:";
                program_time = 1;
 
                // run the computations for a single parameter box
                runSingleComp (boxCoordinates,
                        mapOptFileName, mapOptFileLocal,
                        graphName, shareName,
                        phaseSpaceName, fullPhaseSpace, cubesFilePrefix,
                        morseDecName, procFilePrefix,
                        connOrbits, skipIndices);
 
                // set an appropriate program time message
                program_time = "Total time used:";
 
                // finalize
                return 0;
        }
        catch (const char *msg)
        {
                sout << "ERROR: " << msg << '\n';
                return -1;
        }
        catch (const std::exception &e)
        {
                sout << "ERROR: " << e. what () << '\n';
                return -1;
        }
        catch (...)
        {
                sout << "ABORT: An unknown error occurred.\n";
                return -1;
        }
} /* main */