/////////////////////////////////////////////////////////////////////////////
/// @file indpdemo.cpp
///
/// @author Pawel Pilarczyk
/////////////////////////////////////////////////////////////////////////////

// Copyright (C) 2006-2007 by Pawel Pilarczyk.
//
// This program is free software. It is provided herein and it can be
// redistributed and/or modified 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 program 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.

// Started in May 2006. Last revision: June 17, 2007.


// include relevant parts of the CHomP library
#include "chomp/chomp.h"
#include "chomp/homology/indxpalg.h"
using namespace chomp::homology;

// include relevant parts of the CAPD library
#include "krak/krak.h"
#include "interval/intervalOld.h"
#include "vectalg/vectalglib.h"
#include "dynsys/dynsyslib.h"
#include "dynset/dynsetlib.h"
using namespace capd::dynset;

// include standard libraries
#include <exception>
#include <cstdlib>
#include <new>
#include <iostream>
#include <fstream>
using namespace std;


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

const char *title = "\
INDPDEMO, ver. 0.01. Copyright (C) 2006-2007 by Pawel Pilarczyk.\n\
This is free software. No warranty. Consult 'license.txt' for details.";

const char *helpinfo = "\
Call with: S.cub Q1.cub Q2.cub F.map.\n\
This is a demonstration program which computes an index pair for some\n\
example maps and saves the result to the given files as follows:\n\
Q1 contains P1\\P2, Q2 contains P2, and F contains the entire map\n\
computed while running the program (its domain contains P1).\n\
The set S is an initial guess for an isolated invariant set.\n\
Optional arguments:\n\
-P - use Pilarczyk's algorithm - the entire set S will be contained in Q1,\n\
-M - use Mrozek's definition of an index pair - o(S) must contain some\n\
isolated invariant set, or the constructed index pair may be empty,\n\
Note: The switches -P and -M are complementary, and -P is the default.\n\
-H - use the Henon map,\n\
-V2 - use the time-t map of the Vanderpol flow in the plane,\n\
-R2 - use the reversed time-t map of the Van der Pol flow in the plane,\n\
-V3, -R3 - as above, but for the vector field embedded in R^3,\n\
Note: The switches -H, -V?, -R? are complementary, and -H is the default.\n\
For more information ask the author at http://www.PawelPilarczyk.com/.";


// --------------------------------------------------
// ----------------- PRELIMINARIES ------------------
// --------------------------------------------------

// dynamical systems to chose from
const int HenonMap = 0;
const int VanderpolMap = 1;
const int VanderpolMap3d = 2;
const int ReversedVanderpol = 3;
const int ReversedVanderpol3d = 4;

// the names of the dynamical systems
const char *sysnames [] = {"Henon", "Vanderpol", "Vanderpol 3D",
"Reversed Vanderpol", "Reversed Vanderpol 3D"};

// algorithms to chose from
const int PilarczykPair = 0;
const int MrozekPair = 1;

// the names of the algorithms
const char *algnames [] = {"Pilarczyk's algorithm",
        "Mrozek's index pair definition"};


// --------------------------------------------------
// ------------------- Henon Map --------------------
// --------------------------------------------------

/// The grid size: 1/64.
const int HenonGrid = 64;

/// Computes the image of a cube according to the Henon map
/// h (x, y) = (1 + y/5 - ax^2, 5bx) with a = 1.4 and b = 0.2.
int henon (const Cube::CoordType *coord, double *left, double *right)
{
        const double gridsize = HenonGrid;

        // version with rigorous numerics (interval arithmetic)
        static interval a = interval (7, 7) / 5;
        static interval b = interval (1, 1) / 5;
        interval x = interval (coord [0], coord [0] + 1) / gridsize;
        interval y = interval (coord [1], coord [1] + 1) / gridsize;
        interval result_x = 1 + y / 5 - a * power (x, 2);
        interval result_y = 5 * b * x;
        result_x *= gridsize;
        result_y *= gridsize;
        left [0] = result_x. left_b ();
        right [0] = result_x. right_b ();
        left [1] = result_y. left_b ();
        right [1] = result_y. right_b ();

/*      // version with approximate numerics
        int x2a = coord [0] * coord [0];
        int x2b = (coord [0] + 1) * (coord [0] + 1);
        int x2min = (x2a < x2b) ? x2a : x2b;
        int x2max = (x2a > x2b) ? x2a : x2b;
        left [0] = (gridsize * 5 + coord [1] - 7 * x2max / gridsize) / 5.0;
        right [0] = (gridsize * 5 + (coord [1] + 1) - 7 * x2min / gridsize) / 5.0;
        left [1] = coord [0];
        right [1] = coord [0] + 1;
*/
        return 0;
} /* henon */


// --------------------------------------------------
// ------------- Time-t translation map -------------
// --------------------------------------------------

/// The class that computes the translation by time t in an ODE.
/// NOTE: Both grid and its inverse must be representable!
class Translation
{
public:
        /// The only construction allowed.
        Translation (int _dim, double _grid,
                const char *_mapstring, int _order,
                double _step, int _steps);

        /// Computes the image of a cube given its coordinates.
        int operator () (const Cube::CoordType *coord,
                double *left, double *right);

private:
        /// The dimension of the space.
        int dim;

        /// The inverse of the grid size.
        double invgrid;

        /// A half of the grid size.
        double halfgrid;
       
        /// The dynamical system for the translation by one step.
        Taylor dynam;

        /// The number of steps to apply.
        int steps;

        /// The Lohner cube centered at the middle of each cube.
        intervalVector lohnerCube;

}; /* class Translation */

inline Translation::Translation (int _dim, double _grid,
        const char *_mapstring, int _order, double _step, int _steps):
        dim (_dim), invgrid (1.0 / _grid), halfgrid (_grid / 2),
        dynam (_mapstring, _order, interval (_step, _step)),
        steps (_steps), lohnerCube ()
{
        // create an appropriate Lohner cube
        interval edge = interval (-halfgrid, halfgrid);
        for (int i = 0; i < dim; ++ i)
                lohnerCube [i] = edge;

        return;
} /* Translation::Translation */

inline int Translation::operator () (const Cube::CoordType *coord,
        double *left, double *right)
{
        // compute the middle of this cube in R^n
        intervalVector center;
        for (int i = 0; i < dim; ++ i)
                center [i] = (coord [i] + coord [i] + 1) * halfgrid;

        // create a Lohner set that represents this cube
        rect2Set s (center, lohnerCube);

        // iterate the Lohner set as long as necessary
        for (int j = 0; j < steps; ++ j)
                s. move (dynam);

        // cover the Lohner set with a box (an interval vector)
        intervalVector cover (s);

        // fill in the return arrays
        for (int i = 0; i < dim; ++ i)
        {
                cover [i] *= invgrid;
                left [i] = cover [i]. left_b ();
                right [i] = cover [i]. right_b ();
        }

        return 0;
} /* Translation::operator () */


// --------------------------------------------------
// ----------------- Vanderpol Map ------------------
// --------------------------------------------------

/// The inverse of the grid size.
const int VanderpolGrid = 32;

/// The order of the numerical method for integrating the ODE.
const int VanderpolOrder = 3;

/// The step size for the iterations of the ODE.
const double VpolStep = 1.0 / 128;

/// The number of steps for Vanderpol.
const int VpolSteps = 71;

/// The number of steps for reversed Vanderpol.
const int RevSteps = 35;

/// Computes the image of a cube in the translation by time 1
/// in the Vanderpol flow in the plane.
int vanderpol (const Cube::CoordType *coord, double *left, double *right)
{
        /// Abort if the dimension is different from the compiled one.
        if (DIM != 2)
                throw "The CAPD library must be compiled with DIM=2.";

        /// Translation in the Vanderpol vector field.
        const char *vfield = "var:x,y;fun:-y+(x*((x*x)-1)),x;";
        static Translation vpolfield (2, 1.0 / VanderpolGrid,
                vfield, VanderpolOrder, VpolStep, VpolSteps);

        return vpolfield (coord, left, right);
} /* vanderpol */

/// Computes the image of a cube in the translation by time -1
/// in the Vanderpol flow in the plane.
int reversedvpol (const Cube::CoordType *coord, double *left, double *right)
{
        /// Abort if the dimension is different from the compiled one.
        if (DIM != 2)
                throw "The CAPD library must be compiled with DIM=2.";

        /// The reversed Vanderpol vector field.
        const char *vfield = "var:x,y;fun:y-(x*((x*x)-1)),-x;";
        static Translation revvpolfield (2, 1.0 / VanderpolGrid,
                vfield, VanderpolOrder, VpolStep, RevSteps);

        return revvpolfield (coord, left, right);
} /* reversedvpol */

/// Computes the image of a cube in the translation by time 1
/// in the Vanderpol flow embedded in the space R^3.
int vanderpol3d (const Cube::CoordType *coord, double *left, double *right)
{
        /// Abort if the dimension is different from the compiled one.
        if (DIM != 3)
                throw "The CAPD library must be compiled with DIM=3.";

        /// Translation in the Vanderpol vector field embedded in R^3.
        ostringstream v;
        v << "var:x,y,z;fun:y-(x*((x*x)-1)),-x," <<
                (0.5 / VanderpolGrid) << "-8*(z-" <<
                (0.5 / VanderpolGrid) << ");";
        const char *vfield = v. str (). c_str ();
        static Translation vpolfield3d (3, 1.0 / VanderpolGrid, vfield,
                VanderpolOrder, VpolStep, VpolSteps);

        return vpolfield3d (coord, left, right);
} /* vanderpol3d */

/// Computes the image of a cube in the translation by time -1
/// in the Vanderpol flow embedded in the space R^3.
int reversedvpol3d (const Cube::CoordType *coord, double *left, double *right)
{
        /// Abort if the dimension is different from the compiled one.
        if (DIM != 3)
                throw "The CAPD library must be compiled with DIM=3.";

        /// The reversed Vanderpol vector field embedded in R^3.
        ostringstream v;
        v << "var:x,y,z;fun:-y+(x*((x*x)-1)),x," <<
                (0.5 / VanderpolGrid) << "-8*(z-" <<
                (0.5 / VanderpolGrid) << ");";
        const char *vfield = v. str (). c_str ();
        static Translation revvpolfield3d (3, 1.0 / VanderpolGrid,
                vfield, VanderpolOrder, VpolStep, RevSteps);

        return revvpolfield3d (coord, left, right);
} /* reversedvpol3d */


// --------------------------------------------------
// ------------------- index pair -------------------
// --------------------------------------------------

template <class object>
int savetofile (const char *filename, const object &obj)
{
        if (!filename || !*filename)
                return -1;
        ofstream out (filename);
        out << obj;
        return 0;
} /* savetofile */

template <class object>
int readfromfile (const char *filename, object &obj)
{
        if (!filename || !*filename)
                return -1;
        ifstream in (filename);
        in >> obj;
        return 0;
} /* readfromfile */

void indpdemo (const char *Sname, const char *Q1name, const char *Q2name,
        const char *Fname, int whichSystem, int whichAlgorithm)
{
        // say what you are doing
        sout << "Using " << algnames [whichAlgorithm] << " for " <<
                sysnames [whichSystem] << ".\n";

        // create an initial approximation of S
        SetOfCubes S;
        sout << "Reading S from '" << Sname << "'... ";
        readfromfile (Sname, S);
        sout << static_cast<int> (S) << " cubes.\n";

        // prepare the sets for the index pair
        SetOfCubes Q1, Q2;

        // choose the right map
        BufferedMapClass<>::map f [] = {henon, vanderpol,
                vanderpol3d, reversedvpol, reversedvpol3d};

        // create a map object
        BufferedMapClass<> F (f [whichSystem]);

        // run the index pair construction algorithm
        switch (whichAlgorithm)
        {
        case MrozekPair:
                IndexPairM (F, S, Q1, Q2);
                break;

        case PilarczykPair:
                IndexPairP (F, S, Q1, Q2);
                break;
        }

        // indicate the sizes of the sets
        sout << "There are " << static_cast<int> (Q1) <<
                " cubes in Q1 and " << static_cast<int> (Q2) <<
                " cubes in Q2.\n";
        sout << "The map was computed on " <<
                static_cast<int> (F. F. getdomain ()) << " cubes.\n";

        // save the sets and the map
        savetofile (Q1name, Q1);
        savetofile (Q2name, Q2);
        savetofile (Fname, F);
        return;
} /* HenonPair */


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

/// The main function that runs the program.
/// Returns: 0 = Ok, -1 = Error, 1 = Help displayed, 2 = Wrong arguments.
int main (int argc, char *argv [])
{
        // prepare user-configurable data
        char *Sname = 0, *Q1name = 0, *Q2name = 0, *Fname = 0;
        int whichSystem = HenonMap;
        int whichAlgorithm = PilarczykPair;

        // interprete the command-line arguments
        arguments a;
        arg (a, 0, Sname);
        arg (a, 0, Q1name);
        arg (a, 0, Q2name);
        arg (a, 0, Fname);
        argswitch (a, "M", whichAlgorithm, MrozekPair);
        argswitch (a, "P", whichAlgorithm, PilarczykPair);
        argswitch (a, "H", whichSystem, HenonMap);
        argswitch (a, "V2", whichSystem, VanderpolMap);
        argswitch (a, "V3", whichSystem, VanderpolMap3d);
        argswitch (a, "R2", whichSystem, ReversedVanderpol);
        argswitch (a, "R3", whichSystem, ReversedVanderpol3d);
        argswitch (a, "V", whichSystem, VanderpolMap);
        argswitch (a, "R", whichSystem, ReversedVanderpol);
        arghelp (a);

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

        // show the program's main title
        if (argresult >= 0)
                sout << title << '\n';

        // if there are not enough file names, help should be displayed
        if (!Fname)
                argresult = 1;

        // 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 main procedure
                indpdemo (Sname, Q1name, Q2name, Fname,
                        whichSystem, whichAlgorithm);

                // 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 << "EXCEPTION: " << e. what () << '\n';
                return -1;
        }
        catch (...)
        {
                sout << "ABORT: An unknown error occurred.\n";
                return -1;
        }
} /* main */