package hep.aida.ref.function;

import hep.aida.IFitData;
import hep.aida.dev.IDevFitData;
import hep.aida.dev.IDevFitDataIterator;

/**
 *
 * @author The FreeHEP team @ SLAC
 *
 */
public class NonParametricFunction extends AbstractIFunction {
    
    private boolean mirrorLeft = false;
    private boolean mirrorRight = false;
    private double rho = 1;
    
    private int nPoints = 1000;
    private double[] lookupTable;
    
    private int nEvents;
    private double binWidth;
    
    private double[] dataPts;
    private double[] weights;
    
    double upperLimit = Double.NaN;
    double lowerLimit = Double.NaN;
    
    private static final double sqrt2pi = Math.sqrt(2*Math.PI);
    
    public NonParametricFunction(String title, IFitData data) {
        super(title, data.dimension(), 1);
        if ( data.dimension() != 1 )
            throw new IllegalArgumentException("Only one dimensional non-parametric functions are supported!");
        IDevFitDataIterator dataIter = ((IDevFitData) data).dataIterator();
        
        lookupTable = new double[ nPoints+1 ];
        
        nEvents = dataIter.entries();
        if (mirrorLeft) nEvents += dataIter.entries();
        if (mirrorRight) nEvents += dataIter.entries();
        
        //Calculate upper and lower limits
        dataIter.start();
        while( dataIter.next() ) {
            double x = dataIter.vars()[0];
            if ( Double.isNaN(upperLimit) || x > upperLimit )
                upperLimit = x;
            if ( Double.isNaN(lowerLimit) || x < lowerLimit )
                lowerLimit = x;
        }
        
        binWidth = (upperLimit - lowerLimit)/(nPoints-1);
        
        dataPts = new double[nEvents];
        weights = new double[nEvents];
        
        double entries = 0;
        double m = 0;
        double r = 0;
        
        int dataCount = 0;
        
        dataIter.start();
        while ( dataIter.next() ) {
            double x = dataIter.vars() [0];
            
            dataPts[dataCount] = x;
            
            entries++;
            m += x;
            r += x*x;
            
            dataCount++;
            
            if (mirrorLeft) {
                dataPts[dataCount]= 2*lowerLimit - x;
                dataCount++;
            }
            
            if (mirrorRight) {
                dataPts[dataCount]= 2*upperLimit - x;
                dataCount++;
            }
        }
        
        
        double mean = m/entries;
        double sigma = Math.sqrt(r/entries-mean*mean);
        double h = Math.pow(4./3.,0.2)*Math.pow(nEvents,-0.2)*rho;
        double hmin = h*sigma*Math.sqrt(2.)/10.;
        double norm = h*Math.sqrt(sigma)/(2.0*Math.sqrt(3.0));
        
        for(int j=0; jnPoints-1)
//            throw new IllegalArgumentException("Cannot have point above upper limit.");
            i=nPoints-1;
        
        double dx = (x -(lowerLimit + i*binWidth) )/binWidth;
        
        double val = (lookupTable[i] + dx*(lookupTable[i+1]-lookupTable[i]) );
        return val;
    }
    
    
    private double evaluateFull( double x ) {
        
        double y=0;
        
        for (int i=0;i