package algorithmrepository;

import algorithmrepository.exceptions.NotImplementedException;

/**
 * Class for doing linear interpolation in 2D uniform
 * grids (bilinear interpolation).
 * 
 * Usage:
 * 
 * LinearInterpolation2D ip = new LinearInterpolation2D(xp,yp,fp);
 * ip.eval(3.1, 4.4)
 *
 * From: http://www.algorithmrepository.org
 *
 */
public class LinearInterpolation2D implements Interpolation2D {
    double[] xp, yp;
    double[][] fp;
    int numx, numy;
    double minx, maxx;
    double miny, maxy;
    double dx, dy, dxdy;
    double extrapValue = Double.NaN;


    /**
     * Initialises bilinear interpolation from a given
     * 2D grid of x-values, y-values and corresponding
     * function values. The grid has to be uniformly
     * sampled in both directions, and the arrays
     * monotonically increasing.
     *
     * @param xp The x-coordinates of the grid points
     * @param yp The y-coordinates of the grid points.
     * @param fp The function values at the grid points
     * fp[i][j] is the function value at the ith x-coordinate and jth y-coordinate.
     *
     */
    public LinearInterpolation2D(double[] xp, double[] yp, double[][] fp) {
        this(xp, yp, fp, Double.NaN);
    }
    
    /**
     * Initialises bilinear interpolation from a given
     * 2D grid of x-values, y-values and corresponding
     * function values. The grid has to be uniformly
     * sampled in both directions, and the arrays
     * monotonically increasing.
     *
     * @param xp The x-coordinates of the grid points
     * @param yp The y-coordinates of the grid points.
     * @param fp The function values at the grid points
     * fp[i][j] is the function value at the ith x-coordinate and jth y-coordinate.
     * @param extrapValue The value used for evaluation outside the grid.
     *
     */
    public LinearInterpolation2D(double[] xp, double[] yp, double[][] fp, double extrapValue) {
        
        this.dx = xp[1]-xp[0];
        this.dy = yp[1]-yp[0];
        double eps = 1e-6;
        for(int i=1;i<xp.length;++i) {
            double ratio = (xp[i]-xp[i-1])/dx; 
            if (ratio > 1+eps || ratio < 1-eps) {
                throw new RuntimeException("Interpolation grid in 2d must be monotonically increasing and uniformly spaced.");
            }
        }
        for(int i=1;i<yp.length;++i) {
            double ratio = (yp[i]-yp[i-1])/dy; 
            if (ratio > 1+eps || ratio < 1-eps) {
                throw new RuntimeException("Interpolation grid in 2d must be monotonically increasing and uniformly spaced.");
            }
        }
                
        this.xp = xp;
        this.yp = yp;
        this.fp = fp;
        this.numx = xp.length;
        this.numy = yp.length;
        this.minx = xp[0];
        this.maxx = xp[numx-1];
        this.miny = yp[0];
        this.maxy = yp[numy-1];       
        this.dxdy = dx*dy;
        this.extrapValue = extrapValue;
        
    }

    
    public void setF(double[][] fp) {
        this.fp = fp;
    }
    
    /** Sets the function, given a list of coordinates and function values
     * x[], y[], f[]
     * 
     * The x[]'s and y[]'s are assumed to be equal (or very near equal)
     * to actual grid points set when the class was initialised
     * 
     * Current algorithm is slow, I might upgrade it when I'm bored.
     * 
     * oliford
     * 
     */
    public void setFByList(double xList[], double yList[], double fList[], 
                        double replacementValue, double epsilon){
        setF(makeGridFromList(xp, yp, xList, yList, fList, replacementValue, epsilon));
        
    }
    
    public static double[][] makeGridFromList(double xp[], double yp[], double xList[], double yList[], double fList[], 
        double replacementValue, double epsilon){
        double[][] fmatrix = new double[xp.length][yp.length];
        
        double allX[] = new double[xp.length*yp.length];
        double allY[] = new double[xp.length*yp.length];
        for(int i=0;i<xp.length;++i) {          
            for(int j=0;j<yp.length;++j) {
                int index = i*yp.length+j;              
                allX[index] = xp[i];
                allY[index] = yp[j];                
                boolean foundGridPoint = false;
                for(int k=0;k<xList.length;++k) {
                    if (Math.abs(xList[k]-xp[i]) < epsilon && Math.abs(yList[k]-yp[j]) < epsilon) {
                        fmatrix[i][j] = fList[k]; 
                        foundGridPoint = true;
                        break;
                    }
                }
                if (!foundGridPoint) {
                    fmatrix[i][j] = replacementValue; 
                }
            }
        }
               
        return fmatrix;       
    }


    /**
     *
     * Does a bilinear interpolation on the uniform grid given in the
     * constructor.
     *
     * Extrapolation: Doesn't do extrapolation, returns -Inf if outside grid.
     *
     * @param x The x-coordinates where the function should be interpolated
     * @param y The y-coordinates where the function should be interpolated
     * @return The interpolated values of the function.
     */
    public double[] eval(double[] x, double[] y) {
        double[] ret = new double[x.length];
        for(int i=0;i<x.length;++i) {
            ret[i] = eval(x[i],y[i]);
        }
        return ret;
    }
    
    /**
     * Evaluates the function on a 2D grid with given knots in the x and y-dimension.
     * 
     * @param xKnots The knots along the x dimension.
     * @param yKnots The knots along the y dimension.
     * @return ret[i][j] is the function value for know xi, yj.
     */
    public double[][] evalGrid(double[] xKnots, double[] yKnots) {
        double[][] ret = new double[xKnots.length][yKnots.length];
        
        for(int i=0;i<ret.length;++i) {
            for(int j=0;j<ret[0].length;++j) {
                ret[i][j] = eval(xKnots[i], yKnots[j]);
            }
        }
        
        return ret;
    }
    

    @Override
    public double eval(double x, double y, double[] der) {
       throw new NotImplementedException();
    }

    @Override
    public double[] eval(double[] x, double[] y, double[][] der) {
        throw new NotImplementedException();
    }

    /**
    *
    * Does a single bilinear interpolation on the uniform grid given in the
    * constructor.
    *
    * Extrapolation: Doesn't do extrapolation, returns -Inf if outside grid.
    *
    * @param x The x-coordinate where the function should be interpolated
    * @param y The y-coordinate where the function should be interpolated
    * @return The interpolated value of the function.
    */    
    public double eval(double x, double y) {
        int ix, iy;
        
        //ix = (int) Math.floor((numx-1)*(x-minx)/(maxx-minx));
        //iy = (int) Math.floor((numy-1)*(y-miny)/(maxy-miny));

        
        //maybe it was faster to do this after the ix= and iy=, in
        // order to do the tests in integer arithmetic
        // but then we have to use Math.floor() rather than just (int)
        if(x < minx || y < miny || x > maxx || y > maxy)
            return extrapValue;
        
        ix = (int)((numx-1)*(x-minx)/(maxx-minx));
        iy = (int)((numy-1)*(y-miny)/(maxy-miny));

        //we should be able to hit the end exactly, but no futher
        if (x == maxx) ix -= 1;
        if (y == maxy) iy -= 1;

        
         return 1/dxdy*(fp[ix][iy]*(xp[ix+1]-x)*(yp[iy+1]-y) +
                fp[ix+1][iy]*(x-xp[ix])*(yp[iy+1]-y) +
                fp[ix+1][iy+1]*(x-xp[ix])*(y-yp[iy]) +
                fp[ix][iy+1]*(xp[ix+1]-x)*(y-yp[iy]));

    }
    
    
    /**
     * Static method exposing the core functionality of bilinear interpolation.
     * 
     * Implements a single 2D bilinear interpolation. The bilinear interpolation
     * uses the function value at 2x2 grid points.
     *
     * @param x1 left x-coordinate
     * @param x2 right x-coordinate
     * @param y1 lower y-coordinate
     * @param y2 upper y-coordinate
     * @param fp The four function values corresponding to the given
     * grid points. The order is anticlockwise: lower left, lower right,
     * upper right, upper left.
     * @param x The x-coordinate of the point that should be interpolated.
     * @param y The y-coordinate of the point that should be interpolated.
     *
     * @return The interpolated value of the function at the given point.
     */
    public static double eval(double x1, double x2, double y1, double y2, double[] fp, double x, double y) {
        double dxdy = (x2-x1)*(y2-y1);
        double f = 1/dxdy*(fp[0]*(x2-x)*(y2-y) + fp[1]*(x-x1)*(y2-y) + fp[2]*(x-x1)*(y-y1) + fp[3]*(x2-x)*(y-y1));

        return f;
    }
    
    /**
     * Returns the vector of x-coordinates from which the grid is 
     * defined.
     * 
     * @return The vector of x-coordinates defining the grid division
     * in the x-direction.
     */
    public double[] getxp() {
        return xp;
    }

    
    /**
     * Returns the vector of y-coordinates from which the grid is 
     * defined.
     * 
     * @return The vector of y-coordinates defining the grid division
     * in the y-direction.
     */    
    public double[] getyp() {
        return yp;
    }
    
    
    /**
     * Returns the grid values. 
     * 
     * @return The matrix of grid values given in the constructor.
     * ret[i][j] is the function value corresponding to the grid cell with
     * lower left coordinates x[i],y[j].
     */    
    public double[][] getfp() {
        return fp;
    }


}