/*
 * #%L
 * IsisFish data
 * %%
 * Copyright (C) 2006 - 2011 Ifremer, CodeLutin
 * %%
 * 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 2 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/gpl-2.0.html>.
 * #L%
 */
package scripts;

import java.util.List;

import org.apache.commons.logging.Log;
import org.apache.commons.logging.LogFactory;

import org.nuiton.math.matrix.*;

import fr.ifremer.isisfish.entities.Population;
import fr.ifremer.isisfish.entities.PopulationGroup;
import fr.ifremer.isisfish.entities.Species;
import fr.ifremer.isisfish.simulator.SimulationContext;
import fr.ifremer.isisfish.types.TimeStep;


/**
 * RuleUtil.
 *
 * Created: 6 septembre 2006
 *
 * @author anonymous <anonymous@labs.libre-entreprise.org>
 * @version $Revision: 1.3 $
 *
 * Last update: $Date: 2007-11-20 15:51:05 $
 * by : $Author: bpoussin $
 */
public class RuleUtil {

    /** to use log facility, just put in your code: log.info("..."); */
    static private Log log = LogFactory.getLog(RuleUtil.class);

    /**
     * @param context le context de simulation
     * @param species l'espece sur lequel on souhaite le total
     * @param date la date pour laquel on veut le calcul, cet argument est
     * passe seulement pour que le cache ne retourne pas toujours la meme valeur
     * @return total catch in tons
     */
    public static double getTotalCatchTons(SimulationContext context, Species species, TimeStep step) {
        double result = 0;
        for (Population pop : species.getPopulation()) {
            MatrixND mat = context.getPopulationMonitor().getHoldCatch(pop);
            if (mat != null) {
                mat = mat.copy();
                mat = mat.sumOverDim(0); // sum over Strategies
                mat = mat.sumOverDim(1); // sum over metiers
                mat = mat.sumOverDim(3); // sum over zones

                List<PopulationGroup> groups = pop.getPopulationGroup();

                for (int c = 0; c < groups.size(); c++) {
                    PopulationGroup group = groups.get(c);
                    double weight = group.getMeanWeight();
                    result += mat.getValue(0, 0, c, 0) * weight / 1000.0;
                }
            }
        }
        return result;
    }
	
	public static double getTotalCatchTonsPop(SimulationContext context, Population pop, TimeStep step) {
        double result = 0;
            MatrixND mat = context.getPopulationMonitor().getHoldCatch(pop);
            if (mat != null) {
                mat = mat.copy();
                mat = mat.sumOverDim(0); // sum over Strategies
                mat = mat.sumOverDim(1); // sum over metiers
                mat = mat.sumOverDim(3); // sum over zones

                List<PopulationGroup> groups = pop.getPopulationGroup();

                for (int c = 0; c < groups.size(); c++) {
                    PopulationGroup group = groups.get(c);
                    double weight = group.getMeanWeight();
                    result += mat.getValue(0, 0, c, 0) * weight / 1000.0;
                }
            }
        return result;
    }
	
	/*
	// Le probleme est qu'en l'etat deriveF n'est pas general puisque ff est une equation liee a notre besoin
	public static double deriveF(TimeStep step, double xx, double C, double M, double N) { // Derives Fsq;
	//step:
	//xx:
	//C:
	//ffFBaranov : calcul de la difference entre les captures ISIS (annuelles) et les captures deduites de l'equation de Baranov pour une valuer de F=xx, une abondance en debut d'annee de N et une mortalite naturelle annuelle de M

		double ff;
		ff = Math.pow((C - (xx/(xx+M))*(1-Math.exp(-(xx+M)))*N),2);
		//log.info("ff= " + ff);
		return ff;  
		}
*/	
	
	// Idem : je doute qu'en l'etat fmin soit general	
	//public static double fmin (TimeStep step, double a, double b, double tol, double C, double M, double N) {
	public static double fmin (double a, double b, double tol, FonctionObjectif f) {
	// Algo de minimisation de derivF(step, x,C,M,N) sur x 
	// en prenant une valeur initiale de x dans [a,b]
	// et une tolerance = tol
	// Reference de l'algo : http://www1.fpl.fs.fed.us/optimization.html 
	// Algorithme de Brent qui approxime la fonction a minimiser par
	// une parabole si possible et fait du "golden step" sinon
	
		double c,d,e,eps,xm,p,q,r,tol1,t2,u,v,w,fu,fv,fw,fx,x,tol3;

		c = .5*(3.0 - Math.sqrt(5.0));
		d = 0.0;
		eps = 1.2e-16; // 1.1102e-16 is machine precision
		tol1 = eps + 1.0;
		eps = Math.sqrt(eps);
		v = a + c*(b-a);
		//v = 0.2;
		log.info("Finit=" + v);
		w = v;
		x = v;
		e = 0.0;
		fx = 0.0;
		fu = 0.0;
		
		//fx = deriveF(step, x, C, M, N);
		fx= f.optimize(new Object(), x);
		log.info("fx= " + fx);

		
		fv = fx;
		fw = fx;
		tol3 = tol/3.0;
		xm = .5*(a + b);
		tol1 = eps*Math.abs(x) + tol3;
		t2 = 2.0*tol1;

		while (Math.abs(x-xm) > (t2 - .5*(b-a))) { // main loop;
		/////// De l�...
			p = q = r = 0.0;
			if (Math.abs(e) > tol1) { // fit the parabola;
				r = (x-w)*(fx-fv);
				q = (x-v)*(fx-fw);
				p = (x-v)*q - (x-w)*r;
				q = 2.0*(q-r);
				if (q > 0.0) {
					p = -p;
				}
				else {
					q = -q;
				}
				r = e;
				e = d;         
			}
			if ((Math.abs(p) < Math.abs(.5*q*r)) && (p > q*(a-x)) && (p < q*(b-x))) { // parabolic interpolation step;
				d = p/q;
				u = x+d;
				if (((u-a) < t2) || ((b-u) < t2)) { // f must not be evaluated too close to a or b;
					d = tol1;
					if (x >= xm) d = -d;
				}
			}
			else { // a golden-section step;
				if (x < xm) {
					e = b-x;
				}
				else {
					e = a-x;
				}
				d = c*e;
			}
			if (Math.abs(d) >= tol1) { // f must not be evaluated too close to x;
				u = x+d;
			}
			else {
				if (d > 0.0) {
					u = x + tol1;
				}
				else {
					u = x - tol1;
				}
			}
			////// ... � ici : ne sert qu'� calculer un u (selon 2 mani�res diff�rentes) afin de comparer fx a fu 
			//////--> On peut tr�s bien prendre u � partir d'une autre simu a la place (?)
			
			//fu = deriveF(step, u, C, M, N);
			fu= f.optimize(new Object(), u);
			log.info("fu= " + fu);


			// Update a, b, v, w, and x
			if (fx <= fu) {
				if (u < x) {
					a = u;
				}
				else {
					b = u;
				}
			}
			if (fu <= fx) {
				if (u < x) {
					b = x;
				}
				else {
					a = x;
				}
				v = w;
				fv = fw;
				w = x;
				fw = fx;
				x = u;
				fx = fu;
				xm = .5*(a + b);
				tol1 = eps*Math.abs(x) + tol3;
				t2 = 2.0*tol1;
			}
			else {
				if ((fu <= fw) || (w == x)) {
					v = w;
					fv = fw;
					w = u;
					fw = fu;
					xm = .5*(a + b);
					tol1 = eps*Math.abs(x) + tol3;
					t2 = 2.0*tol1;
				}
				else if ((fu > fv) && (v != x) && (v != w)) {
					xm = .5*(a + b);
					tol1 = eps*Math.abs(x) + tol3;t2 = 2.0*tol1;
				}
				else {
					v = u;
					fv = fu;
					xm = .5*(a + b);
					tol1 = eps*Math.abs(x) + tol3;
					t2 = 2.0*tol1;
				}
			}
		}
		return x;
	}	
}