/*
 * #%L
 * IsisFish data
 * %%
 * Copyright (C) 2012, 2014 Ifremer, CodeLutin, Chatellier Eric
 * %%
 * 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/gpl-3.0.html>.
 * #L%
 */

package scripts;

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

import fr.ifremer.isisfish.annotations.Nocache;

public class MinimisationUtil {

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

    /**
     * Algo de minimisation de fmin(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
     * 
     * @param a
     * @param b
     * @param tol tolérance
     * @param f fonction d'objectif
     * @return
     */
    @Nocache
    public static double fmin(double a, double b, double tol, ObjectiveFunction f) {

        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);
        // log.info("Finit=" + v);
        w = v;
        x = v;
        e = 0.0;
        fx = 0.0;
        fu = 0.0;

        fx = f.compute(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 = f.compute(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;
    }
}