package optimizacion;
import operaciones.Matriz;

/**
 * <p>Title: Metodo de Punto interior</p>
 * <p>Description: Minimiza una funcion utilizando el algoritmo de punto interior</p>
 * <p>Copyright: Copyright (c) 2005</p>
 * <p>Company: UMSNH</p>
 * @author Dr. Felix Calderon Solorio
 * @version 1.0
 */

public class Restricciones {

        /**
         * n => numero de variables
         * m => numero de restricciones
         */

    int n, m;
    double T = 1e-6;
    final double UMBRALZERO = 1e-10;

    /**
     * Metodo de Barrera. Calcula el minimo de una funcion sujeta a una
     * a restricciones del tipo Ax - b > 0
     * @param x Matriz Punto inicial
     * @param mu double Parametro de Barrera
     * @return Matriz Soluci�n
     */

    public Matriz Barrera(Matriz x, double mu) {
        int i, j, iter;
        double tol, val;
        double Aux[][] = Restricciones();
        Matriz g, H, A, At, b, C, C1, C2, dx;

        n = Aux[0].length-1;
        m = Aux.length;

        A = new Matriz(m, n);
        b = new Matriz(m, 1);
        C1 = new Matriz(m,1);
        C2 = new Matriz(m,m);

        for (i = 0; i < m; i++) {
            b.inserta(i, 0, Aux[i][n]);
            for (j = 0; j < n; j++)
                A.inserta(i, j, Aux[i][j]);
        }

        At = A.T(); // calcula la transpuesta

        for(iter = 1; iter <1000; iter++)
        {
            C = (A.por(x)).menos(b);

            for (i = 0; i < m; i++) {
                val = 1.0 / C.obten(i, 0);
                C1.inserta(i, 0, val);
                C2.inserta(i, i, -val * val);
            }

            g = Gradiente(x).menos((A.T().por(C1)).por(mu));
            H = Hessiano(x).menos(((A.T().por(C2)).por(A)).por(mu));

            dx = g.entre(H);
            x = x.menos(dx);
            System.out.print("iteracion " + iter + " mu = " + mu +" X = ");
            x.T().imprime();
            if(Magnitud(g) < T) mu /= 1.5;
            if(mu < T) break;
        }
        return x;
    }

    /**
     * M�todo de Punto interior. Calcula el minino de una funci�n cuadr�tica con
     * restricciones del tipo Ax - b > 0
     * @param X Matriz Punto inicial
     * @param lo double[] Vector de l's inicial
     * @param sigma double tama�o de paso ruta central
     * @param alpha double tama�o de paso
     */

    public void interior(Matriz X, double lo[], double sigma, double alpha) {
        double Aux[][] = Restricciones();
        Matriz H = Hessiano(X), g;
        Matriz M, A, b, y, l, rd, rb, F, Delta;
        int i, j;
        double mu;

        n = H.nren;
        m = Aux.length;

        A = new Matriz(m, n);
        b = new Matriz(m, 1);

        for (i = 0; i < m; i++) {
            b.inserta(i, 0, Aux[i][n]);
            for (j = 0; j < n; j++)
                A.inserta(i, j, Aux[i][j]);
        }

        y = A.por(X).menos(b);
        l = new Matriz(lo);

        for (int k = 1; k < 1000; k++) {
            g = Gradiente(X);
            rd = g.menos(A.T().por(l));
            rb = (A.por(X)).menos(y).menos(b);
            mu = (y.T().por(l)).obten(0, 0) / m;

            M = llena(Aux, H, g, y, l);
            F = llena(rd, rb, l, y, sigma, mu);
            Delta = F.entre(M);

            //M.imprime();
            //F.imprime();
            //Delta.imprime();

            for (i = 0; i < n; i++)
                X.inserta(i, 0, X.obten(i, 0) + alpha * Delta.obten(i, 0));

            for (i = 0; i < m; i++) {
                l.inserta(i, 0, l.obten(i, 0) + alpha * Delta.obten(i + n, 0));
                y.inserta(i, 0,
                          y.obten(i, 0) + alpha * Delta.obten(i + n + m, 0));
            }

            System.out.println("Iteracion " + k);
            System.out.println("X =");
            X.T().imprime();
            System.out.println("Lambda =");
            l.T().imprime();
            System.out.println("y =");
            y.T().imprime();
            if (Mayor(Delta) < T)break;
        }
    }

    /**
     * Llena la matriz caracteristica para el m�todo de punto interior
     * @param rd Matriz
     * @param rb Matriz
     * @param l Matriz
     * @param y Matriz
     * @param sigma double
     * @param mu double
     * @return Matriz
     */

    private Matriz llena(Matriz rd, Matriz rb, Matriz l, Matriz y, double sigma,
                        double mu) {

        Matriz F = new Matriz(n + 2 * m, 1);
        int i;

        for (i = 0; i < n; i++)
            F.inserta(i, 0, -rd.obten(i, 0));

        for (i = 0; i < m; i++) {
            F.inserta(i + n, 0, -rb.obten(i, 0));
            F.inserta(i + n + m, 0, -l.obten(i, 0) * y.obten(i, 0) + sigma * mu);
        }
        return F;
    }

    /**
     * Llena la matriz caracteristica para el m�todo de Barrera
     * @param A double[][]
     * @param H Matriz
     * @param g Matriz
     * @param y Matriz
     * @param l Matriz
     * @return Matriz
     */

    private Matriz llena(double A[][], Matriz H, Matriz g, Matriz y, Matriz l) {
        int i, j;
        Matriz M = new Matriz(n + 2 * m, n + 2 * m);

        for (i = 0; i < n; i++)
            for (j = 0; j < n; j++) {
                //System.out.println(i + " " + j + " " + n);
                M.inserta(i, j, H.obten(i, j));
            }

        for (i = 0; i < m; i++)
            for (j = 0; j < n; j++) {
                M.inserta(i + n, j, A[i][j]);
                M.inserta(j, i + n, -A[i][j]);
            }

        for (i = 0; i < m; i++) {
            M.inserta(n + m + i, n + i, y.obten(i, 0));
            M.inserta(n + m + i, n + m + i, l.obten(i, 0));
            M.inserta(n + i, n + m + i, -1.0);
        }

        for (i = 0; i < m; i++)
            for (j = 0; j < m; j++)
                M.inserta(n + i, n + m + j, -1.0);

        return M;
    }

    /**
     * Imprime un arreglo en dos dimensiones
     * @param texto String
     * @param fuente double[][]
     */

    private static void  imprime(String texto, double fuente[][]) {
        System.out.println(texto);

        for (int i = 0; i < fuente.length; i++) {
            for (int j = 0; j < fuente[0].length; j++) {
                System.out.print(fuente[i][j]);
                if (j < fuente[0].length - 1) System.out.print(", ");
            }

            System.out.println("");
        }
    }

    /**
     * Regresa el mayor de los valores contenido en el vector X
     * @param x Matriz
     * @return double
     */

    private double Mayor(Matriz x) {
        double max = 0;
        for (int i = 0; i < x.nren; i++)
            if (max < Math.abs(x.obten(i, 0)))
                max = Math.abs(x.obten(i, 0));

        return max;
    }

    /**
     * Calcula la magnitud de un vector
     * @param x Matriz
     * @return double
     */

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

    /**
     * Funci�n de error
     * @param texto_error String
     */

    private void error( String texto_error )
    {
       System.err.println( "\n\nError en Clase Restricciones.java: " + texto_error + "." );
       System.exit( 1 );
    }

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

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

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

    public double[][] Restricciones() {
        return null;
    }
}