//#include "header_cpp.h"

//#include "/home/jpei/bin/grishin/random.c"

void nrerror(char error_text[]){

  fprintf(stderr, "%s\n", error_text);
  exit(0);

}

double erff(double x)
/* erff(x)=2/(sqrt(Pi))\Int_0^x exp(-t^2) dt */
{
        double gammp(double a, double x);

        return x < 0.0 ? -gammp(0.5,x*x) : gammp(0.5,x*x);
}

double gammp(double a, double x)

/* returns incomplete gamma-function gammp(a,x)=1/Gamma(a)*\Int_0^{+\infinity}exp(-t)*t^(a-1)dt */
{
        void gcf(double *gammcf, double a, double x, double *gln);
        void gser(double *gamser, double a, double x, double *gln);
        void nrerror(char error_text[]);
        double gamser,gammcf,gln;

        if (x < 0.0 || a <= 0.0) nrerror("Invalid arguments in routine gammp");        if (x < (a+1.0)) {
                gser(&gamser,a,x,&gln);
                return gamser;
        } else {
                gcf(&gammcf,a,x,&gln);
                return 1.0-gammcf;
        }
}

#define ITMAX 200
#define EPS 3.0e-10

void gser(double *gamser, double a, double x, double *gln)
/* used for inclomplete gamma function calculation */
{
        double gammln(double xx);
        void nrerror(char error_text[]);
        int n;
        double sum,del,ap;

        *gln=gammln(a);
        if (x <= 0.0) {
                if (x < 0.0) nrerror("x less than 0 in routine gser");
                *gamser=0.0;
                return;
        } else {
                ap=a;
                del=sum=1.0/a;
                for (n=1;n<=ITMAX;n++) {
                        ++ap;
                        del *= x/ap;
                        sum += del;
                        if (fabs(del) < fabs(sum)*EPS) {
                                *gamser=sum*exp(-x+a*log(x)-(*gln));
                                return;
                        }
                }
                nrerror("a too large, ITMAX too small in routine gser");
                return;
        }
}
#undef ITMAX
#undef EPS

#define ITMAX 200
#define EPS 3.0e-10
#define FPMIN 1.0e-30

void gcf(double *gammcf, double a, double x, double *gln)
/* used for incomplete gamma function calculation */
{
        double gammln(double xx);
        void nrerror(char error_text[]);
        int i;
        double an,b,c,d,del,h;

        *gln=gammln(a);
        b=x+1.0-a;
        c=1.0/FPMIN;
        d=1.0/b;
        h=d;
        for (i=1;i<=ITMAX;i++) {
                an = -i*(i-a);
                b += 2.0;
                d=an*d+b;
                if (fabs(d) < FPMIN) d=FPMIN;
                c=b+an/c;
                if (fabs(c) < FPMIN) c=FPMIN;
                d=1.0/d;
                del=d*c;
                h *= del;
                if (fabs(del-1.0) < EPS) break;
        }
        if (i > ITMAX) nrerror("a too large, ITMAX too small in gcf");
        *gammcf=exp(-x+a*log(x)-(*gln))*h;
}
#undef ITMAX
#undef EPS
#undef FPMIN




double KS_integral(double *x, int N, int sflag) {

/* Kolmogorov-Smirnov Integral.
   Given an array of x-values of N elements, first sort the array
   and then calculate the integral of cumulative distribution function
   x: array of N values, x[0] to x[N-1]
   N: number of elements
   sflag: sorted array (0), unsorted (1)
*/

  void sort(unsigned long n, double arr[]);

  int i,j,k;
  double y_interval = 1.0/N;
  double delta_x;
  double value_y;
  double integral = 0;

  if( sflag) {
	/* sort the array */
	sort(N, x-1);
	/*for(i=0;i<N;i++) {
	  fprintf(stderr, "%f\n", x[i]);
	}*/
	/* to be added */
  }

  for(i=1;i<N;i++) {
    delta_x = x[i] - x[i-1];
    if(delta_x <0) {
	fprintf(stderr, "Not an sorted increasing array \n");
	exit(0);
    }
    value_y = y_interval * i;
    integral += value_y * delta_x;
  } 

  return integral;

}

double KS_integral_weighted(double *x, double *w, int N, int sflag) {

/* Kolmogorov-Smirnov Integral: weighted array.
   Given an array of x-values of N elements and the array of their weight, 
   first sort the array and then calculate the integral of cumulative 
   distribution function
   x: array of N values, x[0] to x[N-1]
   w: array of N weights associated with x
   N: number of elements
   sflag: sorted array (0), unsorted (1)
*/

  void sort2(unsigned long n, double arr[], double brr[]);

  int i,j,k;
  double y_interval = 1.0/N;
  double delta_x;
  double value_y;
  double integral = 0;
  double nw[N]; //normalized weights
  double sum_of_w = 0;

  if( sflag) {
        /* sort the array */
        sort2(N, x-1, w-1);
        /*for(i=0;i<N;i++) {
          fprintf(stderr, "%f\n", x[i]);
        }*/
        /* to be added */
  }

  /* normalize the weights */

  for(i=0;i<N;i++) {nw[i] = 0;sum_of_w += w[i];
	if(sum_of_w==0) {fprintf(stderr,"sum_of_w: 0\n");exit(0); }
  }
  for(i=0;i<N;i++) 
    nw[i] = w[i]/sum_of_w;
  
  
  value_y = 0;
  for(i=1;i<N;i++) {
    delta_x = x[i] - x[i-1];
    if(delta_x <0) {
        fprintf(stderr, "Not an sorted increasing array \n");
        exit(0);
    }
    value_y += nw[i-1];
    integral += value_y * delta_x;
  }

  return integral;

}

double Gaussian_integral_x1_x2 (double x1, double x2, double mean, double sigma){

/* Given a Gaussian function with mean and variance (mean and sigma), and two boundary values x1 and x2, calculate the cumulative integral of the distribution function.
*/

double Gaussian_distribution(double x, double mean, double sigma);
double qsimp3(double (*func)(double,double, double), double param1, double param2, double a, double b);

double   Z = (x2 - mean)/sigma;
if(Z < -5) return 0;
return(qsimp3(&Gaussian_distribution, mean, sigma, x1, x2));

}


double Gaussian_distribution(double x, double mean, double sigma) {

  double x_normalized;

  x_normalized = (x-mean)/sigma;
  //fprintf(stderr, "x: %5.3f;  %5.3f\n", x, 0.5*(1+erff(x_normalized)) );
  return(0.5*(1+erff(x_normalized)));

}

#define EPS 1.0e-7
#define JMAX 30

double qsimp3(double (*func)(double,double, double), double param1, double param2, double a, double b)
/* integrates function func(param,x) by x from a to b, presision EPS */
{
double trapzd3(double (*func)(double,double, double),double param1, double param2, double a, double b,
int n);
        void nrerror(char error_text[]);
        int j;
        double s,st,ost,os;

        ost = os = -1.0e30;
        for (j=1;j<=JMAX;j++) {
                st=trapzd3(func,param1, param2,a,b,j);
                s=(4.0*st-ost)/3.0;
                if (fabs(s-os) < EPS*fabs(os)) return s;
                os=s;
                ost=st;
        }
        nrerror("Too many steps in routine qsimp");
        return 0.0;
}
#undef EPS
#undef JMAX

double trapzd3(double (*func)(double,double, double),double param1, double param2,double a, double b,
int n)
/* used for integration */
{
        double x,tnm,sum,del;
        static double s;
        int it,j;

        if (n == 1) {
                return (s=0.5*(b-a)*((*func)(a,param1, param2)+(*func)(b,param1, param2)));
        } else {
                for (it=1,j=1;j<n-1;j++) it <<= 1;
                tnm=it;
                del=(b-a)/tnm;
                x=a+0.5*del;
                for (sum=0.0,j=1;j<=it;j++,x+=del) sum += (*func)(x,param1,param2);
                s=0.5*(s+(b-a)*sum/tnm);
                return s;
        }
}

double gammln(double xx)

/* returns ln(Gamma(xx)) xx>0 */
{
        double x,y,tmp,ser;
        static double cof[6]={76.18009172947146,-86.50532032941677,
                24.01409824083091,-1.231739572450155,
                0.1208650973866179e-2,-0.5395239384953e-5};
        int j;

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

#define M 7
#define NSTACK 50
#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}

void sort(unsigned long n, double arr[])
/* sorts array arr[1.n] into ascending order */
{
        unsigned long i,ir=n,j,k,l=1;
        int jstack=0;
        double a,temp;

        int istack[NSTACK+1];
        for (;;) {
                if (ir-l < M) {
                        for (j=l+1;j<=ir;j++) {
                                a=arr[j];
                                for (i=j-1;i>=1;i--) {
                                        if (arr[i] <= a) break;
                                        arr[i+1]=arr[i];
                                }
                                arr[i+1]=a;
                        }
                        if (jstack == 0) break;
                        ir=istack[jstack--];
                        l=istack[jstack--];
                } else {
                        k=(l+ir) >> 1;
                        SWAP(arr[k],arr[l+1])
                        if (arr[l+1] > arr[ir]) {
                                SWAP(arr[l+1],arr[ir])
                        }
                        if (arr[l] > arr[ir]) {
                                SWAP(arr[l],arr[ir])
                        }
                        if (arr[l+1] > arr[l]) {
                                SWAP(arr[l+1],arr[l])
                        }
                        i=l+1;
                        j=ir;
                        a=arr[l];
                        for (;;) {
                                do i++; while (arr[i] < a);
                                do j--; while (arr[j] > a);
                                if (j < i) break;
                                SWAP(arr[i],arr[j]);
                        }
                        arr[l]=arr[j];
                        arr[j]=a;
                        jstack += 2;
                        if (jstack > NSTACK) {fprintf(stderr, "NSTACK too small in sort."); exit(0); }
                        if (ir-i+1 >= j-l) {
                                istack[jstack]=ir;
                                istack[jstack-1]=i;
                                ir=j-1;
                        } else {
                                istack[jstack]=j-1;
                                istack[jstack-1]=l;
                                l=i;
                        }
                }
        }
        /*free_ivector(istack,1,NSTACK);*/
}
#undef M
#undef NSTACK
#undef SWAP


#define NRANSI
#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
#define M 7
#define NSTACK 50

void sort2(unsigned long n, double arr[], double brr[])
{
        unsigned long i,ir=n,j,k,l=1;
        int jstack=0;
        double a,b,temp;

	int istack[NSTACK+1];
        for (;;) {
                if (ir-l < M) {
                        for (j=l+1;j<=ir;j++) {
                                a=arr[j];
                                b=brr[j];
                                for (i=j-1;i>=1;i--) {
                                        if (arr[i] <= a) break;
                                        arr[i+1]=arr[i];
                                        brr[i+1]=brr[i];
                                }
                                arr[i+1]=a;
                                brr[i+1]=b;
                        }
                        if (!jstack) {
                                //free_ivector(istack,1,NSTACK);
                                return;
                        }
                        ir=istack[jstack];
                        l=istack[jstack-1];
                        jstack -= 2;
                } else {
                        k=(l+ir) >> 1;
                        SWAP(arr[k],arr[l+1])
                        SWAP(brr[k],brr[l+1])
                        if (arr[l+1] > arr[ir]) {
                                SWAP(arr[l+1],arr[ir])
                                SWAP(brr[l+1],brr[ir])
                        }
                        if (arr[l] > arr[ir]) {
                                SWAP(arr[l],arr[ir])
                                SWAP(brr[l],brr[ir])
                        }
                        if (arr[l+1] > arr[l]) {
                                SWAP(arr[l+1],arr[l])
                                SWAP(brr[l+1],brr[l])
                        }
                        i=l+1;
                        j=ir;
                        a=arr[l];
                        b=brr[l];
                        for (;;) {
                                do i++; while (arr[i] < a);
                                do j--; while (arr[j] > a);
                                if (j < i) break;
                                SWAP(arr[i],arr[j])
                                SWAP(brr[i],brr[j])
                        }
                        arr[l]=arr[j];
                        arr[j]=a;
                        brr[l]=brr[j];
                        brr[j]=b;
                        jstack += 2;
                        if (jstack > NSTACK) {fprintf(stderr, "NSTACK too small in sort2."); exit(0); }
                        if (ir-i+1 >= j-l) {
                                istack[jstack]=ir;
                                istack[jstack-1]=i;
                                ir=j-1;
                        } else {
                                istack[jstack]=j-1;
                                istack[jstack-1]=l;
                                l=i;
                        }
                }
        }
}
#undef M
#undef NSTACK
#undef SWAP
#undef NRANSI

long timeseed(){

/* used to generate a negative long int as a seed for ran1 and ran2 from current time */
long seed;
time_t ttime;
time(&ttime);
seed=921303200-ttime;
if(seed>0)seed=-seed;if(seed==0){fprintf(stderr, "Something is wrong with time");exit(0);}
//fprintf(stderr,"The seed is %ld\n",seed);
return seed;
}


#define MBIG 1000000000
#define MSEED 161803398
#define MZ 0
#define FAC (1.0/MBIG)

double ran3(long *idum)
/* returns uniform random deviate between 0.0 and 1.0
idum - is a negative integer, do not alter it for successive calls*/

{
        static int inext,inextp;
        static long ma[56];
        static int iff=0;
        long mj,mk;
        int i,ii,k;

        if (*idum < 0 || iff == 0) {
                iff=1;
                mj=MSEED-(*idum < 0 ? -*idum : *idum);
                mj %= MBIG;
                ma[55]=mj;
                mk=1;
                for (i=1;i<=54;i++) {
                        ii=(21*i) % 55;
                        ma[ii]=mk;
                        mk=mj-mk;
                        if (mk < MZ) mk += MBIG;
                        mj=ma[ii];
                }
                for (k=1;k<=4;k++)
                        for (i=1;i<=55;i++) {
                                ma[i] -= ma[1+(i+30) % 55];
                                if (ma[i] < MZ) ma[i] += MBIG;
                        }
                inext=0;
                inextp=31;
                *idum=1;
        }
        if (++inext == 56) inext=1;
        if (++inextp == 56) inextp=1;
        mj=ma[inext]-ma[inextp];
        if (mj < MZ) mj += MBIG;
        ma[inext]=mj;
        return mj*FAC;
}
#undef MBIG
#undef MSEED
#undef MZ
#undef FAC

void avevar(double data[], unsigned long n, double *ave, double *var)
{
        unsigned long j;
        double s,ep;

        for (*ave=0.0,j=1;j<=n;j++) *ave += data[j];
        *ave /= n;
        *var=ep=0.0;
        for (j=1;j<=n;j++) {
                s=data[j]-(*ave);
                ep += s;
                *var += s*s;
        }
        *var=(*var-ep*ep/n)/(n-1);
}


double betai(double a, double b, double x)
{
        double betacf(double a, double b, double x);
        double gammln(double xx);
        double bt;

        if (x < 0.0 || x > 1.0) { fprintf(stderr, "Bad x in routine betai"); exit(0); }
        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;
}


#define MAXIT 100
#define EPS 3.0e-7
#define FPMIN 1.0e-30

double betacf(double a, double b, double x)
{
        int m,m2;
        double aa,c,d,del,h,qab,qam,qap;

        qab=a+b;
        qap=a+1.0;
        qam=a-1.0;
        c=1.0;
        d=1.0-qab*x/qap;
        if (fabs(d) < FPMIN) d=FPMIN;
        d=1.0/d;
        h=d;
        for (m=1;m<=MAXIT;m++) {
                m2=2*m;
                aa=m*(b-m)*x/((qam+m2)*(a+m2));
                d=1.0+aa*d;
                if (fabs(d) < FPMIN) d=FPMIN;
                c=1.0+aa/c;
                if (fabs(c) < FPMIN) c=FPMIN;
                d=1.0/d;
                h *= d*c;
                aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
                d=1.0+aa*d;
                if (fabs(d) < FPMIN) d=FPMIN;
                c=1.0+aa/c;
                if (fabs(c) < FPMIN) c=FPMIN;
                d=1.0/d;
                del=d*c;
                h *= del;
                if (fabs(del-1.0) < EPS) break;
        }
        if (m > MAXIT) nrerror("a or b too big, or MAXIT too small in betacf");
        return h;
}
#undef MAXIT
#undef EPS
#undef FPMIN


static double sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
void tutest(double data1[], unsigned long n1, double data2[], unsigned long n2,
        double *t, double *prob)
{
        void avevar(double data[], unsigned long n, double *ave, double *var);
        double betai(double a, double b, double x);
        double var1,var2,df,ave1,ave2;

        avevar(data1,n1,&ave1,&var1);
        avevar(data2,n2,&ave2,&var2);
        *t=(ave1-ave2)/sqrt(var1/n1+var2/n2);
        df=SQR(var1/n1+var2/n2)/(SQR(var1/n1)/(n1-1)+SQR(var2/n2)/(n2-1));
        *prob=betai(0.5*df,0.5,df/(df+SQR(*t)));
}

void tutest1(double data1[], unsigned long n1, double Ave2, double Var2, unsigned long n2, double *t, double *prob) {

/* given an array and the mean and variance of another array */

   	void avevar(double data[], unsigned long n, double *ave, double *var);
        double betai(double a, double b, double x);
	double var1,var2,df,ave1,ave2;
	
	avevar(data1,n1,&ave1,&var1);
	ave2 = Ave2;
	var2 = Var2;
	*t=(ave1-ave2)/sqrt(var1/n1+var2/n2);
	df=SQR(var1/n1+var2/n2)/(SQR(var1/n1)/(n1-1)+SQR(var2/n2)/(n2-1));
	*prob=betai(0.5*df,0.5,df/(df+SQR(*t)));
}


void moment(double data[], int n, double *ave, double *adev, double *sdev,
        double *var, double *skew, double *curt)
{
        void nrerror(char error_text[]);
        int j;
        double ep=0.0,s,p;

        if (n <= 1) nrerror("n must be at least 2 in moment");
        s=0.0;
        for (j=1;j<=n;j++) s += data[j];
        *ave=s/n;
        *adev=(*var)=(*skew)=(*curt)=0.0;
        for (j=1;j<=n;j++) {
                *adev += fabs(s=data[j]-(*ave));
                *var += (p=s*s);
                *skew += (p *= s);
                *curt += (p *= s);
        }
        *adev /= n;
        *var=(*var-ep*ep/n)/(n-1);
        *sdev=sqrt(*var);
        if (*var) {
                *skew /= (n*(*var)*(*sdev));
                *curt=(*curt)/(n*(*var)*(*var))-3.0;
        } else nrerror("No skew/kurtosis when variance = 0 (in moment)");
}

