/*
************************************************************************
*
*   imp.C -
*
*   Copyright (c) 1995
*
*   ETH Zuerich
*   Institut fuer Molekularbiologie und Biophysik
*   ETH-Hoenggerberg
*   CH-8093 Zuerich
*
*   All Rights Reserved
*
*   Date of last modification : 95/09/15
*   Pathname of SCCS file     : /export/home3/cb/garant-1.0/src/SCCS/s.imp.C
*   SCCS identification       : 1.2
*
************************************************************************
*/
/**************************************************************************/
/*                             imp.C                                      */
/*                                                                        */
/* determine significance of deviation of two mean values                 */
/**************************************************************************/

#include <stdio.h>
#include <math.h>
#include "imp.h"

#define ITMAX 100
#define EPS 2.0e-4

/* Error function from numerical recipes in C; Press et al. 1988,      */
/* Cambridge University Press                                          */

float erfcc(float x)
{
   float t,z,ans;
   z=fabs(x);
   t=1.0/(1.0+0.5*z);
   ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
       t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
       t*(-0.82215223+t*0.17087277)))))))));
   return x>=0.0 ? ans : 2.0-ans;
}

float gammln(float xx)
{
   float x,tmp,ser;
   static float cof[6]={76.18009173,-86.50532033,24.01409822,
                        -1.231739516,0.120858003,-0.536382e-5};
   int j;

   x=xx-1.0;
   tmp=x+5.5;
   tmp -= (x+0.5)*log(tmp);
   ser=1.0;
   for(j=0;j<=5;j++) {
      x+= 1.0;
      ser += cof[j]/x;
   }
   return -tmp+log(2.50662827469*ser);
}

float betacf(float a,float b, float x)
{
   float qap,qam,qab,em,tem,d;
   float bz,bm=1.0,bp,bpp;
   float az=1.0,am=1.0,ap,app,aold;
   int m;

   qab=a+b;
   qap=a+1.0;
   qam=a-1.0;
   bz=1.0-qab*x/qap;
   for(m=1;m<=ITMAX;m++) {
      em=(float)m;
      tem=em+em;
      d=em*(b-em)*x/((qam+tem)*(a+tem));
      ap=az+d*am;
      bp=bz+d*bm;
      d= -(a+em)*(qab+em)*x/((qap+tem)*(a+tem));
      app=ap+d*az;
      bpp=bp+d*bz;
      aold=az;
      am=ap/bpp;
      bm=bp/bpp;
      az=app/bpp;
      bz=1.0;
      if( fabs(az-aold)<EPS*fabs(az) ) return az;
   }
   printf("ERROR when calculating students t distribution\n");
   return 0.0;
}


float betai(float a,float b,float x)
{
   float bt;

   if(x==0.0 || x==1.0) bt=0.0;
   else
      bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x));

   if (x<(a+1.0)/(a+b+2.0))
      return bt*betacf(a,b,x)/a;
   else
      return 1.0-bt*betacf(b,a,1.0-x)/b;
}

Imp::Imp()
{
   x1=0.0;
   x2=0.0;
   q1=0.0;
   q2=0.0;
   n1=0;
   n2=0;
   sig=1.0;
}

void Imp::update(float v)
{
   x1 += v;
   q1 += v*v;
   n1++;
}

void Imp::calcsignif()
{
   float s1,m1,s2,m2;
   float t,svar,df;

   if(n1>0) {
      m1=x1/n1;
      s1=q1/n1-m1*m1;
      if(s1<0.0) s1=0.0; 	/* correct numerical problems */
   }
   if(n2>0) {
      m2=x2/n2;
      s2=q2/n2-m2*m2;
      if(s2<0.0) s2=0.0; 	/* correct numerical problems */
   }
   if(n1>1 && n2>1 && s2>0.0 && s1>0.0) {
      df=n1+n2-2;
      svar=((n1-1)*s1 + (n2-1)*s2)/df;
      t=(m1-m2)/sqrt(svar*(1.0/n1+1.0/n2));
      sig=1.0 - betai(0.5*df,0.5,df/(df+t*t));
   }
 
   if(sig<0.0 || sig>1.0)
      printf("%7.3f %7.3f %7.3f %7.3f --> %7.3f %7.3f\n",
             m1,sqrt(s1),m2,sqrt(s2),sqrt(svar*(1.0/n1+1.0/n2)),sig);

   x2=0.9*x2+x1;  q2=0.9*q2+q1;  n2=0.9*n2+n1;
   x1=0.0; q1=0.0; n1=0;
}