package operaciones;


import java.io.IOException;
import javax.swing.JOptionPane;
import java.util.StringTokenizer;
import java.io.RandomAccessFile;

/**
 * <p>Title: Funciones para hacer operaciones con matrices en Arreglos
 * version 1.01</p>
 * <p>Description: Define las operaciones basicas con matrices</p>
 * <p>Copyright: Copyright (c) 2005</p>
 * <p>Company: UMSNH</p>
 * @author Dr. Felix Calderon Solorio
 * @version 1.01
 */

public class Matrices {

        /**
         * Metodo de Cholesky Incompleto
         * @param A double[][]
         */

        public static void CholeskyI(double A[][]) {
        int i, j, k, n, s;
        double fact, suma = 0;
        n = A.length;

        for (i = 0; i < n; i++) {
            for (j = 0; j <= i - 1; j++) {
                suma = 0;
                for (k = 0; k <= j - 1; k++)
                    suma += A[i][k] * A[j][k];
                if (A[i][j] != 0)
                    A[i][j] = (A[i][j] - suma) / A[j][j];
            }
            suma = 0;
            for (k = 0; k <= i - 1; k++)
                suma += A[i][k] * A[i][k];
            A[i][i] = Math.sqrt(A[i][i] - suma);
        }

        for (j = 0; j < n; j++) {
            for (k = 1 + j; k < n; k++)
                A[j][k] = 0;
        }
    }

    /**
     * Metodo de Cholesky completo
     * @param A double[][]
     */

    public static void Cholesky(double A[][]) {
        int i, j, k, n, s;
        double fact, suma = 0;
        n = A.length;

        for (i = 0; i < n; i++) {
            for (j = 0; j <= i - 1; j++) {
                suma = 0;
                for (k = 0; k <= j - 1; k++)
                    suma += A[i][k] * A[j][k];

                A[i][j] = (A[i][j] - suma) / A[j][j];
            }
            suma = 0;
            for (k = 0; k <= i - 1; k++)
                suma += A[i][k] * A[i][k];
            A[i][i] = Math.sqrt(A[i][i] - suma);
        }

        for (j = 0; j < n; j++) {
            for (k = 1 + j; k < n; k++)
                A[j][k] = 0;
        }
    }

    /**
     * Copia un arreglo fuente en un arreglo destino
     * @param destino int[]
     * @param fuente int[]
     */

    public static void copia(int destino[], int fuente[]) {
        for (int i = 0; i < fuente.length; i++)
            destino[i] = fuente[i];
    }

    /**
     * Copia un arreglo bidimensional fuente en un arreglo destino
     * @param destino double[][]
     * @param fuente double[][]
     */

    public static void copia(double destino[][], double fuente[][]) {
        for (int i = 0; i < fuente.length; i++)
            for (int j = 0; j < fuente[0].length; j++)
                destino[i][j] = fuente[i][j];
    }

    /**
     * Lee las dimensiones de un arreglo bidimensional en un archivo de texto
     * @param archivo String
     * @return double[][]
     */

    static private double[][] Dimensiones(String archivo) {
        StringTokenizer st;
        String linea;
        int num_ren = 0, num_col = 0;
        double salida[][] = null;

        try {
            RandomAccessFile DIS = new RandomAccessFile(archivo, "r");

            while ((linea = DIS.readLine()) != null) {
                st = new StringTokenizer(linea, " ");
                if (st.countTokens() != 0) num_ren++;

                if (num_col == 0) {
                    while (st.hasMoreTokens()) {
                        num_col++;
                        st.nextToken(" ");
                    }
                }
            }
            salida = new double[num_ren][num_col];
            DIS.close();
        } catch (IOException e) {
            JOptionPane.showMessageDialog(null, "Error" + e.toString(), "ERROR",
                                          JOptionPane.ERROR_MESSAGE);
        }
        return salida;
    }

    /**
     * Calcula la inversa de la matriz A y deja el resultado en Ainv
     * @param A double[][]
     * @param Ainv double[][]
     */

    static public void inversa(double A[][], double Ainv[][]) {
        int n = A.length;
        int k, j, i;

        for (i = 0; i < n; i++)
            for (j = 0; j < n; j++)
                Ainv[i][j] = A[i][j];

        for (k = 0; k < n; k++) {
            for (i = 0; i < n; i++)
                for (j = 0; j < n; j++) {
                    if ((i != k) && (j != k))
                        Ainv[i][j] -= (Ainv[i][k] * Ainv[k][j]) /
                                Ainv[k][k];
                }

            for (j = 0; j < n; j++) {
                if (j != k)
                    Ainv[k][j] = -Ainv[k][j] / Ainv[k][k];
            }

            for (i = 0; i < n; i++) {
                if (i != k)
                    Ainv[i][k] = Ainv[i][k] / Ainv[k][k];
            }

            Ainv[k][k] = 1 / Ainv[k][k];
        }
    }

    /**
     * Imprime un arreglo
     * @param texto String Mensaje
     * @param fuente double[] Datos
     */

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

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

    }

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

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

    /**
     * Redondea un numero a dos cuatro cifras significativas
     * @param x double
     * @return double
     */

    public static double redondea(double x)
    {
        int aux;
        double resul;
        aux = (int)(x*10000.0 + 0.5);
        resul = aux/10000.0;
        return resul;
    }

    /**
     * Imprime un arreglo bidimensional
     * @param texto String Mensaje
     * @param fuente double[][] Datos a imprimir
     */

    public 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(redondea(fuente[i][j]));
                if(j < fuente[0].length-1) System.out.print(", ");
            }

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

    /**
     * Imprime un arreglo bidimensional
     * @param texto String Mensaje
     * @param fuente double[][] Datos a imprimir
     */

    public static void Guarda(String name, double fuente[][])
    {
        String aux;
        try
        {
          RandomAccessFile DIS = new RandomAccessFile(name, "rw");
          for (int i = 0; i < fuente.length; i++) {
              for (int j = 0; j < fuente[0].length; j++) {
                  aux = fuente[i][j] + " ";
                  DIS.writeBytes(aux);
              }

              DIS.writeBytes("\n");
          }
          DIS.close();
        }
        catch(IOException e)
        {
          System.out.println("no pude abrir el archivo");
        }
    }


    /**
     * Lee un archivo de texto con la informacion de una matriz
     * @param archivo String Nombre y ruta de los datos
     * @return double[][] Datos leidos del archivo
     */

    static public double[][] Lee(String archivo) {
        String linea = "", dato = "";
        double entrada[][] = null;
        int i = 0, j = 0;
        StringTokenizer st;
        int nren, ncol;

        entrada = Dimensiones(archivo);
        nren = entrada.length;
        ncol = entrada[0].length;
        //System.out.println(nren + " " + ncol);

        try {

            RandomAccessFile DIS = new RandomAccessFile(archivo, "r");

            i = 0;
            while (((linea = DIS.readLine()) != null) && i < nren) {
                st = new StringTokenizer(linea, " ");
                if (st.countTokens() != 0) {

                    st = new StringTokenizer(linea, " ");
                    for (j = 0; j < ncol; j++) {
                        dato = st.nextToken(" ");
                        entrada[i][j] = Double.valueOf(dato).doubleValue();
                    }
                    i++;
                }
            }
            DIS.close();
        } catch (IOException e) {
            JOptionPane.showMessageDialog(null, "Error" + e.toString(), "ERROR",
                                          JOptionPane.ERROR_MESSAGE);
        }
        return entrada;
    }

    /**
     * Norma de un vector
     * @param v double[] Vector de datos
     * @return double resultado
     */

    public static double norma(double[] v) {
        double sum = 0;
        for (int i = 0; i < v.length; i++)
            sum += v[i] * v[i];
        return Math.sqrt(sum);
    }

    /**
     * Realiza la muntiplicacion de dos matrices
     * @param A double[][]
     * @param B double[][]
     * @return double[][] producto
     */

    public static double[][] multiplica(double A[][], double B[][]) {

        int nr = A.length, nc = B[0].length, m = B.length;
        int i, j, k;

        double res[][] = new double[nr][nc];

        for (i = 0; i < nr; i++) {
            for (j = 0; j < nc; j++) {
                res[i][j] = 0;

                for (k = 0; k < m; k++)
                    res[i][j] += A[i][k] * B[k][j];
            }
        }
        return res;
    }

    /**
     * Multiplica un arreglo bidimensional por un vector
     * @param A double[][] arreglo bidimensional
     * @param b double[] vecor
     * @return double[] resultado
     */

    public static double[] multiplica(double A[][], double b[]) {

        int nr = A.length, nc = A[0].length;

        double res[] = new double[nr];

        for (int i = 0; i < nr; i++) {
            res[i] = 0;

            for (int j = 0; j < nc; j++)
                res[i] += A[i][j] * b[j];
        }
        return res;
    }

    /**
     * Multiplica un vector por otro vector
     * @param v1 double[]
     * @param v2 double[]
     * @return double
     */


    public static double multiplica(double[] v1, double[] v2) {
        double resultado = 0;

        for (int i = 0; i < v1.length; i++)
            resultado += v1[i] * v2[i];
        return resultado;
    }

    public static double[] multiplica(double[] v, double e) {
        double[] res = new double[v.length];

        for (int i = 0; i < res.length; i++)
            res[i] = e * v[i];

        return res;
    }

    /**
     * Calcula la solucion de un sistema de ecuaciones Ax = b
     * @param a double[][]
     * @param x double[]
     * @param b double[]
     */

    static public void SolucionSL(double a[][], double x[], double b[]) {
        int n = a.length, i, j, k;
        double suma;

        for (k = 0; k <= n - 2; k++) {
            for (i = k + 1; i <= (n - 1); i++) {
                b[i] -= a[i][k] * b[k] / a[k][k];
                for (j = n - 1; j >= k; j--)
                    a[i][j] -= a[i][k] * a[k][j] / a[k][k];
            }
        }
        for (k = n - 1; k >= 0; k--) {
            suma = 0;
            for (j = k + 1; j < n; j++)
                suma += a[k][j] * x[j];
            x[k] = (b[k] - suma) / a[k][k];
        }
    }
    /**
     * Calcula el determinante de una matriz
     * @param entrada double[][] matriz
     * @return double regresa el resultado
     */
    static public double Determinante(double entrada[][])
    {
        int n = entrada.length, i, j, k;
        double suma;
        double a[][] = new double[n][n];
        for(i=0; i<n; i++)
            for(j=0; j<n; j++)
                a[i][j] = entrada[i][j];

        for (k = 0; k <= n - 2; k++) {
            for (i = k + 1; i <= (n - 1); i++) {
                for (j = n - 1; j >= k; j--)
                    a[i][j] -= a[i][k] * a[k][j] / a[k][k];
            }
        }
        suma = a[0][0];
        for(k=1; k<n; k++)
            suma *= a[k][k];
        return suma;
    }
    /**
     * Metodo sustitucion hacia adelante para solucionar un sistema
     * triangular inferior Ay = b
     * @param A double[][]
     * @param y double[]
     * @param b double[]
     */

    public static void susA(double A[][], double y[], double b[]) {
        int n, i, j;

        n = A.length;
        double s;
        i = n;
        j = 0;

        for (i = 0; i < n; i++) {
            s = 0;
            for (j = 0; j < i; j++) {
                s += A[i][j] * y[j];
            }
            y[i] = (b[i] - s) / A[i][i];
        }
    }

    /**
     * Metodo de sustitucion hacia atras para solucionar un sistema de
     * ecuaciones triangula superior Ay= b
     * @param A double[][]
     * @param y double[]
     * @param b double[]
     */

    public static void susAt(double A[][], double y[], double b[]) {
        int n, i, j;

        n = A.length;
        double s;
        i = n;
        j = 0;

        for (i = n - 1; i >= 0; i--) {
            s = 0;
            for (j = n - 1; j > i; j--) {
                s += A[j][i] * y[j];
            }
            y[i] = (b[i] - s) / A[i][i];
        }
    }

    /**
     * Realiza la resta de dos vectores
     * @param minuendo double[]
     * @param sustraendo double[]
     * @return double[]
     */

    public static double[] resta(double[] minuendo, double[] sustraendo) {

        double[] res = new double[minuendo.length];

        for (int i = 0; i < res.length; i++)
            res[i] = minuendo[i] - sustraendo[i];
        return res;
    }

    /**
     * Realiza la suma de dos vectores
     * @param v1 double[]
     * @param v2 double[]
     * @return double[]
     */

    public static double[] suma(double[] v1, double[] v2) {

        double[] res = new double[v1.length];

        for (int i = 0; i < res.length; i++)
            res[i] = v1[i] + v2[i];
        return res;
    }

    static public void quiksort(int lo,int ho, double valor[], int indice[])
    {
        int l = lo, h = ho;
        double aux, mid;

        if (ho > lo) {
            mid = valor[indice[(lo + ho) / 2]];
            while (l < h) {
                while ((l < ho) && (valor[indice[l]] > mid))++l;
                while ((h > lo) && (valor[indice[h]] < mid))--h;
                if (l <= h) {
                    aux = indice[l];
                    indice[l] = indice[h];
                    indice[h] = (int) aux;

                    ++l;
                    --h;
                }
            }

            if (lo < h)
                quiksort(lo, h, valor, indice);
            if (l < ho)
                quiksort(l, ho, valor, indice);
        }
    }

    /**
     * Calcula la Seudo-inversa de una matriz definida positiva
     * @param A double[][] matriz inicial
     * @param A_inv double[][] matriz resultado A+
     * @param N int tamano de la matriz
     * @return double regresa el determinante de A
     */


    public static double svdInversa(double A[][], double A_inv[][], int N)
    {
        double U[][] = new double[N][N], W[] = new double[N];
        int indice[] = new int[N];
        double V[][] = new double[N][N];
        double temp, det, suma, sumaT;
        int i, j, k;

        for(i=0; i<N; i++)
        {
            indice[i] = i;
            for (j = 0; j < N; j++)
                U[i][j] = A[i][j];
        }

        svdcmp(U, N, N, W, V);

        quiksort(0, N-1, W, indice);

        sumaT = 0;
        for(i=0; i<N; i++)
        {
            if(W[indice[i]] > 0)
                sumaT += W[indice[i]];
            else W[indice[i]] = 0;
        }

        suma =0; det = 1;
        for(i=0; i<N; i++)
        {
            if (suma/sumaT < 0.99)
            {
                suma += W[indice[i]];
                det *= W[indice[i]];
                W[indice[i]] = 1.0 / W[indice[i]];
            }
            else W[indice[i]] = 0;
        }

        // Calcula W*U'

        for(i=0; i<N; i++)
            for(j=0; j<N; j++)
            {
                if(j >= i)
                {
                    temp = U[i][j];
                    U[i][j] = U[j][i];
                    U[j][i] = temp;
                }
                U[i][j] *= W[i];
            }
        // Calcula finalmente VWU'

        for (i = 0; i < N; i++) {
            for (j = 0; j < N; j++) {
                A_inv[i][j] = 0;

                for (k = 0; k < N; k++)
                    A_inv[i][j] += V[i][k] * U[k][j];
            }
        }
        return det;
    }

    /**
     * Descomposicion en valores singulares de para una matriz
     * Dada la matriz a esta rutina calcula los valores singulares
     * a = uwvT. La matriz u remplaza a la matriz a
     * @param a double[][] matriz de entrada a y salida u
     * @param m int
     * @param n int
     * @param w double[] diagonal
     * @param v double[][] matriz v
     */

    public static void svdcmp(double a[][], int m, int n, double w[], double v[][]) {
        int flag, i, its, j, jj, k, l = 0, nm = 0;
        double anorm, c, f, g, h, s, scale, x, y, z, rv1[];

        rv1 = new double[n];
        g = scale = anorm = 0.0; /* Householder reduction to bidiagonal form */
        for (i = 0; i < n; i++) {
            l = i + 1;
            rv1[i] = scale * g;
            g = s = scale = 0.0;
            if (i < m) {
                for (k = i; k < m; k++)
                    scale += Math.abs(a[k][i]);
                if (scale != 0) {
                    for (k = i; k < m; k++) {
                        a[k][i] /= scale;
                        s += a[k][i] * a[k][i];
                    }
                    f = a[i][i];
                    g = -SIGN(Math.sqrt(s), f);
                    h = f * g - s;
                    a[i][i] = f - g;
                    for (j = l; j < n; j++) {
                        for (s = 0.0, k = i; k < m; k++)
                            s += a[k][i] * a[k][j];
                        f = s / h;
                        for (k = i; k < m; k++)
                            a[k][j] += f * a[k][i];
                    }
                    for (k = i; k < m; k++)
                        a[k][i] *= scale;
                }
            }
            w[i] = scale * g;
            g = s = scale = 0.0;
            if (i <= m && i != n) {
                for (k = l; k < n; k++)
                    scale += Math.abs(a[i][k]);
                if (scale != 0) {
                    for (k = l; k < n; k++) {
                        a[i][k] /= scale;
                        s += a[i][k] * a[i][k];
                    }
                    f = a[i][l];
                    g = -SIGN(Math.sqrt(s), f);
                    h = f * g - s;
                    a[i][l] = f - g;
                    for (k = l; k < n; k++)
                        rv1[k] = a[i][k] / h;
                    for (j = l; j < m; j++) {
                        for (s = 0.0, k = l; k < n; k++)
                            s += a[j][k] * a[i][k];
                        for (k = l; k < n; k++)
                            a[j][k] += s * rv1[k];
                    }
                    for (k = l; k < n; k++)
                        a[i][k] *= scale;
                }
            }
            anorm = DMAX(anorm, (Math.abs(w[i]) + Math.abs(rv1[i])));
        }
        for (i = n - 1; i >= 0; i--) {
            /* Accumulation of right-hand transformations. */
            if (i < n - 1) {
                if (g != 0) {
                    for (j = l; j < n; j++)
                        /* Double division to avoid possible underflow. */
                        v[j][i] = (a[i][j] / a[i][l]) / g;
                    for (j = l; j < n; j++) {
                        for (s = 0.0, k = l; k < n; k++)
                            s += a[i][k] * v[k][j];
                        for (k = l; k < n; k++)
                            v[k][j] += s * v[k][i];
                    }
                }
                for (j = l; j < n; j++)
                    v[i][j] = v[j][i] = 0.0;
            }
            v[i][i] = 1.0;
            g = rv1[i];
            l = i;
        }
        for (i = (int) IMIN(m - 1, n - 1); i >= 0; i--) {
            /* Accumulation of left-hand transformations. */
            l = i + 1;
            g = w[i];
            for (j = l; j < n; j++)
                a[i][j] = 0.0;
            if (g != 0) {
                g = 1.0 / g;
                for (j = l; j < n; j++) {
                    for (s = 0.0, k = l; k < m; k++)
                        s += a[k][i] * a[k][j];
                    f = (s / a[i][i]) * g;
                    for (k = i; k < m; k++)
                        a[k][j] += f * a[k][i];
                }
                for (j = i; j < m; j++)
                    a[j][i] *= g;
            } else
                for (j = i; j < m; j++)
                    a[j][i] = 0.0;
            ++a[i][i];
        }
        for (k = n - 1; k > 0; k--) {
            /* Diagonalization of the bidiagonal form. */
            for (its = 1; its <= 30; its++) {
                flag = 1;
                for (l = k; l >= 0; l--) { /* Test for splitting. */
                    nm = l - 1; /* Note that rv1[1] is always zero. */
                    if ((double) (Math.abs(rv1[l]) + anorm) == anorm) {
                        flag = 0;
                        break;
                    }
                    if ((double) (Math.abs(w[nm]) + anorm) == anorm)
                        break;
                }
                if (flag != 0) {
                    c = 0.0; /* Cancellation of rv1[l], if l > 1. */
                    s = 1.0;
                    for (i = l; i <= k; i++) {
                        f = s * rv1[i];
                        rv1[i] = c * rv1[i];
                        if ((double) (Math.abs(f) + anorm) == anorm)
                            break;
                        g = w[i];
                        h = pythag(f, g);
                        w[i] = h;
                        h = 1.0 / h;
                        c = g * h;
                        s = -f * h;
                        for (j = 0; j < m; j++) {
                            y = a[j][nm];
                            z = a[j][i];
                            a[j][nm] = y * c + z * s;
                            a[j][i] = z * c - y * s;
                        }
                    }
                }
                z = w[k];
                if (l == k) { /* Convergence. */
                    if (z < 0.0) { /* Singular value is made nonnegative. */
                        w[k] = -z;
                        for (j = 0; j < n; j++)
                            v[j][k] = -v[j][k];
                    }
                    break;
                }
                if (its == 30)
                    System.out.println("no convergence in 30 svdcmp iterations");
                x = w[l]; /* Shift from bottom 2-by-2 minor. */
                nm = k - 1;
                y = w[nm];
                g = rv1[nm];
                h = rv1[k];
                f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
                g = pythag(f, 1.0);
                f = ((x - z) * (x + z) + h * ((y / (f + SIGN(g, f))) - h)) / x;
                c = s = 1.0; /* Next QR transformation: */
                for (j = l; j <= nm; j++) {
                    i = j + 1;
                    g = rv1[i];
                    y = w[i];
                    h = s * g;
                    g = c * g;
                    z = pythag(f, h);
                    rv1[j] = z;
                    c = f / z;
                    s = h / z;
                    f = x * c + g * s;
                    g = g * c - x * s;
                    h = y * s;
                    y *= c;
                    for (jj = 0; jj < n; jj++) {
                        x = v[jj][j];
                        z = v[jj][i];
                        v[jj][j] = x * c + z * s;
                        v[jj][i] = z * c - x * s;
                    }
                    z = pythag(f, h);
                    w[j] = z; /* Rotation can be arbitrary if z = 0. */
                    if (z != 0) {
                        z = 1.0 / z;
                        c = f * z;
                        s = h * z;
                    }
                    f = c * g + s * y;
                    x = c * y - s * g;
                    for (jj = 0; jj < m; jj++) {
                        y = a[jj][j];
                        z = a[jj][i];
                        a[jj][j] = y * c + z * s;
                        a[jj][i] = z * c - y * s;
                    }
                }
                rv1[l] = 0.0;
                rv1[k] = f;
                w[k] = x;
            }
        }
    }

    static public double SIGN(double a, double b) {
        return ((b) >= 0.0 ? Math.abs(a) : -Math.abs(a));
    }

    static public double DMAX(double a, double b) {
        return (a > b ? a : b);
    }

    static public double IMIN(double a, double b) {
        return (a < b ? a : b);
    }

    static public double pythag(double a, double b) {
        double absa, absb;
        absa = Math.abs(a);
        absb = Math.abs(b);
        if (absa > absb)
            return absa * Math.sqrt(1.0 + (absb / absa) * (absb / absa));
        else
            return (absb == 0.0 ? 0.0 :
                    absb * Math.sqrt(1.0 + (absa / absb) * (absa / absb)));
    }

    static public double [][] transpuesta(double origen[][])
    {
        int nr = origen.length, nc = origen[0].length, i, j;

        double res[][] = new double[nc][nr];

        for(i=0; i<nr; i++)
            for(j=0; j<nc; j++)
                res[j][i] = origen[i][j];

        return res;

    }
}