/*
 *  Draco Econometrics
 *  Copyright (C) 2008-2009 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 abstract class SimpleLikelihood implements LikelihoodFunction {

	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;
	}

	public double[] initializeParameters(double[] dependent,
			double[][] independent) {
		double[] results = RunLS.SimpleLS(independent, dependent);
		if(results == null) return null;
		return results;
	}

	protected 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));
	}
	
	protected abstract double evaluatePDF(double[] parameters, double[] variables);
	
	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];

					return_val[i][j] -= p1*p2*evaluatePDF(parameters,indep_row);
				}
			}
		}
		
		// 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;
	}
}
