//
// Whereas SimpsonDemo demonstrates adaptive uniform integration 
// by Simpson's rule, here I explore nonuniform sampling, with 
// integrand samples concentrated near the lower limit of the interval.
//
// My remapper is x=u^2 where u is uniformly sampled and x is integration axis.
// This puts ten percent of the samples into the bottom 1% of the interval.
// Here, I use a fixed integration inverval (0, 1). 
//
// M.Lampton UCB SSL 2009

public class NonuniformDemo
{
    static int ncalls = 0; 

    public static void main(String[] args)
    {
        double funcparms[] = {10.0};
        double ans = doNonuniform(funcparms);
        double err = Math.abs(1 - ans); 
        System.out.println(" Nonuniform err = "+err+";   ncalls = "+ncalls); 
    }



    static double func(double x, double parms[])
    // This is a negative exponential, normalized to unit integral.
    // Its parameter controls the concentration near the origin.
    {
        ncalls++; 
        double coef = parms[0]/(1 - Math.exp(-parms[0])); 
        return coef*Math.exp(-x*parms[0]);
    } 



    //-------Nonunform Simpson package starts here----------------

    static int    NMAX=22; 
    static int    NMIN=5;
    static double NTOL=1E-9;

    static double doNonuniform(double fp[])
    {
        int nnext = 1; 
        // The following initial guess is safer than zero:
        // when one end is big, it forbids an early exit.  
        double sumprev = 0.5*(func(0,fp)+func(1,fp));  

        double p=0, pprev=0;                // predictions
        for (int j=1; j<NMAX; j++)
        {
            double uh = 0.5/nnext;          // half u interval
            double uc = uh;                 // initial ucenter starting at 0 

            double sum = 0.0; 
            for (int k=0; k<nnext; k++)
            {
                double xc = uc*uc;          // xcenter
                double deltax = (uc+uh)*(uc+uh)-(uc-uh)*(uc-uh); 
                sum += func(xc,fp)*deltax; 
                uc += 2*uh; 
            }

            p = (4*sum-sumprev)/3; // prediction based on two most recent sums

            // System.out.println("sum = "+sum); 
            // System.out.println("  sumprev = "+sumprev); 
            // System.out.println("  predict = "+p); 
            // System.out.println("  preprev = "+pprev); 
            if ((j>NMIN) && (Math.abs(p-pprev) < NTOL*(1+Math.abs(pprev))))
              return p;
            sumprev = sum;    // could be subtler!
            pprev = p; 
            nnext *= 2; 
        }
        return -999; // failed to converge in JMAX iterations
    }        
}