/*
 *  Draco Econometrics
 *  Copyright (C) 2006-2007 Approximatrix, LLC
 *  http://www.approximatrix.com
 *  
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package com.approximatrix.draco.fit;


import java.util.List;

import org.netlib.blas.DDOT;
import org.netlib.blas.DGEMM;
import org.netlib.blas.DSCAL;
import org.netlib.blas.DAXPY;
import org.netlib.lapack.DGETRF;
import org.netlib.lapack.DGETRI;
import org.netlib.util.intW;

import com.approximatrix.draco.datamodel.StatTableModel;
import com.approximatrix.draco.fit.functions.*;
import com.approximatrix.draco.fit.statistics.ExtendedStatistic;
import com.approximatrix.draco.io.XMLio.AnalysisLoadStore;
import com.approximatrix.draco.math.MatrixSupport;
import com.approximatrix.draco.math.VectorSupport;
import com.approximatrix.draco.utility.ResultsUtilities;

public class SimpleMLEFit implements LSIterative {

	public static final int TECHNIQUE_SECANT = 1;
	
	private static final double INITIAL_SLOPE_PERCENTAGE = 0.05;
	
	private LikelihoodFunction evaluator = null;
	
	private StatTableModel model = null;
	private String dependent_var = null;
	private String[] independent_var = null;
	
	private double relaxation = 1.0;
	private String failure_message = null;
	private int iteration = 0;
	private double change = Double.POSITIVE_INFINITY;
	
	private double[] parameters = null;
	private double[] old_parameters = null;
	private double[][] covariance_matrix = null;
	
	private int technique = TECHNIQUE_SECANT;
	
	private boolean success = true;
	
	private int[] use_row = null;
	
	private String results = null;
	
	/** The list of extended statistics to compute */
	private List<ExtendedStatistic> estats = null;
	
	/** Sets all data in the fit to the information loaded from
	 * a file
	 * 
	 * @param data a valid AnalysisLoadStore object that has been properly filled
	 */
	public void setLoadedAnalysis(AnalysisLoadStore data, StatTableModel model) {
		independent_var = data.getIndependentVariables();
		dependent_var = data.getDependentVariable();
		results = data.getResults();
		parameters = data.getCoef();
		covariance_matrix = data.getCovarianceMatrix();
		
		if(data.getAnalysisSubType().compareTo(LogitLikelihood.class.getName())==0)
			evaluator = new LogitLikelihood();
		else if(data.getAnalysisSubType().compareTo(ProbitLikelihood.class.getName())==0)
			evaluator = new ProbitLikelihood();
		
		success = true;
	}
	
	public SimpleMLEFit() {
		// Nothing is done here - should only be called before a call to
		// setLoadedAnalysis.
	}
	
	public SimpleMLEFit(StatTableModel model,
			      String dependent,
			      String[] independent,
			      LikelihoodFunction function,
			      int technique) {
		
		this.model = model;
		this.evaluator = function;
		this.dependent_var = dependent;
		this.independent_var = independent;
		this.technique = technique;
	
		use_row = determineUsable();
		
		parameters = function.initializeParameters(model.getDataDouble(use_row,model.getColumnNumber(dependent_var)),
				                                   model.getDataDouble(use_row,model.getColumnNumber(independent_var)));
		
		if(parameters == null) {
			success = false;
			failure_message = "Failed to generate initial parameter estimates.";
		} 
		
	}
	
	public SimpleMLEFit(StatTableModel model,
			String dependent,
			String[] independent,
			LikelihoodFunction function) {
		
		this.model = model;
		this.evaluator = function;
		this.dependent_var = dependent;
		this.independent_var = independent;
		
		use_row = determineUsable();
		
		parameters = function.initializeParameters(model.getDataDouble(use_row,model.getColumnNumber(dependent_var)),
				model.getDataDouble(use_row,model.getColumnNumber(independent_var)));
		
		if(parameters == null) {
			success = false;
			failure_message = "Failed to generate initial parameter estimates.";
		}
		
	}
	
	public SimpleMLEFit(SimpleMLEFit last_iteration) {
		model = last_iteration.getModel();
		dependent_var = last_iteration.getDependentVariable();
		independent_var = last_iteration.getIndependentVariables();
		
		iteration = last_iteration.getCurrentIteration()+1;
		
		technique = last_iteration.getOptimizationTechnique();
		
		use_row = determineUsable();
		
		old_parameters = last_iteration.getOldParameters();
		parameters = last_iteration.getCoef();
		
		evaluator = last_iteration.getFunction();
		
		success = last_iteration.isSuccessful();
		
		this.setExtendedStatistics(last_iteration.getExtendedStatisticsWithoutCompute());
	}
	
	public int getOptimizationTechnique() {
		return technique;
	}
	
	public int getCurrentIteration() {
		// TODO Auto-generated method stub
		return iteration;
	}

	public double getIterationChange() {
		return change;
	}

	public LikelihoodFunction getFunction() {
		return evaluator;
	}
	
	public LSIterative getNextIteration() {
		return new SimpleMLEFit(this);
	}

	public String getResults() {
		if(results == null) {
			StringBuffer text = null;
			String results_header = "Results of the "+evaluator.getLikelihoodType()+" MLE Regression Model\n";
			if(!success) {
				text = LSFitResultsGen.GenerateFailureResults(this,getFailureMessage(),results_header);
			} else {
				try {
					text = LSFitResultsGen.GenerateResults(this,results_header);
				} catch(Exception e) {
					
					this.failure_message = "Unintelligible results: "+e.getMessage()+"<br>Convergence probably did not occur.";
					if(!ResultsUtilities.USE_HTML)
						this.failure_message = this.failure_message.replaceAll("<br>", "\n");
					
					text = LSFitResultsGen.GenerateFailureResults(this,getFailureMessage(),results_header);
					success = false;
				}
				//print_special_options(text);
				//checkExtendedOutput(text);
			}
			results = text.toString();
		}
		return results;
	}

	public double[][] getCovarianceMatrix() {
		if(covariance_matrix != null) return covariance_matrix;
		if(!success) return null;

		double[][] ls_coef = generate_Hessian_matrix();
		
		ls_coef = MatrixSupport.Negate(ls_coef);
		
		// Check if this can be invertedhessian matrix logit
		double rcond = MatrixSupport.GetRecipricolCondition(ls_coef);
		
		if(rcond > RunLS.MINIMUM_R_CONDITION) {

			int[] ipiv = new int[parameters.length];
			double[] work = new double[parameters.length];
			intW result = new intW(0);

			// Together, these invert a matrix (temp)
			DGETRF.DGETRF(parameters.length,parameters.length,ls_coef,ipiv,result);
			DGETRI.DGETRI(parameters.length,ls_coef,ipiv,work,work.length,result);

			// Multiply by residual variance
			//MatrixSupport.ScalarTimesMatrix(getResidualVariance(),ls_coef,parameters.length,parameters.length);

			// Do we need white estimates???
			//if(variance_type == VAR_OLS) {
			//	MatrixSupport.ScalarTimesMatrix(residual_variance(),temp,indep_size[1],indep_size[1]);
			//} else if(variance_type == VAR_WHITE) {				
			//	temp = calculateWhiteVCVMatrix(temp);
			//}
			covariance_matrix = ls_coef;
			
		} else {
			covariance_matrix = MatrixSupport.MatrixofZeros(parameters.length, parameters.length);
		}
		
		
		return covariance_matrix;
	}

	public double[] getResiduals() {
		double[][] data = model.getDataDouble(use_row,model.getColumnNumber(independent_var));
		double[] results = model.getDataDouble(use_row,model.getColumnNumber(dependent_var));
		for(int i=0;i<results.length;i++) {
			results[i] = results[i] - evaluator.evaluateFunction(parameters, data[i]);
		}
		return results;
	}

	public double getResidualSumSqared() {
		double[] residuals = getResiduals();
		if(residuals == null) return 0;
		
		return DDOT.DDOT(residuals.length,residuals,1,residuals,1);
	}

	public double getTotalSumSqared() {
		double[] dep = model.getDataDouble(use_row,model.getColumnNumber(dependent_var));
		if(dep == null) return 0;
		
		return DDOT.DDOT(dep.length,dep,1,dep,1);
	}

	public double[] getCoef() {
		return parameters;
	}

	public String[] getCoefLabels() {
		
		if(parameters.length == independent_var.length) {
			return getIndependentVariables();
		}
		
		String[] names = new String[independent_var.length+1];
		names[0] = "Constant";
		for(int i=0;i<independent_var.length;i++) names[i+1] = independent_var[i];
		
		return names;
	}

	public double[] getCoefVariances() {
		double[] variances = MatrixSupport.GetDiagonal(getCovarianceMatrix());
		return variances;
	}

	public String[] getIndependentVariables() {
		return independent_var;
	}

	public String getDependentVariable() {
		return dependent_var;
	}

	public boolean isSuccessful() {
		return success;
	}

	public Double[] getFittedData() {
		Double[] result = new Double[model.getRowCount()];
		
		double[][] data = model.getDataDouble(use_row,model.getColumnNumber(independent_var));
		
		
		for(int i=0;i<model.getRowCount();i++)
			result[i] = null;
		
		for(int i=0;i<use_row.length;i++) {
			result[use_row[i]] = evaluator.evaluateFunction(parameters, data[i]);
		}
		return result;
	}

	public Double[] getResidualData() {
		Double[] results = new Double[model.getRowCount()];
		for(int i=0;i<results.length;i++)
			results[i] = null;
		
		double[] residuals = getResiduals();
		for(int i=0;i<use_row.length;i++) {
			results[use_row[i]] = new Double(residuals[i]);
		}
		return results;
	}

	public StatTableModel getModel() {
		return model;
	}
	
	protected double[] getOldParameters() {
		return old_parameters;
	}

	private int[] determineUsable() {
		String[] all_columns = new String[independent_var.length+1];
		all_columns[0] = dependent_var;
		for(int i=0;i<independent_var.length;i++)
			all_columns[i+1] = independent_var[i];
		return model.getNonNullRows(all_columns);
	}
	
	
	private double[] secant_rule(double[][] independent, 
			                     double[] dependent, 
			                     double[] derivatives,
			                     double[] old_derivatives) {
		
		double[] updates = new double[parameters.length];
		for(int i=0;i<updates.length;i++) {
			double denom = derivatives[i] - old_derivatives[i];
			if(denom != 0.0) 
				updates[i] = -1.0*derivatives[i]*(parameters[i]-old_parameters[i])/denom;
			else
				updates[i] = 0.0;
		}
		
		return updates;
	}
	
	private boolean regress() {
		
		double[][] independent_data = model.getDataDouble(use_row,model.getColumnNumber(independent_var));
		double[] dependent_data = model.getDataDouble(use_row,model.getColumnNumber(dependent_var));
		double[] derivatives = evaluator.evaluateLikelihoodDerivative(parameters,independent_data,dependent_data);
		
		double[] updates = null;
		
		if(old_parameters == null) {
			old_parameters = (double[])parameters.clone();
			DSCAL.DSCAL(old_parameters.length,
					    1.0+INITIAL_SLOPE_PERCENTAGE,
					    old_parameters,1);
		}
		
		switch(technique) {
			case TECHNIQUE_SECANT:
			default:
				double[] old_derivatives = evaluator.evaluateLikelihoodDerivative(old_parameters,independent_data,dependent_data);
				updates = secant_rule(independent_data,dependent_data,derivatives,old_derivatives);
				break;
		}
		
		success = (updates != null);
		for(int i=0;i<updates.length && success;i++) {
			if(Double.isNaN(updates[i]))
				success = false;
		}
		
		
		if(success) {
			old_parameters = (double[])parameters.clone();
			DAXPY.DAXPY(parameters.length,relaxation,updates,1,parameters,1);
			change = VectorSupport.Magnitude(updates);
		}
		
		return success;
	}

	private String getFailureMessage() {
		if(failure_message != null)
			return failure_message;
		else
			return "Error in SimpleMLEFit algorithm";
	}

	/** Returns the variance of the residuals
	 * 
	 * @return double representing the residual variance
	 */
	private double getResidualVariance() {
		double denom = new Integer(use_row.length).doubleValue() -
		    new Integer(parameters.length).doubleValue();
		return getResidualSumSqared()/denom;
	}

	/** Generates an X matrix for determining the covariance matrix
	 * 
	 * @return the X matrix
	 */
	private double[][] generate_X_matrix() {
		double[][] indep_data = model.getDataDouble(use_row,model.getColumnNumber(independent_var));
		
		
		double[][] solver_data = new double[indep_data.length][parameters.length];
		
		for(int i=0;i<indep_data.length;i++) {
			int addon = 0;
			if(indep_data[i].length+1 == parameters.length) {
				solver_data[i][0] = 1.0;
				addon = 1;
			}
			for(int j=addon;j<parameters.length;j++) {
				solver_data[i][j] = indep_data[i][j-addon];
			}			
		}			
		return solver_data;
	}
	
	private double[][] generate_Hessian_matrix() {
		double[][] independent_data = model.getDataDouble(use_row,model.getColumnNumber(independent_var));
		double[] dependent_data = model.getDataDouble(use_row,model.getColumnNumber(dependent_var));
		//double[] derivatives = evaluator.evaluateLikelihoodDerivative(parameters,independent_data,dependent_data);
		
		return evaluator.getHessian(parameters, independent_data, dependent_data);
	}

	public void perform() {
		if(this.success) {
			this.success = regress();
		}
	}
	
	/** This method computes and returns the extended statistics
	 * 
	 * @see com.approximatrix.draco.fit.LSFit#getExtendedStatistics()
	 */
	public List<ExtendedStatistic> getExtendedStatistics() {
		if(estats == null) return null;
		for(ExtendedStatistic stat : estats) {
			stat.setAnalysis(this);
			stat.compute();
		}
		return estats;
	}

	/** Sets the list of extended statistics
	 * @see com.approximatrix.draco.fit.LSFit#setExtendedStatistics(java.util.List)
	 */
	public void setExtendedStatistics(List<ExtendedStatistic> statistics) {
		estats = statistics;
	}
	
	/** Returns the current list of extended statistics without performing computations
	 * @see com.approximatrix.draco.fit.LSIterative#getExtendedStatisticsWithoutCompute()
	 */
	public List<ExtendedStatistic> getExtendedStatisticsWithoutCompute() {
		return estats;
	}
}
