runhenon.cpp

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file runhenon.cpp
///
/// The main computations for the Henon attractor.
/// This program runs rigorous numerical computations aimed at proving
/// that an attractor for the given dynamical system is mixing
/// at some finite resolution.
/// This program accompanies the paper "Finite resolution dynamics"
/// by S. Luzzatto and P. Pilarczyk, and follows the algorithms
/// introduced there.
///
/// @author Pawel Pilarczyk
///
/////////////////////////////////////////////////////////////////////////////
///
/// @mainpage The Finite Resolution Dynamics Software
///
/// \section intro Introduction
///
/// This website contains the documentation of software designed and
/// programmed by Pawel Pilarczyk which accompanies the paper
/// <em>Finite resolution dynamics</em> by S. Luzzatto and P. Pilarczyk.
///
/// \section License
///
/// In order to make this software freely available for wide audience,
/// it is provided here under the terms of the GNU General Public License,
/// version 2 or any later version. Please, see the
/// <a href="http://www.gnu.org/copyleft/gpl.html">GNU website</a>
/// for details.
///
/// \section browsing Browsing
///
/// Although a lot of effort was put into making this software as transparent
/// and clear as possible, and to keep this documentation complete,
/// I am sure that there is a lot of important information
/// that might have been included here but was not.
/// Therefore, this documentation was generated in such a way
/// that it contains the source code which can be browsed and checked
/// in case of any doubt.
///
/// So, what are you waiting for? Dive and read. Have a happy browsing!
///
/// Pawel Pilarczyk
///
/////////////////////////////////////////////////////////////////////////////

// Copyright (C) 1997-2010 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 December 22, 2008. Last revision: April 12, 2010.


// BOOST (this header must be included first)
#include "boost/numeric/interval.hpp"

// include some standard C++ header files
#include <iomanip>

// include some headers 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 "maphenon.h"
#include "covboxes.h"
#include "attractor.h"
#include "graph.h"
#include "streams.h"
#include "plotcover.h"
#include "savecover.h"


// --------------------------------------------------
// -------------------- OVERTURE --------------------
// --------------------------------------------------

/// The title of the program with brief licensing information.
const char *title = "\
Finite Resolution Dynamics - computations for the Henon attractor.\n\
Ver. 0.01, April 12, 2010. Copyright (C) 1997-2010 by Pawel Pilarczyk.\n\
This is free software. No warranty. Consult 'license.txt' for details.";

/// Short help information about the program's usage.
const char *helpinfo = "\
This program accompanies the paper \"Finite resolution dynamics\"\n\
by S. Luzzatto and P. Pilarczyk. It runs computations aiming at proving\n\
that the Henon attractor is mixing at all resolutions > some epsilon.\n\
Call this program with the following arguments:\n\
-p N - the number of cover parts in each direction (must be given),\n\
-b filename.bmp - a file to save a plot of the cover to (default: none),\n\
-w N - the width of the BMP file (default: 1000 pixels),\n\
-s filename.txt - a text file to save the actual cover to (def: none),\n\
-h - to display this brief help information only and exit.\n\
For more information consult the accompanying documentation (if available)\n\
or ask the program's author at http://www.PawelPilarczyk.com/.";


// --------------------------------------------------
// ---------------- HENON ATTRACTOR -----------------
// --------------------------------------------------

/// Runs the computations for the Henon attractor.
template <class NumType>
void henonAttractor (int parts, const char *bitmapFilename,
        int bitmapWidth, const char *saveFilename)
{
        using chomp::homology::timeused;

        sout << "--- Running computations for the Henon map. ---\n";
        sout << std::setprecision (12);

        // define the type of the map to be used for the computations
        typedef tMapHenon<NumType> theMapType;

        // define the type of the map to be used for the computations
        typedef tCoverBoxes<NumType> theCoverType;

        // set the dimension of the phase space
        const int dim = 2;

        // prepare the parameters of the Henon map
        NumType aNumerator = 14;
        NumType aDenominator = 10;
        NumType bNumerator = 3;
        NumType bDenominator = 10;
        sout << "Parameters of the Henon map: a = " <<
                aNumerator / aDenominator <<
                ", b = " << bNumerator / bDenominator << ".\n";

        // create a dynamical system generator
        tMapHenon<NumType> function (aNumerator, aDenominator,
                bNumerator, bDenominator);

        // prepare the phase space box to cover
        NumType leftCorner [dim] = {-1.4, -0.5}; //{-1.3, -0.4};
        NumType boxWidths [dim] = {2.8, 1.0}; //{2.6, 0.8};
        for (int i = 0; i < dim; ++ i)
        {
                sout << (i ? ") x (" : "Phase space box: (") <<
                        leftCorner [i] << "," <<
                        (leftCorner [i] + boxWidths [i]);
        }
        sout << ").\n";

        // prepare the number of elements in each direction
        int partCounts [dim];
        partCounts [0] = parts;
        for (int i = 1; i < dim; ++ i)
        {
                partCounts [i] = static_cast<int>
                        (parts * boxWidths [i] / boxWidths [0] + 0.5);
        }
        NumType *overlap = 0;
        double partsCount = 1;
        for (int i = 0; i < dim; ++ i)
        {
                sout << (i ? " x " : "Number of cover elements: ") <<
                        partCounts [i];
                partsCount *= partCounts [i];
        }
        sout << " = " << partsCount << ".\n";
        double area1 = 1;
        for (int i = 0; i < dim; ++ i)
        {
                double size = static_cast<double> (boxWidths [i]) /
                        partCounts [i];
                sout << (i ? " x " : "Size of each box: ") << size;
                area1 *= size;
        }
        sout << ".\n";
        sout << "Area of each box: " << area1 << ".\n";

        // create an initially empty cover
        tCoverBoxes<NumType> cover (dim, leftCorner, boxWidths,
                partCounts, overlap);

        // prepare an empty directed graph to store the map on the attractor
        Graph graph;

        // set up an initial point that is supposedly in the attraction
        // basin of the expected attractor
        NumType initialPoint [dim] = {0.61989426930989, 0.17586130934794};

        // set the number of iterations that need to be done
        // until the point gets close enough to the attractor
        int iterCount = 1000;

        // define the maximal allowed cover size
        int maxCoverSize = 0;

        // construct a cover of the attractor
        sout << "Computing a cover of the attractor...\n";
        timeused coverTime;
        double maxDiameter2 = constructAttractor (function, cover, graph,
                initialPoint, iterCount, maxCoverSize);
        double percentage = static_cast<double> (cover. size ()) * 100 /
                partsCount;
        sout << "A cover consisting of " << cover. size () <<
                " elements (" << percentage << "%) constructed.\n";
        sout << "Computed in " << coverTime << ".\n";
        double area = 1;
        for (int i = 0; i < dim; ++ i)
                area *= boxWidths [i] / partCounts [i];
        area *= cover. size ();
        sout << "Approximate area of the cover: " << area << ".\n";
        sout << "Max square of the diameter of map images: " <<
                maxDiameter2 << ".\n";
        sout << "Approx. max diameter of map images: " <<
                std::sqrt (maxDiameter2) << ".\n";
        sout << "The map graph has " << graph. countEdges () <<
                " edges.\n";

        // plot the resulting cover in a bitmap image file
        if (bitmapFilename)
        {
                sout << "Plotting the computed cover in a "
                        "black-and-white bitmap image...\n";
                plotCover (cover, bitmapFilename, bitmapWidth);
        }

        // save the computed cover to a text file
        if (saveFilename)
        {
                sout << "Saving the computed cover to a text file...\n";
                saveCover (cover, saveFilename);
        }

        // clean up the cover and cubical sets from the memory
        sout << "Releasing the cover from the memory...\n";
        timeused forgetTime;
        cover. forget ();
        chomp::homology::PointBase::forget ();
        sout << "Memory released in " << forgetTime << ".\n";

        // verify if the constructed graph is strongly connected
        sout << "Checking if the graph is strongly connected...\n";
        timeused strConnTime;
        if (graph. checkStronglyConnected ())
                sout << "Yes! The graph is strongly connected.\n";
        else
                throw "The constructed graph is not strongly connected.\n";
        sout << "Verified in " << strConnTime << ".\n";

        // verify if the constructed graph is aperiodic
        sout << "Checking if the graph is aperiodic...\n";
        timeused aperTime;
        if (graph. checkAperiodic ())
                sout << "Good! The graph is aperiodic.\n";
        else
                throw "The constructed graph is not aperiodic.\n";
        sout << "Verified in " << aperTime << ".\n";

        sout << "--- The computations completed successfully. ---\n";
        return;
} /* henonAttractor */


// --------------------------------------------------
// ---------------------- MAIN ----------------------
// --------------------------------------------------

/// The main function of the program.
/// Returns: 0 = Ok, -1 = Error, 1 = Help displayed, 2 = Wrong arguments.
int main (int argc, char *argv [])
{
        using namespace chomp::homology;

        // turn on a message that will appear if the program does not finish
        program_time = "Aborted after";

        // prepare user-configurable data
        int parts = 0;
        char *bitmapFilename = 0;
        int bitmapWidth = 1000;
        char *saveFilename = 0;

        // analyze the command line
        arguments a;
        arg (a, "p", parts);
        arg (a, "b", bitmapFilename);
        arg (a, "w", bitmapWidth);
        arg (a, "s", saveFilename);
        arghelp (a);

        argstreamprepare (a);
        int argresult = a. analyze (argc, argv);
        argstreamset ();

        // show the program's 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 or no part number defined, show help information
        if ((argresult > 0) || (parts <= 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:";

                // run the main procedure
                henonAttractor<double> (parts, bitmapFilename, bitmapWidth,
                        saveFilename);

                // set an appropriate program time message
                program_time = "Total time used:";
                program_time = 1;

                // 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;
        }
        return 0;
} /* main */