/* ARIADNE V.1.3

Copyright 

Richard Mott 2000

Wellcome Trust Centre For Human Genetics
Univeristy of Oxford
Roosevelt Drive
Oxford OX3 7AD
UK

The software package ARIADNE is distributed in the hope that it will be
useful, but in order that the University as a charitable foundation
protects its assets for the benefit of its educational and research
purposes, the University makes clear that no condition is made or to
be implied, nor is any warranty given or to be implied, as to the
accuracy of ARIADNE, or that it will be suitable for
any particular purpose or for use under any specific conditions, or that
the content or use of ARIADNE will not constitute
or result in infringement of third-party rights.  Furthermore, the
University disclaims all responsibility for the use which is made of 
ARIADNE.


*/
#include<stdio.h>
#include<limits.h>
#include<math.h>
#include <stdlib.h>
#include"seq_util.h"
#include"profile.h"
#include"cl.h"
#include"gapstat.h"
    
/* functions for Karlin-Altschul stats on profiles */

int KarlinAltschulProfileStatistics( PROFILE *psi, double *freq, double *lambda, double *Kminus, double *Kplus, double *H, double *r, double *s ) {

/* compute all necessary quantities for ungapped statistics 

   inputs are the matrix and symbol frequencies freq1, freq2

   outputs are 
   lambda            the exponential rate 
   Kminus, Kplus     bounds on K
   H                 the entropy (for edge corrections)
   r, s              quantities needed for GEM gapped statistics

*/

  double *h;
  int hmin, hmax;
  double R1, R2, R3;
  double mean;

  h = GetPsiH( psi, freq, &hmin, &hmax, &mean );
  if ( mean < 0.0 && (*lambda = solve_for_lambda( h, hmin, hmax )) > 0 ) {
    
    *H = entropy_H(h,hmin,hmax,*lambda);
    
    iglehart( h, hmin, hmax, *lambda, -1, 1, &R1, &R2, &R3, Kplus, Kminus );
    
    *Kplus /= R3;
    *Kminus /= R3;
    
    *s = R1;
    *r = R2;
    
    free(h+hmin);

    /*    fprintf(stderr, "lambda=%8.6f K-=%8.6f K+=%8.6f H=%8.4f r=%8.4f s=%8.4f\n", *lambda,  *Kplus, *Kminus, *H, *r, *s);  */
    return 1;
  }
  else {
    int n;
    /*    fprintf( stderr, "mean : %e h distribution:\n", mean);
        for(n=hmin;n<=hmax;n++){
      fprintf(stderr, "%5d %e\n", n, h[n]); 
    } */
    free(h+hmin);
    return 0;
  }

}

double *GetPsiH( PROFILE *psi, double *freq, int *hmin, int *hmax, double *mean ) {
  
  int i, j, n;
  double *h;
  *hmin = INT_MAX;
  *hmax = INT_MIN;
  *mean = 0.0;

  for(i=0;i<psi->length;i++) {
    for(j='A';j<='Z';j++) {
      if ( (n=psi->profile[i][j]) > *hmax ) 
	*hmax = n;
      if ( (n=psi->profile[i][j]) < *hmin ) 
	*hmin = n;
    }
  }

  h = (double*)calloc(*hmax-*hmin+1,sizeof(double))-*hmin;

  for(i=0;i<psi->length;i++) 
    for(j='A';j<='Z';j++) 
      h[psi->profile[i][j]]+=freq[j];

  for(j=*hmin;j<=*hmax;j++)
    h[j] /= (double)psi->length;


  for(j=*hmin;j<=*hmax;j++)
    *mean += h[j]*j;

  return h;
}

int SetDefaultParameters( PROFILE *psi ) {

  double *freq = RobinsonResidueFrequencies( );
  double Kplus, r;
  int status;

  status = KarlinAltschulProfileStatistics( psi, freq, &(psi->L), &(psi->K), &Kplus, &(psi->H), &r, &(psi->s) );
  free(freq);
  return status;
}