#include "nr3.h"
#include "gamma.h"  // for Gaussian quadrature
#include "gauss_wgts.h" // for Gaussian quadrature
#include "quadrature.h" // for Romberg quadrature
#include "interp_1d.h" // for Romberg quadrature
#include "romberg.h" // for Romberg
#include "chebyshev.h" // for Chebyshev

double PI=4.0*atan(1); // PI to machine precision
Doub E; // for landau integration, must be set before calling LandauIntegrand
Doub eps; // to vary integration precision, test for convergence
Doub normalization; // to set the Landau normalization.
Chebyshev *cheb1signalptr, *cheb1partptr; // used to be able to call previously-determined chebyshev structs in the folding functions A1part and A2part below 
int NORDER=100; // chebyshev order - can be increased if necessary
int NORDERFOLD=200; // chebyshev order for folded distributions - structure is not at -5 but at 0 - more internal points needed

inline double MINIMUM(double a, double b) { return (b<a) ? b : a; }; // MIN and MAX from NR3 returned floats, do never use them
inline double MAXIMUM(double a, double b) { return (b>a) ? b : a; };

Doub LandauIntegrand(Doub x) {
   // gives integrand in the Landau function
   if (x<1e-16) return sin(x*PI);
   double y=x*log(x)+(E-2.4)*x;
   if (y>100) return 0; // f(x) << e(-100)
   return exp(-y)*sin(x*PI);
}

Doub LandauIntegrand1overx(Doub x) {
   // gives integrand in the Landau function with x replaced by 1/x
   if (x<1e-6) return 0;
   x=1/x;
   double y=x*log(x)+(E-2.4)*x;
   if (y>100) return 0; // f(x) << e(-100)
   return exp(-y)*sin(x*PI)*x*x;
}

Doub Landau(Doub x) {
   // calculates the unnormalized landau function
   // integral runs from 0 to infinity. make coordinate transformation: y= y for x = 0 to Ep, 1/y for Ep to inf
   if (x<0) return 0;
   E=x;
   if (x<5) return qromb(LandauIntegrand,0,100,eps);
   else return qromb(LandauIntegrand,0,1/E,eps) + qromb(LandauIntegrand1overx,0.01,E,eps); // 
}

Doub Asignal(Doub x) {
   return 0.2*normalization*Landau(x/5); // A=5E, pdf(A)=1/5pdf(E=A/5);
}

Doub Anoise(Doub x) {
// gaussian-distributed with mean 0 and sigma=0.5
    if (fabs(x)>20) return 0; // Guard against underflow
    return sqrt(2/PI)*exp(-2*x*x); // sqrt(1/2pisigma^2)e^-(x^2/2sigm^2)
}

Doub NoiseInt(Doub x) {
// integral from -inf to x of normal distribution with mean 0 and sigma 0.5
   if (x<-10) return 0;
   if (x>10) return 1;
   return (1+erf(x*sqrt(2)))/2; //
}

Doub A1part(Doub Atot) {
// fold the signal with the noise. The integral over the noise runs
// from -10sigma = -5, to +10sigma = +5. this is within machine accuracy
// the full precision needed 
// 20-point Gauss-legendre is accurate enough, but you may check that 
// by hand by doing a 30-point integration and comparing the difference
   Doub sum=0;
   if (Atot<-5) return 0; // Asignal is always positive, and Anoise>-5
   if (Atot>255) return 0; // Asignal can be calculated to 250
   double lowlim=MAXIMUM(-5,Atot-249.9999999); //can calculate Asignal until 250
   double uplim=MINIMUM(Atot-1e-8,5); // integral from y=-5 to +5 when Atot>5, from -5 to Atot else (because Asignal(Atot-y) should be positive) 
   VecDoub xx(20),ww(20);
   gauleg(lowlim,uplim,xx,ww); // integrate noise from -5 to uplim and Asignal from Atot-uplim to Atot+5 
   for (int i=0; i<20; i++) sum+=ww[i]*Anoise(xx[i])* cheb1signalptr->eval(Atot-xx[i],NORDER);
   return sum;
}

Doub Fold2Part(Doub Atot) {
  // Fold 1 particle(signal+noise) with signal of 1 particle, to get response of 2 particles + noise
   // Atot = A1part+Asignal
   // Asignal runs from 0 to 250
   // Apart runs from -5 to 255
   // folded distribution is correct from -5 to 245 and slightly deviates
   // for the part above 245 
   // Folding: with Asignal(y)Aparticle(Atot-y)
   double upperlim=MINIMUM(Atot+4.9999999,250);
   if (upperlim<1e-8) return 0;
   VecDoub xx(100),ww(100);
   VecDoub x2(110),w2(110); // to test whether the precision is high enough
   gauleg(0,upperlim,xx,ww);
   gauleg(0,upperlim,x2,w2);
   double sum100=0; // calculate folded distributions with 100 and 110-point Gaussian quadrature
   double sum110=0;
   for (int i=0; i<100; i++) sum100+=ww[i]* cheb1signalptr->eval(xx[i],NORDER)* cheb1partptr->eval(Atot-xx[i],NORDERFOLD);
   for (int i=0; i<110; i++) sum110+=w2[i]* cheb1signalptr->eval(x2[i],NORDER)* cheb1partptr->eval(Atot-x2[i],NORDERFOLD);
   if (fabs(sum100/sum110-1)>1e-9) cerr << "Fold2part : precision below 1e-9 , " << Atot << " " << sum100 << " " << sum110 <<  " " << sum100/sum110-1 << endl;
   return sum110;
}        

int main() {
   eps=5e-15; // can be set to lower accuracy too
   // E can run towards infinity. you may cut it at a high value, e.g. 100,
   // but the integral is not saturated then.
   // we will make a chebyshev object describing the signal response in A from
   // 0 to 200, but we will check the integral and the normalization from 0 to
   // infinity
   // to check the integral, I use Gaussian quadrature and verify it from
   // 0 to 10 and from 10 to infinity.
   // I use Gaussian quadrature to calculate the integral from 0 to 10 and
   // from 10 to infinity.
   int NGAUSS=40;
   double previntegral=0;
   double totintegral=1;
   while (fabs(previntegral/totintegral-1)>200*eps) { // the integral over all landau values E has a precision less than 50 times the precision with which we do the landau integral for 1 value of E
	previntegral=totintegral;
        totintegral=0;
	NGAUSS+=10;
	VecDoub xx(NGAUSS),weight(NGAUSS);
  	gauleg(0,50,xx,weight);
        for (int i=0; i<NGAUSS; i++) totintegral+=weight[i]*Landau(xx[i]);
        cout << "NGRID " << NGAUSS << " int[0,50] " << totintegral;
  	gauleg(0,0.02,xx,weight); // make replacement t=1/x, dx/dt=-1/t^2
        for (int i=0; i<NGAUSS; i++) totintegral+=weight[i]/xx[i]/xx[i]*Landau(1/xx[i]);
	cout << "  total integral [0,inf] " << totintegral << " " << previntegral-totintegral << endl;
    }
    cout << " after integral part " << endl;
    normalization=1/totintegral;
        
    Chebyshev cheb1signal(Asignal,0,250,NORDER); // make a 200-order chebyshev in the range  A=5E, E=0,50, A=0-250
    Chebyshev cheb1int=cheb1signal.integral();
    cheb1signalptr = &cheb1signal; // for folding with noise in function A1part
    // verify accuracy, print it out
    cout << endl;
    cout << "test as a function of A signal the precision of 80, 90, 100-order chebyshev and print integral " << endl;
    for (double A=0; A<250.1; A+=1) {
        double trueval=Asignal(A);
        if (trueval<1e-16) trueval=1e-16;
	cout << A << " " << cheb1signal.eval(A,NORDER-20)/trueval-1 << " ";
	cout << cheb1signal.eval(A,NORDER-10)/trueval-1 << " ";
	cout << cheb1signal.eval(A,NORDER)/trueval-1 << " " << trueval << " " << cheb1int.eval(A,NORDER) << endl;
    }
    cout << endl << " next: folding " << endl;
    
    // create the chebyshev object where the signal is folded with the noise
    // since I want to use the cheb1signal here to speed up the calculation,
    // I define the function at this place
   // create a chebyshev object for it
   Chebyshev cheb1part(A1part,-5,250,NORDERFOLD);
   cheb1partptr = &cheb1part;
   cheb1int=cheb1part.integral(); // reset cheb1int
   
   // determine the treshold for A such that P(A_noise)<1e-8 : only 1e-8 percent   // of the noise distribution falls above it
   // use a simple bisection
   double xmin=2;
   double xmax=5;
   double xmid=(xmin+xmax)/2;
   double y=NoiseInt(xmid)+1e-8-1;
   while (fabs(y)>1e-13) {
      if (y>0) xmax=xmid;
      else xmin=xmid;
      xmid=(xmin+xmax)/2;
      y=NoiseInt(xmid)+1e-8-1;
   }
   cout.precision(12);
   cout << endl << "probability noise below A = " <<xmid << " equals " << NoiseInt(xmid) << " ; False detection rate equals " << 1e8*(1-NoiseInt(xmid)) << endl;
   cout << " False dismissal: probability particle below threshold A " << xmid << " equals " << cheb1int.eval(xmid,NORDERFOLD)-cheb1int.eval(-5,NORDERFOLD) << endl;
   

   // now, create a chebyshev object for 2-particle passage. The total signal is   // Atot = A_1signal + A_1particle ( draw twice from the signal and once from the noise contribution)
   
 
  Chebyshev cheb2part(Fold2Part,-5,250,NORDERFOLD);
  Chebyshev cheb2int = cheb2part.integral();
  
  // plot results
  for (double x=-5; x<250; x+=0.1) {
     cout << x << " " << Anoise(x) << " ";
     if (x>0) cout << cheb1signal.eval(x,NORDER);  else  cout << "0.0";
     cout << " " << cheb1part.eval(x,NORDERFOLD) << " " << cheb2part.eval(x,NORDERFOLD) << " " <<  NoiseInt(x) << " " << cheb1int.eval(x,NORDERFOLD) << " " << cheb2int.eval(x,NORDERFOLD) << endl;
  }
  cout << endl << endl;
  // check, where 1 particle and 2 particle distributions are equal
  xmin=10;
  xmax=50;
  xmid=(xmin+xmax)/2;
  y=cheb1part.eval(xmid,NORDERFOLD)/cheb2part.eval(xmid,NORDERFOLD)-1;
  while (fabs(y)>1e-8) {
     if (y>0) xmin=xmid; else xmax=xmid;
     xmid=(xmin+xmax)/2;
     y=cheb1part.eval(xmid,NORDERFOLD)/cheb2part.eval(xmid,NORDERFOLD)-1;
  }  
  cout << " Particle distributions pdfs 1 and 2 particles equal at " << xmid << endl;
  cout << " Amount of 1 particle below that threshold " << cheb1int.eval(xmid,NORDERFOLD) << endl;
  cout << " Amount of 2 particles below that threshold " << cheb2int.eval(xmid,NORDERFOLD) << endl;
       
  return 0;
}