/*
 *  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.analysis;

import org.netlib.blas.DASUM;

import com.approximatrix.draco.analysis.statistics.FDistribution;
import com.approximatrix.draco.datamodel.StatTableModel;
import com.approximatrix.draco.math.VectorSupport;
import com.approximatrix.draco.utility.ResultsUtilities;

/** Class supporting Analaysis of Variance runs
 * 
 * @author Jeffrey Armstrong
 * @see http://davidmlane.com/hyperstat/intro_ANOVA.html
 *
 */ 
public class RunANOVA implements StandardAnalysis {

	//private static final double[] CERTAINTY_TESTS = { 90.0, 95.0, 97.5, 99.0 };
	
	private StatTableModel model = null;
	private String[] variables = null;
	
	private boolean success = true;
	private StringBuffer failure_message = null;
	
	private int df_numerator = 1;
	private int df_denominator = 1;
	private double F_statistic = 0.0;
	
	private String results = null;
	
	private Double grand_mean = null;
	private Double SSbg = null;
	private Double SSt = null;
	
	private Integer totalSamples = null;
	
	public RunANOVA() {
		this.success = false;
		failure_message = new StringBuffer("No analysis performed...");
	}
	
	public RunANOVA(StatTableModel model, String[] variables) {
		this.model = model;
		this.setVariables(variables);
		
		// Set up the failure_message buffer
		failure_message = new StringBuffer();
		
		// Run analysis
		success = doANOVA();
	}
	
	public String[] getVariables() {
		return variables;
	}
	
	public boolean isSuccessful() {
		return success;
	}

	public void setVariables(String[] variables) {
		this.variables = variables;
	}
	
	public void setSuccess(boolean value) {
		this.success = value;
	}
	
	public void setResults(String results) {
		this.results = results;
	}
	
	public String getFailureMessage() {
		return failure_message.toString();
	}
	
	public int getNumberSteps() {
		return 1;
	}
	
	public int getCurrentStep() {
		return 2;
	}
	
	public String getResults() {
		if(this.results == null) {
			StringBuffer results = new StringBuffer();
			
			if(ResultsUtilities.USE_HTML)
				results.append(ResultsUtilities.GenHTMLHeader());
			
			results.append(ResultsUtilities.GetTitle("Results of the One-Factor Analysis of Variance (ANOVA)"));

			results.append(variables.length);
			results.append(" Variables Considered:");
			
			String[] list = new String[variables.length];
			
			for(int i=0;i<variables.length;i++) {
				StringBuffer builder = new StringBuffer();
				builder.append(variables[i]);
				builder.append("\t");
				builder.append(model.getObeservationCount(variables[i]));
				builder.append(" Observations");
				list[i] = builder.toString();
			}
			
			results.append(ResultsUtilities.GenerateList(list));
			
			results.append("\n\n");
			results.append(ResultsUtilities.SeparatorLine());
			if(this.success) {
				results.append("Resulting F-Statistic: ");
				results.append(ResultsUtilities.OutNumber(F_statistic));
				
				if(ResultsUtilities.USE_HTML)
					results.append("<BR>");
				
				results.append("\n\nDegrees of Freedom: \n ");
				
				String[][] table = new String[2][2];
				table[0][0] = "Numerator: "; 
				table[0][1] = Integer.toString(this.df_numerator);
				table[1][0] = "Denominator:";
				table[1][1] = Integer.toString(this.df_denominator);
				
				results.append(ResultsUtilities.GenerateTable(table, false, true));
				
				results.append("\nConfidence Interval: ");
				results.append(ResultsUtilities.OutNumber(100.0*FDistribution.getCDF(F_statistic, this.df_numerator, this.df_denominator)));
				results.append("%\n");
				results.append(ResultsUtilities.SeparatorLine());
				
				// Sum-of-Squares partitioning
				results.append("Sum-of-Squares Partioning\n");				
				String[][] part = new String[4][4];
				
				part[0][0] = "Source";
				part[0][1] = "df";
				part[0][2] = "Ssq";
				part[0][3] = "Msq";
				
				double bgssq = getBGSumSquares();
				double tssq = getTotalSumSquares();
				double errssq = tssq-bgssq;
				
				part[1][0] = "Groups";
				part[2][0] = "Error";
				part[3][0] = "Total";
				
				part[1][1] = Integer.toString(this.df_numerator);
				part[2][1] = Integer.toString(this.df_denominator);
				part[3][1] = Integer.toString(this.df_numerator+this.df_denominator);

				part[1][2] = ResultsUtilities.OutNumber(bgssq);
				part[1][3] = ResultsUtilities.OutNumber(calculateMSB());

				part[2][2] = ResultsUtilities.OutNumber(errssq);
				part[2][3] = ResultsUtilities.OutNumber(calculateMSE());
				
				part[3][2] = ResultsUtilities.OutNumber(tssq);
				part[3][3] = "  ";
				//part[3][3] = ResultsUtilities.OutNumber(tssq/((double)this.df_numerator+(double)this.df_denominator));
				
				results.append(ResultsUtilities.GenerateTable(part, true, false));

				// Certainty tests
				//results.append("Hypothesis Test (Means are equal):\n\n");
				//results.append("\tCertainty\tHypothesis Status\n");
				//for(int i=0;i<CERTAINTY_TESTS.length;i++){
				//	results.append("\t");
				//	results.append(CERTAINTY_TESTS[i]);
				//	results.append("%    \t");
				//	if(F_statistic<=FDistribution.FValue(CERTAINTY_TESTS[i]/100.0,df_numerator,df_denominator)) 
				//		results.append("Cannot Reject");
				//	else
				//		results.append("Rejected");
				//	
				//	results.append(" (");
				//	results.append(FDistribution.FValue(CERTAINTY_TESTS[i]/100.0,df_numerator,df_denominator));
				//	results.append(")");
				//	results.append("\n");

				//}
				//results.append("\n");
				//results.append(LSFitResultsGen.SeparatorLine());

			} else {
				results.append("Analysis Failed:\n");
				results.append(this.getFailureMessage());
			}
			
			if(ResultsUtilities.USE_HTML)
				results.append(ResultsUtilities.GenHTMLFooter());
			
			this.results = results.toString();
		}
		return this.results;
	}
	
	private boolean doANOVA() {
		
		// Default to success
		success = true;
		
		double mse = calculateMSE();
		double msb = calculateMSB();
		
		F_statistic = getFStatistic(msb,mse);
		
		// Degrees of freedom are necessary for calculating the F value
		df_numerator = variables.length-1;
		df_denominator = calculateTotalSamples() - variables.length;
		
		if(new Double(F_statistic).isNaN()) { 
			failure_message.append("Error calculating F-statistic: Divide by zero occurred\n");
			success = false;
		}
		
		if(df_numerator == 0) {
			failure_message.append("Error in analysis: No degrees of freedom in MSB (numerator)\n");
			success = false;
		}
		
		if(df_denominator == 0) {
			failure_message.append("Error in analysis: No degrees of freedom in MSE (denominator)\n");
			success = false;
		}
		
		return success;
	}
	
	private double calculateMSE() {
//		double sum = 0.0;
//		for(int i=0;i<variables.length;i++) {
//			double sigma = model.getColumnStandardDeviation(variables[i]);
//			sum += sigma*sigma;
//		}
//		return sum/((double)variables.length);
		return getErrorSumSquares()/(double)(calculateTotalSamples() - variables.length);
	}
	
	private double calculateMSB() {
//		double means[] = new double[variables.length];
//		for(int i=0;i<variables.length;i++) {
//			means[i] = model.getColumnMean(variables[i]);
//		}
//		double std = VectorSupport.StandardDeviation(means);
//		
//		return std*std*((double)variables.length);
		return getBGSumSquares()/((double)variables.length-1.0);
	}
	
	private double getFStatistic(double msb, double mse) {
		return msb/mse;
	}
	
	private int calculateTotalSamples() {
		if(totalSamples == null) {
			int sum = 0;
			for(int i=0;i<variables.length;i++)
				sum+=model.getObeservationCount(variables[i]);
			totalSamples = new Integer(sum);
		} 
		return totalSamples.intValue();
	}
	
	/** 
	 * @return Sum squares between groups
	 */
	private double getBGSumSquares() {
		
		if(SSbg == null) {

			double[] means = new double[variables.length];
			double[] samples = new double[variables.length];
			for(int i=0;i<variables.length;i++) {
				means[i] = model.getColumnMean(variables[i]);
				samples[i] = (double)model.getObeservationCount(variables[i]);
			}

			// Square the elements
			means = VectorSupport.MultiplyElements(means, means);

			// Divide by each sample size
			//means = VectorSupport.DivideElements(means, samples);
			means = VectorSupport.MultiplyElements(means, samples);

			//DAXPY.DAXPY(means.length, -1.0*getGrandMean()*getGrandMean()/VectorSupport.VectorElementSum(samples), 
			//				    VectorSupport.VectorOfOnes(means.length), 1, 
			//				    means, 1);
			
			// Sum the squares
			double ret = DASUM.DASUM(means.length,means,1);

			// Subtract the grand mean squared over the total number of samples
			ret = ret - getGrandMean()*getGrandMean()*VectorSupport.VectorElementSum(samples);

			SSbg = new Double(ret);

			// Compute differences from the mean
			//DAXPY.DAXPY(means.length, -1.0*getGrandMean()/VectorSupport.VectorElementSum(samples), 
			//		    VectorSupport.VectorOfOnes(means.length), 1, 
			//		    means, 1);

			// Calculate the sum of the resulting vector (DASUM is safe
			// since all elements are positive)
			//return DASUM.DASUM(means.length,means,1);

		}
		return SSbg.doubleValue();
	}
	
	private double getErrorSumSquares() {
		return getTotalSumSquares() - getBGSumSquares();
	}
	
	private double getTotalSumSquares() {
		if(SSt == null) {
			double ret = 0.0;
			int total_vars = 0;
			for(int i=0;i<variables.length;i++) {
				// Get one variable
				double[] one_var = model.getDataDouble(model.getNonNullRows(variables[i]), model.getColumnNumber(variables[i]));
				
				// Store the observation count
				total_vars += one_var.length;
				
				// Square all values
				one_var = VectorSupport.MultiplyElements(one_var, one_var);
				
				// Sum the elements and add them to the total sum
				ret = ret + DASUM.DASUM(one_var.length,one_var,1);
			}

			// Subtract the grand mean squared divided by the total number of observations
			ret = ret - getGrandMean()*getGrandMean()*((double)total_vars);
			
			// Store the result
			SSt = new Double(ret);
		}
		return SSt.doubleValue();
	}
	
	private double getGrandMean() {
		if(grand_mean == null) {
			double return_val = 0.0;
			for(int i=0;i<variables.length;i++) 
				return_val += model.getColumnMean(variables[i])*((double)model.getObeservationCount(variables[i]));
		
			grand_mean = new Double(return_val/((double)calculateTotalSamples()));
		}
		return grand_mean.doubleValue();
	}
}
