#include "header_cpp.h"
#include "myreadpdb.h"
#include "ematrix.h"
#include "program_options.h"
#include "glob_var.h"
#include "ReadProgramOption.c"
#include "stat.h"


#define ABS(x)   ( ((x)<0) ? (-(x)) : (x))

long sseed;
int char2int(char a);
double CompoEntrSeq(char *seq, int n);

void DCLError(char *msg){
  cerr<<msg<<endl<<flush;
  exit(1);
}

char mjaa[20] = {'C', 'M', 'F', 'I', 'L', 'V','W','Y','A','G','T','S','Q','N','E','D','H','R','K','P' };

ifstream cf,cff,sf,af,pf;
ofstream out;


int main(int argc, char **argv){

//READ PROGRAM OPTIONS

  ReadProgramOptions(argc, argv);

//READING INPUT PARAMETERS

  /*if(argc<13){
    cout<<"Simulating conservatism (evolution) using Z-score.\n"<<flush;
    cout<<"DesignS-acr.linux  \n"<<flush;
    cout<<"               (1) pdb-file\n"<<flush;
    cout<<"               (2) threshold for defining a contact\n"<<flush;
    cout<<"               (3) number of types of amino acids\n"<<flush;
    cout<<"               (4) temperature\n"<<flush;
    cout<<"               (5) average profile scores\n"<<flush;
    cout<<"               (6) number of sequences generated\n"<<flush;
    cout<<"               (7) number of mutations per sequence\n"<<flush;
    cout<<"               (8) energy matrix file name\n"<<flush;
    cout<<"               (9) file with typical concentrations of aa\n"<<flush;
    cout<<"               (10) the error bar in aa concentrations\n"<<flush;
    cout<<"               (11) file with the frequencies of contacts\n"<<flush;
    cout<<"               (12) output file name\n"<<flush;
    cout<<"               [(13) profile matrix]\n"<<flush;
    cout<<"		  [(14) relative weight between mj energy and profile scores]\n"<<flush;
    exit(1);
  }*/
 
  //sseed =   timeseed();
  sseed = -1;

  char pfile[100];
  strcpy(pfile,prog_opt.pfile); //strcpy(pfile,argv[1]);
  double thr = prog_opt.thr; //double thr=atof(argv[2]);

  DCLPDB *p =new DCLPDB[1](pfile);
  if(!p){cerr<<"Cannot allocate memory for p\n"<<flush;exit(1);}
  p->Read();
  //p->GetCenterContMap(thr);
  p->GetCBContMap(thr);
  //p->GetCAContMap(thr);
  p->GetSeq();

  int N=p->N;
  int Ntypes = prog_opt.Ntypes; //int Ntypes=atoi(argv[3]);
  double T = prog_opt.T; //double T=atof(argv[4]);
  double beta=-1/T; 
  long nrun = prog_opt.nrun; //long nrun=atol(argv[6]);  // number of sequences
  long step = prog_opt.step; //long step=atol(argv[7]);  // number of mutations per sequence

  char emfname[100];    // energy matrix
  char ofile[100];      // output file
  char cfile[100];      // concentration file
  char ffile[100];      // frequency file
  char afile[100];       // average score file
  char prfile[100];      // profile score file
  strcpy(afile, prog_opt.afile); //strcpy (afile, argv[5]);
  strcpy(emfname, prog_opt.emfname); //strcpy(emfname,argv[8]);
  strcpy(cfile, prog_opt.cfile); //strcpy(cfile,argv[9]);
  double eps = prog_opt.eps; //double eps=atof(argv[10]);
  //fprintf(stderr, "eps: %f\n", eps);
  strcpy(ffile, prog_opt.ffile); //strcpy(ffile,argv[11]);
  strcpy(ofile, prog_opt.ofile); //strcpy(ofile,argv[12]);
  strcpy(prfile, prog_opt.prfile); //strcpy(prfile,argv[13]);
  double rw = prog_opt.rw;  //double rw = atof(argv[14]);
  if( (rw>1) || (rw<0) ) {
    cout<<"relative weight should be in between 0 and 1"<<endl;
    exit(0);
  }
  cout << "pfile: " << pfile << endl;
  cout << "thr: " << thr << endl;
  cout << "Ntypes: " << Ntypes << endl;
  cout << "T: " << T << endl;
  cout << "nrun: " << nrun << endl;
  cout << "step: " << step << endl;
  cout << "emfname: " << emfname << endl;
  cout << "ofile: " << ofile << endl;
  cout << "cfile: " << cfile << endl;
  cout << "ffile: " << ffile << endl;
  cout << "afile: " << afile << endl;
  cout << "prfile: " << prfile << endl;
  cout << "rw: " << rw << endl;
  cout << "nalign: " << prog_opt.nalign<< endl;
  cout << "calign: " << prog_opt.calign<< endl;
  cout << "Use CM file: " << prog_opt.CM << endl;
  cout << "CM file: " << prog_opt.CMfile << endl;

  
// <--END READING INPUT PARAMETERS


  int i,j,ij;
  /*int Nc=0;
  for(i=0;i<N;i++){
    for(j=0;j<N;j++){
      if(p->CM[i][j])Nc++;
    }
  }*/

// READING ENERGY MATRIX
  ematrix em(emfname,20);
// END READING ENERGY MATRIX

//     for(i=0;i<20;i++) {for(j=0;j<20;j++) cout << em.e[i][j] << " "; cout << endl; }

// OUTPUT FILE
  out.open(ofile,ios::out);
  if(!out)DCLError("Cannot open output file");
// END OUTPUT FILE


// INITIALIZING THE RANGE FOR THE RANDOM NUMBERS
  double factm=Ntypes/(RAND_MAX+1.0);
  double factn=N/(RAND_MAX+1.0);
// END INITIALIZING THE RANGE FOR THE RANDOM NUMBERS


// READING CONCENTRATIONS AND DISPERSION EPS
  cf.open(cfile,ios::in);
  if(!cf)DCLError("Cannot open file with concentrations");

  double *conc  = new double[Ntypes];
  double *concA = new double[20];
  double *concl = new double[Ntypes];
  double *conch = new double[Ntypes];
  double *rconc = new double[Ntypes];
  char *aa      = new char[20]; 

  i=0;
  while((!cf.eof())&&(i<20)){
    cf>>aa[i]>>concA[i];
    //fprintf(stderr, "%f\n", concA[i]);
    i++;
  }
  cf.close();

  switch(Ntypes){
  case 20:
    for(i=0;i<Ntypes;i++)conc[i]=concA[i];
    break;
  case 6:
    for(i=0;i<Ntypes;i++)conc[i]=0;
    for(i=0;i<20;i++)conc[em.ind6(aa[i])]+=concA[i];
    break;
  }

  for(i=0;i<Ntypes;i++){
    concl[i]=conc[i]*(1-eps)*((double)N)/100.0;
    conch[i]=conc[i]*(1+eps)*((double)N)/100.0;
    //fprintf(stderr, "%c\t%d\t%f\t%f\t%f\n", aa[i], i, conc[i]*N/100, concl[i], conch[i]);
  }

  //exit(0);

  delete [] conc;
  delete [] concA;
// END READING CONCENTRATIONS AND DISPERSION EPS



// READING FREQUENCY OF CONTACTS
  cff.open(ffile,ios::in);
  if(!cff)DCLError("Cannot open file with frequency of contacts");

  double *freq = new double[N+1];
  double tmp;
  int itmp;
  for(i=0;i<N+1;i++)freq[i]=0;

  for(i=0;i<N;i++){
    cff>>itmp;
    cff>>freq[itmp];
    if(itmp==1) freq[itmp] = 0;
  }
  cff.close();
// END READING FREQUENCY OF CONTACTS

//READING AVERAGE PROFILE SCORES AND DISPERSION
  af.open(afile, ios::in);  
  if(!af) DCLError("Cannot open file with the average profile scores");

  double *avprof = new double[Ntypes+1];
  double *sigmaprof = new double[Ntypes+1];
  char letter;
  for(i=0;i<Ntypes;i++) {
    af>>letter;
    for(j=0;j<Ntypes;j++) {
      if(aa[j]==letter){af>>avprof[j]>>sigmaprof[j];break;}
    }

  }  

// READ THE PSSM
  pf.open(prfile, ios::in);
  if(!pf) DCLError("Cannot open file with the PSSM");
  
  int *letcast = new int[Ntypes+1]; // facciliate the conversion
  int **pssm = new int *[N+1];  // position specific scoring matrix
  for(i=0;i<N;i++) pssm[i] = new int[Ntypes+1];
  for(i=0;i<Ntypes;i++) { // first line is the letter codes
    pf>>letter;
    for(j=0;j<Ntypes;j++) {
      if(aa[j]==letter) {letcast[i] = j;break;}
    }
  }
  for(i=0;i<N;i++){
    pf>>j; pf>>letter;
    for(j=0;j<Ntypes;j++){
      pf>>pssm[i][letcast[j]];
      //cout<<pssm[i][letcast[j]]<<" ";
    }
  }
//for(i=0;i<N;i++){
//  for(j=0;j<Ntypes;j++){
//    cout<<pssm[i][j]<<" ";
//  }
//  cout<<endl;
//}
  
// INITIALIZING PROBABILITIES
  double **P = new double*[N];
  for(i=0;i<N;i++)P[i] = new double[Ntypes];

  for(i=0;i<N;i++)for(j=0;j<Ntypes;j++)P[i][j]=0;
// END INITIALIZING PROBABILITIES


//--> SIMULATIONS START HERE

  int irod,s,m=0,k;
  double delta=0, delta1=0, delta2=0;
  double en_new,en1_new,en2_new,z_new;
  double enp_new, enp1_new, enp2_new;
  short f1;

  double en=0,gk=0,gk1=0,gk2=0,en1=0,en2=0;
  double sigma,z,minz;
 
  short *f = new short[N];   // the sequence
  short *nn = new short[N];  // number of the neighbours
  double *g = new double[N]; // energy of an aa
  double *g1 = new double[N]; // new energy of an aa
  double *g2 = new double[N]; // new energy of an aa

  short **cc = new short*[N];            // list of neighbours
  for(i=0;i<N;i++)cc[i] = new short[N];

  short *f0 = new short[N];   // the sequence original sequence

  j=0;

  switch(Ntypes){
  case 20:
    for(i=0;i<N;i++){
      j=em.ind(p->Seq[i]);
      //cout<<i<< " " <<j<<"  "<<p->Seq[i]<<endl;
      f0[i]=j;                 // initializing initial sequence
    }
    break;
  case 6:
    for(i=0;i<N;i++){
      j=em.ind6(p->Seq[i]);
      f0[i]=j;                 // initializing initial sequence
    }
    break;
  }


  double averageZP, sigmaZP;
  double averageZN, sigmaZN;
  unsigned long NP;
  double *zs0, *zs1;
  int seqnum = 0;
  int seqnum0 = 0;
  int seqnum1 = 0;
  char outputzs[100];
  for(irod=0;irod<2;irod++){

    if(irod==0) {sprintf(outputzs, "%s.%g.z", prog_opt.calign, rw);}
    else {sprintf(outputzs,"%s.%g.z", prog_opt.nalign, rw);}

    //cout << prog_opt.nalign << endl;

    ofstream zsoutfile(outputzs, ios::out);
    if(!zsoutfile) {
	fprintf(stderr, "native alignment cannot be opened for writing\n");
	exit(0);
    }

    //cout << outputzs << endl;


// CALCULATING REAL CONCENTRATIONS
// INITIALIZING RANDOM SEQUENCE
    for(i=0;i<Ntypes;i++)rconc[i]=0;
/*

    j=0;
    for(i=0;i<N;i++){
      j=int(factm*random());
      f[i]=j;                 // initializing initial sequence
      rconc[j]+=1.0;
    }
*/
// END INITIALIZING RANDOM SEQUENCE

// INITIALIZING ORIGINAL SEQUENCE

    switch(Ntypes){
    case 20:
      for(i=0;i<N;i++){
	f[i]=f0[i];                 // initializing initial sequence
	//cout << i << "\t" << f[i] << endl;
	rconc[f[i]]+=1.0;
      }
      break;
    case 6:
      for(i=0;i<N;i++){
	f[i]=f0[i];                 // initializing initial sequence
	rconc[f[i]]+=1.0;
      }
      break;
    }

    //for(i=0;i<Ntypes;i++) fprintf(stderr, "i: %d;  rconc: %f\n", i, rconc[i]);
    //exit(0);


// END INITIALIZING ORIGINAL SEQUENCE
// END CALCULATING REAL CONCENTRATIONS


// CREATING THE LIST OF NEIGHBOURS
    //for(i=0;i<N;i++)for(j=0;j<N;j++) cc[i][j] = -1;

    ifstream cmf; 
    if(prog_opt.CM) {
  	cmf.open(prog_opt.CMfile, ios::in);
	if(!cmf){cout<<"cannot open the contact matrix file"<<endl;exit(0);}
	for(i=0;i<N;i++) {
	   for(j=0;j<N;j++) {
		cmf >> p->CM[i][j];
		if(ABS(i-j)==1) p->CM[i][j] = 0; // nearest neighbor
	   }
	}
    }

  for(j=0;j<N;j++){
      for(i=0;i<N;i++) {
	if(ABS(i-j)==1) p->CM[j][i]=0;
      }
  }

  for(j=0;j<N;j++){
      for(i=0;i<N;i++) {
      cout<< p->CM[j][i] <<" ";
      }
      cout<<endl;
    }

   //exit(0); 
    k=0;
    for(i=0;i<N;i++){
      nn[i]=0;
      cc[i]=cc[0]+k;
      for(j=0;j<N;j++){
	if(p->CM[i][j]>0){
	  cc[0][k]=j;
	  k++;
	  nn[i]++;
	}
	else p->CM[i][j]=0;
      }
    }
//for(j=0;j<N;j++){
//      for(i=0;i<N;i++) {
//	cout<< p->CM[j][i] <<" ";
//      }
//      cout<<endl;
//    }
//    for(j=0;j<N;j++){
//      for(i=0;i<nn[j];i++) {
//	cout<< cc[j][i] <<" ";
//      }
//      cout<<endl;
//    }
// END CREATING THE LIST OF NEIGHBOURS
   //for(i=0;i<N;i++) fprintf(stderr, "%d\t%d\n", i+1, nn[i]);
   //exit(0);


// CALCULATING INITIAL ENERGY
    en=gk=gk1=gk2=0;
    for(i=0;i<N;i++){
      gk=0;
      switch(Ntypes){
      case 20:
	for(j=0;j<nn[i];j++)gk+=em.e[f[i]][f[cc[i][j]]];
	break;
      case 6:
	for(j=0;j<nn[i];j++)gk+=em.e6[f[i]][f[cc[i][j]]];
	break;
      }
      g[i]=gk;
      en+=gk;
    }
    en/=2;
// END CALCULATING INITIAL ENERGY


// CALCULATING INITIAL AVERAGE ENERGY, DISPERSION AND Z-SCORE
    en1=en2=0;
    for(i=0;i<N;i++){ 
      gk1=0;
      gk2=0;
      for(j=0;j<N;j++){
	if(i!=j){
	  ij=ABS(i-j);
	  switch(Ntypes){
	  case 20:
	    gk1+=(em.e[f[i]][f[j]]*(freq[ij]));
	    gk2+=(em.e2[f[i]][f[j]]*(freq[ij])*(1-freq[ij]));
	    break;
	  case 6:
	    gk1+=(em.e6[f[i]][f[j]]*(freq[ij]));
	    gk2+=(em.e62[f[i]][f[j]]*(freq[ij])*(1-freq[ij]));
	    break;
	  }
	}
      }
      g1[i]=gk1;
      g2[i]=gk2;
      en1+=gk1;
      en2+=gk2;
    }
    en1/=2;
    en2/=2;
    sigma=sqrt(en2);
    if(!sigma)DCLError("Init: sigma = 0");
    z=(en-en1)/sigma;
    minz=z;

// END CALCULATING INITIAL AVERAGE ENERGY, DISPERSION AND Z-SCORE

// CALCULATE INITIAL ENERGY-2: THE PROFILE SCORE
    double enp=0, enp1=0, enp2=0;
    for(i=0;i<N;i++){ 
      enp -= pssm[i][f[i]];}
    for(i=0;i<N;i++){
      enp1 -= avprof[f[i]];}
    for(i=0;i<N;i++){
      enp2 += sigmaprof[f[i]]*sigmaprof[f[i]];}
    if(enp2==0) DCLError("Init: sigma = 0");
    z = (rw*(en-en1)+(1-rw)*(enp-enp1)) /sqrt(rw*rw*en2+(1-rw)*(1-rw)*enp2);
    minz = z;
    cout << "INITIAL Z: " << z<<endl;
    cout << "INITIAL en: " << en << " "<< en1 << " " << en2 <<endl;



/*******************************************/    
// CALCULATE Z-SCORES FOR THE MULTIPLE SEQUENCE ALIGNMENT
    char msa[100];
    ifstream msafp;
    //cout << "Please give the name of the multiple sequence alignment: " <<endl;
    //cin >> msa;
    strcpy(msa, prog_opt.calign); // 100000 profile sequences
    if(irod==1) strcpy(msa, prog_opt.nalign); // native sequences
    msafp.open(msa, ios::in);
    if(!msafp) {
      cerr << "File of multiple alignment not readable"<<endl;
      exit(0);
    }


    char seqname[100];
    char sequen[N+1];
    int sequence[N];
    //char *sequen = new char[N+1];
    //int *sequence = new int[N];
    seqnum = 0; seqnum0 = 0; seqnum1 = 0;

    // print the energy matrix
    /*for(i=0;i<20;i++) {
      for(j=0;j<20;j++) {
	cout << em.e[i][j]<< " ";
      }
      cout << endl;
    }
    for(i=0;i<N;i++) {
      for(j=0;j<N;j++) {
	cout << aa[f[i]] << " "<< aa[f[j]]<< " "<< em.e[f[i]][f[j]] << " ";
      }
      cout << endl;
    }*/

// COUNT THE SEQUENCE NUMBER
while (msafp >> seqname >>sequen   ) {
   seqnum ++;
}

if(irod==0) NP = seqnum;
if(irod==0) seqnum0 = seqnum;
if(irod==1) seqnum1 = seqnum;
cout << "NP " << NP << endl;

/*delete zs;*/
if(irod==0) {zs0 = new double [seqnum0];}  // initialize the zs array
if(irod==1) {zs1 = new double [seqnum1];} 
double averageZ =0;
double sigmaZ = 0;
int seqcount = 0;
msafp.clear();
msafp.seekg( 0 );


while(msafp >> seqname >>sequen   ) {
 
  //cout << "seqname: " << seqname << endl;
  //cout << "sequen:  " << sequen << endl;
  //seqnum++;
  for(i=0;i<N;i++) {
    sequence[i] = char2int(sequen[i]);
    //cout << "I: " << i << " " << sequen[i] << " " << sequence[i] << endl;
  }

  //cout << "here1" << endl;

// CALCULATING INITIAL ENERGY
     en=gk=gk1=gk2=0;
     for(i=0;i<N;i++){

       if(sequence[i]==20) continue;
       gk=0;
       switch(Ntypes){
       case 20:
	//cout << "nni: " << nn[i] << endl;
	 for(j=0;j<nn[i];j++){
	   if( sequence[cc[i][j]]==20) continue;
	//cout << "i: " << i << " " << j << " " << sequence[i] << " " << sequen[i] << " " << sequence[cc[i][j]] << " " << gk << endl;
	   // excluding the adjacent contacts
	   //if((i==cc[i][j]+1)||(i==cc[i][j]-1)) continue;
	   gk+=em.e[sequence[i]][sequence[cc[i][j]]];
	//cout << "i: " << i << " " << j << " " << sequence[i] << " " << sequence[cc[i][j]] << " " << gk << endl;
	 }
	 break;
       case 6:
	 for(j=0;j<nn[i];j++)gk+=em.e6[sequence[i]][sequence[cc[i][j]]];
	 break;
       }
       g[i]=gk;
       en+=gk;
     }
     en/=2;
// END CALCULATING INITIAL ENERGY

   //cout << "here2" << endl;

// CALCULATING INITIAL AVERAGE ENERGY, DISPERSION AND Z-SCORE
    en1=en2=0;
    for(i=0;i<N;i++){
      if(sequence[i]==20) continue;
      gk1=0;
      gk2=0;
      for(j=0;j<N;j++){
	if(sequence[j]==20) continue;
	if(i!=j){
	  ij=ABS(i-j);
	  //if( (i==j+1) ||( i==j-1) ) continue;
	  switch(Ntypes){
	  case 20:
	    gk1+=(em.e[sequence[i]][sequence[j]]*(freq[ij]));
	    gk2+=(em.e2[sequence[i]][sequence[j]]*(freq[ij])*(1-freq[ij]));
	    break;
	  case 6:
	    gk1+=(em.e6[sequence[i]][sequence[j]]*(freq[ij]));
	    gk2+=(em.e62[f[i]][f[j]]*(freq[ij])*(1-freq[ij]));
	    break;
	  }
	}
      }
      g1[i]=gk1;
      g2[i]=gk2;
      en1+=gk1;
      en2+=gk2;
    }
    en1/=2;
    en2/=2;
    sigma=sqrt(en2);
    if(!sigma)DCLError("Init 1: sigma = 0");
    //cout << "en: " << en << endl;
    //cout << "en1: " << en1 <<endl;
    //cout << "sigma: " << sigma << endl;
    z=(en-en1)/sigma;
    //cout << "z: " << z << endl;
    minz=z;

// END CALCULATING INITIAL AVERAGE ENERGY, DISPERSION AND Z-SCORE

   //cout << "here3" << endl;

// CALCULATE INITIAL ENERGY-2: THE PROFILE SCORE
    double enp=0, enp1=0, enp2=0;
    for(i=0;i<N;i++){ 
      if(sequence[i]==20) continue;
      enp -= pssm[i][sequence[i]];}
    for(i=0;i<N;i++){
      if(sequence[i]==20) continue;
      enp1 -= avprof[sequence[i]];}
    for(i=0;i<N;i++){
      if(sequence[i]==20) continue;
      enp2 += sigmaprof[sequence[i]]*sigmaprof[sequence[i]];}
    if(enp2==0) DCLError("Init 2: sigma = 0");
    z = (rw*(en-en1)+(1-rw)*(enp-enp1))/sqrt(rw*rw*en2+(1-rw)*(1-rw)*enp2);
    if(irod==0)zs0[seqcount] = z; 
    if(irod==1)zs1[seqcount] = z;seqcount ++; averageZ += z;

    //cout << "here4" << endl;

// calculate entropy of the sequence
    char designSeq[N];
    for(i=0;i<N;i++) designSeq[i] = mjaa[sequence[i]];
    double designSeqEntr;
    designSeqEntr = CompoEntrSeq(designSeq, N);

    zsoutfile << z << endl;
    //cout << z << endl;
    //cout << seqname << " "<< "z: " << z << endl; 
    //cout << "en: " << en << "enp: " << enp << " entropy: " << designSeqEntr <<endl;
    //  " " << en1 << " " << en2 << endl;
}

   //cout << "here: " << endl;


  //CALCULATE Z-SCORE AVERAGE AND VARIANCE OF THE PROFILE ALIGNMENT
  if(irod == 0) {
    averageZ /= seqnum0;
    for(i=0;i<seqnum0;i++) {
      sigmaZ += (zs0[i]-averageZ) * (zs0[i]-averageZ);
    }
    sigmaZ /= (seqnum0 - 1);
    sigmaZ = sqrt(sigmaZ);
  
    averageZP = averageZ; sigmaZP = sigmaZ;
    cout << "averagezp: " << averageZP << endl;
  }

  if(irod==1) {
    averageZ /= seqnum1;
    for(i=0;i<seqnum1;i++) {
	sigmaZ += (zs1[i]-averageZ) * (zs1[i]-averageZ);
    }
    sigmaZ /= (seqnum1 - 1);
    sigmaZ = sqrt(sigmaZ);
    
    averageZN = averageZ; sigmaZN = sigmaZ;
    cout << "averagezn: " << averageZN << endl;
  }

} // End of 2 irod cycles

////CALCULATE THE HENIKOFF WEIGHT OF THE NATIVE SEQUENCES
//  align_info *nativeAlign = new align_info(prog_opt.nalign);
//  int native_nal = nativeAlign->getNal();
//  double *native_hw;
//  native_hw = nativeAlign->hw(1); // index [1..N], not [0..N-1]
//  /*for(i=0;i<=native_nal;i++) {
//	cout << i << ":\t" << native_hw[i] << endl;
//  }*/
//
//
//
//
////CALCULATE KS INTEGRAL
//  double ks_value = KS_integral(zs1, seqnum1,  1);
//  cout <<"***"<< endl;
//  double ks_value_w = KS_integral_weighted(zs1, native_hw+1, seqnum1, 1);
//  cout << "ks_intergral: " << ks_value << endl;
//  cout << "ks_intergral_w: " << ks_value_w << endl;
//  cout << "averageZP: " << averageZP << " sigmaZP: " << sigmaZP << endl;
//  //cout << zs1[0] << " " << zs1[seqnum1-1] << " " << averageZP << " " << sigmaZP << endl;
//  double gaussian_i = Gaussian_integral_x1_x2(zs1[0], zs1[seqnum1-1], averageZP, sigmaZP);
//  //cout << gaussian_i << endl;
//  double gaussian_r = 10 - Gaussian_integral_x1_x2(zs1[0], zs1[seqnum1-1], averageZP, sigmaZP);
//  //double gaussian_r = 10 - Gaussian_integral_x1_x2(zs1[seqnum1-1],zs1[seqnum1-1]+10,averageZP,sigmaZP);
//  cout << "gaussian integral: " << gaussian_i << endl;
//  //cout << "gaussian residual: " << gaussian_r << endl;
//  cout << rw << " Difference of the integral: " << ks_value - gaussian_i << endl;
//
////CALCULATE AVERAGE WEIGHTED AVERAGE Z-SCORE
//  double sum_of_weights = 0;
//  double average_ZN = 0;
//  double average_ZNUNWEIGHTED  = 0;
//  for(i=0;i<native_nal;i++) {
//	sum_of_weights += native_hw[i+1];
//	average_ZN += native_hw[i+1] * zs1[i];
//	average_ZNUNWEIGHTED += zs1[i];
//	//cout << native_hw[i+1] << "\t" << zs1[i]<< endl;
//  }
//  average_ZN /= sum_of_weights;
//  average_ZNUNWEIGHTED /= native_nal;
//  double sigma_ZNUNWEIGHTED = 0;
//  for(i=0;i<native_nal;i++) {
//	sigma_ZNUNWEIGHTED += (zs1[i] - sigma_ZNUNWEIGHTED) * (zs1[i] - sigma_ZNUNWEIGHTED) ;
//  }
//  sigma_ZNUNWEIGHTED = sqrt(sigma_ZNUNWEIGHTED/(native_nal-1));
//  cout <<  rw << " Average native z score: " << average_ZN << endl;
//
////CALCULATE UNWEIGHTED MEAN DIFFERENCE STATISTICS
  double mean_diff;
  double p_mean_diff;
  tutest1(zs1-1, seqnum1,averageZP, sigmaZP*sigmaZP, NP, &mean_diff, &p_mean_diff);

//  cout << rw << " student t statistic for mean difference: "<< mean_diff << "\t" << p_mean_diff << endl;

  cout << "RW: " << rw << " " << "0" << " " << mean_diff << " " << (averageZN - averageZP) /(sigmaZN + sigmaZP) << endl;


  return 1;
}


int char2int(char a) {
  switch(a) {
    case 'c': case 'C': return 0; break;
    case 'm': case 'M': return 1; break;
    case 'f': case 'F': return 2;break;
    case 'i': case 'I': return 3; break;
    case 'l': case 'L': return 4; break;
    case 'v': case 'V': return 5; break;
    case 'w': case 'W': return 6; break;
    case 'y': case 'Y': return 7; break;
    case 'a': case 'A': return 8; break;
    case 'g': case 'G': return 9; break;
    case 't': case 'T': return 10; break;
    case 's': case 'S': return 11; break;
    case 'q': case 'Q': return 12; break;
    case 'n': case 'N': return 13; break;
    case 'e': case 'E': return 14; break;
    case 'd': case 'D': return 15; break;
    case 'h': case 'H': return 16; break;
    case 'r': case 'R': return 17; break;
    case 'k': case 'K': return 18; break;
    case 'p': case 'P': return 19; break;
    case 'x': case 'X': case '.': case '-': return 20; break;
  }
}






//Compositional entropy for a sequence 
double CompoEntrSeq(char *seq, int n) 
{

  int i, j;
  int aaCount[20];
  double aaFreq[20];
  double entr = 0;;
  
  for(i=0;i<20;i++) aaCount[i] = 0;
  for(i=0;i<n;i++) {
    aaCount[char2int(seq[i])]++;
  }

  for(i=0;i<20;i++) {
    aaFreq[i] = 1.0*aaCount[i]/n;
    if(aaFreq[i]) entr+=aaFreq[i]*log(aaFreq[i]);
  }

  return entr;
}
  



