/////////////////////////////////////////////////////////////////////////////
/// @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 */