package montecarlo;

import java.util.Random;

/**
 * Collection of methods for performing numerical quadrature on a given multivariate function
 */

public class MultiVariateQuadrature {
    
    private MultiVariateFunction integrand; //function being integrated
    private double[] xMin, xMax;  //lower and upper limit of integration
    private double[][] x;         //x values used in quadrature  [no of points][dimension]
    private double[] X;           //x values for one term in quadrature sum
    private int nPoints;        //number of quadrature points (function evaluations)
    private int nPointsD;       //number of points for each dimension of integration
    private int jPoint;
    
    public MultiVariateQuadrature(MultiVariateFunction i) {
        setFunction(i);
        setN(20);
        setXMin(0.0);
        setXMax(1.0);
    }
    
    public void setFunction(MultiVariateFunction i) {integrand = i; setN(nPoints);}
    
    public void setN(int n) {
        nPointsD = (int)Math.pow((double)n,1.0/(double)integrand.D());
        nPoints = (int)Math.pow((double)nPointsD,(double)integrand.D());
        System.out.println(n + " " + nPointsD + " " + nPoints);
        x = new double[nPoints][integrand.D()];
        X = new double[integrand.D()];
    }
    public int getN() {return nPoints;}
    
    public double[][] getX() {return x;}
    
    public void setXMin(double xmin) {
        xMin = new double[integrand.D()];
        for(int i=0; i<x[0].length; i++) {xMin[i]=xmin;}
    }
    public double getXMin() {return xMin[0];}
    public void setXMax(double xmax) {
        xMax = new double[x[0].length];
        for(int i=0; i<x[0].length; i++) {xMax[i]=xmax;}
    }
    public double getXMax() {return xMax[0];}
    
    public double givenPoints(double[][] points) {  //computes integral using given points (assumes unweighted points in xMin, xMax values)
        double sum = 0.0;
        for(int i=0; i<points.length; i++) {
            sum += integrand.f(points[i]);
        }
        double volume = 1.0;
        for(int k=0; k<integrand.D(); k++) {volume *= (xMax[k]-xMin[k]);}
        return sum*volume/(double)points.length;
    }    
    public double rectangle(int n) {setN(n); return rectangle();}
    public double rectangle() {jPoint = 0; return rectangleRecurse(0);}
    private double rectangleRecurse(int k) {  //rectangle-rule quadrature formula
        double deltaX = (xMax[k]-xMin[k])/(double)nPointsD;
        double sum = 0.0;
        for(int i=0; i<nPointsD; i++) {
            X[k] = xMin[k] + (i+0.5)*deltaX;
            if(k == integrand.D()-1) { 
                sum += integrand.f(X);
                for(int m=0; m<=k; m++) x[jPoint][m] = X[m];
                jPoint++;
            }
            else {
                sum += rectangleRecurse(k+1);
            }
        }
        return sum*deltaX;
    }
 /*   public double trapezoid(int n) {setN(n); return trapezoid();}
    public double trapezoid() {  //trapezoid-rule quadrature formula
        double deltaX = (xMax-xMin)/((double)nPoints-1.0);
        double sum = 0.0;
        for(int i=0; i<nPoints; i++) {
            x[i] = xMin + i*deltaX;
            sum += (i==0 || i==nPoints-1) ? 0.5*integrand.f(x[i]) : integrand.f(x[i]);
        }
        return sum*deltaX;
    }
 */   
    public double simpleMC(int n) {setN(n); return simpleMC();}
    public double simpleMC() {  //unbiased Monte Carlo quadrature
        Random rand = new MyRandom();
        double sum = 0.0;
        for(int i=0; i<nPoints; i++) {
            for(int k=0; k<integrand.D(); k++) {
                x[i][k] = xMin[k] + (xMax[k]-xMin[k])*rand.nextDouble();
            }
            sum += integrand.f(x[i]);
        }
        double volume = 1.0;
        for(int k=0; k<integrand.D(); k++) {volume *= (xMax[k]-xMin[k]);}
        return sum*volume/(double)nPoints;
    }
/*    
    public double importanceMC(int n, Random rand, Function w) {setN(n); return importanceMC(rand, w);}
    public double importanceMC(Random rand, Function w) {
        double deltaX = (xMax-xMin);
        double sum = 0.0;
        for(int i=0; i<nPoints; i++) {
            x[i] = xMin + deltaX*rand.nextDouble();
            sum += integrand.f(x[i])/w.f(x[i]);
        }
        return sum*deltaX/(double)nPoints;
    }
    */
}