#include "hdim.hh"
#include <iostream.h>
#include <fstream.h>
#include "Vector.hh"
#include "Matrix.hh"
#include <string.h>
#include <stdio.h>
#include "lorenz1.hh"
#include <math.h>
#include <algorithm>

static int LARGENUM=100000;
/* 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;
}

//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;
}



/*** Implementation of poData and poList ***/
void poData::readData(ifstream& ifile){
  ifile>>s;
  ifile>>T;
  ifile>>err;
  ifile>>lyap;
  if (ifile.eof()){
    T = -1;
    err = -1;
    lyap = -1;
  }
}


poList::poList(char *fname){
  list = new poData[200000];
  ifstream ifile(fname);
  count = 0;
  while(!ifile.fail()){
    (list[count]).readData(ifile);
    if((list[count]).dataOK())
      count++;
    else
      break;
  }
  //cout<<"COUNT: "<<count<<endl;
  ifile.close();
}
 
void poList::getSymStr(int i, char *ss){
  if((0<=i) && (i < count))
    list[i].getSymStr(ss);
}

double poList::returnT(int i){
  if((0<=i) && (i < count))
    return list[i].returnT();
}

double poList::returnErr(int i){
  if((0<=i) && (i < count))
    return list[i].returnErr();
}

double poList::returnLyap(int i){
  if((0<=i) && (i < count))
    return list[i].returnLyap();
}

/*** end implementation of poData and poList ***/

void makeLambdaList(int m, poList& poinfo, 
		    double *llist, int * counts, int& num){
  int maxnum = poinfo.returnCount();
  int i;
  num = 0;
  for(i=0; i < maxnum; i++){
    char s[100]; int slen; int quot;
    double T, lyap;
    poinfo.getSymStr(i, s);
    slen = strlen(s);
    if(slen > m)
      break;
    if(m%slen != 0)
      continue;
    else
      quot = m/slen;
    T = poinfo.returnT(i);
    lyap = poinfo.returnLyap(i);
    lyap = -11.0 - 8.0/3.0 - lyap;
    llist[num] = pow(2.718281828459046e+00, T*lyap*quot);
    char ss[100];
    int j;
    for(j=0; j < slen; j++)
      ss[j] = (s[j]=='A')?'B':'A';
    ss[slen] = '\0';
    if (isequal(s, ss))
      counts[num] = slen;
    else
      counts[num] = 2*slen;
    num++;
  }
}
		    
void makeExLambdaList(int m, poList& poinfo, 
		    double *ulist, double *slist, int * counts, int& num){
  int maxnum = poinfo.returnCount();
  int i;
  num = 0;
  for(i=0; i < maxnum; i++){
    char s[100]; int slen; int quot;
    double T, lyap;
    poinfo.getSymStr(i, s);
    slen = strlen(s);
    if(slen > m)
      break;
    if(m%slen != 0)
      continue;
    else
      quot = m/slen;
    T = poinfo.returnT(i);
    lyap = poinfo.returnLyap(i);
    ulist[num] = pow(2.718281828459045e+00, T*lyap*quot);
    lyap = -11.0 - 8.0/3.0 - lyap;
    slist[num] = pow(2.718281828459045e+00, T*lyap*quot);
    char ss[100];
    int j;
    for(j=0; j < slen; j++)
      ss[j] = (s[j]=='A')?'B':'A';
    ss[slen] = '\0';
    if (isequal(s, ss))
      counts[num] = slen;
    else
      counts[num] = 2*slen;
    num++;
  }
}

void makeEXXLambdaList(int m, poList& poinfo, 
		       double *ulist, double *slist, double *plist,
		       int * counts, int& num){
  int maxnum = poinfo.returnCount();
  int i;
  num = 0;
  for(i=0; i < maxnum; i++){
    char s[100]; int slen; int quot;
    double T, lyap;
    poinfo.getSymStr(i, s);
    slen = strlen(s);
    if(slen > m)
      break;
    if(m%slen != 0)
      continue;
    else
      quot = m/slen;
    T = poinfo.returnT(i);
    lyap = poinfo.returnLyap(i);
    ulist[num] = pow(2.718281828459045e+00, T*lyap*quot);
    lyap = -11.0 - 8.0/3.0 - lyap;
    slist[num] = pow(2.718281828459045e+00, T*lyap*quot);
    plist[num] = T*quot;
    char ss[100];
    int j;
    for(j=0; j < slen; j++)
      ss[j] = (s[j]=='A')?'B':'A';
    ss[slen] = '\0';
    if (isequal(s, ss))
      counts[num] = slen;
    else
      counts[num] = 2*slen;
    num++;
  }
}


double hdimsum(double g, double *llist, int *counts, int num){
  double sum = 0;
  int i;
  for(i=0; i < num; i++){
    sum += counts[i]*pow(llist[i], g);
  }
  return sum;
}

double hdim(poList& poinfo, int m){
  int counts[55000];
  double llist[55000];
  int num=0;
  makeLambdaList(m, poinfo, llist, counts, num);
#ifdef DEBUG
  int totalnum = 0;
  for(int i=0; i < num; i++)
    totalnum += counts[i];
  //totalnum must be 2^m-2
  cout<<"m:"<<m<<"  counts:"<<totalnum<<endl;
#endif
  double l = .055;
  double r = .065;
  double ls = hdimsum(l, llist, counts, num);
  double rs = hdimsum(r, llist, counts, num);
  while(fabs(l-r) > 1e-8){
    double m = (l+r)/2;
    double ms = hdimsum(m, llist, counts, num);
    if (ms > 1)
      l = m;
    else
      r = m;
  }
  return (l+r)/2;
}

double trace(poList& poinfo, int m){
  int counts[100000];
  double ulist[100000];
  double slist[100000];
  int num=0;
  makeExLambdaList(m, poinfo, ulist, slist, counts, num);
#ifdef DEBUG
  int totalnum = 0;
  for(int i=0; i < num; i++)
    totalnum += counts[i];
  //totalnum must be 2^m-2
  cout<<"m:"<<m<<"  counts:"<<totalnum<<endl;
#endif
  double ans = 0;
  for(int i=0; i < num; i++)
    ans += counts[i]/fabs((1-slist[i])*(1-ulist[i]));
  return ans;
}

double traceXX(poList& poinfo, double s, double beta, int m){
  int counts[100000];
  double ulist[100000];
  double slist[100000];
  double plist[100000];
  int num=0;
  makeEXXLambdaList(m, poinfo, ulist, slist, plist, counts, num);
  assert(num<100000);
  double ans = 0;
  for(int i=0; i < num; i++){
    double fac=((pow(ulist[i],2-beta)+pow(slist[i],2-beta))/fabs(ulist[i]-slist[i]));
    ans += (counts[i]/fabs((1-slist[i])*(1-ulist[i])))
      *(pow(2.718281828459045e+00, s*plist[i])*fac
	/m);
  }
  return -ans;
}


void polyInit(int degree, double *poly){
  int i;
  for(i=0; i<=degree; i++)
    poly[i] = 0.0;
}

//poly1 = poly1*poly2
void polyMult(int degree, double *poly1, double *poly2){
  double *poly = new double[degree+1];
  polyInit(degree, poly);
  int i,j, k;
  for(i=0; i<=degree; i++)
    for(j=i; j>=0; j--){
      k = i-j;
      poly[i] += poly1[j]*poly2[k];
    }
  for(i=0; i<=degree; i++)
    poly1[i] = poly[i];
  delete[] poly;
}
    
//poly1 = poly1+poly2
void polyAdd(int degree, double *poly1, double *poly2){
  int i;
  for(i=0; i<=degree; i++)
    poly1[i] += poly2[i];
}

void polyExp(int degree, double *poly){
  double *ppoly = new double[degree+1];
  double *spoly = new double[degree+1];
  polyInit(degree, ppoly);
  polyInit(degree, spoly);
  ppoly[0] = 1;
  spoly[0] = 1;
  int i;
  for(i=1; i<= degree; i++){
    polyMult(degree, ppoly, poly);
    int j;
    for(j=0; j<=degree; j++)
      ppoly[j] /= i;
    polyAdd(degree, spoly, ppoly);
  }
  for(i=0; i<= degree; i++)
    poly[i] = spoly[i];
  delete[] spoly;
  delete[] ppoly;
}

void tracePoly(poList& poinfo, double s, double beta, int degree,
	       double *poly){
  poly[0] = 0;
  int i;
  for(i=1; i <= degree; i++)
    poly[i] = traceXX(poinfo, s, beta, i);
}


class fredholm{
private:
  poList *poinfo;
  static int fredDegree;
public:
  fredholm();
  ~fredholm();
  double fredValue(double s, double beta);
  void printPoly(double s, double beta);
};

fredholm::fredDegree = 20;

fredholm::fredholm(){
  char fname[200] = "DATA2/WinnerInfo.dat";
  poinfo = new poList(fname);
}

fredholm::~fredholm(){
  delete poinfo;
}

double fredholm::fredValue(double s, double beta){
  int degree = fredDegree;
  double *poly = new double[degree+1];
  tracePoly(*poinfo, s, beta, degree, poly);
  polyExp(degree, poly);
  double ans = 0.0;
  int i;
  for(i=0; i<= degree; i++)
    ans += poly[i];
  delete[] poly;
  return ans;
}

void fredholm::printPoly(double s, double beta){
  int degree = fredDegree;
  double *poly = new double[degree+1];
  tracePoly(*poinfo, s, beta, degree, poly);
  polyExp(degree, poly);
  int i;
  for(i=0; i<= degree; i++){
    cout<<poly[i]<<" * z^"<<i;
    if(i<degree)
      cout<<" + ";
    cout<<endl;
  }
}

double solveFORs(double beta){
  double dbeta = beta - 1.0;
  double ls = dbeta * 0.88-.0001;
  double rs = dbeta * 0.92+.0001;
  fredholm det;
  double lval = det.fredValue(ls, beta);
  double rval = det.fredValue(rs, beta);
  while(fabs(ls-rs) > 1e-14){
    assert(lval*rval < 0);
    double ms = (ls+rs)/2;
    double mval = det.fredValue(ms, beta);
    cout<<lval<<" "<<mval<<" "<<rval<<endl;
    if(mval*lval > 0){
      ls = ms;
      lval = mval;
    }
    else{
      rs = ms;
      rval = mval;
    }
  }
  return((ls+rs)/2.0);
}

int issymm(char *s){
  char ss[100];
  int j;
  int slen = strlen(s);
  for(j=0; j < slen; j++)
    ss[j] = (s[j]=='A')?'B':'A';
  ss[slen] = '\0';
  return isequal(s, ss);
}

void ParryPollicott(){
  char fname[100] = "DATA2/WinnerInfo.dat";
  poList poinfo(fname);
  int count = poinfo.returnCount();
  double *tlist = new double[3*count];
  char s[50];
  int i;
  int j=0;
  int num = 0;
  for(i=0; i < count; i++){
    tlist[j] = poinfo.returnT(i);
    poinfo.getSymStr(i, s);
    j++; num++;
    if(!issymm(s)){
      tlist[j] = tlist[j-1];
      j++;
      num++;
    }
  }
  std::sort(tlist, tlist+num);
  Vector entropy(num);
  for(i=0; i < num; i++)
    entropy(i) = log(i+1)/tlist[i];
  entropy.output("DBG/entropy.dat");
  double h = 0.745;
  Vector C(num);
  for(i=0; i < num; i++)
    C(i) = (i+1)/pow(EEXP,h*tlist[i])/(h*tlist[i]);
  C.output("DBG/C.dat");
  Vector T(num);
  for(i=0; i < num; i++)
    T(i) = tlist[i];
  T.output("DBG/T.dat");
  delete[] tlist;
}

/*int main(){
  ParryPollicott();
  }*/

int main(){//computes Lyapunov exponents
  fredholm det;
  double s = 0;
  double beta = 1;
  s = solveFORs(beta);
  cout<<"s = "<<s<<endl;
  cout<<det.fredValue(s, beta)<<endl;
  det.printPoly(s,beta);
  double dbeta  = 1e-7;
  beta = 1+dbeta;
  double ss = solveFORs(beta);
  cout<<"s = "<<ss<<endl;
  cout<<det.fredValue(ss, beta)<<endl;
  det.printPoly(ss,beta);
  cout.precision(8);
  cout<<"lyap = "<<(ss-s)/dbeta<<endl;
}
 
 

/*
int main(){//compute escape rate
  char fname[200] = "DATA2/WinnerInfo.dat";
  poList poinfo(fname);
  int i;
  ofstream ofile("temp.math");
  ofile.precision(15);
  for(i=1; i<= 20; i++){
    double tr = trace(poinfo, i);
    cout<<" i:" <<i<<" trace:"<<tr<<endl;
    ofile<<tr<<endl;
  }    
}
*/

/*
int main(){//computes Hausdorff dimension
  char fname[200] = "DATA2/WinnerInfo.dat";
  poList poinfo(fname);
  int i;
  sprintf(fname, "DATA9/HDimension.dat");
  ofstream ofile(fname);
  ofile.setf(ios::scientific, ios::floatfield);
  ofile.precision(8);
  for(i=4; i<=20; i++)
    ofile<<i<<"  "<<2.0+hdim(poinfo, i)<<endl;
  ofile.close();
}   
*/

/*
int main(){//outputs list of prime sequences upto length 20
  char fname[200] = "DATA2/WinnerInfo.dat";
  poList poinfo(fname);
  int maxnum = poinfo.returnCount();
  int i;
  ofstream ofile("PrimeSeq.dat");
  ofile<<"A"<<endl;
  ofile<<"B"<<endl;
  int count = 0;
  for(i=0; i < maxnum; i++){
    char s[50], ss[50];
    poinfo.getSymStr(i, s);
    ofile<<s<<endl;
    count++;
    int slen = strlen(s);
    int j;
    for(j=0; j < slen; j++)
      ss[j] = (s[j]=='A')?'B':'A';
    ss[slen] = '\0';
    if(!isequal(s, ss)){
      ofile<<ss<<endl;
      count++;
    }
  }
  ofile.close();
  cout<<"count: "<<count<<endl;
}
*/