/*   File: validate.cc 
 
     A non-rigorous validation program for
     the rigorous RODES computations.
 
     Latest edit: Tue Jan 16, 2001

     Compilation: g++ -O3 -D__I386__ -Wall -o validate validate.cc 
                  -I/PROFIL/incl -IBIAS -L/PROFIL/lib -lProfil -LBIAS -lBias

     In the compilation line, PROFIL is the path to PROFIL's incl/ and lib/ 
     directories. If you don't have PROFIL installed, you must create classes
     that correspond to INTERVAL etc. Directed rounding is NOT neccessary.

*/

////////////////////////////////////////////////////////////////////////////
// Standard header files.

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

////////////////////////////////////////////////////////////////////////////
// PROFIL/BIAS header files.

#include "Constants.h"
#include "Functions.h"
#include "Interval.h"
#include "IntervalMatrix.h"
#include "IntervalVector.h"
#include "Functions.h"
#include "Matrix.h"
#include "Vector.h"

////////////////////////////////////////////////////////////////////////////
// User-defined header files.

#include "list.h"

////////////////////////////////////////////////////////////////////////////
// Short hand.

#define BOX INTERVAL_VECTOR
#define I_M INTERVAL_MATRIX
#define CMDSHELL "sh"
#define CMDOPT "-c"
#define LoadsOfOutput // Gives lots of output.
#undef  LoadsOfOutput

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

class InData
{
public:
  BOX in_box;
  BOX ret_box;
  INTERVAL angle;
  double pre_exp;
  double min_exp;
};

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

// Classical constants for the 'lorenz' system.
const double S = 10.0;
const double R = 28.0;
const double B = 8./3;

// Parameters in the Jordan normal form.
const double TEMP = sqrt((S + 1) * (S + 1) + 4 * S * (R - 1));
const double K1 = S / TEMP;
const double K2 = (S - 1 + TEMP) / (2 * S);
const double K3 = (S - 1 - TEMP) / (2 * S);

// Eigenvalues at the origin.
const double E1 = (- (S + 1) + TEMP) / 2;
const double E2 = (- (S + 1) - TEMP) / 2;
const double E3 = - B;

// Shortcuts for DF.
const double K2_PLUS_K3 = K2 + K3; 
const double TWO_K2     = 2 * K2; 
const double TWO_K3     = 2 * K3; 

// The step size for the RK4 integrator.
const double STEP_SIZE = 0.01;

// x3-tolerance for a return to be valid.
const double TOLERANCE = 1e-8;

// Scale factor for the input data.
const double SCALE = pow(0.5, 8);

// Number of points the sides of the initial
// rectangles are split into. Must be greater than one.
const int NUMBER_OF_POINTS = 10;

// The size of the points to be plotted.
const double POINT_SIZE = pow(0.5, 10);

// Trigonometric constants.
const double PI = 2 * asin(1);
const double DEG_TO_RAD = PI / 180.0;
const double RAD_TO_DEG = 180.0 / PI;

////////////////////////////////////////////////////////////////////////////
// Function declarations

VECTOR SmallF            (const VECTOR &x);
MATRIX DF                (const VECTOR &x);
VECTOR F                 (const VECTOR &x);
VECTOR RK4               (const VECTOR &x, double h);
VECTOR ReturnMap         (const VECTOR &x);
void   GetTheFlags       (const int &argc, char *argv[], char *SourceFile);
void   GetAllInData      (const char *SourceFile, List<InData> &InDataList);
void   ValidateAllInData (List<InData> &InDataList);
void   HardCopy          (List<InData> &idList, List<InData> &InDataList);
void   GetReturnList     (List<InData> &ReturnList, const InData &id, 
			  List<InData> &idList);
bool   ValidateOneInData (const InData &id, List<InData> &ReturnList, 
			  const int &NumberOfPoints);
void   CheckC1Info       (const InData &id, const InData &RetBox, 
			  const VECTOR &x_ret);

////////////////////////////////////////////////////////////////////
//
//                     ODE-DEPENDENT FUNCTIONS
//
////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////
// The Lorenz vector field in normal form. Here x is 3-dimensional.

VECTOR SmallF(const VECTOR &x)
{
  VECTOR result(3);
  double C1 = x(1) + x(2);
  double C2 = K1 * C1 * x(3);

  result(1) = E1 * x(1) - C2;
  result(2) = E2 * x(2) + C2;
  result(3) = E3 * x(3) + C1 * (K2 * x(1) + K3 * x(2));

  return result;
}

////////////////////////////////////////////////////////////////////
// The Lorenz partial derivatives.

MATRIX DF(const VECTOR &x)
{
  MATRIX result(3, 3);
  double C1 = K1 * (x(1) + x(2));
  double C2 = K1 * x(3);

  result(1, 1) = E1 - C2;
  result(1, 2) = - C2;
  result(1, 3) = - C1;
  result(2, 1) = + C2;
  result(2, 2) = E2 + C2;
  result(2, 3) = + C1;
  result(3, 1) = TWO_K2 * x(1) + K2_PLUS_K3 * x(2);
  result(3, 2) = K2_PLUS_K3 * x(1) + TWO_K3 * x(2);
  result(3, 3) = E3;

  return result;
}
////////////////////////////////////////////////////////////////////
// The extended Lorenz system. Here x is 12-dimensional.

VECTOR F(const VECTOR &x)
{
  VECTOR result(12);
  double C1 = x(1) + x(2);
  double C2 = K1 * C1 * x(3);
 
  result(1) = E1 * x(1) - C2;
  result(2) = E2 * x(2) + C2;
  result(3) = E3 * x(3) + C1 * (K2 * x(1) + K3 * x(2));

  MATRIX M   = DF(x); // Only uses the three first elements of x.
  result(4)  = M(1, 1) * x(4)  + M(1, 2) * x(5)  + M(1, 3) * x(6);
  result(5)  = M(2, 1) * x(4)  + M(2, 2) * x(5)  + M(2, 3) * x(6);
  result(6)  = M(3, 1) * x(4)  + M(3, 2) * x(5)  + M(3, 3) * x(6);

  result(7)  = M(1, 1) * x(7)  + M(1, 2) * x(8)  + M(1, 3) * x(9);
  result(8)  = M(2, 1) * x(7)  + M(2, 2) * x(8)  + M(2, 3) * x(9);
  result(9)  = M(3, 1) * x(7)  + M(3, 2) * x(8)  + M(3, 3) * x(9);

  result(10) = M(1, 1) * x(10) + M(1, 2) * x(11) + M(1, 3) * x(12);
  result(11) = M(2, 1) * x(10) + M(2, 2) * x(11) + M(2, 3) * x(12);
  result(12) = M(3, 1) * x(10) + M(3, 2) * x(11) + M(3, 3) * x(12);

  return result;
}

////////////////////////////////////////////////////////////////////
// The basic Runge-Kutta integrator for F.

VECTOR RK4(const VECTOR &x, double h)
{
  VECTOR k1 = h * F(x);
  VECTOR k2 = h * F(x + 0.5 * k1);
  VECTOR k3 = h * F(x + 0.5 * k2);
  VECTOR k4 = h * F(x + k3);
 
  return x + (k1 + 2 * (k2 + k3) + k4) / 6;
}

////////////////////////////////////////////////////////////////////
// The Poincare' return map for the plane x(3) = R - 1.

VECTOR ReturnMap(const VECTOR &x)
{
  bool   FULL_RETURN = false;
  double StepSize = STEP_SIZE;
  VECTOR x_old = x;
  VECTOR x_new;

  while ( !FULL_RETURN )
    {
      x_new = RK4(x_old, StepSize);

      if ( x_old(3) > R - 1 && x_new(3) <= R - 1 )
	{    
	  if ( fabs(x_new(3) - (R - 1)) < TOLERANCE )
	    FULL_RETURN = true;
	  else
	    StepSize *= 0.5;
	}
      else
	x_old = x_new;
    }

  return x_new;
}

////////////////////////////////////////////////////////////////////
//
//                     THE PROGRAM PROPER
//
////////////////////////////////////////////////////////////////////

int main(int argc, char *argv[])
{
  char SourceFile[99] = "";
  List<InData> InDataList;

  GetTheFlags(argc, argv, SourceFile);
  GetAllInData(SourceFile, InDataList);
  ValidateAllInData(InDataList);

  return 0;
}

////////////////////////////////////////////////////////////////////
//
//                     ODE-INDEPENDENT FUNCTIONS
//
////////////////////////////////////////////////////////////////////

void GetTheFlags(const int &argc, char *argv[], char *SourceFile)
{
  if ( argc != 2 )
    { 
      cout << endl << "Usage: " << endl;
      cout << "\t" << argv[0] << " <source_file>" << endl << endl;
      exit(0);
    }  
  strcpy(SourceFile, argv[1]); 
  cout << "Reading data from " << SourceFile << endl;
}

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

void GetAllInData(const char *SourceFile, List<InData> &InDataList)
{
  FILE *InputFile;
  InputFile = fopen(SourceFile,"r");      

  if ( !InputFile )
    {
      cout << "File " << SourceFile << " could not be opened!" << endl;
      exit(1);
    }
 
  char   InputString[200];
  int    x_in, y_in;
  int    x_ret_lo, y_ret_lo, x_ret_hi, y_ret_hi;
  int    d1, d2, d3, d4, d5; // Dummy variables.
  double f1, f2, f3, f4;
  BOX    InBox(3);  InBox(3)  = R - 1;
  BOX    RetBox(3); RetBox(3) = R - 1;
  InData id;
  
  while ( !feof(InputFile) )	     	
    {
      fgets(InputString, sizeof(InputString), InputFile);
      sscanf(InputString, "%d %d %d %d %d %lf %lf %lf %lf %d %d %d %d %d %d", 
	     &x_in, &y_in, &d1, &d2, &d3, &f1, &f2, &f3, &f4, 
	     &x_ret_lo, &y_ret_lo, &d4, &x_ret_hi, &y_ret_hi, &d5);
      InBox(1)  = SCALE * Hull(x_in - 1, x_in + 1);
      InBox(2)  = SCALE * Hull(y_in - 1, y_in + 1);
      RetBox(1) = SCALE * Hull(x_ret_lo - 1, x_ret_hi + 1);  
      RetBox(2) = SCALE * Hull(y_ret_lo - 1, y_ret_hi + 1);
      id.in_box  = InBox; 
      id.ret_box = RetBox; 
      id.angle = Hull(f1, f2);
      id.pre_exp = f3;
      id.min_exp = f4;
      InDataList += id;
    }
  
  cout << "Read " << Length(InDataList) << " lines of indata." << endl;
  fclose(InputFile);
}

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

void ValidateAllInData(List<InData> &InDataList)
{
  unsigned int counter = 1;
  unsigned int errors  = 0;
  bool         valid_data;
  List<InData> idList;

  HardCopy(idList, InDataList);
  First(InDataList);
  while( !Finished(InDataList) ) 
    { 
      List<InData> ReturnList;
      InData id  = Current(InDataList);

      GetReturnList(ReturnList, id, idList);
      valid_data = ValidateOneInData(id, ReturnList, NUMBER_OF_POINTS);

      if ( valid_data )
	cout << "Inclusion ";
      else
	{
	  cout << "Leakage ";	
	  errors++;
	}
      cout << "for rectangle # " << counter << endl;

#ifdef LoadsOfOutput
      cout << "\t x = (" << Mid(id.in_box(1)) << "," 
	   << Mid(id.in_box(2)) << "); "
	   << "angle = " << id.angle << endl;
      cout << "Number of registerd returns: " 
	   << Length(ReturnList) << endl;
#endif

      Next(InDataList);
      counter++;
    }
  cout << endl << endl << errors << " rectangles" 
       << (errors != 1 ? "s " : " ") << "had leakage." << endl;   
}

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

void HardCopy(List<InData> &idList, List<InData> &InDataList)
{
  First(InDataList);
  while( !Finished(InDataList) ) 
    { 
      idList += Current(InDataList);
      Next(InDataList);
    }
  cout << "HardCopy complete." << endl;
}

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

void GetReturnList(List<InData> &ReturnList, const InData &id, 
		   List<InData> &idList)
{
  BOX    GlobalRetBox(2);
  BOX    LocalRetBox(2);
  InData ID;

  GlobalRetBox(1) = id.ret_box(1);
  GlobalRetBox(2) = id.ret_box(2);

  First(idList);
  while( !Finished(idList) ) 
    { 
      BOX dummy(2);
      ID = Current(idList);
      LocalRetBox(1) = ID.in_box(1);
      LocalRetBox(2) = ID.in_box(2);
      if( Intersection(dummy, LocalRetBox, GlobalRetBox) )
	ReturnList += ID;
      if( Intersection(dummy, - LocalRetBox, GlobalRetBox) )
	{
	  ID.in_box(1)  = - ID.in_box(1);
	  ID.in_box(2)  = - ID.in_box(2);
	  ID.ret_box(1) = - ID.ret_box(1);
	  ID.ret_box(2) = - ID.ret_box(2);
	  ReturnList += ID;
	}
      Next(idList);
    }  
}

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

bool ValidateOneInData(const InData &id, List<InData> &ReturnList, 
		       const int &NumberOfPoints)
{
  bool         total_inclusion = true;
  char         answer;
  INTERVAL     x1 = id.in_box(1);
  INTERVAL     x2 = id.in_box(2);
  double       dx1 = Diam(x1) / (NumberOfPoints - 1);
  double       dx2 = Diam(x2) / (NumberOfPoints - 1);
  VECTOR       x_in(12), x_ret(12);
  List<VECTOR> VectorList;
#ifdef LoadsOfOutput
  double       MaxPreExp = 0.0;
#endif

  x_in(3) = Mid(id.in_box(3)); // Always x3 = R - 1. 
  
  // The identity matrix.
  x_in(4)  = 1; x_in(5)  = 0; x_in(6)  = 0; 
  x_in(7)  = 0; x_in(8)  = 1; x_in(9)  = 0; 
  x_in(10) = 0; x_in(11) = 0; x_in(12) = 1; 

  for ( int i = 0; i < NumberOfPoints; i++ )
    {
      x_in(1) = Inf(x1) + i * dx1;
      for ( int j = 0; j < NumberOfPoints; j++ ) 
	{
	  x_in(2) = Inf(x2) + j * dx2;

	  // Compute the (i,j):th return.
	  x_ret = ReturnMap(x_in);
	  VectorList += x_ret;

	  // Validate inclusion.
	  bool inclusion = false;
	  First(ReturnList);
	  InData RetBox;
	  while( !Finished(ReturnList) ) 
	    { 
	      RetBox = Current(ReturnList);
	      if ( (x_ret(1) <= RetBox.in_box(1)) && 
		   (x_ret(2) <= RetBox.in_box(2)) )
		{
		  inclusion = true;
		  CheckC1Info(id, RetBox, x_ret);	
#ifdef LoadsOfOutput
		  MaxPreExp = Max(MaxPreExp, RetBox.pre_exp);
#endif	    
		  break;
		}
	      Next(ReturnList);
	    }  
	  if ( !inclusion ) // No inclusion in ReturnList.
	    {
	      total_inclusion = false;
	      cout << "No inclusion for subdata (" 
		   << i << "," << j << "):" << endl;
	      cout << "\t The point  " << x_in << " from " 
		   << endl << "\t " << id.in_box << endl 
		   << "\t returns as " << x_ret 
		   << " which is NOT in " << endl
		   << "\t " << id.ret_box << endl << endl; 
	      cin >> answer;
	    }
	} // Done with the (i,j):th initial point.
    }

#ifdef LoadsOfOutput
  cout << "MaxPreExp = " << MaxPreExp << endl;
#endif

  return total_inclusion;
}

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

void CheckC1Info(const InData &id, const InData &RetBox, 
		 const VECTOR &x_ret)
{
  bool   error         = false;
  double InitialUAngle = DEG_TO_RAD * Inf(id.angle);
  double InitialVAngle = DEG_TO_RAD * Sup(id.angle);

  VECTOR UVector(2); 
  VECTOR VVector(2); 

  UVector(1) = cos(InitialUAngle); 
  UVector(2) = sin(InitialUAngle);  
  VVector(1) = cos(InitialVAngle); 
  VVector(2) = sin(InitialVAngle);  

  MATRIX DPhi(3, 3);
  DPhi(1, 1) = x_ret(4);  DPhi(1, 2) = x_ret(7);  DPhi(1, 3) = x_ret(10);
  DPhi(2, 1) = x_ret(5);  DPhi(2, 2) = x_ret(8);  DPhi(2, 3) = x_ret(11);
  DPhi(3, 1) = x_ret(6);  DPhi(3, 2) = x_ret(9);  DPhi(3, 3) = x_ret(12);

  VECTOR vf = SmallF(x_ret);
  MATRIX DP(2, 2);
  DP(1, 1) = DPhi(1, 1) - vf(1) / vf(3) * DPhi(3, 1);
  DP(1, 2) = DPhi(1, 2) - vf(1) / vf(3) * DPhi(3, 2);
  DP(2, 1) = DPhi(2, 1) - vf(2) / vf(3) * DPhi(3, 1);
  DP(2, 2) = DPhi(2, 2) - vf(2) / vf(3) * DPhi(3, 2);

  VECTOR RetUVector = DP * UVector;
  VECTOR RetVVector = DP * VVector;

  double FinalUAngle = RAD_TO_DEG * atan(RetUVector(2) / RetUVector(1));
  double FinalVAngle = RAD_TO_DEG * atan(RetVVector(2) / RetVVector(1));
  INTERVAL FinalAngle = Hull(FinalUAngle, FinalVAngle);

  if ( !(FinalAngle <= RetBox.angle) )
    {
      error = true;
      cout << "Cone leakage: " << endl 
	   << FinalAngle << " not in " << RetBox.angle << endl;
    }
#ifdef LoadsOfOutput
  else
    cout << "\t Cone inclusion: " 
   	 << FinalAngle << " in " << RetBox.angle << endl;
#endif

  double Expansion = Min(Norm(RetUVector), Norm(RetUVector));
  if ( Expansion < id.min_exp )
    {
      error = true;
      cout << "Expansion leakage: " << endl 
	   << Expansion << " less than " << id.min_exp << endl;
    } 
#ifdef LoadsOfOutput 
  else
    cout << "\t Expansion bound good: "
	 << Expansion << " greater than " << id.min_exp << endl; 
#endif

  if ( Expansion < RetBox.pre_exp )
    {
      error = true;
      cout << "Pre-expansion leakage: " << endl 
	   << Expansion << " less than " << RetBox.pre_exp << endl;
    }
#ifdef LoadsOfOutput
  else
    cout << "\t Pre-expansion bound good: "
	 << Expansion << " greater than " << RetBox.pre_exp << endl; 
#endif

  if ( error )
    {
      char answer;
      cout << "Initial angles: "
	   << "U = " << Inf(id.angle) << " V = " << Sup(id.angle) << endl;
      cout << "Initial vectors: "
	   << "U = " << UVector << " V = " << VVector << endl;
      cout << "DPhi = " << endl << DPhi << endl;
      cout << "vf = " << vf << endl;
      cout << "DP = " << endl << DP << endl;
      cout << "Return vectors: "
	   << "R(U) = " << RetUVector << " R(V) = " << RetVVector << endl; 
      cout << "Return angles: "
	   << "R(U) = " << FinalUAngle << " R(V) = " << FinalVAngle << endl;
      cin >> answer;
    }
}

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