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


import org.netlib.blas.DAXPY;

import com.approximatrix.draco.fit.RunLS;
import com.approximatrix.draco.math.VectorSupport;

public class LogitLikelihood implements LikelihoodFunction {

	private static final String FUNCTION_NAME = "Logit";
	
	public LogitLikelihood() {
		// Dummy
	}
	
	public double[] initializeParameters(double[] dependent,
			double[][] independent) {
		double[] results = RunLS.SimpleLS(independent, dependent);
		if(results == null) return null;
		return results;
		
		// The following excludes the constant, but I think we need it.
		//double[] answer = new double[results.length-1];
		//for(int i=0;i<answer.length;i++)
		//	answer[i] = results[i+1];
		//return answer;
	}

	public double evaluateFunction(double[] parameters, double[] variables) {
		double exponential = calculateExponential(parameters,variables);
		return exponential/(1.0+exponential);
	}
	
	public double evaluateLikelihood(double[] parameters, 
									 double[][] independent, 
									 double[] dependent) {
		double L = 0.0;
		
		for(int i=0;i<dependent.length;i++) {
			double[] row = new double[independent[i].length];
			for(int j=0;j<independent[i].length;j++)
				row[j] = independent[i][j];
			
			L += (dependent[i]*calculateLinearTerm(parameters,row) - 
			     Math.log(1.0+calculateExponential(parameters,row)));
		}
		
		return L;
	}

	public double[] evaluateLikelihoodDerivative(double[] parameters,
												 double[][] independent, 
												 double[] dependent) {
		
		double[] derivatives = new double[parameters.length];
		for(int i=0;i<derivatives.length;i++)
			derivatives[i] = 0;
		
		for(int i=0;i<dependent.length;i++) {
			double[] row = new double[independent[i].length];
			double[] complete_vars = new double[row.length+1];
			complete_vars[0] = 1.0;
			for(int j=0;j<independent[i].length;j++) {
				row[j] = independent[i][j];
				complete_vars[j+1] = row[j];
			}
			
			double difference = dependent[i] - evaluateFunction(parameters,row);
			
			DAXPY.DAXPY(complete_vars.length,difference,complete_vars,1,derivatives,1);
			
		}
		
		return derivatives;
	}
	
	private double calculateExponential(double[] parameters, double[] variables) {
		return Math.exp(calculateLinearTerm(parameters,variables));
	}

	private double calculateLinearTerm(double[] parameters, double[] variables) {
		// For the constant to be included
		double[] complete_vars = new double[variables.length+1];
		complete_vars[0] = 1.0;
		for(int i=0;i<variables.length;i++)
			complete_vars[i+1] = variables[i];
		
		return VectorSupport.VectorElementSum(VectorSupport.MultiplyElements(parameters,complete_vars));
	}
	
	public String getLikelihoodType() {
		return FUNCTION_NAME;
	}
	
	public double[][] getHessian(double[] parameters, double[][] independent, double[] dependent) {
		double[][] return_val = new double[parameters.length][parameters.length];
		
		for(int i = 0;i<parameters.length;i++) {
			
			for(int j = i;j<parameters.length;j++) {
				
				return_val[i][j] = 0.0;
				for(int indep = 0;indep < independent.length;indep++) {
					double p1=0.0; double p2 = 0.0;
					
					double[] indep_row = new double[independent[indep].length];
					for(int row_build = 0;row_build < independent[indep].length; row_build++)
						indep_row[row_build] = independent[indep][row_build];
					
					// Determine the x1, x2 variables
					if(i==0)
						p1 = 1.0;
					else
						p1 = indep_row[i-1];

					if(j==0)
						p2 = 1.0;
					else
						p2 = indep_row[j-1];
					
					double exponential = calculateExponential(parameters,indep_row);
					
					return_val[i][j] -= p1*p2*exponential/((1.0+exponential)*(1.0+exponential));
				}
			}
		}
		
		// Only the top was filled in, copy the top to the bottom
		for(int i=0;i<parameters.length;i++) {
			for(int j=i+1;j<parameters.length;j++) {
				return_val[j][i] = return_val[i][j];
			}
		}
		
		return return_val;
	}
	
}
