package org.seqcode.math.stats;

import Jama.Matrix;

/**
 * Lowess: This implementation of Lowess interpolation is taken from MIDAS (v2.22).
 * http://www.tm4.org/midas.html
 * 
 * References:
 * Quackenbush, J. Microarray data normalization and transformation. Nature Genetics. Vol.32 supplement pp496-501 (2002).
 * Yang, I.V. et al. Within the fold: assessing differential expression measures and reproducibility in microarray assays. Genome Biol. 3, research0062.1-0062.12 (2002).
 * Cleveland, W.S. Robust locally weighted regression and smoothing scatterplots. J. Amer. Stat. Assoc. 74, 829-836 (1979).
 *  
 * Copyright @ 2001-2003, The Institute for Genomic Research (TIGR).
 * All rights reserved.
 *
 * Lowess.java
 *
 * Created on December 10, 2001, 10:39 AM
 * @author  wliang
 * @version
 */
public class Lowess {
	private double bandwidthPct;
	private int windowSize;
	private static double[] xArray;
    private static double[] yArray;
    private static double[] yEst;
    private double[] xSub;
    private double[] ySub;
    private int subQueryPointIndex;
    private double[] weights;
    
    /**
     * Constructor: Note that xInput,yInput pairs should be sorted along the x-axis (increasing).
     * 
     * @param xInput
     * @param yInput
     * @param bandwidthPct : Reasonable values are 0.25 to 0.5
     */
    public Lowess(double[] xInput, double[] yInput, double bandwidthPct) {
        this.bandwidthPct = bandwidthPct;
        xArray = xInput;
        yArray = yInput;
        int totalDataCount = xArray.length;
        yEst = new double[totalDataCount];
        
        windowSize = new Double(bandwidthPct * totalDataCount).intValue();
        
        if (windowSize > 1){    //Issue Locfit2.0, 08-01-2002
            xSub = new double[windowSize];
            ySub = new double[windowSize];
            weights = new double[windowSize];
            
            //System.out.println("    ---- Lowess Bandwidth = " + windowSize);
            
            //Construct matrix
            Matrix X = new Matrix(windowSize, 2);
            Matrix WX = new Matrix(windowSize, 2);
            Matrix WY = new Matrix(windowSize, 1);
            
            for (int q = 0; q < totalDataCount; q++){
                
                getEstWindowSubArrayAndWeights(q, windowSize);
                
                //Construct matrix
                //Matrix X = new Matrix(windowSize, 2);
                //Matrix WX = new Matrix(windowSize, 2);
                //Matrix WY = new Matrix(windowSize, 1);
                
                for ( int row = 0; row < windowSize; row++ ) {
                    WX.set(row, 0, Math.sqrt(weights[row]));
                    WX.set(row, 1, Math.sqrt(weights[row]) * xSub[row]);
                    WY.set(row, 0, Math.sqrt(weights[row]) * ySub[row]);
                }
                
                //Apply first order estimation: local linear model
                Matrix WXT = WX.transpose();
                Matrix B = (WXT.times(WX)).inverse().times((WXT).times(WY));
                
                double B0 = (double)B.get(0, 0);
                double B1 = (double)B.get(1, 0);
                yEst[q] = B0 + B1 * xSub[subQueryPointIndex];
            }
        }else{ //windowSize = 0 or 1 //Issue Locfit2.0, 08-01-2002
            for (int q = 0; q < totalDataCount; q++){
                yEst[q] = 0;
            }
        }
    }
    
    /**
     * Estimate the y fit values at each x value
     * @param x
     * @return
     */
    public double[] estimateValues(double[] x){
    	double[] yVal = new double[x.length];
    	int totalDataCount = x.length;
    	if (windowSize > 1){ 
            xSub = new double[windowSize];
            ySub = new double[windowSize];
            weights = new double[windowSize];
            //Construct matrix
            Matrix WX = new Matrix(windowSize, 2);
            Matrix WY = new Matrix(windowSize, 1);
           
            for (int q = 0; q < totalDataCount; q++){
                getEstWindowSubArrayAndWeightsEV(x[q], windowSize);
                for ( int row = 0; row < windowSize; row++ ) {
                    WX.set(row, 0, Math.sqrt(weights[row]));
                    WX.set(row, 1, Math.sqrt(weights[row]) * xSub[row]);
                    WY.set(row, 0, Math.sqrt(weights[row]) * ySub[row]);
                }
                
                //Apply first order estimation: local linear model
                Matrix WXT = WX.transpose();
                Matrix B = (WXT.times(WX)).inverse().times((WXT).times(WY));
                
                double B0 = (double)B.get(0, 0);
                double B1 = (double)B.get(1, 0);
                yVal[q] = B0 + B1 * x[q];
            }
        }else{ //windowSize = 0 or 1 //Issue Locfit2.0, 08-01-2002
            for (int q = 0; q < totalDataCount; q++){
                yVal[q] = 0;
            }
        }
    	return yVal;
    }
    
    public double[] getXarray(){
        return xArray;
    }
    
    public double[] getYEst(){
        return yEst;
    }
    
    private void getEstWindowSubArrayAndWeights(int queryPointIndex, int window){
        
        int fullLength = xArray.length;
        int subStartIndex = 0;
        int windowSize = window;
        double maxDistance = 0;
        
        //Issue Locfit1.0, 08-01-2002
        //if ((queryPointIndex >= windowSize) && (queryPointIndex <= fullLength - windowSize)){//middle region data
        if ((queryPointIndex >= windowSize) && (queryPointIndex < fullLength - windowSize)){//middle region data
            for (int i = 0; i < windowSize; i++){
                if ( (xArray[queryPointIndex] - (xArray[queryPointIndex - windowSize + 1 + i])) <= (xArray[queryPointIndex + 1 + i] - xArray[queryPointIndex]) ){
                    subStartIndex = queryPointIndex - windowSize + 1 + i;
                    subQueryPointIndex = queryPointIndex - subStartIndex;
                    maxDistance = ( ( (xArray[queryPointIndex] - xArray[subStartIndex]) > (xArray[subStartIndex + windowSize - 1]) - xArray[queryPointIndex] ) ? (xArray[queryPointIndex] - xArray[subStartIndex]) : (xArray[subStartIndex + windowSize - 1] - xArray[queryPointIndex] ) );
                    break;
                }
            }
        }else if (queryPointIndex < windowSize){
            for (int i = 0; i < windowSize; i++){
                if ( (xArray[queryPointIndex] - xArray[i]) <= (xArray[windowSize + i] - xArray[queryPointIndex]) ){
                    subStartIndex = i;
                    subQueryPointIndex = queryPointIndex - subStartIndex;
                    maxDistance = ( (xArray[queryPointIndex] - xArray[subStartIndex]) > (xArray[subStartIndex + windowSize - 1] - xArray[queryPointIndex] ) ) ? (xArray[queryPointIndex] - xArray[subStartIndex]) : (xArray[subStartIndex + windowSize - 1] - xArray[queryPointIndex] );
                    break;
                }
            }
        }else{//right margin data
            for (int i = 0; i < windowSize; i++){
                if ( (xArray[fullLength - 1 - i] - xArray[queryPointIndex]) <= (xArray[queryPointIndex] - xArray[fullLength - windowSize - 1 - i]) ){
                    subStartIndex = fullLength - windowSize - i;
                    subQueryPointIndex = queryPointIndex - subStartIndex;
                    maxDistance = ( (xArray[queryPointIndex] - xArray[subStartIndex]) > (xArray[subStartIndex + windowSize - 1] - xArray[queryPointIndex] ) ) ? (xArray[queryPointIndex] - xArray[subStartIndex]) : (xArray[subStartIndex + windowSize - 1] - xArray[queryPointIndex] );
                    break;
                }
            }
        }
        
        System.arraycopy(xArray, subStartIndex, xSub, 0, windowSize);
        System.arraycopy(yArray, subStartIndex, ySub, 0, windowSize);
        
        //Calculate distances and find tricube weights
        for (int j = 0; j < windowSize; j ++) {
            double distance = Math.abs(xSub[subQueryPointIndex] - xSub[j]);
            weights[j] = (double)Math.pow(1.0 - Math.pow((distance / maxDistance), 3.0), 3.0);
        }
    }
    
    private void getEstWindowSubArrayAndWeightsEV(double queryX, int window){
        
        int fullLength = xArray.length;
        int subStartIndex = 0;
        int windowSize = window;
        double maxDistance = 0;
        
        int closestXPoint=0; double closestXDist=Double.MAX_VALUE; 
        for(int i=0; i<fullLength; i++){
        	if(Math.abs(xArray[i]-queryX)<closestXDist){
        		closestXDist = Math.abs(xArray[i]-queryX);
        		closestXPoint=i;
        	}
        }

        if ((closestXPoint >= windowSize) && (closestXPoint < fullLength - windowSize)){//middle region data
            for (int i = 0; i < windowSize; i++){
                if ( (xArray[closestXPoint] - (xArray[closestXPoint - windowSize + 1 + i])) <= (xArray[closestXPoint + 1 + i] - xArray[closestXPoint]) ){
                    subStartIndex = closestXPoint - windowSize + 1 + i;
                    subQueryPointIndex = closestXPoint - subStartIndex;
                    maxDistance = ( ( (queryX - xArray[subStartIndex]) > (xArray[subStartIndex + windowSize - 1]) - queryX ) ? (queryX - xArray[subStartIndex]) : (xArray[subStartIndex + windowSize - 1] - queryX ) );
                    break;
                }
            }
        }else if (closestXPoint < windowSize){
            for (int i = 0; i < windowSize; i++){
                if ( (xArray[closestXPoint] - xArray[i]) <= (xArray[windowSize + i] - xArray[closestXPoint]) ){
                    subStartIndex = i;
                    subQueryPointIndex = closestXPoint - subStartIndex;
                    maxDistance = ( (queryX - xArray[subStartIndex]) > (xArray[subStartIndex + windowSize - 1] - queryX ) ) ? (queryX - xArray[subStartIndex]) : (xArray[subStartIndex + windowSize - 1] - queryX );
                    break;
                }
            }
        }else{//right margin data
            for (int i = 0; i < windowSize; i++){
                if ( (xArray[fullLength - 1 - i] - xArray[closestXPoint]) <= (xArray[closestXPoint] - xArray[fullLength - windowSize - 1 - i]) ){
                    subStartIndex = fullLength - windowSize - i;
                    subQueryPointIndex = closestXPoint - subStartIndex;
                    maxDistance = ( (queryX - xArray[subStartIndex]) > (xArray[subStartIndex + windowSize - 1] - queryX ) ) ? (queryX - xArray[subStartIndex]) : (xArray[subStartIndex + windowSize - 1] - queryX );
                    break;
                }
            }
        }
        
        System.arraycopy(xArray, subStartIndex, xSub, 0, windowSize);
        System.arraycopy(yArray, subStartIndex, ySub, 0, windowSize);
        
        //Calculate distances and find tricube weights
        for (int j = 0; j < windowSize; j ++) {
            double distance = Math.abs(queryX - xSub[j]);
            weights[j] = (double)Math.pow(1.0 - Math.pow((distance / maxDistance), 3.0), 3.0);
        }
    }
}

