package optimizacion;

import operaciones.Matriz;

/**
 * <p>Title: Algoritmos para encontrar minimos de una funcion</p>
 * <p>Description: Se dan un conjunto de algoritmos para
 * minimizacion de funciones basados en gradiente</p>
 * <p>Copyright: Copyright (c) 2005</p>
 * <p>Company: UMSNH</p>
 * @author Dr Felix Calderon Solorio
 * @version 1.0
 */

public class minimos {
    public double c1 = 1e-10, c2 = 1, da = 0.1;
    double T = 1e-3;

    /**
     * M�todo de Maximo descenso de Gradiente
     * @param x0 Matriz
     * @return Matriz
     */

    public Matriz Maximo_Descenso_Gradiente(Matriz x0)
    {
        Matriz g, h, x = null, p;
        double alpha;
        int iter = 1;

        x = Matriz.igual_a(x0);

        do
        {
            // calcula el gradiente de la funcion
            g = Gradiente(x);
            // calcula el Hessiano
            h = Hessiano(x);

            // Direcci�n opuesta al gradiente
            p = Direccion(g);
            // alfa de maximo descenso
            alpha = alpha_MD(g, h);
            System.out.println("Alfa = " + alpha);
            // Actualizaci�n x = x + alpha*p
            x = x.mas(p.por(alpha));

            System.out.println("Iteracion " + iter++);
            x.imprime();
        }while (Magnitud(g) > T);
        return x;
    }

    /**
     * Minimiza una funci�n utilizando Razon dorada
     * @param x0 Matriz
     * @return Matriz
     */

    public Matriz Minimiza_Wolfe_RD(Matriz x0)
    {
        Matriz g, h, x = null, p;
        double alpha, alpha_max;
        int iter = 1;

        x = Matriz.igual_a(x0);

        do
        {
            // calcula el gradiente de la funcion
            g = Gradiente(x);

            // Direcci�n opuesta al gradiente
            p = Direccion(g);

            // Calcula la maxima alpha que cumple condiciones de wolfe
            alpha_max = alpha_wolfe(x, Direccion(p));

            // alfa de maximo descenso
            alpha = alpha_RD(x, p, 0, alpha_max);

            // Actualizaci�n x = x + alpha*p
            x = x.mas(p.por(alpha));

            System.out.println("Iteracion " + iter++);
            x.imprime();
        }while (Magnitud(g) > T);
        return x;
    }

    public Matriz Forsite(Matriz X)
    {
        Matriz X1 = new Matriz(), X2 = new Matriz();
        Matriz Y = new Matriz();
        Matriz g = new Matriz(), p = new Matriz();
        int i, N =2, k =1;
        double alpha;

        X1 = Matriz.igual_a(X);
        X2 = Matriz.igual_a(X);

        do
        {
            for (i = 0; i < N; i++) {
                g = Gradiente(X2);
                p = Direccion(g);
                alpha = alpha_Direcccion(X2, p);
                X2 = X2.mas(p.por(alpha));
            }

            Y = X2.menos(X1);
            alpha = alpha_Direcccion(X2, Y);

            X2 = X2.mas(p.por(alpha));
            X1 = Matriz.igual_a(X2);

            System.out.println("Iteracion " + k++);
            X2.imprime();
        } while(Magnitud(Gradiente(X2)) > T);

        return X2;
    }

    /**
     * Calcula el tama�o de paso para una funci�n en la direcci�n d
     * @param x Matriz
     * @param d Matriz
     * @return double
     */

    public double alpha_Direcccion(Matriz x, Matriz d)
    {
        double f1, f2, da = 0.1, alpha = 0;
        Matriz p = Unitario(d);
        f2 = funcion(x);
        do
        {
          alpha += da;
          f1 = f2;
          f2 = funcion(x.mas(p.por(alpha)));
        } while (f2 < f1);

        return (alpha_RD(x, p, 0, alpha));
    }

    /**
     * Calcula el m�ximo tama�o de paso que cumple con las condiciones de
     * Wolfe
     * @param x Matriz punto inicial
     * @param p Matriz direcci�n
     * @return double
     */

    public double alpha_wolfe(Matriz x, Matriz p)
    {
        double alpha =0, gp1, gp;
        Matriz x1, g;
        boolean armijo, curvatura;

        g = Gradiente(x);

        gp = (p.T().por(g)).obten(0,0);

        do
        {
            alpha += da;
            x1 = x.mas(p.por(alpha));

            g = Gradiente(x1).T();
            gp1 = g.por(p).obten(0,0);

            curvatura =  gp1 > c2 *gp;
            armijo = phi(x,alpha, p) <= phi(x, 0, p) + c1 *gp;

        }while((armijo && curvatura) == true) ;

        return alpha;
    }

    /**
     * Tama�o de paso de maximo descenso
     * @param g Matriz Vector gradiente
     * @param H Matriz Matriz Hessiana
     * @return double
     */

    public double alpha_MD(Matriz g, Matriz H)
    {
        Matriz p = Direccion(g);
        Matriz alfa = (p.T().por(g)).entre(p.T().por(H.por(p)));
        return (-alfa.obten(0,0));
    }

    /**
     * Calcula el tamano de paso en un intervalo [alow, aupper] utilizando
     * razon dorada
     * @param x Matriz punto de evaluacion
     * @param p Matriz direccion de busqueda
     * @param alow double limite inferior de alfa
     * @param aupper double limite superior de alfa
     * @return double tamano de paso calculado
     */

    public double alpha_RD(Matriz x, Matriz p, double alow, double aupper)
    {
        double d, a1, a2, R = (Math.sqrt(5.0) -1.0)/2.0;

        do
        {
            d = (aupper - alow)*R;

            a1 = alow   + d;
            a2 = aupper - d;

            //System.out.println(alow + " " + a1 + " " + a2 + " " + aupper);

            if(phi(x, a1, p) < phi(x, a2, p))
                alow = a2;
             else
                aupper = a1;

        } while (aupper - alow > 1e-6);

        return(phi(x,a1, p) < phi(x,a2,p) ? a1 : a2);
    }

    /**
     * Calcula el tama�o de paso utilizando el metodo de interpolaci�n
     * @param x Matriz punto de evaluaci�n
     * @param p Matriz direcci�n
     * @param alpha_0 double alfa inicial
     * @return double
     */

    public double alpha_Int(Matriz x, Matriz p, double alpha_0)
    {
        double a;
        a = -phi_p(x) * alpha_0 * alpha_0/ 2.0 /
            (phi(x, alpha_0, p) - phi_p(x) * alpha_0 - phi(x, 0, p));
        return a;
    }

    /**
     * Calcula un vector unitario en la direccion x
     * @param x Matriz
     * @return Matriz
     */

    public Matriz Direccion(Matriz x)
    {
        double xx = -1.0/Magnitud(x);

        Matriz p = x.por(xx);
        return p;
    }

    /**
     * Calcula un vector unitario del vector x
     * @param x Matriz
     * @return Matriz
     */

    public Matriz Unitario(Matriz x)
    {
        double xx = 1.0/Magnitud(x);

        Matriz p = x.por(xx);
        return p;
    }

    /**
     * Encuentra el m�nimo de un funci�n por el m�todo de Newton
     * @param x Matriz
     * @return Matriz
     */

    public Matriz Newton(Matriz x) {
      int k = 0;
      double tol;

      do {
          System.out.print("iteracion \t" + k  + " \t f(x) = \t " + funcion(x) + "\t");
          x.T().imprime();

          x = x.menos(Gradiente(x).entre(Hessiano(x)));
          k++;
          tol = Magnitud(Gradiente(x));
          //System.out.println("Tol = " + tol);
      } while ( tol > T);
      System.out.println("La solucion en "+ k + " iteraciones es x* = " );
      x.imprime();
      return x;
    }

    /**
     * Minimiza una funci�n utilizando el m�todo de Broyden
     *
     * @param x Matriz
     * @return Matriz
     */

    public Matriz Broyden(Matriz x) {
      int k = 0;
      double tol, ss;
      Matriz A = Hessiano(x);
      Matriz s = null, y = null, g1 = null, g0 = Gradiente(x);
      Matriz aux = null;

      do {
          System.out.print("iteracion \t" + k  + " \t f(x) = \t " + funcion(x) + "\t");
          x.T().imprime();
          s = (g0.entre(A)).por(-1);
          x = x.mas(s);
          g1 = Gradiente(x);
          y = g1.menos(g0);

          ss = 1.0/(s.T().por(s)).obten(0,0);

          aux = (y.menos(A.por(s))).por(s.T());

          A = A.mas(aux.por(ss));
          k++;
          tol = Magnitud(g1);
          g0 = Matriz.igual_a(g1);
          System.out.println("Tol = " + tol);
      } while ( tol > T);
      System.out.println("La solucion en "+ k + " iteraciones es x* = " );
      x.imprime();
      return x;
    }

    /**
     * Metodo BFGS para minimiza una funci�n
     * @param x Matriz Valor inicial
     * @return Matriz Soluci�n
     */

    public Matriz BFGS(Matriz x) {
      int k = 0;
      double tol, ss;
      Matriz A = Hessiano(x);
      Matriz s = null, y = null, g1 = null, g0 = Gradiente(x);
      Matriz aux = null;

      do {
          System.out.print("iteracion \t" + k  + " \t f(x) = \t " + funcion(x) + "\t");
          x.T().imprime();
          s = (g0.entre(A));
          x = x.menos(s);
          g1 = Gradiente(x);
          y = g1.menos(g0);

          ss = 1.0/(s.T().por(s)).obten(0,0);

          aux = s.T().por(A).por(s);

          A = A.menos(((A.por(s)).por(s.T())).por(A).por(1.0/aux.obten(0,0)));

          aux = y.T().por(s);
          A = A.menos((y.por(y.T())).por(1.0/aux.obten(0,0)));

          k++;
          tol = Magnitud(g1);
          g0 = Matriz.igual_a(g1);
          System.out.println("Tol = " + tol);
      } while ( tol > T);
      System.out.println("La solucion en "+ k + " iteraciones es x* = " );
      x.imprime();
      return x;
    }

    /**
     * Minimiza una funci�n utilizando el m�todo de GC.
     *
     * @param x0 Matriz Valor inicial
     * @return Matriz Solucion encontrada
     */

    public Matriz GC_lineal(Matriz x0)
    {
        Matriz x = null, dx = null;
        x = Matriz.igual_a(x0);
        double tol;
        int k =0;
        do
        {
            System.out.print("iteracion \t" + k  + " \t f(x) = \t " + funcion(x) + "\t");
            x.T().imprime();
            dx = GC(x);
            dx.T().imprime();
            x = x.menos(dx);
            tol = Magnitud(Gradiente(x));
            k++;
        }   while (tol > T);
        x.T().imprime();
        return x;
    }

    /**
     * M�todo de Gradiente conjugado de Fletcher-Reeves
     * @param x Matriz valor inicial
     * @return Matriz soluci�n
     */

    public Matriz GC_FR(Matriz x)
    {
        double alpha, alpha_max, beta, rr0, tol;
        int k =0;
        Matriz A = Hessiano(x), b = Gradiente(x), gg =null;

        Matriz g = new Matriz(A.nren, 1);
        Matriz p = new Matriz(A.nren, 1);

        g = Gradiente(x);
        gg = (g.T()).por(g);
        rr0 = gg.obten(0,0);

        p = g.por( -1);

        do {
            alpha_max = alpha_wolfe(x, p);

            alpha = alpha_RD(x, p, 0, alpha_max);

            x = x.mas(p.por(alpha));

            g = Gradiente(x);
            gg = (g.T()).por(g);

            beta = gg.obten(0, 0) / rr0;

            p = (p.por(beta)).menos(g);

            rr0 = gg.obten(0, 0);

            System.out.print("Iteracion \t" + (k) + "\t f(x) = \t" + funcion(x) + "\t");
            x.T().imprime();
            tol = Magnitud(g);
            k++;
        } while(tol > T);
        return x;
    }

    /**
     * M�todo de Gradiente Conjugado Polak-Ribiere
     * @param x Matriz
     * @return Matriz
     */

    public Matriz GC_PR(Matriz x)
    {
        double alpha, alpha_max, beta, rr0, tol;
        int k =0;

        Matriz g0 = null, g1 = null;
        Matriz p = null;

        g0 = Gradiente(x);

        p = g0.por( -1);

        do {
            alpha_max = alpha_wolfe(x, p);
            alpha = alpha_RD(x, p, 0, alpha_max);

            x = x.mas(p.por(alpha));

            g1 = Gradiente(x);

            beta = ((g1.T()).por(g1.menos(g0))).obten(0,0) / (g0.T().por(g0)).obten(0,0);
            //if(beta < 0) beta = 0;

            p = (p.por(beta)).menos(g1);

            System.out.print("Iteracion \t" + (k) + "\t f(x) = \t" + funcion(x) + "\t");
            x.T().imprime();
            tol = Magnitud(g1);
            g0 = Matriz.igual_a(g1);
            k++;
        } while(tol > T);
        return x;
    }

    /**
     * Metodo de gradiente conjugado para solucionar el sistema de ecuacione
     * Ax = b
     *
     * @param b Vector de terminos independientes
     * @return x Solucion del sistema de ecuaciones
     */

    public Matriz GC(Matriz x)
    {
        double alpha, beta, rr0, rr1;
        int Nmax = x.nren + 1;
        Matriz A = Hessiano(x), b = Gradiente(x);

        Matriz rr = new Matriz(1, 1);
        Matriz pAp = new Matriz(1, 1);

        Matriz Ap = new Matriz(A.nren, 1);
        Matriz r = new Matriz(A.nren, 1);
        Matriz p = new Matriz(A.nren, 1);
        Matriz aux1 = new Matriz(1, A.nren);
        Matriz aux2 = new Matriz(A.nren, 1);

        Ap = A.por(x);
        r = Ap.menos(b);
        rr = (r.T()).por(r);
        rr0 = rr.obten(0, 0);
        p = r.por( -1);

        for (int k = 0; k < Nmax; k++) {

            System.out.println("r = ");
            r.imprime();

            System.out.println("p = ");
            p.imprime();

            System.out.println("rr = ");
            rr.imprime();


            Ap = A.por(p);
            pAp = (p.T()).por(Ap);

            alpha = rr0 / pAp.obten(0, 0);


            x = x.mas(p.por(alpha));
            r = r.mas(Ap.por(alpha));
            rr = (r.T()).por(r);

            //System.out.println("alpha = " + alpha);
            beta = rr.obten(0, 0) / rr0;

            p = (p.por(beta)).menos(r);

            System.out.println("beta = " + beta);

            System.out.println("p = ");
            p.imprime();

            System.out.println("rr = ");
            rr.imprime();

            rr0 = rr.obten(0, 0);
            System.out.println("Iteracion \t" + (k) + "\t f(x) = \t" + funcion(x));
            x.imprime();
        }

        return x;
    }

    /**
     * Calcula la magnitud de un vector
     * @param x Matriz vector a calcular su magnitud
     * @return double magnitud calculada
     */

    public double Magnitud(Matriz x)
    {
      double mag = (x.T().por(x)).obten(0,0);
      return(Math.sqrt(mag));
    }

    /**
     * funci�n f(x+alpha p) en la direcci�n p
     * @param x Matriz Valor inicial
     * @param alfa double tama�o de paso
     * @param p Matriz Direcci�n
     * @return double
     */

    public double phi(Matriz x, double alfa, Matriz p)
    {
        Matriz x1 = x.mas(p.por(alfa));
        return(funcion(x1));
    }

    /**
     * Derivada de la funci�n phi
     * @param x Matriz
     * @return double
     */

    public double phi_p(Matriz x)
    {
        Matriz g = Gradiente(x);
        Matriz p = Direccion(g);

        return ((p.T().por(g)).obten(0,0));
    }

    /**
     * Funci�n a minimizar
     * @param x Matriz
     * @return double
     */

    public double funcion(Matriz x)
    {
        return 0;
    }

    /**
     * Gradiente de la funci�n
     * @param x Matriz Valor de evaluaci�n
     * @return Matriz vector gradiente
     */

    public Matriz Gradiente(Matriz x)
    {
        return null;
    }

    /**
     * Calcula la matriz Hessiana
     * @param x Matriz punto de evalucci�n
     * @return Matriz Matriz calculada
     */

    public Matriz Hessiano(Matriz x)
    {
        return null;
    }

}