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

#include "genlimit.h"

static double csint0(double xlo,double logscale,
		     int nchan,const int nobs[],const EB chan[],
		     int ngl,const double xgl[],const double lwgl[],
		     PRIOR prior);

/*  Joel Heinrich  8 April 2005

Function returns integral from xlo to infinity.

*uncertainty (if not null pointer) returned with uncertainty of
integral due to Monte Carlo statistical fluctuations of the prior
ensemble.

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


*/


double csint(double xlo,int nchan,int nens,const int nobs[],const EB* ens,
	     int* ngl,double xgl[],double lwgl[],
	     PRIOR prior,double* uncertainty) {
  int i,ntot=0;
  double sum=0, sum2=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) {
    const double t = csint0(xlo,logscale,nchan,nobs,p,*ngl,xgl,lwgl,prior);
    sum += t;
    sum2 += t*t;
    p+=nchan;
  }
  sum /= nens;
  sum2 /= nens;
  if(uncertainty)
    *uncertainty = (nens>1) ? sqrt((sum2-sum*sum)/(nens-1)) : 1.0 ;
  return sum;
}

static double csint0(double xlo,double logscale,
		     int nchan,const int nobs[],const EB chan[],
		     int ngl,const double xgl[],const double lwgl[],
		     PRIOR prior) {
  int i,k;
  double sum=0, esum=0, bsum=0, resum;

  for(i=0;i<nchan;++i) {
    esum += chan[i].e;
    bsum += chan[i].b + xlo*chan[i].e;
  }
  assert(esum>0);
  resum=1/esum;

  for(k=0;k<ngl;++k) {
    const double xr = xgl[k]*resum + xlo;
    double t = logscale-bsum, v;
    for(i=0;i<nchan;++i)
      if(nobs[i]>0)
	t += nobs[i] * log( xr*chan[i].e + chan[i].b );
    sum += v = exp(lwgl[k]+t);
    if(v<DBL_EPSILON*sum) break;
  }
  if(prior==flat)
    sum *= resum;
  return sum;
}