#include <fstream>
#include <iostream>
#include "Vector.hh"
#include "Matrix.hh"
#include <string.h>
#include <stdio.h>
#include "lorenz1.hh"
#include "lorenz3.hh"
const double ZP = 27.0;
/**************************************************/
/* global variables */
static const int MAXNUM = 2000;//Maximum size of library
static const int LARGENUM = 1000000;
static Matrix *poList[MAXNUM];
static char *symList[MAXNUM];
static double periodList[MAXNUM];
static size_t numOfPos = 0;



/**************************************************/
/* symbol string related functions */
void stringfit(char *str1, int indx1, char *str2, int indx2, 
	       int &lengthOfFit){
  int len1 = strlen(str1);
  int len2 = strlen(str2);
  assert((0<=indx1) && (indx1 < len1));
  assert((0<=indx2) && (indx2 < len2));
  int fwdfit = 0;
  int i1 = indx1; 
  int i2 = indx2;
  while(str1[i1] == str2[i2]){
    fwdfit++;
    i1 = (i1 + 1)%len1;
    i2 = (i2 + 1)%len2;
    if(fwdfit >= len1*len2){
      fwdfit = LARGENUM;
      break;
    }
  }
  int bwdfit = 0;
  i1 = indx1;
  i2 = indx2;
  while(str1[i1]==str2[i2]){
   bwdfit++;
   i1 = (i1 + (len1-1))%len1;
   i2 = (i2 + (len2-1))%len2;
   if(bwdfit >= len1*len2){
     bwdfit = LARGENUM;
     break;
   }
  }
  if(bwdfit > 0)
    bwdfit = 10 + (bwdfit-1)*14;
  lengthOfFit = (bwdfit < fwdfit)? bwdfit:fwdfit;
  //lengthOfFit = fwdfit;
}
  
void stringbestfit(char *str1, int indx1, char *str2,
		   int& indx2, int& lengthOfFit){
  int len1 = strlen(str1);
  int len2 = strlen(str2);
  assert((0<=indx1)&&(indx1<len1));
  assert((len1 > 0) && (len2 > 0));
  int i, maxi, fiti, fitmaxi;
  i = 0;
  maxi = 0;
  stringfit(str1, indx1, str2, i, fiti);
  fitmaxi = fiti;
  for(i = 1; i < len2; i++){
    stringfit(str1, indx1, str2, i, fiti);
    if(fiti > fitmaxi){
      maxi = i;
      fitmaxi = fiti;
    }
  }
  indx2 = maxi;
  lengthOfFit = fitmaxi;
}

//function below touches global data
void stringbestfit(char *str1, int indx1,
		   int& poindx, int& symindx, int& lengthOfFit){
  int i=0;
  int si, lfit;
  lengthOfFit = -100;
  for(i=0; i < numOfPos; i++){
    stringbestfit(str1, indx1, symList[i], si, lfit);
    if(lfit > lengthOfFit){
      lengthOfFit = lfit;
      poindx = i;
      symindx = si;
    }
  }
}

//are the two strings equal as pos?
int isequal(char *str1, char* str2){
  int lfit;
  int indx2;
  stringbestfit(str1, 0, str2, indx2, lfit);
  if(lfit == LARGENUM)
    return 1;
  else
    return 0;
}

//is one string a flip of the other as pos?
int isflip(char *str1, char *str2){
  int len1 = strlen(str1);
  int len2 = strlen(str2);
  char *str3 = (char *)malloc((len2+1)*sizeof(char));
  int i;
  for(i=0; i < len2; i++)
    str3[i] = (str2[i]=='A')? 'B':'A';
  str3[len2] = '\0';
  int ans;
  ans = isequal(str1, str3);
  free(str3);
  return(ans);
}


/**************************************************/
/* functions related to pos and symbol strings */
//po = input 3xn matrix
//s = its sym string after return
void findSymString(Matrix& po, char *s){
  assert(po.getm()==3);
  int n = po.getn();
  Matrix x(po, 0, 0, 0, n-1);  
  Matrix y(po, 1, 1, 0, n-1);
  Matrix z(po, 2, 2, 0, n-1);
  int i, si;
  si = 0;
  for(i=0; i < n ; i++){
    int j = (i+1)%n;
    if((z(i) <= ZP) && (z(j) > ZP)){
      if(y(i) < 0)
	s[si] = 'A';
      else
	s[si] = 'B';
      si++;
    }
  }
  s[si] = '\0';
}

//po(input) =  3xn matrix representing po
//symIndx(input) = index of symbol in its symbol string
//transIndx(output) = column index into po just preceding 
//                   occurrence of indexed symbol
void findTransIndx(Matrix& po, int symIndx, 
		   int& transIndx){
  assert(po.getm()==3);
  assert(symIndx >= 0);
  int n = po.getn();
  Matrix z(po, 2, 2, 0, n-1);
  int i, si;
  si = -1;
  for(i=0; i < n; i++){
    int j = (i+1)%n;
    if((z(i) <= ZP) && (z(j) > ZP))
      si++;
    if(si == symIndx){
      transIndx = i;
      return;
    }
  }
  cout<<"symIndx is too large!"<<endl;
}

//po(input) = 3xn Matrix
//transIndx(input) = column index just preceding a symbol
//h(output) = shift (taking period=2*PI) so that column index
//            now points to a Poincare section
//flag(input)(optional) = if flag!=0 shift is more accurate.
void findShift(Matrix& po, int transIndx,
	       double &h, int flag){
  int n = po.getn();
  Matrix ztemp(po, 2, 2, 0, n-1);
  Vector z(n), zz(n);
  int j;
  for(j=0; j < n; j++)
    z(j) = ztemp(j);
  int i = transIndx;
  assert((0<=i)&&(i<n));
  double t1 = 0;
  double x1 = z(i)-ZP;
  double t2 = 2*PI/n;
  double x2 = z((i+1)%n)-ZP;
  double tt;
  double dx1;
  double tol = (flag==0)? 1e-9:1e-12;
  int numItns = (flag==0)? 10:20;
  

  for(j=0; j < numItns; j++){//upto 10 or 20 iterations of the secant method
    if(fabs(x1) < tol)
      break;
    if(fabs(x2) < tol){
      t1 = t2;
      break;
    }
    dx1 = (x2 - x1)/ (t2 - t1);
    tt = t1 - x1/dx1;
    t2 = t1;
    x2 = x1;
    t1 = tt;
    shift(z, t1, zz);
    x1 = zz(i)-ZP;
  }
  assert((t1 >= 0) && (t1 <= 2*PI/n));
  h = t1;
}

//po(input) = 3xn matrix
//transIndx(input) = column index of po
//transIndx(output) = colum index > transIndx1 just preceding a symbol
void findNextTransIndx(Matrix& po, int transIndx1,
		       int& transIndx2){
  int n = po.getn();
  assert((0<=transIndx1)&&(transIndx1 < n));
  Matrix z(po, 2, 2, 0, n-1);
  int j = (transIndx1+1)%n;
  while(!((z(j)<=ZP)&&(z((j+1)%n)>ZP)))
    j = (j+1)%n;
  transIndx2 = j;
}
 
//po(input) = 3xn Matrix
//T(input) = period of po
//transIndx(input) = column index of po preceding a symbol
//deltat(output) = time to hit another symbol from far end of slice
//returned Matrix = 4xn Matrix that slices po from one symbol to next
Matrix *returnSlice(Matrix& po, double T, int transIndx,
		    double& deltat){
  int n = po.getn();
  int i1 = transIndx;
  int i2;
  findNextTransIndx(po, i1, i2);
  double h1, h2;
  findShift(po, i1, h1);
  findShift(po, i2, h2);
  deltat = h2*T/(2*PI);
  int slicen;
  if(i1 < i2)
    slicen = (i2-i1)+1;
  else
    slicen = ((i2+n)-i1)+1;
  Matrix* sliceptr = new Matrix(4, slicen);
  double h = T/n;
  h1 = h1*T/(2*PI);
  int j;
  for(j=0; j < slicen; j++){
    (*sliceptr)(0,j) = po(0,(i1+j)%n);
    (*sliceptr)(1,j) = po(1,(i1+j)%n);
    (*sliceptr)(2,j) = po(2,(i1+j)%n);
    (*sliceptr)(3,j) = j*h-h1;
  }

  Matrix tempo(3,n);
  shift(po, h1*(2*PI)/T, tempo);
  (*sliceptr)(0,0) = tempo(0,i1);
  (*sliceptr)(1,0) = tempo(1,i1);
  (*sliceptr)(2,0) = tempo(2,i1);
  (*sliceptr)(3,0) = 0.0;
  return sliceptr;
}

//slice1(input) = ptr to 4xn slice (deallocated inside)
//slice2(input) = ptr to 4xn slice (deallocated inside)
//deltat1(input) = time from end of slice1 to beginning of slice2
//return = pointer to another 4xn slice  allocated inside
Matrix* returnMergedSlice(Matrix *slice1, double deltat1, Matrix *slice2){
  int n1 = (*slice1).getn();
  int n2 = (*slice2).getn();
  Matrix *slice = new Matrix(4, n1+n2);
  int j;
  double T = (*slice1)(3,n1-1)+deltat1;
  for(j=0; j < n2; j++)
    (*slice2)(3,j) += T;
  Matrix shadowm;
  shadowm.shadow((*slice), 0, 3, 0, n1-1);
  shadowm = (*slice1);
  shadowm.shadow((*slice), 0, 3, n1,n1+n2-1);
  shadowm = (*slice2);
  delete slice1;
  delete slice2;
  return slice;
}

//po(input) = 3xn matrix
//T(input) = period of po
//transIndx(input) = column index of po preceding a symbol
//vd(output) = time derivative of po at left end of slice
void findLeftDervOfSlice(Matrix& po, double T, int transIndx,
			 Vector& vd){
  assert(vd.getSize()==3);
  int i1 = transIndx;
  double h1;
  findShift(po, i1, h1);
  int n = po.getn();
  Matrix temp(3,n);
  shift(po, h1, temp);
  Matrix dtemp(3,n);
  diff(temp, dtemp);
  dtemp.scale(2*PI/T);
  vd = dtemp.getColumn(i1);
}

//po(input) = 3xn matrix
//T(input) = period of po
//transIndx(input) = column index of po preceding slice
//vd(output) = time derivative at right end of slice 
void findRightDervOfSlice(Matrix& po, double T, int transIndx,
			  Vector& vd){
  assert(vd.getSize()==3);
  int i1 = transIndx;
  int i2;
  findNextTransIndx(po, i1, i2);
  int n = po.getn();
  Matrix dtemp(3,n);
  diff(po, dtemp);
  dtemp.scale(2*PI/T);
  vd = dtemp.getColumn(i2);
}

/**************************************************/
/* routines for loading pos and a driver */

void initLorenz1(){
}

//dir(input) = directory
//indx(input) = index of orbit in that directory
//T(output) = period of loaded orbit
//err(output) = error of loaded orbit
//return = ptr of Matrix holding orbit allocated inside
Matrix *loadANDreturnOrbit(char *dir, int indx,
			   double& T, double& err){
  char s1[100] = "orbit";
  char s2[100] = "T";
  char s3[100] = "n";
  char s4[100] = "err";
  char tempstr[100];
  strcpy(tempstr, dir); strcat(tempstr, s1); strcpy(s1, tempstr);
  strcpy(tempstr, dir); strcat(tempstr, s2); strcpy(s2, tempstr);
  strcpy(tempstr, dir); strcat(tempstr, s3); strcpy(s3, tempstr);
  strcpy(tempstr, dir); strcat(tempstr, s4); strcpy(s4, tempstr);
  char s[50];
  sprintf(s, "%d.dat", indx);
  strcat(s1, s);
  strcat(s2, s);
  strcat(s3, s);
  strcat(s4, s);
  Matrix tt(1,1);
  tt.input(s2);
  T = tt(0,0);
  tt.input(s3);
  int n = tt(0,0);
  tt.input(s4);
  err = tt(0,0);
  Matrix *po = new Matrix(3, n);
  po->input(s1);
  return po;
}


  
//N(input) = index of orbit
//dir(input) = directory
//loads into library maintained using global variables  
void loadOrbit(int N, char *dir){
  //cout<<"loadOrbit: point A"<<endl;
  double T;
  double err;
  int index = numOfPos;
  numOfPos++;
  assert(numOfPos < MAXNUM);
  poList[index] = loadANDreturnOrbit(dir, N, T, err);
  //cout<<"loadOrbit: point B"<<endl;
  Matrix po(*(poList[index]));
  int n = po.getn();
  periodList[index] = T;
  double h = 0.8*(2*PI)/n;
  Matrix temp(3,n);
  shift(po, h, temp);//shift the loaded orbit by less than an interval
  //cout<<"loadOrbit: point C"<<endl;
  po = temp;
  int len;
  len = (n/256+1)*4+30;
  symList[index] = (char *)malloc(len*sizeof(char));
  findSymString(po, symList[index]);
  //cout<<"loadOrbit: point D"<<endl;
}

//s(input) = name of directory
//loads all orbits from m(input) to n(input) into library
void loadOrbits(char *s, int m, int n){
  int i;
  for(i=m; i <= n; i++){
    cout<<"loading orbit #"<<i<<" from "<<s<<" ..."<<endl;
    loadOrbit(i, s);
  }
}

//po(input) = 3xn matrix
//T(input) = its period
//err(input) = its error
//dir(input) = directory to save the orbit
//indx(input) = index to use for saving
//saves orbit in directory with desired index
void saveOrbit(Matrix& po, double T, double err, const char *dir, int indx){
  char fname[200];
  sprintf(fname, "%sorbit%d.dat", dir, indx);
  po.output(fname);
  sprintf(fname, "%sT%d.dat", dir, indx);
  Matrix tempm(1,1);
  tempm(0,0) = T;
  tempm.output(fname);
  int n = po.getn();
  sprintf(fname, "%sn%d.dat",dir,indx);
  tempm(0,0) = n;
  tempm.output(fname);
  sprintf(fname,"%serr%d.dat",dir, indx);
  tempm(0,0) = err;
  tempm.output(fname);
}

//symStr(input) = symbol string to approximate
//T(output) = approx period of slice returned
//leftD(output) = derv at left of slice
//rightD(output) = derv at right of slice
//return = 4xn slice allocated inside
//uses library to construct approximate slice
Matrix *returnApproxSlice(char *symStr, 
			  double& T, Vector& leftD, Vector& rightD){
  int len = strlen(symStr);
  Matrix *slice1, *slice2;
  double deltat1, deltat2;
  int i;
  int minlft=LARGENUM;
  for(i=0; i < len; i++){
    //cout<<"slice no: "<<i<<endl;
    int poindx;
    int symindx;
    int lengthOfFit;
    stringbestfit(symStr, i, 
		  poindx, symindx, lengthOfFit);
    if(lengthOfFit < minlft)
      minlft = lengthOfFit;
    int transIndx;
    findTransIndx(*(poList[poindx]), symindx, 
		  transIndx);
    slice2 = returnSlice(*(poList[poindx]), periodList[poindx],
			 transIndx, deltat2);
    if(i==0)
      findLeftDervOfSlice(*(poList[poindx]), periodList[poindx],
			  transIndx, leftD);
    else if(i==(len-1))
      findRightDervOfSlice(*(poList[poindx]), periodList[poindx],
			  transIndx, rightD);
    if(i==0){
      slice1 = slice2;
      deltat1 = deltat2;
    }
    else{
      slice1 = returnMergedSlice(slice1, deltat1, slice2);
      deltat1 = deltat2;
    }
  }
  T = (*slice1)(3,(*slice1).getn()-1)+deltat1;
#ifdef DEBUG
  cout<<"Goodness of fit: "<< minlft<<endl;
#endif
  return slice1;
}

static int FourierFactor = 1;

void setFourierFactor(int ff){
  FourierFactor = ff;
}
  
//symStr = symbol string to be approximated
//T = approximate period
//return = 3xn matrix that is an approx po allocated inside
//uses slices from library and cubic spline interpolation to 
//construct approximate po
//interpolation may or may not use right and left derivatives
Matrix *returnApprox(char *symStr, double& T){
  int len = strlen(symStr);
  Vector leftD(3), rightD(3);
  Matrix *slice = returnApproxSlice(symStr, T, leftD, rightD);
  //int factor = (len/4);
  //  if(!(len%4==0))
  //  factor++;
  int factor = len;
  int rfactor = 1;
  while(rfactor < factor)
    rfactor *= 2;
  //int n = 256*rfactor;
  int n = 64*rfactor;
  n = n*FourierFactor;
  Matrix *poptr = new Matrix(3, n);
  Vector tgrid(n);
  double h = T/n;
  int i;
  for(i=0; i < n; i++)
    tgrid(i) = h*i;
  Matrix slicepo;
  Vector slicet(slice->getn());
  slicepo.shadow(*slice, 0, 2, 0, slice->getn()-1);
  for(i=0; i < slice->getn(); i++)
    slicet(i) = (*slice)(3,i);
  //cspline(slicet, slicepo, tgrid, *poptr);
  spline(slicet, slicepo, leftD, rightD, tgrid, *poptr);
  delete slice;
  return poptr;
}



//deallocates all pos and symbol string in library
void freeLorenz1(){
  int i;
  for(i=0; i < numOfPos; i++){
    delete poList[i];
    free(symList[i]);
  }
}


void findClosestString(Vector& v, int len, char *s){
  int i;
  int index=-1;
  double dist_index = 10000000;
  if((len > strlen(symList[numOfPos-1]))||(len <= 1)){
    cout<<"no string of length "<< len <<" in library"<<endl;
    return;
  }
  for(i=0; i < numOfPos; i++){
    if(strlen(symList[i])!=len)
      continue;
    double disti;
    disti = returnMinDist(*(poList[i]), v);
    if(disti < dist_index){
      index = i;
      dist_index = disti;
    }
  }
  strcpy(s, symList[index]);
}