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

#define	Debug		0
#define	Static		static

/********************************************************
	Trapezoid Integral
********************************************************/

/*							*/
/*	integral of function f over [a,b]		*/
/*	function f should have the following form:	*/
/*		dbl f(dbl x)				*/
/*							*/

extern dbl TrapezoidIntegral(f, a, b, n)
dbl	(*f)();
dbl	a, b;
int	n;
{
  int	i;
  dbl	delta, x, sum;
  
  delta = (b-a)/n;
  sum = 0;
  for (i=0, x=a; i<n; i++, x+=delta) {
    sum += (*f)(x)*delta;
  }
  return sum;
}

/********************************************************
	Gauss Integral
********************************************************/

#define MaxPoints	1000

Static int	numpoints = 15;
Static dbl	point[MaxPoints] = {
	
	-0.987992518020485428489565718566,
	-0.937273392400705904307758947710,
	-0.848206583410427216200648320774,
	-0.724417731360170047416186054613,
	-0.570972172608538847537226737253,
	-0.394151347077563369897207370981,
	-0.201194093997434522300628303394,
	0.000000000000000000000000000000,
	0.201194093997434522300628303394,
	0.394151347077563369897207370981,
	0.570972172608538847537226737253,
	0.724417731360170047416186054613,
	0.848206583410427216200648320774,
	0.937273392400705904307758947710,
	0.987992518020485428489565718566
};

Static dbl	weight[MaxPoints] = {
	
	0.030753241996117268354628393577,
	0.070366047488108124709267416450,
	0.107159220467171935011869546685,
	0.139570677926154314447804794511,
	0.166269205816993933553200860481,
	0.186161000015562211026800561866,
	0.198431485327111576456118326443,
	0.202578241925561272880620199967,
	0.198431485327111576456118326443,
	0.186161000015562211026800561866,
	0.166269205816993933553200860481,
	0.139570677926154314447804794511,
	0.107159220467171935011869546685,
	0.070366047488108124709267416450,
	0.030753241996117268354628393577
};

/*						*/
/*	integral of function f over [-1,1]	*/
/*						*/

extern dbl GaussIntegralBase(f)
dbl	(*f)();
{
	int	i;
	dbl	sum;
	
	sum = 0;
	for (i=0; i<numpoints; i++) {
		sum += weight[i] * (*f)(point[i]);
	}
	return(sum);
}

/*							*/
/*	integral of function f over [a,b]		*/
/*	function f should have the following form:	*/
/*		dbl f(dbl x)				*/
/*							*/

extern dbl GaussIntegral(f, a, b)
dbl	(*f)();
dbl	a, b;
{
	int	i;
	dbl	p, q, sum;
	
	p = (b-a)/2.00;
	q = (b+a)/2.00;
	sum = 0;
	for (i=0; i<numpoints; i++) {
		sum += weight[i] * (*f)(p*point[i]+q);
	}
	return(p*sum);
}

/*							*/
/*	integral of function f with parameters		*/
/*		parameters: params[0...npar-1]		*/
/*		function f: f(x, npar, params)		*/
/*			dbl x, int npar, dbl params[]	*/
/*							*/

extern dbl GaussIntegralParams(f, a, b, npar, params)
dbl	(*f)();
dbl	a, b;
int	npar;
dbl	*params;
{
	int	i;
	dbl	p, q, sum;
	
	p = (b-a)/2.00;
	q = (b+a)/2.00;
	sum = 0;
	for (i=0; i<numpoints; i++) {
		sum += weight[i] * (*f)(p*point[i]+q, npar, params);
	}
	return(p*sum);
}

/*******************************
	multiple integral
********************************/

/*							*/
/*	integral of function f over [-1,1]\times[-1,1]	*/
/*							*/

extern dbl GaussIntegralBase2D(f)
dbl	(*f)();
{
	int	i, j;
	dbl	sum;
	
	sum = 0;
	for (i=0; i<numpoints; i++) {
		for (j=0; j<numpoints; j++) {
			sum += weight[i] * weight[j] * (*f)(point[i],point[j]);
		}
	}
	return(sum);
}

/*								*/
/*	integral of function f over [xa,xb]\times[ya,yb]	*/
/*	function f should have the following form:		*/
/*		dbl f(dbl x, dbl y)				*/
/*								*/

extern dbl GaussIntegral2D(f, xa, xb, ya, yb)
dbl	(*f)();
dbl	xa, xb, ya, yb;
{
	int	i, j;
	dbl	px, qx, py, qy, sum;
	
	px = (xb-xa)/2.00;	qx = (xb+xa)/2.00;
	py = (yb-ya)/2.00;	qy = (yb+ya)/2.00;
	sum = 0;
	for (i=0; i<numpoints; i++) {
		for (j=0; j<numpoints; j++) {
			sum += weight[i] * weight[j] * (*f)(px*point[i]+qx,py*point[j]+qy);
		}
	}
	return(px*py*sum);
}

extern dbl GaussIntegral2Dmarginal1(f, xa, xb, y)
dbl	(*f)();
dbl	xa, xb, y;
{
	int	i;
	dbl	px, qx, sum;
	
	px = (xb-xa)/2.00;	qx = (xb+xa)/2.00;
	sum = 0;
	for (i=0; i<numpoints; i++) {
		sum += weight[i] * (*f)(px*point[i]+qx, y);
	}
	return(px*sum);
}

extern dbl GaussIntegral2Dmarginal2(f, x, ya, yb)
dbl	(*f)();
dbl	x, ya, yb;
{
	int	j;
	dbl	py, qy, sum;
	
	py = (yb-ya)/2.00;	qy = (yb+ya)/2.00;
	sum = 0;
	for (j=0; j<numpoints; j++) {
		sum += weight[j] * (*f)(x, py*point[j]+qy);
	}
	return(py*sum);
}

/*								*/
/*	integral of function f with parameters			*/
/*		parameters: params[0...npar-1]			*/
/*		function f: f(x, y, npar, params)		*/
/*			dbl x, dbl y, int npar, dbl params[]	*/
/*								*/

extern dbl GaussIntegral2DParams(f, xa, xb, ya, yb, npar, params)
dbl	(*f)();
dbl	xa, xb, ya, yb;
int	npar;
dbl	*params;
{
	int	i, j;
	dbl	px, qx, py, qy, sum;
	
	px = (xb-xa)/2.00;	qx = (xb+xa)/2.00;
	py = (yb-ya)/2.00;	qy = (yb+ya)/2.00;
	sum = 0;
	for (i=0; i<numpoints; i++) {
		for (j=0; j<numpoints; j++) {
			sum += weight[i] * weight[j] *
			(*f)(px*point[i]+qx, py*point[j]+qy, npar, params);
		}
	}
	return(px*py*sum);
}

extern dbl GaussIntegral2Dmarginal1Params(f, xa, xb, y, npar, params)
dbl	(*f)();
dbl	xa, xb, y;
int	npar;
dbl	*params;
{
	int	i;
	dbl	px, qx, sum;
	
	px = (xb-xa)/2.00;	qx = (xb+xa)/2.00;
	sum = 0;
	for (i=0; i<numpoints; i++) {
		sum += weight[i] * (*f)(px*point[i]+qx, y, npar, params);
	}
	return(px*sum);
}

extern dbl GaussIntegral2Dmarginal2Params(f, x, ya, yb, npar, params)
dbl	(*f)();
dbl	x, ya, yb;
int	npar;
dbl	*params;
{
	int	j;
	dbl	py, qy, sum;
	
	py = (yb-ya)/2.00;	qy = (yb+ya)/2.00;
	sum = 0;
	for (j=0; j<numpoints; j++) {
		sum += weight[j] * (*f)(x, py*point[j]+qy, npar, params);
	}
	return(py*sum);
}



/*******************************
	Utilities
********************************/

/*							*/
/*	set parameters for Gauss integral from file	*/
/*							*/

extern void GaussIntegralInitialize(filename)
char	*filename;
{
	FILE	*fd;
	int	i, n;
	
	if ((fd = fopen(filename,"r")) == NULL) {
		fprintf(stderr, "cannot find %s\n", filename);
		exit(1);
	}
	
	fscanf(fd, "%d", &numpoints);
#if Debug
	fprintf(stderr, "Number of points = %d\n", numpoints);
#endif
	
	for (i=0; i<numpoints; i++) {
		fscanf(fd, "%d %lf", &n, point+i);
	}
	for (i=0; i<numpoints; i++) {
		fscanf(fd, "%d %lf", &n, weight+i);
	}
	
	fclose(fd);
#if Debug
	matrixprintformat("%l20.18f\n");
	vectorfprint(stderr, numpoints, point);
	vectorfprint(stderr, numpoints, weight);
#endif
}