/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package org.apache.commons.math.linear;

import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.util.MathUtils;

/**
 * Calculates the Singular Value Decomposition of a matrix.
 * <p>
 * The Singular Value Decomposition of matrix A is a set of three matrices: U,
 * &Sigma; and V such that A = U &times; &Sigma; &times; V<sup>T</sup>. Let A be
 * an m &times; n matrix, then U is an m &times; m orthogonal matrix, &Sigma; is
 * a m &times; n diagonal matrix with positive diagonal elements, and V is an n
 * &times; n orthogonal matrix.
 * </p>
 *
 * @version $Revision$ $Date$
 * @since 2.0
 */
public class SingularValueDecompositionImpl implements SingularValueDecomposition {

	/** Number of rows of the initial matrix. */
	private int m;

	/** Number of columns of the initial matrix. */
	private int n;

	/** Transformer to bidiagonal. */
	private BiDiagonalTransformer transformer;

	/** Main diagonal of the bidiagonal matrix. */
	private double[] mainBidiagonal;

	/** Secondary diagonal of the bidiagonal matrix. */
	private double[] secondaryBidiagonal;

	/** Main diagonal of the tridiagonal matrix. */
	private double[] mainTridiagonal;

	/** Secondary diagonal of the tridiagonal matrix. */
	private double[] secondaryTridiagonal;

	/** Eigen decomposition of the tridiagonal matrix. */
	private EigenDecomposition eigenDecomposition;

	/** Singular values. */
	private double[] singularValues;

	/** Cached value of U. */
	private RealMatrix cachedU;

	/** Cached value of U<sup>T</sup>. */
	private RealMatrix cachedUt;

	/** Cached value of S. */
	private RealMatrix cachedS;

	/** Cached value of V. */
	private RealMatrix cachedV;

	/** Cached value of V<sup>T</sup>. */
	private RealMatrix cachedVt;

	/**
	 * Calculates the Singular Value Decomposition of the given matrix.
	 * 
	 * @param matrix
	 *            The matrix to decompose.
	 * @exception InvalidMatrixException
	 *                (wrapping a {@link ConvergenceException} if algorithm
	 *                fails to converge
	 */
	public SingularValueDecompositionImpl(RealMatrix matrix) throws InvalidMatrixException {

		m = matrix.getRowDimension();
		n = matrix.getColumnDimension();

		cachedU = null;
		cachedS = null;
		cachedV = null;
		cachedVt = null;

		// transform the matrix to bidiagonal
		transformer = new BiDiagonalTransformer(matrix);
		mainBidiagonal = transformer.getMainDiagonalRef();
		secondaryBidiagonal = transformer.getSecondaryDiagonalRef();

		// compute Bt.B (if upper diagonal) or B.Bt (if lower diagonal)
		mainTridiagonal = new double[mainBidiagonal.length];
		secondaryTridiagonal = new double[mainBidiagonal.length - 1];
		double a = mainBidiagonal[0];
		mainTridiagonal[0] = a * a;
		for (int i = 1; i < mainBidiagonal.length; ++i) {
			final double b = secondaryBidiagonal[i - 1];
			secondaryTridiagonal[i - 1] = a * b;
			a = mainBidiagonal[i];
			mainTridiagonal[i] = a * a + b * b;
		}

		// compute singular values
		eigenDecomposition = new EigenDecompositionImpl(mainTridiagonal, secondaryTridiagonal, MathUtils.SAFE_MIN);
		singularValues = eigenDecomposition.getRealEigenvalues();
		for (int i = 0; i < singularValues.length; ++i) {
			singularValues[i] = Math.sqrt(singularValues[i]);
		}

	}

	/** {@inheritDoc} */
	public RealMatrix getU() throws InvalidMatrixException {

		if (cachedU == null) {

			if (m >= n) {
				// the tridiagonal matrix is Bt.B, where B is upper bidiagonal
				final double[][] eData = eigenDecomposition.getV().getData();
				final double[][] iData = new double[m][];
				double[] ei1 = eData[0];
				iData[0] = ei1;
				for (int i = 0; i < n - 1; ++i) {
					// compute Bt.E.S^(-1) where E is the eigenvectors matrix
					// we reuse the array from matrix E to store the result
					final double[] ei0 = ei1;
					ei1 = eData[i + 1];
					iData[i + 1] = ei1;
					for (int j = 0; j < n; ++j) {
						ei0[j] = (mainBidiagonal[i] * ei0[j] + secondaryBidiagonal[i] * ei1[j]) / singularValues[j];
					}
				}
				// last row
				final double lastMain = mainBidiagonal[n - 1];
				for (int j = 0; j < n; ++j) {
					ei1[j] *= lastMain / singularValues[j];
				}
				for (int i = n; i < m; ++i) {
					iData[i] = new double[n];
				}
				cachedU = transformer.getU().multiply(MatrixUtils.createRealMatrix(iData));
			} else {
				// the tridiagonal matrix is B.Bt, where B is lower bidiagonal
				cachedU = transformer.getU().multiply(eigenDecomposition.getV());
			}

		}

		// return the cached matrix
		return cachedU;

	}

	/** {@inheritDoc} */
	public RealMatrix getUT() throws InvalidMatrixException {

		if (cachedUt == null) {
			cachedUt = getU().transpose();
		}

		// return the cached matrix
		return cachedUt;

	}

	/** {@inheritDoc} */
	public RealMatrix getS() throws InvalidMatrixException {

		if (cachedS == null) {

			// cache the matrix for subsequent calls
			cachedS = MatrixUtils.createRealDiagonalMatrix(singularValues);

		}
		return cachedS;
	}

	/** {@inheritDoc} */
	public double[] getSingularValues() throws InvalidMatrixException {
		return singularValues.clone();
	}

	/** {@inheritDoc} */
	public RealMatrix getV() throws InvalidMatrixException {

		if (cachedV == null) {

			if (m >= n) {
				// the tridiagonal matrix is Bt.B, where B is upper bidiagonal
				cachedV = transformer.getV().multiply(eigenDecomposition.getV());
			} else {
				// the tridiagonal matrix is B.Bt, where B is lower bidiagonal
				final double[][] eData = eigenDecomposition.getV().getData();
				final double[][] iData = new double[n][];
				double[] ei1 = eData[0];
				iData[0] = ei1;
				for (int i = 0; i < m - 1; ++i) {
					// compute Bt.E.S^(-1) where E is the eigenvectors matrix
					// we reuse the array from matrix E to store the result
					final double[] ei0 = ei1;
					ei1 = eData[i + 1];
					iData[i + 1] = ei1;
					for (int j = 0; j < m; ++j) {
						ei0[j] = (mainBidiagonal[i] * ei0[j] + secondaryBidiagonal[i] * ei1[j]) / singularValues[j];
					}
				}
				// last row
				final double lastMain = mainBidiagonal[m - 1];
				for (int j = 0; j < m; ++j) {
					ei1[j] *= lastMain / singularValues[j];
				}
				for (int i = m; i < n; ++i) {
					iData[i] = new double[m];
				}
				cachedV = transformer.getV().multiply(MatrixUtils.createRealMatrix(iData));
			}

		}

		// return the cached matrix
		return cachedV;

	}

	/** {@inheritDoc} */
	public RealMatrix getVT() throws InvalidMatrixException {

		if (cachedVt == null) {
			cachedVt = getV().transpose();
		}

		// return the cached matrix
		return cachedVt;

	}

	/** {@inheritDoc} */
	public RealMatrix getCovariance(final double minSingularValue) {

		// get the number of singular values to consider
		int dimension = 0;
		while ((dimension < n) && (singularValues[dimension] >= minSingularValue)) {
			++dimension;
		}

		if (dimension == 0) {
			throw MathRuntimeException.createIllegalArgumentException(
					"cutoff singular value is {0}, should be at most {1}", minSingularValue, singularValues[0]);
		}

		final double[][] data = new double[dimension][n];
		getVT().walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
			/** {@inheritDoc} */
			@Override
			public void visit(final int row, final int column, final double value) {
				data[row][column] = value / singularValues[row];
			}
		}, 0, dimension - 1, 0, n - 1);

		RealMatrix jv = new Array2DRowRealMatrix(data, false);
		return jv.transpose().multiply(jv);

	}

	/** {@inheritDoc} */
	public double getNorm() throws InvalidMatrixException {
		return singularValues[0];
	}

	/** {@inheritDoc} */
	public double getConditionNumber() throws InvalidMatrixException {
		return singularValues[0] / singularValues[singularValues.length - 1];
	}

	/** {@inheritDoc} */
	public int getRank() throws IllegalStateException {

		final double threshold = Math.max(m, n) * Math.ulp(singularValues[0]);

		for (int i = singularValues.length - 1; i >= 0; --i) {
			if (singularValues[i] > threshold) {
				return i + 1;
			}
		}
		return 0;

	}

	/** {@inheritDoc} */
	public DecompositionSolver getSolver() {
		return new Solver(singularValues, getUT(), getV(), getRank() == singularValues.length);
	}

	/** Specialized solver. */
	private static class Solver implements DecompositionSolver {

		/** Singular values. */
		private final double[] singularValues;

		/** U<sup>T</sup> matrix of the decomposition. */
		private final RealMatrix uT;

		/** V matrix of the decomposition. */
		private final RealMatrix v;

		/** Singularity indicator. */
		private boolean nonSingular;

		/**
		 * Build a solver from decomposed matrix.
		 * 
		 * @param singularValues
		 *            singularValues
		 * @param uT
		 *            U<sup>T</sup> matrix of the decomposition
		 * @param v
		 *            V matrix of the decomposition
		 * @param nonSingular
		 *            singularity indicator
		 */
		private Solver(final double[] singularValues, final RealMatrix uT, final RealMatrix v,
				final boolean nonSingular) {
			this.singularValues = singularValues;
			this.uT = uT;
			this.v = v;
			this.nonSingular = nonSingular;
		}

		/**
		 * Solve the linear equation A &times; X = B in least square sense.
		 * <p>
		 * The m&times;n matrix A may not be square, the solution X is such that
		 * ||A &times; X - B|| is minimal.
		 * </p>
		 * 
		 * @param b
		 *            right-hand side of the equation A &times; X = B
		 * @return a vector X that minimizes the two norm of A &times; X - B
		 * @exception IllegalArgumentException
		 *                if matrices dimensions don't match
		 * @exception InvalidMatrixException
		 *                if decomposed matrix is singular
		 */
		public double[] solve(final double[] b) throws IllegalArgumentException, InvalidMatrixException {

			if (b.length != singularValues.length) {
				throw MathRuntimeException.createIllegalArgumentException(
						"vector length mismatch: got {0} but expected {1}", b.length, singularValues.length);
			}

			final double[] w = uT.operate(b);
			for (int i = 0; i < singularValues.length; ++i) {
				final double si = singularValues[i];
				if (si == 0) {
					throw new SingularMatrixException();
				}
				w[i] /= si;
			}
			return v.operate(w);

		}

		/**
		 * Solve the linear equation A &times; X = B in least square sense.
		 * <p>
		 * The m&times;n matrix A may not be square, the solution X is such that
		 * ||A &times; X - B|| is minimal.
		 * </p>
		 * 
		 * @param b
		 *            right-hand side of the equation A &times; X = B
		 * @return a vector X that minimizes the two norm of A &times; X - B
		 * @exception IllegalArgumentException
		 *                if matrices dimensions don't match
		 * @exception InvalidMatrixException
		 *                if decomposed matrix is singular
		 */
		public RealVector solve(final RealVector b) throws IllegalArgumentException, InvalidMatrixException {

			if (b.getDimension() != singularValues.length) {
				throw MathRuntimeException.createIllegalArgumentException(
						"vector length mismatch: got {0} but expected {1}", b.getDimension(), singularValues.length);
			}

			final RealVector w = uT.operate(b);
			for (int i = 0; i < singularValues.length; ++i) {
				final double si = singularValues[i];
				if (si == 0) {
					throw new SingularMatrixException();
				}
				w.setEntry(i, w.getEntry(i) / si);
			}
			return v.operate(w);

		}

		/**
		 * Solve the linear equation A &times; X = B in least square sense.
		 * <p>
		 * The m&times;n matrix A may not be square, the solution X is such that
		 * ||A &times; X - B|| is minimal.
		 * </p>
		 * 
		 * @param b
		 *            right-hand side of the equation A &times; X = B
		 * @return a matrix X that minimizes the two norm of A &times; X - B
		 * @exception IllegalArgumentException
		 *                if matrices dimensions don't match
		 * @exception InvalidMatrixException
		 *                if decomposed matrix is singular
		 */
		public RealMatrix solve(final RealMatrix b) throws IllegalArgumentException, InvalidMatrixException {

			if (b.getRowDimension() != singularValues.length) {
				throw MathRuntimeException.createIllegalArgumentException(
						"dimensions mismatch: got {0}x{1} but expected {2}x{3}", b.getRowDimension(),
						b.getColumnDimension(), singularValues.length, "n");
			}

			final RealMatrix w = uT.multiply(b);
			for (int i = 0; i < singularValues.length; ++i) {
				final double si = singularValues[i];
				if (si == 0) {
					throw new SingularMatrixException();
				}
				final double inv = 1.0 / si;
				for (int j = 0; j < b.getColumnDimension(); ++j) {
					w.multiplyEntry(i, j, inv);
				}
			}
			return v.multiply(w);

		}

		/**
		 * Check if the decomposed matrix is non-singular.
		 * 
		 * @return true if the decomposed matrix is non-singular
		 */
		public boolean isNonSingular() {
			return nonSingular;
		}

		/**
		 * Get the pseudo-inverse of the decomposed matrix.
		 * 
		 * @return inverse matrix
		 * @throws InvalidMatrixException
		 *             if decomposed matrix is singular
		 */
		public RealMatrix getInverse() throws InvalidMatrixException {

			if (!isNonSingular()) {
				throw new SingularMatrixException();
			}

			return solve(MatrixUtils.createRealIdentityMatrix(singularValues.length));

		}

	}

}
