#include <fstream>
#include <iostream>
#include <iomanip>
#include "Vector.hh"
#include "Matrix.hh"
#include "Solve.hh"
#include "Clock.hh"
#include "Random.hh"
#include "PeriodicOrbit.hh"
#include <string.h>
#include "lorenz.hh"
#include "lorenz1.hh"
#include "lorenz2.hh"
#include "lorenz3.hh"
#include "lorenz4.hh"
#include <stdlib.h>

static int NUMITNS = 2000;

/* vfield and jacobian */
void lorenz(const Vector& x, Vector& y){
  double s = 10.0;
  double r = 28.0;
  double b = 8.0/3.0;
  y(0) = s*(x(1)-x(0));
  y(1) = x(0)*(r-x(2)) - x(1);
  y(2) = x(0)*x(1) - b*x(2);
}
 
 
void lorenzJ(const Vector& in, Matrix & M){
  double s = 10.0;
  double r = 28.0;
  double b = 8.0/3.0;
  double x = in(0), y = in(1), z = in(2);
  M(0,0) = -s;
  M(0,1) = s;
  M(0,2) = 0;
  M(1,0) = r-z;
  M(1,1) = -1;
  M(1,2) = -x;
  M(2,0) = y;
  M(2,1) = x;
  M(2,2) = -b;
}

//accurately shifts po so that column 0 is on the Poincare section
void adjustPO(Matrix * po){
  int n = po->getn();
  int ind;
  findTransIndx(*po, 0, ind);
  double h;
  int flag = 1;
  findShift(*po, ind, h, flag);
  Matrix temp(3,n);
  shift(*po, h, temp);
  int i;
  for(i=0; i < n; i++)
    (*po).getColumn(i) = temp.getColumn((i+ind)%n);
}


int checkSymString(Matrix& po, char *s){
  char* ss = new char[strlen(s)+1];//verify symbol sequence
  findSymString(po, ss);
  int result = isequal(s, ss);
  if(result==0){
    cout<<"checkSymString failed: "<<endl;
    cout<<"s = "<<s<<endl;
    cout<<"ss = "<<ss<<endl;
  }
  else{
    cout<<"s1 = "<<s<<endl;
    cout<<"s2 = "<<ss<<endl;
  }
  delete[] ss;
  return result;
}

int checkInputString(char *s){
  int i;
  for(i=0; i < strlen(s); i++)
    if((s[i] != 'A')&&(s[i]!='B'))
      return 0;
  return 1;
}


//s(input) = desired symbol string
//T(output) = period of po
//err(output) = error of po
//sintvl = shooting interval for PO
//return = ptr to po allocated by routine called from inside
Matrix *returnPO(char *s, double &T, double& err, int sintvl){
  Matrix *po;
  if(!checkInputString(s)){
    cout<<"input string s has characters other than A and B: "<<s<<endl;
    err = 10000000.0;
    return NULL;
  }
  po = returnApprox(s, T);
  if(!checkSymString(*po, s))
    cout<<"checkSymString: returnApprox gave bad guess for "<<s<<endl;
#ifdef DEBUG
  else
    cout<<"returnApprox OK"<<endl;
#endif
  PeriodicOrbit PO(*po, T, lorenz, lorenzJ);
  PO.setSintvl(sintvl);
  int n = po->getn();
  int i;
  Matrix step(3,n);
  double wstep;
  double delta = -1;
  int icount = 0;
  for(i=0; i < 20; i++){
    PO.smooth();
    err = PO.relativeResidualError();
    cout<<"RELATIVE ERROR: "<<err<<endl;
#ifdef DEBUG    
    cout<<"delta: "<<delta<<endl;
#endif
    if(err < 4.0e-14)
      if(icount > 0)
	break;
      else
	icount++;
    ComputeTrustStep(PO, delta, step, wstep);
    PO.addStep(step, wstep);
  }
  err = PO.relativeResidualError();
  if(err < 1e-13){
    PO.getPeriodicOrbit(*po, T);
    if(!checkSymString(*po, s))
       err = 100.0;
    adjustPO(po);
    return po;
  }
#ifdef DEBUG
  cout<<"MAKING SECOND PASS"<<endl;
#endif
  delta = -1;
  int numItns=NUMITNS;
  int j=0;
  for(i=0; i < numItns; i++){
    PO.smooth();
    err = PO.relativeResidualError();
    cout<<"RELATIVE ERROR: "<<err<<" (i="<<i<<")"<<endl;
    /*if(i%100==0){
      cout<<"refinePO: Saving ..."<<endl;
      PO.getPeriodicOrbit(*po, T);
      saveOrbit(*po, T, err, "DATAX\\", i/100);
      }*/
#ifdef DEBUG
    cout<<"delta: "<<delta<<endl;
#endif
    if(err < 1e-10)
      j++;
    if(j > 3)
      break;
    ComputeTrustStep(PO, delta, step, wstep);
    PO.addStep(step, wstep);
    
  }
  PO.getPeriodicOrbit(*po, T);
  err = PO.relativeResidualError();
  if(!checkSymString(*po, s))
     err = 100.0;
  adjustPO(po);
  return po;
}

//default sintvl is 64;
Matrix *returnPO(char *s, double& T, double& err){
  return returnPO(s,T,err,64);
}
  
void  loadOrbits(){
  static int status = 0;
  if(status == 0){
    char dir[50] = "DATA1/";
    loadOrbits(dir, 1, 1375);
    status = 1;
  }
}