#include "header_cpp.h"
#include "myreadpdb.h"
#include "ematrix.h"
//#include "align_info.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);
double ran3(long *idum);

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

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

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

//READING INPUT PARAMETERS

  ReadProgramOptions(argc, argv);

  // sseed =   timeseed(); 
  sseed = -1;
	 //cout << sseed << endl;

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

  if(prog_opt.sseed) {
	if(prog_opt.sseed>0) {
	   cout << "seed must be less than zero\n";
 	   cout << prog_opt.sseed <<endl;
	   exit(0);
	}
	else sseed = prog_opt.sseed;
  }

  
  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 << "seed: " << prog_opt.sseed << endl;
  if(prog_opt.seqlength)cout << "sequence length: "<< prog_opt.seqlength << endl;
  if(prog_opt.initialseq)cout << "initial sequence: " << prog_opt.initialseq << endl;

  cout << "Use CM file: " << prog_opt.CM << endl;
  if(prog_opt.CM) cout << "CM file: " << prog_opt.CMfile << endl;
  cout << 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


// 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();
//  for(i=0;i<Ntypes;i++){
//    cout<<i<<" "<<aa[i]<<endl;
//  }
  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];
  }
  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;}
    }

  }  
//  for(i=0;i<Ntypes;i++){
//        cout<<aa[i]<<" "<<avprof[i]<<" "<<sigmaprof[i]<<endl;
//  }

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


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

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



  double Z_endrun[nrun]; // install the zscore of the endpoint at each run

  for(irod=0;irod<nrun;irod++){


// 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
	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;
    
    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");
    sigma = sqrt(rw*rw*en2+(1-rw)*(1-rw)*enp2);
    z = (rw*(en-en1)+(1-rw)*(enp-enp1))/sqrt(rw*rw*en2+(1-rw)*(1-rw)*enp2);
    minz = z;
    /*cout << "MJ terms: " << en << "\t" << en1 << "\t" << en2 << endl;
    cout << "Profile terms: " << enp << "\t" << enp1 << "\t" << enp2 << endl;
    cout << "Z: " << z<<endl;
	*/
    


// MONTE-CARLO SIMULATIONS

    for(s=0;s<step;s++){
      en_new=en1_new=en2_new=z_new=0;
      enp_new = enp1_new = enp2_new = 0;
      gk=gk1=gk2=0;
      k=int(N*ran3(&sseed));
      m=int(Ntypes*ran3(&sseed));
      //cout << factn<<" "<<factm<<" "<<k << " " << m << endl;

// Check for the concentration limitations
      f1=f[k];
      if(((rconc[m]+1.0)<conch[m])&&((rconc[f1]-1.0)>concl[f1])){
// Check for the concentration limitations is ok
        f[k]=m;

	switch(Ntypes){
	case 20:
	  for(j=0;j<nn[k];j++)gk+=em.e[f[k]][f[cc[k][j]]];
	  break;
	case 6:
	  for(j=0;j<nn[k];j++)gk+=em.e6[f[k]][f[cc[k][j]]];
	  break;
	}

        for(j=0;j<N;j++){
          if(j!=k){
	    ij=ABS(k-j);
	    switch(Ntypes){
	    case 20:
	      gk1+=(em.e[f[k]][f[j]]*(freq[ij]));
	      gk2+=(em.e2[f[k]][f[j]]*(freq[ij])*(1-freq[ij]));
	      break;
	    case 6:
	      gk1+=(em.e6[f[k]][f[j]]*(freq[ij]));
	      gk2+=(em.e62[f[k]][f[j]]*(freq[ij])*(1-freq[ij]));
	      break;
	    }

          }
        }
        delta=gk-g[k];
        delta1=gk1-g1[k];
        delta2=gk2-g2[k];
        en_new=en+delta;
        en1_new=en1+delta1;
        en2_new=en2+delta2;
        sigma=sqrt(en2_new);

// profile score of the mutation sequence
        enp_new = enp-pssm[k][f[k]] + pssm[k][f1];
        enp1_new = enp1 - avprof[f[k]] + avprof[f1];
	enp2_new = enp2 + sigmaprof[f[k]]*sigmaprof[f[k]]-sigmaprof[f1]*sigmaprof[f1];
	if(enp2_new==0) DCLError("Init: sigma = 0");
	sigma = sqrt( rw*rw*en2_new+(1-rw)*(1-rw)*enp2_new );
	z_new =  ( rw*(en_new-en1_new)+(1-rw)*(enp_new-enp1_new) )/sqrt(rw*rw*en2_new+(1-rw)*(1-rw)*enp2_new);
//	cout << "Z_new: " << z_new <<endl;

        delta=z_new-z;
        if(delta>0){
          if( ran3(&sseed) > exp(delta*beta)){
            f[k]=f1;
	    //cout <<s+1<< " " << z<< endl; /* individual zscores output */
            continue;
          }
        }
// modify conc list
        if(f[k]!=f1){
          rconc[m]+=1.0;
          rconc[f1]-=1.0;
        }
// end of modifying of conc list

// update the pairwise energy items of each position
        g[k]=gk;
        g1[k]=gk1;
        g2[k]=gk2;
        for(i=0;i<N;i++){
	  if(k!=i){
	    ij=ABS(i-k);

	    switch(Ntypes){
	    case 20:
	      if(p->CM[i][k])g[i]+=(em.e[f[i]][f[k]]-em.e[f[i]][f1]);
	      g1[i]+=((em.e[f[i]][f[k]]-em.e[f[i]][f1])*(freq[ij]));
	      g2[i]+=((em.e2[f[i]][f[k]]-em.e2[f[i]][f1])*(freq[ij])*(1-freq[ij]));
	      break;
	    case 6:
	      if(p->CM[i][k])g[i]+=(em.e6[f[i]][f[k]]-em.e6[f[i]][f1]);
	      g1[i]+=((em.e6[f[i]][f[k]]-em.e6[f[i]][f1])*(freq[ij]));
	      g2[i]+=((em.e62[f[i]][f[k]]-em.e62[f[i]][f1])*(freq[ij])*(1-freq[ij]));
	      break;
	    }
          } 
        }
        z=z_new;
	//fprintf(stderr, "%d %5.3f\n", s+1,z); /* individual zscore output */
        en=en_new;
        en1=en1_new;
        en2=en2_new;
	enp = enp_new;
	enp1 = enp1_new;
	enp2 = enp2_new;

        if(minz>z)minz=z;
//	for(i=0;i<N;i++)P[i][f[i]]+=1;

      }
    }

// END MONTE-CARLO SIMULATIONS

// BEGIN SEQUENCE OUTPUT  
    // calculate complexity of the designed sequence
    char designSeq[N+1];
    for(i=0;i<N;i++) designSeq[i] = mjaa[f[i]];
    designSeq[N] = '\0';
    double designSeqEntr;
    designSeqEntr = CompoEntrSeq(designSeq, N);

    // print out the designed sequence at the end of each run
    cout << "design" << z << "\t";
    cout << designSeq << endl;
    Z_endrun[irod] = z;
    //cout<<irod+1<<"  E= "<<rw*en+(1-rw)*enp<< " " <<rw*en1+(1-rw)*enp1 << " +- "<<sigma<<" Z=  "<<z<< " Entropy: " << designSeqEntr << endl<<flush;

    // output to the fasta file
    /*if(Ntypes==20) {
	fprintf(seqfile, ">gii|%d|t%3.2f|cshr%2.1f|mstep%d|%s|%3.1f|%3.1f\n", irod, T, thr, step, cfile, eps,rw);
        for(i=0;i<N;i++)cout<<mjaa[f[i]]<<flush;
	for(i=0;i<N;i++) fprintf(seqfile, "%c", mjaa[f[i]]);
	fprintf(seqfile, "\n");
    }
    cout<<endl<<flush;
    */
// END SEQUENCE OUTPUT  

   for(i=0;i<N;i++)P[i][f[i]]+=1;

  }


  // calculate the average energy of the nrun zscores
  //double Z_AVE;
  //double Z_SIGMA;
  //avevar(Z_endrun, nrun-1, &Z_AVE, &Z_SIGMA);

  //cout << "T " << T << "\t" << Z_AVE << endl;

//--> END SIMULATIONS

// CLEAN 
  delete [] concl;
  delete [] conch;
  delete [] rconc;
  delete [] aa;

  delete [] f;
  delete [] f0;
  delete [] nn;

  delete [] g1;
  delete [] g2;
  delete [] g;
  delete [] freq;
// END CLEAN


// COMPUTING THE ENTROPY
  double *S = new double[N];
  
/*  fprintf(stderr, "nrun: %d\n", nrun);
  for(i=0;i<N;i++){
    S[i]=0;
    for(j=0;j<Ntypes;j++){
      P[i][j]/=((double)nrun);
      //fprintf(stderr, "%5.3lf ", P[i][j]);
      if(P[i][j])S[i]-=(P[i][j]*log(P[i][j]));
    }
    
//    cout<<(i+1)<<"\t"<<S[i]<<endl<<flush;
  fprintf(stderr, "%d\t%lf\n", i+1, S[i]);
  }
  */
// END COMPUTING THE ENTROPY

// COMPUTING CORRELATION COEEFICIENT
  /*if(argc==14){
    double *Sacr = new double[N];
    sf.open(sfile,ios::in);
    for(i=0;i<N;i++){
      switch(Ntypes){
      case 20:
	sf>>tmp>>tmp>>tmp>>tmp>>tmp>>tmp>>tmp>>tmp;
	sf>>Sacr[i];
	break;
      case 6:
	sf>>tmp>>tmp>>tmp>>tmp>>tmp>>tmp>>tmp;
	sf>>Sacr[i]>>tmp;
	break;
      }
    }
    delete [] Sacr;
  }*/
// END COMPUTING CORRELATION COEEFICIENT  


// CLEANUP
  out.close();

  for(i=0;i<N;i++)delete [] P[i];
  delete [] P;
  delete [] S;


  //for(i=0;i<N;i++)delete [] cc[i];
  //delete [] cc;

  //delete p;
  
  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;
}
  
#define MBIG 1000000000
#define MSEED 161803398
#define MZ 0
#define FAC (1.0/MBIG)

double ran3(long *idum)
/* returns uniform random deviate between 0.0 and 1.0
idum - is a negative integer, do not alter it for successive calls*/

{
        static int inext,inextp;
        static long ma[56];
        static int iff=0;
        long mj,mk;
        int i,ii,k;

        if (*idum < 0 || iff == 0) {
                iff=1;
                mj=MSEED-(*idum < 0 ? -*idum : *idum);
                mj %= MBIG;
                ma[55]=mj;
                mk=1;
                for (i=1;i<=54;i++) {
                        ii=(21*i) % 55;
                        ma[ii]=mk;
                        mk=mj-mk;
                        if (mk < MZ) mk += MBIG;
                        mj=ma[ii];
                }
                for (k=1;k<=4;k++)
                        for (i=1;i<=55;i++) {
                                ma[i] -= ma[1+(i+30) % 55];
                                if (ma[i] < MZ) ma[i] += MBIG;
                        }
                inext=0;
                inextp=31;
                *idum=1;
        }
        if (++inext == 56) inext=1;
        if (++inextp == 56) inextp=1;
        mj=ma[inext]-ma[inextp];
        if (mj < MZ) mj += MBIG;
        ma[inext]=mj;
        return mj*FAC;
}
#undef MBIG
#undef MSEED
#undef MZ
#undef FAC






