compmap.cpp

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file compmap.cpp
///
/// Computing the map on a given set of cubes at the final subdivision depth.
/// This program may be useful for the computation of the index pairs
/// and the index map for combinatorial Morse sets saved previously.
///
/// @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: March 19, 2013.
 
 
// include some standard C++ header files
#include <exception>
#include <sstream>
#include <iostream>
#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 "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 combinatorial cubical multivalued map.\n\
Ver. 0.01, March 26, 2013. Copyright (C) 2005-2012 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 combinatorial\n\
cubical multivalued map on a cubical set at the final subdivision depth\n\
for a single parameter box. It is a supplementary program to 'cmgraphs',\n\
and may be used for the computation of the index pairs and the index map\n\
for combinatorial Morse sets saved previously.\n\
Command line arguments (the argument '-b' is mandatory):\n\
-b x1,x2,...,xN - the integer coordinates of a parameter box to use,\n\
file.cub - a file with the set of cubes on which to compute the map,\n\
file.map - a file to which the computed map will be saved,\n\
-d disjoint.cub - a set of cubes that will be removed from the images,\n\
-t filename - file name with optimization information (read only),\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/.";
 
 
// --------------------------------------------------
// ---------------- Map Computation -----------------
// --------------------------------------------------
 
int runMapComp (const char *boxCoordinates, const char *mapOptFileName,
        const char *cubName, const char *mapName, const char *disjointName)
{
        using chomp::homology::sout;
        using chomp::homology::scon;
        using chomp::homology::sbug;
        using chomp::homology::ignorecomments;
        using chomp::homology::ignoreline;
        using chomp::homology::fileerror;
 
        // 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;
 
        // 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
        sout << "Parameter box: " << paramBox << ".\n";
        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;
        if (mapOptFileName && *mapOptFileName)
        {
                std::ifstream in (mapOptFileName);
                if (in)
                {
                        sout << "Reading map optimization data... ";
                        mapOptimization. loadData (in);
                        sout << " Done.\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 finest subdivision depth to be used in the computation
        theCubMapType theCubMap (offset, width, 1 << finalDepth,
                finalDepth, theMap);
        theCubMap. cache = false;
        theCubMap. cropping = false;
 
        // read the set of cubes with which the images will be made disjoint
        spcCubes disjoint;
        if (disjointName && *disjointName)
        {
                sout << "Reading '" << disjointName << "'... ";
                std::ifstream disjointFile (disjointName);
                if (!disjointFile)
                        fileerror (disjointName, "open");
                disjointFile >> disjoint;
                sout << disjoint. size () << " cubes read.\n";
        }
        bool usingDisjoint = !disjoint. empty ();
 
        // open the file containing the set of cubes
        std::ifstream cubFile (cubName);
        if (!cubFile)
                fileerror (cubName, "open");
 
        // ignore the line "dimension N" if it is there
        ignorecomments (cubFile);
        if (cubFile. peek () == 'd')
        {
                ignoreline (cubFile);
                ignorecomments (cubFile);
        }
 
        // create the file for saving the map
        std::ofstream mapFile (mapName);
        if (!mapFile)
                fileerror (mapName, "create");
 
        // compute the map on all the cubes and save it to the file
        sout << "Computing the map... ";
        int_t allImages = 0;
        int_t removedCubes = 0;
        int_t emptyImages = 0;
        int_t cur = 0;
        while (1)
        {
                // read a subsequent cube from the file unless EOF
                bool cubeRead = false;
                spcCube q;
                if (!cubFile. eof ())
                {
                        cubFile >> q;
                        ignorecomments (cubFile);
                        cubeRead = true;
                        ++ cur;
                }
 
                // prepare data for a computed image of a cube
                spcCubes img;
                bool imgComputed = false;
 
                // compute the image of the cube
                if (cubeRead)
                {
                        theCubMap (q, &img, 0, 0, 0, false);
                        imgComputed = true;
                }
 
                // if the image of a cube has been computed
                if (imgComputed)
                {
                        // save the map assignment to the output file
                        mapFile << q << " -> {";
                        int_t imgSize = img. size ();
                        allImages += imgSize;
                        bool first = true;
                        for (int_t i = 0; i < imgSize; ++ i)
                        {
                                // if the cube is in the set with which
                                // the image should be disjoint
                                if (usingDisjoint &&
                                        disjoint. check (img [i]))
                                {
                                        ++ removedCubes;
                                        -- allImages;
                                        continue;
                                }
 
                                // write the cube; separate cubes with spaces
                                if (first)
                                        first = false;
                                else
                                        mapFile << ' ';
                                mapFile << img [i];
                        }
                        mapFile << "}\n";
 
                        // if there were no cubes written, take a note of it
                        if (first)
                                ++ emptyImages;
 
                        // show progress indicator if necessary
                        if (cur && !(cur % 1000))
                        {
                                scon << std::setw (8) << (cur / 1000) <<
                                        "k\b\b\b\b\b\b\b\b\b";
                        }
                }
 
                // if there are no more cubes that must be read
                // and no more images that must be computed then quit
                if (cubFile. eof ())
                        break;
        }
 
        // show the number of processed cubes
        sout << "Images of " << cur << " cubes computed.\n";
        sout << "There are " << allImages << " cubes in total "
                "in all the images.\n";
        if (usingDisjoint)
        {
                sout << removedCubes << " cubes were removed "
                        "from the images.\n";
        }
        if (emptyImages)
        {
                sout << "WARNING: There are " << emptyImages <<
                        " empty images.\n";
        }
        else
        {
                sout << "All the images are nonempty.\n";
        }
 
        // close the files
        cubFile. close ();
        mapFile. close ();
 
        // show the average image size
        if (cur > 0)
                sout << "Average image size: " << (allImages / cur) << ".\n";
 
        // show computed max image diameter and volume
        sout << "Max img diam = " << theCubMapType::maxImgDiam <<
                ", max img vol = " << theCubMapType::maxImgVol << ".\n";
 
        // 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;
} /* runMapComp */
 
 
// --------------------------------------------------
// ---------------------- 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
        char *boxCoordinates = 0;
        char *cubName = 0;
        char *mapName = 0;
        char *disjointName = 0;
        char *mapOptFileName = 0;
        const int maxComments = 16;
        int nComments = 0;
        char *comments [maxComments];
 
        // interprete the command-line arguments
        arguments a;
        arg (a, NULL, cubName);
        arg (a, NULL, mapName);
        arg (a, "b", boxCoordinates);
        arg (a, "t", mapOptFileName);
        arg (a, "d", disjointName);
        arg (a, "-comment", comments, nComments, maxComments);
        arghelp (a);
 
        argstreamprepare (a);
        int argresult = a. analyze (argc, argv);
        argstreamset ();
 
        // if no box coordinates defined then show the help message
        if (!boxCoordinates || !cubName || !mapName)
                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
                runMapComp (boxCoordinates, mapOptFileName,
                        cubName, mapName, disjointName);
 
                // 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 */