#include<math.h>
#include<float.h>
#include<assert.h>

#include "genlimit.h"

static void csint02cut(double xlo1,double xlo2,double xhi,double logscale,
	          int nchan,const int nobs[],const EB chan[],
		  int ngl,const double xgl[],const double lwgl[],PRIOR prior,
	          double* int1,double* int2);

/*      Joel Heinrich  8 April 2005

    returns:
      *int1 = integral from xlo1 to xhi
      *int2 = integral from xlo2 to xhi
      *v11  = variance of *int1
      *v12  = covariance between *int1 and *int2
      *v22  = variance of *int2

See CDF note 7587 for details:
http://www-cdf.fnal.gov/publications/cdf7587_genlimit.pdf

*/



void csint2cut(double xlo1,double xlo2,double xhi,
	       int nchan,int nens,const int nobs[],const EB* ens,
	       int* ngl,double xgl[],double lwgl[],PRIOR prior,
	       double* int1,double* int2,
	       double* v11,double* v12,double* v22) {
  int i,ntot=0;
  double sum1=0, sum21=0, sum2=0, sum22=0, sump=0, logscale=0;
  const EB* p=ens;

  assert(nens>0);
  assert(prior==flat || prior==corr);
  for(i=0;i<nchan;++i) {
    ntot += nobs[i];
    logscale -= lgamma(nobs[i]+1);
  }
  if(*ngl<=0)
    gausslaguerre(xgl,lwgl,*ngl=1+ntot/2,0.0);
  
  for(i=0;i<nens;++i) {
    double t1=0, t2=0;
    csint02cut(xlo1,xlo2,xhi,logscale,nchan,nobs,p,*ngl,xgl,lwgl,
	       prior,&t1,&t2);
    sum1 += t1;
    sum21 += t1*t1;
    sum2 += t2;
    sum22 += t2*t2;
    sump += t1*t2;
    p+=nchan;
  }

  {
    const double rnens = 1.0/nens;
    sum1 *= rnens;
    sum21 *= rnens;
    sum2 *= rnens;
    sum22 *= rnens;
    sump *= rnens; 
    if(nens>1) {
      const double rn1 = 1.0/(nens-1);
      *v11 = (sum21-sum1*sum1)*rn1;
      *v22 = (sum22-sum2*sum2)*rn1;
      *v12 = (sump-sum1*sum2)*rn1;
    } else {
      *v11 = 1;
      *v22 = 1;
      *v12 = 0;
    }
  }

  *int1 = sum1;
  *int2 = sum2;
  return;
}

static void csint02cut(double xlo1,double xlo2,double xhi,double logscale,
		    int nchan,const int nobs[],const EB chan[],
		    int ngl,const double xgl[],const double lwgl[],PRIOR prior,
		    double* int1,double* int2) {
  int i,k;
  double sum1=0, sum2=0, sum3=0, esum=0, bsum1=0, bsum2=0, bsum3=0, resum;

  for(i=0;i<nchan;++i) {
    const double ee=chan[i].e, bb=chan[i].b;
    esum += ee;
    bsum1 += bb + xlo1*ee;
    bsum2 += bb + xlo2*ee;
    bsum3 += bb + xhi*ee;
  }

  if(esum==0 && prior==corr) {
    *int1 = *int2 = 0;
    return;
  }

  assert(esum>0);
  resum=1/esum;

  for(k=0;k<ngl;++k) {
    const double xr = xgl[k]*resum + xlo1;
    double t1 = logscale-bsum1, v1;
    for(i=0;i<nchan;++i)
      if(nobs[i]>0)
	t1 += nobs[i] * log( xr*chan[i].e + chan[i].b );
    sum1 += v1 = exp(lwgl[k]+t1);
    if(v1<DBL_EPSILON*sum1) break;
  }

  for(k=0;k<ngl;++k) {
    const double xr = xgl[k]*resum + xlo2;
    double t2 = logscale-bsum2, v2;
    for(i=0;i<nchan;++i)
      if(nobs[i]>0)
	t2 += nobs[i] * log( xr*chan[i].e + chan[i].b );
    sum2 += v2 = exp(lwgl[k]+t2);
    if(v2<DBL_EPSILON*sum2) break;
  }

  for(k=0;k<ngl;++k) {
    const double xr = xgl[k]*resum + xhi;
    double t3 = logscale-bsum3, v3;
    for(i=0;i<nchan;++i)
      if(nobs[i]>0)
	t3 += nobs[i] * log( xr*chan[i].e + chan[i].b );
    sum3 += v3 = exp(lwgl[k]+t3);
    if(v3<DBL_EPSILON*sum3) break;
  }

  if(prior==flat) {
    *int1 = (sum1-sum3)*resum;
    *int2 = (sum2-sum3)*resum;
  } else {
    *int1 = sum1-sum3;
    *int2 = sum2-sum3;
  }
  return;
}