code: package gr.iti.mklab.reveal.forensics.util.dwt;
  
  /**
   * Copyright 2009 Mark Bishop 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: http://www.gnu.org/licenses
   * 
   * The author makes no warranty for the accuracy, completeness, safety, or
   * usefulness of any information provided and does not represent that its use
   * would not infringe privately owned right.
   */
  /**
   * 
   * Reference: Golub, Gene H., Charles F. Van Loan, Matrix Computations 3rd
   * Edition, The Johns Hopkins University Press, Baltimore, MD, 1996.
   * 
   * 
   * Class responsibility: Provide methods for basic vector and matrix
   * computations. Note: only a few of these methods are used in the DWT demo.
   * 
   * @author mark bishop
   *
   */
  public class MatrixOps {
  
  	/**
  	 * 
  	 * @return gives an upper bound on the relative error due to rounding in
  	 *         floating point arithmetic
  	 */
  	public static double machineEpsilonDouble() {
  		double eps = 1.0;
  		do
  			eps /= 2.0;
  		while ((1.0 + (eps / 2.0)) != 1.0);
  		return eps;
  	}
  
  	/**
  	 * 
  	 * @param n To be explained
  	 * @return Identity matrix I[n,n]
  	 */
  	public static double[][] eye(int n) {
  		double[][] eye = new double[n][n];
  		for (int i = 0; i < n; i++) {
  			eye[i][i] = 1.0;
  		}
  		return eye;
  	}
  
  	public static double[] getColumnAsVector(double[][] A, int colIndex) {
  		int m = A.length;
  		double[] col = new double[m];
  		for (int i = 0; i < m; i++) {
  			col[i] = A[i][colIndex];
  		}
  		return col;
  	}
  
  	public static double[] getRowAsVector(double[][] A, int rowIndex) {
  		int n = A[0].length;
  		System.out.print(A.length + "\n");
  		System.out.print(A[0].length + "\n");
  		double[] col = new double[n];
  		for (int i = 0; i < n; i++) {
  			col[i] = A[rowIndex][i];
  		}
  		return col;
  	}
  
  	public static double[][] vectorToRowMatrix(double[] v) {
  		int n = v.length;
  		double[][] row = new double[1][n];
  		for (int i = 0; i < n; i++) {
  			row[0][i] = v[i];
  		}
  		return row;
  	}
  
  	public static double[][] vectorToColumnMatrix(double[] v) {
  		int n = v.length;
  		double[][] row = new double[n][1];
  		for (int i = 0; i < n; i++) {
  			row[i][0] = v[i];
  		}
  		return row;
  	}
  
  	public static double[] deepCopy(double[] v) {
  		int n = v.length;
  		double[] copy = new double[n];
  		for (int j = 0; j < n; j++) {
  			copy[j] = v[j];
  		}
  		return copy;
  	}
  
  	public static double[][] deepCopy(double[][] A) {
  		int m = A.length;
  		int n = A[0].length;
  		double[][] copy = new double[m][n];
  		for (int i = 0; i < m; i++) {
  			for (int j = 0; j < n; j++) {
  				copy[i][j] = A[i][j];
  			}
  		}
  		return copy;
  	}
  
  	public static double[][] transpose(double[][] A) {
  		int m = A.length;
  		int n = A[0].length;
  		double[][] AT = new double[n][m];
  		for (int i = 0; i < n; i++) {
  			for (int j = 0; j < m; j++) {
  				AT[i][j] = A[j][i];
  			}
  		}
  		return AT;
  	}
  
  	public static double[][] scale(double alpha, double[][] A) {
  		int m = A.length;
  		int n = A[0].length;
  		double[][] alphaA = new double[m][n];
  		for (int i = 0; i < m; i++) {
  			for (int j = 0; j < n; j++) {
  				alphaA[i][j] = alpha * A[i][j];
  			}
  		}
  		return alphaA;
  	}
  
  	public static double[] scale(double alpha, double[] v) {
  		int n = v.length;
  		double[] alphaV = new double[n];
  		for (int i = 0; i < n; i++) {
  			alphaV[i] = alpha * v[i];
  		}
  		return alphaV;
  	}
  
  	public static double[][] add(double[][] A, double[][] B) {
  		int m = A.length;
  		int n = A[0].length;
  		double[][] C = new double[m][n];
  		for (int i = 0; i < m; i++) {
  			for (int j = 0; j < n; j++) {
  				C[i][j] = (A[i][j] + B[i][j]);
  			}
  		}
  		return C;
  	}
  
  	public static double[] add(double[] v1, double[] v2) {
  		int n = v1.length;
  		double[] v1Pv2 = new double[n];
  		for (int i = 0; i < n; i++) {
  			v1Pv2[i] = v1[i] + v2[i];
  		}
  		return v1Pv2;
  	}
  
  	public static double[][] multiply(double[][] A, double[][] B) {
  		int m = A.length;
  		int p = A[0].length;
  		int n = B[0].length;
  		double[][] C = new double[m][n];
  		for (int i = 0; i < m; i++) {
  			for (int k = 0; k < p; k++) {
  				for (int j = 0; j < n; j++) {
  					C[i][j] = C[i][j] + A[i][k] * B[k][j];
  				}
  			}
  		}
  		return C;
  	}
  
  	public static double[] multiply(double[][] A, double[] b) {
  		int m = A.length;
  		int n = A[0].length;
  		double[] v = new double[m];
  		for (int k = 0; k < n; k++) {
  			for (int i = 0; i < m; i++) {
  				v[i] += A[i][k] * b[k];
  			}
  		}
  		return v;
  	}
  
  	public static double[][] hadamard(double[][] A, double[][] B) {
  		int m = A.length;
  		int n = A[0].length;
  		double[][] C = new double[m][n];
  		for (int i = 0; i < m; i++) {
  			for (int j = 0; j < n; j++) {
  				C[i][j] = A[i][j] * B[i][j];
  			}
  		}
  		return C;
  	}
  
  	public static double[] hadamard(double[] v1, double[] v2) {
  		int m = v1.length;
  		double[] c = new double[m];
  		for (int i = 0; i < m; i++) {
  			c[i] = v1[i] * v2[i];
  		}
  		return c;
  	}
  
  	public static double vectorInnerProduct(double[] v1, double[] v2) {
  		int n = v1.length;
  		double ab = 0;
  		for (int i = 0; i < n; i++) {
  			ab += v1[i] * v2[i];
  		}
  		return ab;
  	}
  
  	public static double[][] vectorTensorProduct(double[] v1, double[] v2) {
  		int m = v1.length;
  		int n = v2.length;
  		double[][] c = new double[m][n];
  		for (int i = 0; i < m; i++) {
  			for (int j = 0; j < n; j++) {
  				c[i][j] += v1[i] * v2[j];
  			}
  		}
  		return c;
  	}
  
  	/**
  	 * 
  	 * @param A
  	 *            A(m,n)
  	 * @return maximum row sum
  	 */
  	public static double matrixInfinityNorm(double[][] A) {
  		int m = A.length;
  		double lFinity = 0;
  		for (int i = 0; i < m; i++) {
  			double vIN = vectorInfinityNorm(A[i]);
  			if (vIN > lFinity)
  				lFinity = vIN;
  		}
  		return lFinity;
  	}
  
  	public static double vector2Norm(double[] v) {
  		int n = v.length;
  		double norm = 0;
  		for (int i = 0; i < n; i++) {
  			norm += Math.pow(v[i], 2);
  		}
  		norm = Math.pow(norm, 0.5);
  		return norm;
  	}
  
  	public static double matrixFNorm(double[][] A) {
  		int m = A.length;
  		int n = A[0].length;
  		double norm = 0;
  		for (int i = 0; i < m; i++) {
  			for (int j = 0; j < n; j++) {
  				norm += (Math.pow(A[i][j], 2));
  			}
  		}
  		norm = Math.pow(norm, 0.5);
  		return norm;
  	}
  
  	/**
  	 * 
  	 * @param v
  	 *            vector v[n] as double[]
  	 * @return InfinityNorm(v) ... absolute value of highest magnitude element
  	 */
  	public static double vectorInfinityNorm(double[] v) {
  		double lFinity = 0;
  		int n = v.length;
  		lFinity = Math.abs(v[0]);
  		for (int i = 1; i < n; i++) {
  			double abs = Math.abs(v[i]);
  			if (abs > lFinity) {
  				lFinity = abs;
  			}
  		}
  		return lFinity;
  	}
  
  	/**
  	 * Golub, Algorithm 3.1.2
  	 * 
  	 * @param A
  	 *            A[m,n] where A is upper triangular
  	 * @param b
  	 *            b[m]
  	 * @return x[n] where Ax = b
  	 */
  	public static double[] backCalculateX(double[][] upperTriangularMatrix,
  			double[] b) {
  		int n = upperTriangularMatrix[0].length;
  		b[n - 1] = b[n - 1] / upperTriangularMatrix[n - 1][n - 1];
  		for (int i = n - 2; i >= 0; i += -1) {
  			double temp = 0;
  			for (int j = n - 1; j >= i + 1; j += -1) {
  				temp = temp - upperTriangularMatrix[i][j] * b[j];
  			}
  			temp = temp + b[i];
  			b[i] = temp / upperTriangularMatrix[i][i];
  		}
  		return b;
  	}
  
  	/**
  	 * 
  	 * @param x
  	 *            x(m)
  	 * @param order To be explained
  	 * @return V(m, order+1)such that V(i,j) = x(i)^j
  	 */
  	public static double[][] createVandermonde(double[] x, int order) {
  		int m = x.length;
  		int n = order;
  		double[][] V = new double[m][n + 1];
  		for (int i = 0; i < m; i++) {
  			for (int j = 0; j <= n; j++) {
  				V[i][j] = Math.pow(x[i], (double) j);
  			}
  		}
  		return V;
  	}
  }
  