#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 "lorenz5.hh"
#include <stdlib.h>
static int NUMITNS;
//leaves a long trajectory in out.dat
//this is loaded using DATA/scr1.m
//a snapshot is generated using DATA/scr2.m
//a guess is generated and left in DATA/in.dat DATA/T.dat using DATA/scr3.m
void LorenzA(){
  Random rng(clock());
  Solve ode1(lorenz, lorenzJ, 3, 128*10000, RK7);
  // double T = 150;
  Vector v(3);
  v(0) = rng.randn(); v(1) = rng.randn(); v(2) = rng.randn();
  ode1.setStart(v, 1.0/128.0);
  ode1.advance(128*1000);
  ode1.getSolnAtT(1000.0, v);
  cout<<"x = "<<v(0)<<endl;
  cout<<"y = "<<v(1)<<endl;
  cout<<"z = "<<v(2)<<endl;
  ode1.setStart(v, 1.0/128.0);
  int t = 1000;
  ode1.advance(128*t);
  Vector soln(3*128*t);
  ode1.getSoln(0, 128*t-1, soln);
  Matrix orbit;
  orbit.shadow(soln, 3, 128*t);
  orbit.output("DATA/out.dat");
}
 


//loads a starting guess from DATA/in.dat and DATA/T.dat
//the width of starting guess must be set every run
//tries to improve the starting guess
//saves the result in orbitx.dat nx.dat Tx.dat
//x needs to be set every run 
//all orbit in DATA generated using LorenzA() and LorenzB()
void LorenzB(){
  Matrix po(3, 256*4);//needs to be set for every run
  po.input("DATA/in.dat");
  double T;
  ifstream ifile("DATA/T.dat");
  ifile>>T;
  ifile.close();
  cout<<"T = "<<T<<endl;
  PeriodicOrbit PO(po, T, lorenz, lorenzJ);
  cout.precision(10);
  cout<<"Rel res error: "<<PO.relativeResidualError()<<endl;
  Matrix step(po);
  double wstep;
  double delta=-1;
  size_t i;
  for(i=0; i < 20; i++){//num of itns may need to be set for every run
    PO.smooth();
    cout<<"Rel res error: "<<PO.relativeResidualError()<<endl;
    ComputeTrustStep(PO, delta, step, wstep);
    PO.addStep(step, wstep);
    cout<<"Rel res error: "<<PO.relativeResidualError()<<endl;
    cout<<"delta: "<<delta<<endl;
  }
  PO.getPeriodicOrbit(po, T);
  char s1[100] = "DATA/orbit";
  char s2[100] = "DATA/T";
  char s3[100] = "DATA/n";
  char x[3] = "x";//needs to be set every run
  char s4[10] = ".dat";
  strcat(s1,x); strcat(s1, s4);
  strcat(s2,x); strcat(s2, s4);
  strcat(s3,x); strcat(s3, s4);
  cout<<s1<<endl; cout<<s2<<endl; cout<<s3<<endl;
  po.output(s1);
  Matrix TT(1,1);
  TT(0,0) = T;
  TT.output(s2);
  TT(0,0) = po.getn();
  TT.output(s3);
  findSymString(po, s1);
  cout<<"Symbol dynamics: "<<s1<<endl;
}

//computes all symbol strings of length below a limit
//starting set of pos is loaded from dir1
//data for winners (includeing po) is saved in dir2
//all orbits in DATA1 generated using Lorenz1()
//first run: dir1 = DATA/ and dir2 = DATA1/
//second run: dir1 = DATA1/ and dir2 = DATA1x/
//finally: mv DATA1x DATA1
void  Lorenz1(){
  Clock clk;
  clk.wtic();
  clk.ptic();
  initLorenz1();
  char dir1[100] = "DATA1/";
  loadOrbits(dir1,1,1375);
  char dir2[50] = "/opt/data/divakar/DATA1x/";
  initLorenz2(dir2);
  char s[100] = "A";
  addASloser(s, 1.0);
  s[0] = 'A'; s[1]='A'; s[2] = '\0';
  char ss[100];
  Matrix *po;
  double T;
  double err;
  int num;
  cout<<strlen(s)<<endl;
  while(strlen(s) < 14){
    match2winner(s, num);
    if(num != -1){
      advanceString(s);
      continue;
    }
    match2loser(s, num);
    if(num != -1){
      advanceString(s);
      continue;
    }
    flipANDmatch2loser(s, num);
    if(num != -1){
      addASloser(s, 1.0);
      advanceString(s);
      continue;
    }
    flipANDmatch2winner(s, num);
    if(num != -1){
      addASwinner(s);
      advanceString(s);
      continue;
    }
    po = returnPO(s, T, err);
    cout<<s<<" ERR:"<<err<<endl;
    if(err < 1e-10){
      addASwinner(s, po, T, err);
      findSymString(*po, ss);
      if(!isequal(s,ss)){
	cout<<"computed sym string ss different from s"<<endl;
	abort();
      }
    }
    else
      addASloser(s, err);
    delete po;
    advanceString(s);
  }
  printWinners();
  printLosers();
  cout<<"elapsed time: "<<clk.wtoc()<<" secs"<<endl;
  cout<<"processor time: "<<clk.ptoc()<<" microsecs"<<endl;
  freeLorenz1();
  freeLorenz2();
  FreeSpace();
}

//computes all symbol strings of length below a limit
//starting set of pos is loaded from dir1
//the list of winners and losers is saved in dir2
//but the pos themselves are not stored
//information in DATA2 generated using Lorenz2()
void Lorenz2(){
  Clock clk;
  clk.wtic();
  clk.ptic();
  initLorenz1();
  char dir1[100] = "DATA1/";
  loadOrbits(dir1,1,1375);
  char dir2[50] = "/opt/data/divakar/DATA2/";
  ofstream ofile1("/opt/data/divakar/DATA2/WinnerInfo.dat");
  ofile1.setf(ios::scientific, ios::floatfield);
  ofile1.precision(20);
  ofstream ofile2("/opt/data/divakar/DATA2/Flips.dat");
  initLorenz2(dir2);
  char s[100] = "A";
  addASloser(s, 1.0);
  s[0] = 'A'; s[1]='A'; s[2] = '\0';
  char ss[100];
  Matrix *po;
  double T;
  double err;
  int num;
  int i=0;
  while(strlen(s) <= 20){
   i++;
   if(i%1000 == 10){
     cout<<"printing winners and losers ..."<<endl;
     printWinners();
     printLosers();
   }
   match2winner(s, num);
   if(num != -1){
     advanceString(s);
     continue;
   }
   match2loser(s, num);
   if(num != -1){
     advanceString(s);
     continue;
   }
   flipANDmatch2winner(s, num);
   if(num != -1){
     addASwinnerNOsave(s);
     advanceString(s);
     ofile2<<s<<endl;
     continue;
   }
   flipANDmatch2loser(s, num);
   if(num != -1){
     addASloser(s, 1.0);
     advanceString(s);
     continue;
   }
   po = returnPO(s, T, err);
   if(err < 1e-10){
     addASwinnerNOsave(s);
     findSymString(*po, ss);
     if(!isequal(s,ss)){
       cout<<"computed sym string ss different from s"<<endl;
       abort();
     }
     cout<<s<<endl;
     ofile1<<s<<endl;
     ofile1<<T<<endl;
     ofile1<<err<<endl;
     ofile1<<returnLyap(*po, T)<<endl<<endl;
   }
   else
     addASloser(s, err);
   delete po;
   advanceString(s);
  }
  printWinners();
  printLosers();
  cout<<"Elapsed time: "<<clk.wtoc()/86400.0<<" days"<<endl;
  cout<<"Processor time: "<<clk.ptoc()<<" secs"<<endl;
  freeLorenz1();
  freeLorenz2();
  FreeSpace();
  ofile1.close();
  ofile2.close();
}


void advanceSymmString(char *str){
  int len = strlen(str);
  char *tempstr = (char *)malloc((len+3)*sizeof(char));
  tempstr[0] = 'A';
  strcpy(tempstr+1, str);
  tempstr[len+1] = 'B';
  tempstr[len+2] = '\0';
  strcpy(str, tempstr);
  free(tempstr);
}


  
//generates orbits of the form A^nB^n
//data in DATA3 generated using Lorenz3()
//default shooting interval needs to be set to 256
void Lorenz3(){
  Clock clk;
  clk.wtic();
  char dir1[50] = "DATA1\\";
  char dir2[50] = "DATA3\\";
  initLorenz1();
  loadOrbits(dir1, 1, 1375);
  int curindx = 24;
  loadOrbits(dir2, 1, curindx);
  char s[200] = "AB";
  while(strlen(s) <= 2*curindx)
    advanceSymmString(s);
  Matrix *po;
  double T;
  double err;
  int indx = curindx+1;
  int ffactor = 2;
  setFourierFactor(ffactor);//see note about shooting interval above
  int sintvl = 64*4;
  cout<<"sintvl = "<<sintvl<<endl;
  cout<<"fourier factor = "<<ffactor<<endl;
  char ss[200];
  while(strlen(s) <= 2*(curindx+1)){
    cout<<"strlen: "<<strlen(s)<<endl;
    po = returnPO(s, T, err, sintvl);
    // saveOrbit(*po, T, err, dir2, indx);
    findSymString(*po, ss);
    cout<<"The strings ss and s are: "<<isequal(s,ss)<<endl;
    cout<<endl;
    if(err < 1e-10)
      loadOrbits(dir2, indx, indx);
    else
      cout<<"orbit #"<<indx<<" is inaccurate"<<endl;
    indx++;
    advanceSymmString(s);
    delete po;
  }
  freeLorenz1();
  FreeSpace();
  cout<<"Elapsed time: "<<clk.wtoc()<<" secs"<<endl;
}

//saves distances between pos in DATA3 
void Lorenz3b(){
  char dir[200] = "DATA3/";
  saveLyap(dir, 1, 24);
  ofstream ofile("DATA3/DISTANCES");
  int i;
  for(i=1; i <= 23; i++){
    Matrix *po1;
    Matrix *po2;
    double T;
    double err;
    po1 = loadANDreturnOrbit(dir, i, T, err);
    po2 = loadANDreturnOrbit(dir, i+1, T, err);
    double dist = returnMinDist(*po1, *po2);
    ofile<<"distance between ("<<i<<","<<i+1<<"): "<<dist<<endl;
    delete po1;
    delete po2;
  }
  FreeSpace();
}

//saves distances between pos and the saddle at zero in DATA3
void Lorenz3c(){
  char dir[200] = "DATA3/";
  //ofstream ofile("DATA3/DistFromZero.dat");
  //ofile.precision(20);
  //ofile.setf(ios::scientific, ios::floatfield);
  int i;
  Vector zero(3);
  zero(0) = 0; zero(1) = 0; zero(2) = 0;
  for(i=25; i <= 25; i++){
    Matrix *po;
    double T;
    double err;
    po = loadANDreturnOrbit(dir, i, T, err);
    double dist = returnMinDist(*po, zero);
    cout.setf(ios::scientific, ios::floatfield);
    cout.precision(20);
    cout<<dist<<endl;
    //ofile<<dist<<endl;
    delete po;
  }
  //ofile.close();
}
 
//tries very hard to compute orbits of the form A^nB^n after
//Lorenz3() has failed
void  Lorenz3d(int indx){
  char dir[50] = "DATA3/";
  double T;
  double err;
  Matrix *po = loadANDreturnOrbit(dir, indx, T, err);
  PeriodicOrbit PO(*po, T, lorenz, lorenzJ);
  PO.setSintvl(256);
  double delta = -1;
  Matrix step(po->getm(), po->getn());
  double wstep;
  int i;
  for(i=0; i < 10000; i++){
    PO.smooth();
    err = PO.relativeResidualError();
    cout<<"RELATIVE ERROR: "<<err<<endl;
#ifdef DEBUG    
    cout<<"delta: "<<delta<<endl;
#endif
    if(err < 1e-13)
      break;
    ComputeTrustStep(PO, delta, step, wstep);
    PO.addStep(step, wstep);
  }
  PO.getPeriodicOrbit(*po, T);
  saveOrbit(*po, T, err, dir, 100*indx+indx);
  delete po;
} 

double findAmBn(int m, int n){
  char s[300];
  int i;
  for(i=0; i < m; i++)
    s[i] = 'A';
  for(i=m; i < m+n; i++)
    s[i] = 'B';
  s[m+n] = '\0';
  setFourierFactor(4);
  NUMITNS = 20000;
  Matrix *po;
  double T;
  double err;
  po = returnPO(s, T, err, 2048);
  int index;
  index = 100*m + n;
  char dir[50] = "DATA3/";
  saveOrbit(*po, T, err, dir, index);
  delete po;
  return err;
}

void Lorenz3e(){
  char dir1[50] = "DATA1/";
  loadOrbits(dir1, 1, 1375);
  char dir2[50] = "DATA3/";
  loadOrbits(dir2, 1, 25);
  char dir3[50] = "DATA4/";
  loadOrbits(dir3, 1, 23);
  int m = 25;
  double err;
  for(int n=24; n>=1; n--){
    err = findAmBn(m,n);
    if(err>1e-8)
      break;
  }
}

//saves distances between pos and the saddle at zero in DATA3
void Lorenz3f(){
  char dir[200] = "DATA3/";
  ofstream ofile("DATA3/DistFromZeroX.dat");
  ofile.precision(20);
  ofile.setf(ios::scientific, ios::floatfield);
  int i;
  Vector zero(3);
  zero(0) = 0; zero(1) = 0; zero(2) = 0;
  for(i=2507; i <= 2524; i++){
    Matrix *po;
    double T;
    double err;
    po = loadANDreturnOrbit(dir, i, T, err);
    double dist = returnMinDist(*po, zero);
    cout<<dist<<endl;
    ofile<<dist<<endl;
    delete po;
  }
  double T, err;
  Matrix *po = loadANDreturnOrbit(dir, 25, T, err);
  double dist = returnMinDist(*po, zero);
  ofile<<dist<<endl;
  ofile.close();
}
 
//attempts to find the symbol sequence of zero
void Lorenz3g(){
  char dir1[50] = "DATA1/";
  loadOrbits(dir1, 1, 375);
  char dir2[50] = "DATA3/";
  loadOrbits(dir2, 1, 25);
  char dir3[50] = "DATA4/";
  loadOrbits(dir3, 1, 23);
  char s[50] = "AAAAAAAAAAAAAAAAAAAAAAAAABBBBBB";
  double error = .01;
  setRepeatMax(10);
  setFourierFactor(2);
  ZeroFit zf(s, error);
  NUMITNS = 1000;
  for(int i = 0; i < 10; i++){
    zf.updateFit();
    zf.getBestFit(s, error);
    cout<<s<<endl;
    cout<<error<<endl;
    cout<<"dist: "<<zf.returnDistance()<<endl;
  }
}

//generates orbits of the form BA^n
void Lorenz4(){
  Clock clk;
  clk.wtic();
  char dir1[50] = "DATA1/";
  char dir2[50] = "DATA4/";
  initLorenz1();
  loadOrbits(dir1, 1, 1375);
  int indx = 23;
  loadOrbits(dir2, 1, indx-1);
  char s[100];
  s[0] = 'B';
  int i;
  for(i=1; i <=indx; i++)
    s[i] = 'A';
  s[indx+1] = '\0';
  Matrix *po;
  double T;
  double err;
  char ss[200];
  while(strlen(s) <= 25){
    cout<<"strlen: "<<strlen(s)<<endl;
    po = returnPO(s, T, err);
    saveOrbit(*po, T, err, dir2, indx);
    findSymString(*po, ss);
    cout<<"The strings ss and s are: "<<isequal(s,ss)<<endl;
    if(err < 1e-10)
      loadOrbits(dir2, indx, indx);
    else
      cout<<"orbit #"<<indx<<" is inaccurate"<<endl;
    indx++;
    int len = strlen(s);//advance string
    s[len] = 'A';
    s[len+1] = '\0';
    delete po;
  }
  freeLorenz1();
  FreeSpace();
  cout<<"Elapsed time: "<<clk.wtoc()<<" secs"<<endl;
}  

//saves distances between pos in DATA4 and their lyap exponents
void Lorenz4b(){
  char dir[200] = "DATA4/";
  saveLyap(dir, 1, 23);
  ofstream ofile("DATA4/DISTANCES");
  int i;
  for(i=1; i <= 22; i++){
    Matrix *po1;
    Matrix *po2;
    double T;
    double err;
    po1 = loadANDreturnOrbit(dir, i, T, err);
    po2 = loadANDreturnOrbit(dir, i+1, T, err);
    double dist = returnMinDist(*po1, *po2);
    ofile<<"distance between ("<<i<<","<<i+1<<"): "<<dist<<endl;
    delete po1;
    delete po2;
  }
  FreeSpace();
}

//saves distances between pos and the saddle at zero in DATA4
void Lorenz4c(){
  char dir[200] = "DATA4/";
  ofstream ofile("DATA4/DistFromZero.dat");
  ofile.precision(20);
  ofile.setf(ios::scientific, ios::floatfield);
  int i;
  Vector zero(3);
  zero(0) = 0; zero(1) = 0; zero(2) = 0;
  for(i=1; i <= 23; i++){
    Matrix *po;
    double T;
    double err;
    po = loadANDreturnOrbit(dir, i, T, err);
    double dist = returnMinDist(*po, zero);
    cout<<dist<<endl;
    ofile<<dist<<endl;
    delete po;
  }
  ofile.close();
  
}

//finds the po A^{24}B
void Lorenz4d(){
  Clock clk;
  clk.wtic();
  char dir1[50] = "DATA1\\";
  loadOrbits(dir1, 1, 1375);
  char dir2[50] = "DATA3\\";
  loadOrbits(dir2, 1, 25);
  char dir3[50] = "DATA4\\";
  loadOrbits(dir3, 1, 23);
  int ffactor = 1;
  setFourierFactor(ffactor);
  char s[50] = "AAAAAAAAAAAAAAAAAAAAAAAAB";
  Matrix *po;
  //NUMITNS = 10000;
  double err, T;
  int sintvl = 16;
  cout<<"sintvl = "<<sintvl<<endl;
  cout<<"fourier factor = "<<ffactor<<endl;
  po = returnPO(s, T, err, sintvl);
  //saveOrbit(*po, T, err, dir3, 24);
  cout<<"Elapsed time: "<<clk.wtoc()<<" secs"<<endl;
}

//makes a guess for A24B using A25B25
void Lorenz4ex(){
  char s[50] = "AAAAAAAAAAAAAAAAAAAAAAAAB";
  Matrix orbit25(3,4096);
  Matrix orbit24(3, 2048);
  orbit25.input("DATA3\\orbit25.dat");
  for(int i=0; i < 1869; i++)
    orbit24.getColumn(i) = orbit25.getColumn(i);
  for(int i=3918; i < 4096; i++)
    orbit24.getColumn(1869+(i-3918)) = orbit25.getColumn(i);
  orbit24.output("DATA4\\orbit24.dat"); 
}

//tries very hard to compute the orbit of  A^24B after
//Lorenz4d() has failed
void  Lorenz4e(){
  char s[50] = "AAAAAAAAAAAAAAAAAAAAAAAAB";
  int indx = 24;
  char dir[50] = "DATA4/";
  double T;
  double err;
  Matrix *po = loadANDreturnOrbit(dir, indx, T, err);
  PeriodicOrbit PO(*po, T, lorenz, lorenzJ);
  PO.setSintvl(64);
  double delta = -1;
  Matrix step(po->getm(), po->getn());
  double wstep;
  int i;
  for(i=0; i < 4000; i++){
    PO.smooth();
    err = PO.relativeResidualError();
    cout<<"RELATIVE ERROR: "<<err<<endl;
#ifdef DEBUG    
    cout<<"delta: "<<delta<<endl;
#endif
    if(err < 1e-13)
      break;
    ComputeTrustStep(PO, delta, step, wstep);
    PO.addStep(step, wstep);
  }
  PO.getPeriodicOrbit(*po, T);
  checkSymString(*po, s);
  saveOrbit(*po, T, err, dir, indx);
  delete po;
} 


void randomString(char *s, int len){
  static Random rng;
  s[len] = '\0';
  int i;
  for(i=0; i < len; i++)
    s[i] = (rng.randn()>0)?'A':'B';
}

//generates orbits with random symbol sequences
//saves in dir2/ with INDX
//finds distance between po and po with one last symbol missing
//saves data about distances in ofile
void Lorenz5(){
  char s[200];
  int i;
  char dir1[100] = "DATA1/";
  char dir2[100] = "DATA5/";
  ofstream ofile("DATA5/SYMDIST.DAT");
  loadOrbits(dir1, 1, 1375);
  Matrix *po;
  double T;
  double err;
  for(i=5; i <= 50; i++){
    int LENGTH=i;
    int INDX=i;
    randomString(s, LENGTH);
    cout<<s<<endl;
    po = returnPO(s, T, err);
    double dist=100.0;
    if(err < 1e-10){
      char ss[200];
      findSymString(*po, ss);
      if(!(isequal(s,ss))){
	cout<<"s and ss are not equal!"<<endl;
	abort();
      }
      saveOrbit(*po, T, err, dir2, INDX);
      {//find nearby po
	int len = strlen(s);
	strcpy(ss, s);
	ss[len-1] = '\0';
	Matrix *potwin = returnPO(ss, T, err);
	if(err > 1e-10)
	  dist = 100.0;
	else
	  dist = returnMinDist(*po, *potwin);
	cout<<"Distance is: "<<dist<<endl;
	delete potwin;
      }
    }
    ofile<<"orbit #"<<i<<": "<<s<<"  Dist:"<<dist<<endl;
    delete po;
  }
  ofile.close();
  freeLorenz1();
  FreeSpace();
}

//generate and save random po with symbol string of length 347
void Lorenz5a(){
  char s[500];
  int i=347;
  char dir1[100] = "DATA1/";
  char dir2[100] = "DATA5/";
  loadOrbits(dir1, 1, 1375);
  Matrix *po;
  double T;
  double err;
  int LENGTH=i;
  int INDX=i;
  randomString(s, LENGTH);
  cout<<s<<endl;
  po = returnPO(s, T, err);
  cout<<"saving orbit ..."<<endl;
  saveOrbit(*po, T, err, dir2, INDX);
  delete po;
  saveLyap(dir2, INDX, INDX);
}
  
//generate and save ABA^2B^2...A^nB^n
void Lorenz5b(){
  char s[1000] = "AB";
  char ss[100] = "AABB";
  while(strlen(ss) <= 30){
    strcat(s, ss);
    advanceSymmString(ss);
  }
  int i = strlen(s);
  char dir1[100] = "DATA1/";
  char dir2[100] = "DATA5/";
  loadOrbits(dir1, 1, 1375);
  Matrix *po;
  double T;
  double err;
  int LENGTH=i;
  int INDX=i;
  cout<<s<<endl;
  po = returnPO(s, T, err);
  cout<<"saving orbit ..."<<endl;
  saveOrbit(*po, T, err, dir2, INDX);
  delete po;
  saveLyap(dir2, INDX, INDX);
}


//prints a random point close to the Lorenz attractor
void spitRandomPoint(int N, int indx){
  //timeval ctime;
  //gettimeofday(&ctime, NULL);
  //static Random rng(ctime.tv_sec);
  Random rng;
  Solve ode1(lorenz, lorenzJ, 3, N*120, RK7);
  Vector v(3);
  v(0) = rng.randn(15.0); v(1) = rng.randn(15.0); v(2) = rng.randn();
  v(0) += 20; v(1) += 20; v(2) += 27;
  ode1.setStart(v, 1.0/N);
  ode1.advance(N*100);
  int i;
  for(i=0; i < 100; i++){ 
    if(i%10==0)
      cout<<i<<endl;
    ode1.getSoln(N*99, v);
    ode1.setStart(v, 1.0/N);
    ode1.advance(N*100);
  }
  Vector w(3);
  char  fname[100];
  assert(indx > 0);
  sprintf(fname, "DATA6/pt%d.dat", indx);
  ofstream ofile(fname);
  sprintf(fname, "DATA6/N%d.dat", indx);
  ofstream ofilex(fname);
  ofilex<<"reciprocal of timestep is "<<N<<endl;
  ofilex.close();
  ofile.setf(ios::scientific, ios::floatfield);
  ofile.precision(20);
  for(i=1; i < N*99; i++){
    ode1.getSoln(i, v);
    ode1.getSoln(i+1, w);
    if((v(2) < 27) && (w(2) > 27) && (v(1) < 0)){
      ofile<<v(0)<<endl;
      ofile<<v(1)<<endl;
      ofile<<v(2)<<endl;
      break;
    }
  }
  ofile.close();
}

void spitRandomPoints(){
  //spitRandomPoint(128,1);
  //spitRandomPoint(128,2);
  //spitRandomPoint(512, 3);
  //spitRandomPoint(512, 4);
  //spitRandomPoint(2096, 5);
  //spitRandomPoint(2096, 6);
  //spitRandomPoint(4192, 7);
  //spitRandomPoint(4192, 8);
  //spitRandomPoint(2096, 9);
  //spitRandomPoint(2096, 10);
  //spitRandomPoint(2096, 11);
}

//finds closest pos to random points on Lorenz attractors
//spitRandomPoints() must be run before this routine
void Lorenz6(){
  const int maxlen = 32;
  char dir1[50] = "DATA1/";
  loadOrbits(dir1, 1, 1375);
  int indx; 
  for(indx = 1; indx <= 11; indx++){
    Vector v(3);
    char fname[100];
    sprintf(fname, "DATA6/pt%d.dat", indx);
    v.input(fname);
    int i;
    sprintf(fname, "DATA6/trace%d.dat", indx);
    ofstream trfile(fname);
    char *slist[45];
    for(i=0; i < 45; i++)
      slist[i] = (char *)malloc(200*sizeof(char));
    double dist[45];
    for(i=2; i <=13; i++){
      matchPoint2String(v, i, slist[i], dist[i]);
      cout<<"strng: "<<slist[i]<<"("<<i<<")"<<"  dist:"<<dist[i]<<endl;
      trfile<<"strng: "<<slist[i]<<"("<<i<<")"<<"  dist:"<<dist[i]<<endl;
    }
    for(i=14; i <= maxlen; i++){
      if (i <= 24)
	updateFit(v, slist[i-1], slist[i-2], slist[i-3],
		  slist[i], dist[i]);
      if((24 < i) && (i <= maxlen))
	updateFit(v, slist[i-1], slist[i], dist[i]);
      rotateString2Front(v, slist[i]);
      cout<<"strng: "<<slist[i]<<"("<<i<<")"<<"  dist:"<<dist[i]<<endl;
      trfile<<"strng: "<<slist[i]<<"("<<i<<")"<<"  dist:"<<dist[i]<<endl;
    }
    double mindist = dist[2];
    int minindx = 2;
    for(i=2; i <= maxlen; i++)
      if(dist[i] < mindist){
	mindist = dist[i];
	minindx = i;
      }
    
    Matrix *po;
    double T;
    double err;
    po = returnPO(slist[minindx], T, err);
    //char dir[100];
    //sprintf(dir, "DATA6/");
    //saveOrbit(*po, T, err, dir, indx);
    if(err > 1e-10)
      trfile<<"error in po too high"<<endl;
    char symstr[100];
    findSymString(*po, symstr);
    if(!isequal(symstr, slist[minindx]))
      trfile<<symstr<<" does not match "<<slist[minindx]<<endl;
    else
      trfile<<"min distance is "<<mindist<<endl;
    //saveLyap(dir, indx, indx);
    delete po;
    
    trfile.close();
  }  

  freeLorenz1();
  FreeSpace();
}



void fitPoint(int index, int N, char* head, char* tail, int numitns){
  static int LOADFLAG = 0;
  if(LOADFLAG==0){
    char dir[50] = "DATA1/";
    loadOrbits(dir, 1, 1375);
    LOADFLAG = 1;
  }
  Vector pt(3);
  char fname[50];
  sprintf(fname, "DATA6/pt%d.dat", index);
  pt.input(fname);
  advancePoint2Section(pt, 1.0/(double)N);
  sprintf(fname, "DATA6/TRACE%d.dat", index);
  ofstream ofile(fname);
  FwdFit fitpt(pt, head, tail);
  fitpt.printfit();
  fitpt.printfit(ofile);
  fitpt.printallleaves();
  fitpt.printallleaves(ofile);
  int i;
  for(i=0; i < numitns; i++){
    fitpt.updatefit();
    fitpt.printfit();
    fitpt.printallleaves();
    fitpt.printallleaves(ofile);
  }
  char s[100];
  double dist;
  fitpt.getBestFit(s, dist);
  ofile<<endl<<endl;
  ofile<<"Best fit:"<<s<<" Dist:"<<dist<<endl;
  ofile.close();
  if(dist < 1e-2){
    Matrix* po;
    double T;
    double err;
    po = returnPO(s, T, err);
    char dirout[50] = "DATA6/";
    saveOrbit(*po, T, err, dirout, index); 
  }
}

void fitPoint(int index, int N, char* tail, int numitns){
  char emptystring[10] = "";
  fitPoint(index, N, emptystring, tail, numitns);
}

//rerun
void Point1(){
  char tail[10] = "AA"; 
  char head[50] = "ABBBBABABBBBBBAABABBABBBBBBAB";
  fitPoint(1, 128, head, tail, 3);
}

//rerun
void Point2(){
  char tail[10] = "AA";
  char head[50] = "AAABAAABAAAABABBAAAAAAABBBAABBBA";
  fitPoint(2, 128, head, tail, 1);
}

//rerun
void Point3(){
  char tail[10] = "AA";
  char head[50] = "ABAAABABAAAABBABAABBABBAABBA";
  fitPoint(3, 512, head, tail, 6);
}

//ok
void Point4(){
  char tail[10] = "BA";
  char head[50] = "AAAAAAAAAABAAAAABABAABB";
  setRepeatMax(14);
  fitPoint(4, 512, head, tail, 10);
  setRepeatMax(10);
}
 
//rerun
void Point5(){
  char tail[10] = "AA";
  char head[50] = "AABBBBBBBBBAAABBBABAABBBABBABBBBABBAABA";
  fitPoint(5, 2096, head, tail, 3);
}

//rerun
void Point6(){
  char tail[10] = "AA";
  char head[50] = "ABAABBAAABBABBBBBBBBBBBBB";
  setRepeatMax(15);
  fitPoint(6, 2096, head, tail, 9);
  setRepeatMax(10);
}

//ok 
void Point7(){
  char tail[10] = "ABB";
  fitPoint(7, 4192, tail, 10);
}

//rerun
void Point8(){
  char tail[10] = "BB";
  char head[50] = "AAABABABAAAABBABBBABBBAAA";
  fitPoint(8, 4192, head, tail, 5);
}

//rerun
void Point9(){
  char tail[10] = "BAB";
  char head[50] = "AAABBBAAAAABAAABBBBBBBBBBABABAA";
  setRepeatMax(12);
  fitPoint(9, 2096, head, tail, 4);
  setRepeatMax(10);
}

//ok
void Point10(){
  char tail[10] = "AB";
  char head[50] = "ABAAAAAAAAABBBBAAABABBAAABAAAAB";
  fitPoint(10, 2096, head, tail, 6);
}

//rerun
void Point11(){
  char tail[10] = "BB";
  char head[50] = "ABBABAAABAAAAAAAABABBAABAABABA";
  fitPoint(11, 2096, head, tail, 2);
}
  

void findithstring(char *head, int indx, int N, char *tail, char *s){
  char ss[30];
  int i;
  for(i=0; i < N; i++){
    ss[i] = (indx%2==0)? 'A':'B';
    indx = indx/2;
  }
  ss[N] = '\0';
  strcpy(s, head);
  strcat(s, ss);
  strcat(s, tail);
}


void savefirstpoints(char *head, char *tail, int N, char *fname){
  char s[100];
  char fn1[50];
  char fn2[50];
  sprintf(fn1, "DATA7/%s.dat",fname);
  sprintf(fn2, "DATA7/%s.str",fname);
  ofstream ofile1(fn1);
  ofile1.precision(17);
  ofstream ofile2(fn2);
  int NN = pow(2.0, N);
  int i;
  for(i=0; i < NN; i++){
    findithstring(head, i, N, tail, s);
    cout<<s<<endl;
    Vector pt(3);
    findFirstPoint(s, pt);
    ofile1<<pt(0)<<" "<<pt(1)<<" "<<pt(2)<<endl;
    ofile2<<s<<endl;
  }
  ofile1.close();
  ofile2.close();
}

void findLeaf(char *head, char *tail, int code = 1){
  char fname[100];
  if(code != 1)
    sprintf(fname,"%sx%s", head, tail);
  else
    sprintf(fname, "%s", tail);
  int N;
  N = (code==1)?7:5;
  savefirstpoints(head, tail, N, fname);
}

//must load orbits into memory before calling these routines
void Lorenz7a(){
  char head[30] = "ABABABABAB";
  char tail1[10] = "A";
  char tail2[10] = "B";
  findLeaf(head, tail1);
  findLeaf(head, tail2);
}

void Lorenz7b(){
  char head[40] = "ABABABABABABABABABABABAB";
  char tail1[10] = "AA";
  char tail2[10] = "BA";
  findLeaf(head, tail1);
  findLeaf(head, tail2);
}
 


void Lorenz7c(){
  char dir[50];
  sprintf(dir, "DATA1/");
  loadOrbits(dir, 1, 1375);
  sprintf(dir, "DATA3/");
  loadOrbits(dir, 1, 25);
  sprintf(dir, "DATA4/");
  loadOrbits(dir, 1, 23);
  char head[30] = "AAAAAAAAAAAAAAAAAAAAAAAB";
  char tail1[10] = "BBB";
  char fname[50];
  sprintf(fname, "BBX");
  setFourierFactor(4);
  // savefirstpoints(head, tail1, 5, fname);
  char tail2[10] = "BBAB";
  sprintf(fname, "ABX");
  savefirstpoints(head, tail2, 4, fname);
} 

//saves data for intersection with z=27 section
void Lorenz7d(){
  Vector pt(3);
  pt(0) = 0.001; pt(1) = 0.001; pt(2) = 0.001;
  double h = .01;
  ofstream ofile("DATA7/section.dat");
  ofile.precision(8);
  for(int i=0; i < 512*16; i++){
    advancePoint2ANYSection(pt, h);
    ofile<<pt(0)<<" "<<pt(1)<<" "<<pt(2)<<endl;
    cout<<pt(0)<<" "<<pt(1)<<" "<<pt(2)<<endl;
    advancePoint(pt, h);
  }
  ofile.close();
}

void Lorenz10(){
  char fname[100] = "DATA10/AnBnAngles.dat";
  ofstream ofile(fname);
  printAngles(ofile);
}


int main(){
  //Lorenz3();
  //Lorenz4d();
  Lorenz4ex();
  Lorenz4e();
}