/*
 *  Draco Econometrics
 *  Copyright (C) 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.estimation;

import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Set;

import org.apache.commons.math.estimation.EstimatedParameter;
import org.apache.commons.math.estimation.SimpleEstimationProblem;
import org.apache.commons.math.estimation.WeightedMeasurement;

import com.approximatrix.draco.datamodel.Constant;
import com.graphbuilder.math.Expression;
import com.graphbuilder.math.ExpressionTree;
import com.graphbuilder.math.FuncMap;
import com.graphbuilder.math.VarMap;

/** Estimation problem class for solving arbitrary nonlinear least squares
 * problems using MESP for evaluation and Apache Commons Math for minimization
 * algorithms.  Called exclusively by RunNLLS.
 * 
 * @author jeffrey.armstrong@approximatrix.com
 *
 */
public class ArbitraryNonLinearEstimationProblem extends SimpleEstimationProblem {

	/** Absolute step to use in computations */
	private static final double INITIAL_STEP = 1.0E-5;
	
	private Expression exp = null;
	private FuncMap fm = null;
	
	private String[] independent = null;
	
	private Map<String,Double> constants = null;
	
	public ArbitraryNonLinearEstimationProblem(String equation, 
			                                   double[][] A,
			                                   String[] independent,
			                                   double[] b,
			                                   String dependent,
			                                   List<Constant> constants) {
		exp = ExpressionTree.parse(equation);
		
		fm = new FuncMap();
		fm.loadDefaultFunctions();
		
		if(A.length != b.length)
			throw new IllegalArgumentException("Mismatch of indpendent and dependent data sizes");
		
		this.independent = independent;
		//this.dependent = dependent;
		this.buildMeasurements(A, b);
		
		this.constants = this.buildConstantMap(constants);
	}
	
	public void addParameters(String[] labels, double[] initial) {
		if(labels.length != initial.length)
			throw new ArrayIndexOutOfBoundsException("paramter array length must equal initial guess array length");
		
		for(int i=0;i<labels.length;i++)
			this.addParameter(new EstimatedParameter(labels[i],initial[i]));
	}

	private void buildMeasurements(double[][] A, double[] b) {
		for(int i=0;i<A.length;i++)
			this.addMeasurement(new ArbitraryMeasurement(this.getMeasurementData(i, A), b[i]));
	}
	
	/** Updates the object's parser with all the estimates of parameters
	 * 
	 */
	private VarMap setMESPValues(Map<String,Double> estimates, VarMap variables) {
		if(variables == null)
			variables = new VarMap();
		
		Set<String> varNames = estimates.keySet();
		for(String key : varNames)
			variables.setValue(key, estimates.get(key).doubleValue());
		
		return variables;
	}
	
	private Map<String,Double> buildConstantMap(List<Constant> constantList) {
		HashMap<String,Double> map = new HashMap<String,Double>();
		for(Constant c : constantList)
			map.put(c.getName(),c.getValue().doubleValue());
		return map;
	}
	
	public double[] evaluateAll() {
		WeightedMeasurement[] wm = this.getMeasurements();
		double[] ret = new double[this.getMeasurements().length];
		
		for(int i=0;i<ret.length;i++)
			ret[i] = wm[i].getTheoreticalValue();
		
		return ret;
	}
	
	private double evaluate(Map<String,Double> data, Map<String,Double> estimates) {
		
		VarMap variables = setMESPValues(estimates,null);
		variables = setMESPValues(constants,variables);
		variables = setMESPValues(data,variables);

		double value = 0.0;
		try {
			value = exp.eval(variables,fm);
		} catch(java.lang.RuntimeException re) {
			//JOptionPane.showMessageDialog(jContentPane,re.getMessage(), "Equation Error", JOptionPane.ERROR_MESSAGE);
			value = Double.NaN;
		}
		
		return value;
	}
	
	private Map<String,Double> getCurrentEstimates() {
		HashMap<String,Double> map = new HashMap<String,Double>();
		EstimatedParameter[] params = this.getAllParameters();
		for(int i=0;i<params.length;i++)
			map.put(params[i].getName(), new Double(params[i].getEstimate()));
		return map;
	}
	
	private Map<String,Double> getMeasurementData(int count, double[][] A) {
		HashMap<String,Double> map = new HashMap<String,Double>();
		double[] values = A[count];
		if(values.length != independent.length)
			throw new ArrayIndexOutOfBoundsException("Independent data and variable name arrays do not have matching lengths");
		
		for(int i=0;i<values.length;i++)
			map.put(independent[i], values[i]);
		
		return map;
	}
	
	private class ArbitraryMeasurement extends WeightedMeasurement {

		private Map<String,Double> indep = null;
		
		public ArbitraryMeasurement(Map<String,Double> independent, double value) {
			super(1.0,value);
			this.indep = independent;
		}
		
		@Override
		public double getPartial(EstimatedParameter arg0) {
			Map<String,Double> params = getCurrentEstimates();
			
			double base = arg0.getEstimate();
			
			
			params.put(arg0.getName(), new Double(base+INITIAL_STEP));
			double y_plus = evaluate(indep,params);
			
			params.put(arg0.getName(), new Double(base-INITIAL_STEP));
			double y_minus = evaluate(indep,params);
			
			double partial = 0.0;
			if(Double.isNaN(y_plus) && Double.isNaN(y_minus))
				partial = 0.0;
			
			else if(Double.isNaN(y_plus)) {
				params.put(arg0.getName(), new Double(base));
				double y_base = evaluate(indep,params);
				partial = (y_base-y_minus)/INITIAL_STEP;
			
			} else if(Double.isNaN(y_minus)) {
				params.put(arg0.getName(), new Double(base));
				double y_base = evaluate(indep,params);
				partial = (y_plus-y_base)/INITIAL_STEP;
			}  else
				partial = 0.5*(y_plus-y_minus)/INITIAL_STEP;
			
			return partial;
		}

		@Override
		public double getTheoreticalValue() {
			return evaluate(indep,getCurrentEstimates());
		}
		
	}
	
}
