/*   File: poincare.cc 
 
     Computes Poincare maps for x' = f(x) using a distance-d method. 
     By keeping track of the two largest components of the vector field, 
     the program detects (and warns) if the stepsize used is too large.
     The desired stepsize is independent of the order of the numerical 
     method employed.

     Latest edit: Mon April 30, 2001

     Compilation: g++ -D__I386__ -Wall -Wno-deprecated -o poincare poincare.cc \
                      -I/PROFIL_DIR/include -IBIAS -ILR \
                      -L/PROFIL_DIR/lib -lProfil -llr -LBIAS -lBias
*/


//////////////////////////////////////////////////////////////////////////////
// Standard headers

#include <iostream.h>
#include <math.h>
#include <stdio.h>
#include <time.h>

//////////////////////////////////////////////////////////////////////////////
// IA headers

#include "Functions.h"
#include "Matrix.h"
#include "Vector.h"

//////////////////////////////////////////////////////////////////////////////
// Flags:

#define COMPUTE_C1

//#define EULER
//#define MPEULER
#define RUNGEKUTTA4

//#define ROTATION
#define VOLTERRA
//#define LORENZ
//#define ROSSLER

// Set the order of the integrator.
#ifdef EULER
const int ORDER = 1;
#endif
#ifdef MPEULER
const int ORDER = 2;
#endif
#ifdef RUNGEKUTTA4
const int ORDER = 4;
#endif

// Set the dimension of the autonomous vector field.
#if defined LORENZ || defined ROSSLER
const int PRE_DIM = 3;  
#endif
#if defined VOLTERRA || defined ROTATION 
const int PRE_DIM = 2;
#endif

// Set the parameters.
#ifdef LORENZ
const int INTERSECTIONS = 2;
#else
const int INTERSECTIONS = 1;
#endif

#ifdef VOLTERRA
const double A = 0.4;
const double B = 0.1;
const double C = 0.3;
const double D = 0.5;
#endif

#ifdef LORENZ
const double S = 10.0;
const double R = 250.0;
const double B = 8.0/3;
#endif

#ifdef ROSSLER
const double A = 2.2;
//const double A = 5.7;
#endif

//////////////////////////////////////////////////////////////////////////////
// Constants etc.

const int    DIM       = PRE_DIM + 1;// The dimension of the embedded system.
const int    MAX_ITER  = 10000000;   // Max number of iterations.
const double MACH_EPS  = 1.0e-6;     // A very small numer.
const double DELTA_MIN = Power(2, -8);
const double DELTA_MAX = Power(2, -8);// A rough estimate of the global error.
//const double MIN_STEP  = Power(DELTA_MIN, 1.0 / ORDER); // Minimal step size.
//const double MAX_STEP  = Power(DELTA_MAX, 1.0 / ORDER); // Maximal step size.
const double MIN_STEP  = DELTA_MIN;  // Minimal step size.
const double MAX_STEP  = DELTA_MAX;  // Maximal step size.

enum CLOCK_ORDER {START, STOP};

typedef struct
{
  VECTOR firstX;                  
  int    firstMaxI;            
  int    firstSign;
} POINCARE;

//////////////////////////////////////////////////////////////////////////////

int sign(const double &arg)
{ 
  return (arg >= 0 ? + 1 : - 1); 
}

//////////////////////////////////////////////////////////////////////////////

void id(MATRIX &M)
{ 
  int d = RowDimension(M);
  Clear(M);
  for (int i = 1; i <= d; i++)
    for (int j = 1; j <= d; j++)
      if ( i == j )
	M(i, j) = 1.0;
}

//////////////////////////////////////////////////////////////////////////////
// The vector field proper. Modify this one only.

void vectorField(VECTOR &result, const VECTOR &x)
{
#ifdef ROTATION
  result(1) = + x(2);
  result(2) = - x(1); 
#endif
#ifdef VOLTERRA
  result(1) = + A * x(1) - B * x(1) * x(2);
  result(2) = - C * x(2) + D * x(1) * x(2); 
#endif
#ifdef LORENZ
  result(1) = S * (x(2) - x(1));
  result(2) = R * x(1) - x(2) - x(1) * x(3); 
  result(3) = - B * x(3) + x(1) * x(2); 
#endif
#ifdef ROSSLER
  result(1) = - x(2) - x(3);
  result(2) = + x(1) + x(2) / 5; 
  result(3) = + 1.0 / 5 + x(3) * (x(1) - A); 
#endif
}

//////////////////////////////////////////////////////////////////////////////
// The partial derivatives of the vector field proper. Modify this one only.

#ifdef COMPUTE_C1
void partialVectorField(MATRIX &result, const VECTOR &x)
{
#ifdef ROTATION
  result(1, 1) =   0.0;
  result(1, 2) = + 1.0; 
  result(2, 1) = - 1.0;
  result(2, 2) =   0.0; 
#endif
#ifdef VOLTERRA  
  result(1, 1) = + A - B * x(2);
  result(1, 2) = - B * x(1); 
  result(2, 1) = + D * x(2);
  result(2, 2) = - C + D * x(1); 
#endif
#ifdef LORENZ
  result(1, 1) = - S;
  result(1, 2) = + S; 
  result(1, 3) =    0.0;
  result(2, 1) = + R - x(3);
  result(2, 2) = -  1.0; 
  result(2, 3) = - x(1);
  result(3, 1) = + x(2);
  result(3, 2) = + x(1); 
  result(3, 3) = - B; 
#endif
#ifdef ROSSLER
  result(1, 1) = + 0.0;
  result(1, 2) = - 1.0; 
  result(1, 3) = - 1.0;
  result(2, 1) = + 1.0;
  result(2, 2) = + 1.0 / 5; 
  result(2, 3) = + 0.0;
  result(3, 1) = + x(3);
  result(3, 2) = + 0.0; 
  result(3, 3) = + x(1) - A; 
#endif
}
#endif // COMPUTE_C1

//////////////////////////////////////////////////////////////////////////////
// The rescaled partial derivatives of the vector field. 

#ifdef COMPUTE_C1
void rescaledPartialVectorField(MATRIX &result, const VECTOR &x, 
				const MATRIX &m, const int &maxI)
{
  VECTOR f(PRE_DIM);
  vectorField(f, x);
  MATRIX df(PRE_DIM, PRE_DIM);
  partialVectorField(df, x); 
  for (int i = 1; i <= PRE_DIM; i++)
    for (int j = 1; j <= PRE_DIM; j++)
      result(i, j) = df(i, j)*f(maxI) - f(i)*df(maxI, j);

  result = result * m;        // DPhi'(x, t) = Dg(Phi(x, t)) * DPhi(x, t),  
  result /= Sqr(f(maxI));     // where g(x) = f(x) / f_{maxI}(x).
}
#endif // COMPUTE_C1

//////////////////////////////////////////////////////////////////////////////

double vectorFieldComponent(const VECTOR &x, const int &i)
{
  VECTOR result(PRE_DIM);
  vectorField(result, x);

  return result(i);
}

//////////////////////////////////////////////////////////////////////////////

void embeddedVectorField(VECTOR &result, const VECTOR &x)
{
  vectorField(result, x); // Set the dependent variables.
  result(DIM) = + 1.0;    // Set the time variable.
}

//////////////////////////////////////////////////////////////////////////////

void rescaledVectorField(VECTOR &result, const VECTOR &x, const int &maxI)
{
  VECTOR v(DIM);
  
  embeddedVectorField(v, x);
  result =  v / v(maxI);
}

//////////////////////////////////////////////////////////////////////////////

void getMaxIJ(int &maxI, int &maxJ, const VECTOR &x)
{
  VECTOR v(PRE_DIM);
  vectorField(v, x);

  double maxValue = fabs(v(1));
  maxI = 1;
  
  for (int i = 2; i <= PRE_DIM; i++)
    if ( maxValue < fabs(v(i)) ) {
      maxValue = fabs(v(i));
      maxI     = i;
    }

  if ( maxValue < MACH_EPS ) {
    cout << "Warning (getMaxI): "
	 << "The vector field may be too small." << endl
	 << "Do you want to continue? (y/n) ";
    char answer;
    cin >> answer; cout << endl;
    if ( answer == 'n' )
      exit(1);
  }

  // Now find maxJ.
  if (maxI != 1) {
    maxValue = fabs(v(1));
    maxJ     = 1;
  }
  else {
    maxValue = fabs(v(2));
    maxJ     = 2;
  }    
  for (int i = 2; i <= PRE_DIM; i++)
    if ( (i != maxI) && (maxValue < fabs(v(i))) ) {
      maxValue = fabs(v(i));
      maxJ     = i;
    }
}

//////////////////////////////////////////////////////////////////////////////

void Euler(VECTOR &vResult, MATRIX &DPsi, const VECTOR &x, 
	   const double &dx, const int &maxI)
{
  VECTOR s1(DIM);
  rescaledVectorField(s1, x, maxI);
  vResult = x + dx * s1;

#ifdef COMPUTE_C1
  MATRIX m1(PRE_DIM, PRE_DIM);
  MATRIX ID(PRE_DIM, PRE_DIM); id(ID);
  rescaledPartialVectorField(m1, x, ID, maxI);
  DPsi = ID + dx * m1;
#endif // COMPUTE_C1
}

//////////////////////////////////////////////////////////////////////////////

void MPE(VECTOR &vResult, MATRIX &DPsi, const VECTOR &x, 
	 const double &dx, const int &maxI)
{
  VECTOR s1(DIM);
  VECTOR s2(DIM);

  rescaledVectorField(s1, x, maxI);
  VECTOR p1 = x + dx / 2.0 * s1;
  rescaledVectorField(s2, p1, maxI);
  vResult = x + dx * s2;

#ifdef COMPUTE_C1
  MATRIX ds1(PRE_DIM, PRE_DIM);
  MATRIX ds2(PRE_DIM, PRE_DIM);
  MATRIX ID(PRE_DIM, PRE_DIM); id(ID);
  rescaledPartialVectorField(ds1, x, ID, maxI);
  MATRIX m1 = ID + dx / 2.0 * ds1;
  rescaledPartialVectorField(ds2, p1, m1, maxI);
  DPsi = ID + dx * ds2;
#endif // COMPUTE_C1
}

//////////////////////////////////////////////////////////////////////////////

void RK4(VECTOR &vResult, MATRIX &DPsi, const VECTOR &x, 
	 const double &dx, const int &maxI)
{
  VECTOR s1(DIM);
  VECTOR s2(DIM);
  VECTOR s3(DIM);
  VECTOR s4(DIM);

  rescaledVectorField(s1, x, maxI);
  VECTOR p1 = x + dx / 2.0 * s1;
  rescaledVectorField(s2, p1, maxI);
  VECTOR p2 = x + dx / 2.0 * s2;
  rescaledVectorField(s3, p2, maxI);
  VECTOR p3 = x + dx * s3;
  rescaledVectorField(s4, p3, maxI);
  vResult = x + dx * (s1 + 2 * s2 + 2 * s3 + s4) / 6.0;

#ifdef COMPUTE_C1
  MATRIX ds1(PRE_DIM, PRE_DIM);
  MATRIX ds2(PRE_DIM, PRE_DIM);
  MATRIX ds3(PRE_DIM, PRE_DIM);
  MATRIX ds4(PRE_DIM, PRE_DIM);
  MATRIX ID(PRE_DIM, PRE_DIM); id(ID);
  rescaledPartialVectorField(ds1, x, ID, maxI);
  MATRIX m1 = ID + dx / 2.0 * ds1;
  rescaledPartialVectorField(ds2, p1, m1, maxI);
  MATRIX m2 = ID + dx / 2.0 * ds2;
  rescaledPartialVectorField(ds3, p2, m2, maxI);
  MATRIX m3 = ID + dx * ds3;
  rescaledPartialVectorField(ds4, p3, m3, maxI);
  DPsi = ID + dx * (ds1 + 2 * ds2 + 2 * ds3 + ds4) / 6.0; 
#endif // COMPUTE_C1
}

//////////////////////////////////////////////////////////////////////////////

void computeDPi(MATRIX &DPi, const MATRIX &DPsi, const VECTOR &nextX,
		const int &nextMaxI)
{
  VECTOR rvf(DIM);
  rescaledVectorField(rvf, nextX, nextMaxI);
  DPi = DPsi;
  for (int i = 1; i <= PRE_DIM; i++)
    for (int j = 1; j <= PRE_DIM; j++)
      DPi(i, j) -= rvf(i) * DPsi(nextMaxI, j);
}

//////////////////////////////////////////////////////////////////////////////

int checkStepSize(const VECTOR &x)
{
  static bool firstTime = true;
  static int  oldMaxI, oldMaxJ; 
  static int  oldSign;
  int  maxI, maxJ;
  getMaxIJ(maxI, maxJ, x);
  int  thisSign = sign(vectorFieldComponent(x, maxI));
  bool warning  = false;

  if (firstTime)
    firstTime = false;
  else {                       // Sanity check for the step size.        
    if ( maxI != oldMaxI ) {   // Enforce (oldMaxI, oldMaxJ) -> (maxJ, maxI)  
      if ( (maxI != oldMaxJ) || (maxJ != oldMaxI) )
	warning = true;        // when changing dominating component. 
    }
    else if ( oldSign == - thisSign )
      warning = true;          // Warn when oldMaxI = maxI, but sign changes. 
    if ( warning ) {
      cout << endl << "Warning (computeNextStep): "
	   << "You may be using a too large step size." << endl
	   << "(" << oldMaxI << ", " << oldMaxJ << ") -> "
	   << "(" << maxI << ", " << maxJ << ")." << endl
	   << "Do you want to continue? (y/n) ";
      char answer;
      cin >> answer; cout << endl;
      if ( answer == 'n' )
	exit(1);
    }
  }
  oldMaxI = maxI;
  oldMaxJ = maxJ;
  oldSign = thisSign;
  
  return maxI;
}

//////////////////////////////////////////////////////////////////////////////

void computeNextStep(VECTOR &vResult, const VECTOR &x, 
		     MATRIX &mResult, const double &dx)
{
  int maxI = checkStepSize(x);

#ifdef EULER
  Euler(vResult, mResult, x, dx, maxI);
#endif
#ifdef MPEULER
  MPE(vResult, mResult, x, dx, maxI);
#endif
#ifdef RUNGEKUTTA4
  RK4(vResult, mResult, x, dx, maxI);
#endif
}

//////////////////////////////////////////////////////////////////////////////

void displayData(const VECTOR &thisX, const MATRIX &DX, const int &i)
{
  const int WIDTH = 6;

  //cout.setf(ios::fixed);   // Use a fixed number of decimals.
  cout.precision(15);        // Specify the number of decimal places.
  cout.width(WIDTH);         // The first column should have a fixed width.
  cout.setf(ios::scientific);
  cout << "Step " << i << " \t";
  cout.setf(ios::showpos);   // Show the sign.
  cout << "(x,t) = " << thisX << "\t";
  
  int maxI, maxJ; 
  getMaxIJ(maxI, maxJ, thisX);
  cout << "maxI = " << maxI << "; maxJ = " << maxJ << endl; 

#ifdef COMPUTE_C1
  cout << "\t\t[DX] = ";
  for (int i = 1; i <= PRE_DIM; i++)
    if ( i != maxI ) {
      cout << "[";
      int counter = 0;
      for (int j = 1; j <= PRE_DIM; j++)
	if ( j != maxI ) {
	  cout <<  DX(i, j);
	  counter++;
	  if ( counter < PRE_DIM - 1 )
	    cout << "; ";
	}
      cout << "]"<< endl << "\t\t       ";
    }
  cout << endl;
#endif // COMPUTE_C1
  cout.unsetf(ios::scientific);
  cout.unsetf(ios::showpos); // Hide the sign.
}

//////////////////////////////////////////////////////////////////////////////

bool returnDetected(double &dx, const VECTOR &thisX, const int &thisMaxI, 
		    const POINCARE &pInfo) 
{  
  if ( thisMaxI == pInfo.firstMaxI ) {      // If the transversals coincide.
    int currentSign = sign(vectorFieldComponent(thisX, thisMaxI));
    if ( currentSign == pInfo.firstSign ) { // If the flow directions coincide.
      double thisXValue = thisX(thisMaxI);
      double nextXValue = thisXValue + dx;

      // If we are crossing the plane.
      if ( (thisXValue - pInfo.firstX(pInfo.firstMaxI)) * 
	   (nextXValue - pInfo.firstX(pInfo.firstMaxI)) < 0.0 ) {
	dx = currentSign * fabs(thisXValue - pInfo.firstX(pInfo.firstMaxI));
	return true;
      }
    }
  }
  return false;
}

//////////////////////////////////////////////////////////////////////////////

void computeSectionIterates(VECTOR &finalX, const VECTOR &initialX,
			    MATRIX &finalDX)
{
  POINCARE pInfo;   // Get Poincare' map info.
  int initialMaxI, dummyJ;
  getMaxIJ(initialMaxI, dummyJ, initialX);
  pInfo.firstX    = initialX; 
  pInfo.firstMaxI = initialMaxI;            
  pInfo.firstSign = sign(vectorFieldComponent(initialX, initialMaxI));
  
  VECTOR thisX(DIM); // x^(k).
  VECTOR nextX(DIM); // x^(k+1).
  
  for (int i = 1; i <= PRE_DIM; i++)
    thisX(i) = initialX(i);
  thisX(DIM) = 0.0; // Initialize the time variable.
  
  int  step       = 0;
  bool stayInLoop = true;
  
  /*-------------------- BEGIN ALGORITHM 3 ---------------------*/

  int thisMaxI;      // Used to flow in this direction.
  int thisMaxJ;      // Used to determine the step size.
  getMaxIJ(thisMaxI, thisMaxJ, thisX);
  
  MATRIX DX(PRE_DIM, PRE_DIM); id(DX); // Set DX = I.
  VECTOR NextX(DIM);
  displayData(thisX, DX, step);
  
  while ( stayInLoop && (step <= MAX_ITER) )
    {
      // Step 0: Display data.
      //displayData(thisX, DX, step);

      // Step 1: rescale the vector field.
      VECTOR v;
      rescaledVectorField(v, thisX, thisMaxI);

      // Step 2: Compute the step size, and adjust it for a return.
      double S  = fabs(v(thisMaxJ));     // Note that S \in [0, 1].
      double dx = S * MIN_STEP + (1 - S) * MAX_STEP;
      dx *= sign(vectorFieldComponent(thisX, thisMaxI)); 
      static int hits = 0;
      if ( returnDetected(dx, thisX, thisMaxI, pInfo) ) {
	hits++;
	if ( hits == INTERSECTIONS ) {
	  stayInLoop = false;
	  hits = 0;
	}
      }

      // Step 3: Compute the next integration point, nextX.
      // Step 4: Compute the partial derivatives of the rescaled flow, DPsi.
      MATRIX DPsi(PRE_DIM, PRE_DIM);
      computeNextStep(nextX, thisX, DPsi, dx);
      
      // Step 5: Update the dominating coordinates.
      getMaxIJ(thisMaxI, thisMaxJ, nextX);
      
      // Step 6: Compute the partial derivatives of the local Poincare' map.
      MATRIX DPi(PRE_DIM, PRE_DIM);
      computeDPi(DPi, DPsi, nextX, thisMaxI);
      
      // Step 7: Update the accumulated partial derivatives.
      DX = DPi * DX;
      
      thisX = nextX;
      step++;
    } // Step 8: detect a return and break, else continue. 

  /*--------------------- END ALGORITHM 3 -------------------------*/

  displayData(thisX, DX, step);
  finalDX = DX;
  for (int i = 1; i <= PRE_DIM; i++)
    finalX(i) = thisX(i); 
}

//////////////////////////////////////////////////////////////////////////////

void clock(const CLOCK_ORDER &order)
{
  struct tm *ptr;
  static time_t time_start = time(NULL);
  
  if ( order == START )  { // Start the clock      
    ptr = localtime(&time_start);
    cout << "Started computations: " << asctime(ptr) << endl;
  }
  else {                   // Stop the clock
    time_t time_end = time(NULL);
    ptr = localtime(&time_end);
    cout << endl << "Completed computations: " << asctime(ptr)
	 << "Elapsed run time: " << difftime(time_end,time_start) 
	 << " seconds\n"; 
  }
}

//////////////////////////////////////////////////////////////////////////////

int main()
{
  cout << endl << "Computing on the "
#ifdef ROTATION
       << "pure rotation "
#endif
#ifdef VOLTERRA
       << "Volterra "
#endif
#ifdef LORENZ
       << "Lorenz "
#endif
#ifdef ROSSLER
       << "Rossler "
#endif
       << "system." << endl
       << "Using the "  
#ifdef EULER
       << "Euler "
#endif
#ifdef MPEULER
       << "Midpoint Euler "
#endif
#ifdef RUNGEKUTTA4
       << "Runge-Kutta "
#endif
       << "method, order = " << ORDER << "." << endl;
  cout << "The step sizes are taken between "
       << MIN_STEP << " and " << MAX_STEP << "." << endl;

  VECTOR initialX(PRE_DIM);
  initialX(1) = +  1.5;
  initialX(2) = + 14.0;

#ifdef LORENZ  // Stable periodic orbit for R = 250.
  initialX(1) = + 1.621325444114593e+01;
  initialX(2) = - 5.578140243373939e+01;
  initialX(3) = + R - 1;
#endif

#ifdef ROSSLER
  initialX(1) = + 0.0;
  initialX(2) = - 3.9205;
  initialX(3) = + 0.063858;
//   initialX(2) = - 8.38095;
//   initialX(3) = + 0.0295902;
#endif

  clock(START);
  for (int iterate = 0; iterate < 1; iterate++) {
    VECTOR finalX(PRE_DIM);
    MATRIX finalDX(PRE_DIM, PRE_DIM);
    cout << endl << "iterate " << iterate << endl;
    computeSectionIterates(finalX, initialX, finalDX);
    cout.setf(ios::scientific);
    cout << endl << "|initialX - finalX| = " 
	 << Norm(initialX - finalX) << endl;
    if ( Norm(initialX - finalX) < 1e-15 )
      break;
#ifdef LORENZ
#ifdef COMPUTE_C1
    double det = finalDX(1, 1) * finalDX(2, 2) - finalDX(1, 2) * finalDX(2, 1);
    double tr  = finalDX(1, 1) + finalDX(2, 2);
    double lambda1 = tr/2 + sqrt(Sqr(tr/2) - det);
    double lambda2 = tr/2 - sqrt(Sqr(tr/2) - det);
    cout << "lambda = " << tr/2 << " +/- " << sqrt(Sqr(tr/2) - det);
    cout << " = " << lambda1 << "; " << lambda2 << endl;

    double lambda1Star = -6.217323621724878e-03;
    double lambda2Star = -2.989324529670951e-01;  
    cout << "lambda1 - lambda1* = " << lambda1  - lambda1Star << endl;
    cout << "lambda2 - lambda2* = " << lambda2  - lambda2Star << endl;
#endif // COMPUTE_C1
#endif // LORENZ

#ifdef ROSSLER
#ifdef COMPUTE_C1
    double det = finalDX(2, 2) * finalDX(3, 3) - finalDX(2, 3) * finalDX(3, 2);
    double tr  = finalDX(2, 2) + finalDX(3, 3);
    double lambda1 = tr/2 + sqrt(Sqr(tr/2) - det);
    double lambda2 = tr/2 - sqrt(Sqr(tr/2) - det);
    cout << "lambda = " << tr/2 << " +/- " << sqrt(Sqr(tr/2) - det);
    cout << " = " << lambda1 << "; " << lambda2 << endl;
#endif // COMPUTE_C1
#endif // ROSSLER
    cout.unsetf(ios::scientific);
    initialX = finalX;
  }
  clock(STOP);

  return 0;
}

//////////////////////////////////////////////////////////////////////////////