/* This code is compiled using the gcc compiler versio 4.1.2 on the high
 * performance computing cluster, tara, at UMBC. This code was written by
 * Andrew Coates in partnership with Alexey Ilchenko and advisor Matthias K.
 * Gobbert as a part of the REU Site: Interdisciplinary Program in High
 * Performance Computing at UMBC in Summer 2011. For more information see our
 * Tech Report at www.umbc.edu/hpcf > Publications as Technical Report
 * HPCF-2011-13. This code is last updated August 11, 2011.
 *
 * Calculates the expected entropy of a string of length n from a set of size 4
 * as specifically applied to base pair analysis
 *
 * To compile, run " >> make "" at the linux command prompt. This will create
 * an executable named "EHn". To run, enter " >> ./EHn pa pt pc pg n " at the
 * command prompt. The values of pa,pt,pc and pg are the probabilities
 * associated with A,T,C and G. The value n is the number of samples.
 * This will print out n,pa,pt,pc,pg and the expected value of Hn.
 */
 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "EHn.h"


/* The main function which takes in the inputted probabilites and n and computes
 * the expected value of Hn.
 */
int main(int argc, char *argv[]){
  double sum,ps,hs,coef;
  int i,j,*cc, *row,n,rows;
  float pa,pt,pc,pg;

  pa = atof(argv[1]);
  pt = atof(argv[2]);
  pc = atof(argv[3]);
  pg = atof(argv[4]);
  n = atoi(argv[5]);
  
  printf("n=%d pa=%f pt=%f pc=%f pg=%f\n",n,pa,pt,pc,pg);

  /* allocates the memory of the matrix and row arrays*/
  /* the number of rows in cc */
  rows = (n+1)*(n+2)*(n+3)/6;
  /* the compostions matrix */
  cc = malloc(4*rows*sizeof(int));
  /* will act as the current row during iterations through cc */
  row = malloc(4*sizeof(int));
  /* creates the composition matrix, cc */
  Counts(cc,n);

  sum = 0.0;
  /* iterates through the compositions matrix's rows */
  for(i=0;i<rows;i++){
    /* selects the values of the row in cc and assigns them to *row */
    for(j=0;j<4;j++){
      row[j] = cc[j+4*i];
    }
    /* computes p(s)=ps, H(s)=hs and (n choose ki)=coef and takes the product
     * which is then added to the sum
     */
    ps = prob(pa,pt,pc,pg,n,row);
    hs = Hs(n,row);
    coef = Multinomial(n,row);
    sum -= coef*hs*ps;  
  }

  /* prints results */
  printf("E(H_n) = %f\n",sum);
  free(cc);
  free(row);
}

/* computes the log of factorials
 */
long double factorial(int n){
  return lgamma(n+1);
}

/* computes the multinomial coefficients based on the row in the matrix
 */
long double Multinomial(int n, int *row){
  return exp(factorial(n)-((factorial(row[0])+factorial(row[1])
			    +factorial(row[2])+factorial(row[3]))));
}

/* computes the entropy of the row in the matrix
 */
double Hs(int n, int *row){
  int i;
  double hs,p;
  hs = 0;
  for(i=0;i<4;i++){
    if(row[i] != 0){
      p = (double)(row[i])/(double)(n);
      hs += p*log(p)/log(2);
    }
  }
  return hs;
}

/* Computes the probability of a certain string by raising each probability
 * to its associated value in the row and multiplying them all together 
 */
double prob(float pa, float pt, float pc,float pg, int n, int *row){
  double val;
  val = pow((double)pa,row[0])*pow((double)pt,row[1])
    *pow((double)pc,row[2])*pow((double)pg,row[3]);
  return val;
}

/* Creates the array cc of length 4*rows based on n.
 * cc is the Compositions Matrix
 */
void Counts(int *cc, int n){
  int m,i,j,k,l;
  m = 0;
  for(i=0;i<=n;i++){
    for(j=0;j<=n;j++){
      for(k=0;k<=n;k++){
	l = n - (i + j + k);
	if(l>=0){
	  cc[0+4*m] = i;
	  cc[1+4*m] = j;
	  cc[2+4*m] = k;
	  cc[3+4*m] = l;
	  m++;
	}
      }
    }
  }
}