/*
 * 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.math3.linear;

import org.apache.commons.math3.util.FastMath;


Calculates the rank-revealing QR-decomposition of a matrix, with column pivoting.
The rank-revealing QR-decomposition of a matrix A consists of three matrices Q,
R and P such that AP=QR.  Q is orthogonal (QTQ = I), and R is upper triangular.
If A is m×n, Q is m×m and R is m×n and P is n×n.
QR decomposition with column pivoting produces a rank-revealing QR decomposition and the getRank(double) method may be used to return the rank of the input matrix A.
This class compute the decomposition using Householder reflectors.
For efficiency purposes, the decomposition in packed form is transposed.
This allows inner loop to iterate inside rows, which is much more cache-efficient
in Java.
This class is based on the class with similar name from the
JAMA library, with the
following changes:

  
a getQT method has been added,
  
the solve and isFullRank methods have been replaced by a getSolver method and the equivalent methods provided by the returned DecompositionSolver.

See Also:
MathWorld
Wikipedia
Since:
3.2/**
 * Calculates the rank-revealing QR-decomposition of a matrix, with column pivoting.
 * <p>The rank-revealing QR-decomposition of a matrix A consists of three matrices Q,
 * R and P such that AP=QR.  Q is orthogonal (Q<sup>T</sup>Q = I), and R is upper triangular.
 * If A is m&times;n, Q is m&times;m and R is m&times;n and P is n&times;n.</p>
 * <p>QR decomposition with column pivoting produces a rank-revealing QR
 * decomposition and the {@link #getRank(double)} method may be used to return the rank of the
 * input matrix A.</p>
 * <p>This class compute the decomposition using Householder reflectors.</p>
 * <p>For efficiency purposes, the decomposition in packed form is transposed.
 * This allows inner loop to iterate inside rows, which is much more cache-efficient
 * in Java.</p>
 * <p>This class is based on the class with similar name from the
 * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library, with the
 * following changes:</p>
 * <ul>
 *   <li>a {@link #getQT() getQT} method has been added,</li>
 *   <li>the {@code solve} and {@code isFullRank} methods have been replaced
 *   by a {@link #getSolver() getSolver} method and the equivalent methods
 *   provided by the returned {@link DecompositionSolver}.</li>
 * </ul>
 *
 * @see <a href="http://mathworld.wolfram.com/QRDecomposition.html">MathWorld</a>
 * @see <a href="http://en.wikipedia.org/wiki/QR_decomposition">Wikipedia</a>
 *
 * @since 3.2
 */
public class RRQRDecomposition extends QRDecomposition {

    
An array to record the column pivoting for later creation of P. /** An array to record the column pivoting for later creation of P. */
    private int[] p;

    
Cached value of P. /** Cached value of P. */
    private RealMatrix cachedP;


    
Calculates the QR-decomposition of the given matrix.
The singularity threshold defaults to zero.
Params:
matrix – The matrix to decompose.
See Also:
RRQRDecomposition(RealMatrix, double)/**
     * Calculates the QR-decomposition of the given matrix.
     * The singularity threshold defaults to zero.
     *
     * @param matrix The matrix to decompose.
     *
     * @see #RRQRDecomposition(RealMatrix, double)
     */
    public RRQRDecomposition(RealMatrix matrix) {
        this(matrix, 0d);
    }

   
Calculates the QR-decomposition of the given matrix.
Params:
matrix – The matrix to decompose.
threshold – Singularity threshold.
See Also:
RRQRDecomposition(RealMatrix)/**
     * Calculates the QR-decomposition of the given matrix.
     *
     * @param matrix The matrix to decompose.
     * @param threshold Singularity threshold.
     * @see #RRQRDecomposition(RealMatrix)
     */
    public RRQRDecomposition(RealMatrix matrix,  double threshold) {
        super(matrix, threshold);
    }

    
Decompose matrix.
Params:
qrt – transposed matrix/** Decompose matrix.
     * @param qrt transposed matrix
     */
    @Override
    protected void decompose(double[][] qrt) {
        p = new int[qrt.length];
        for (int i = 0; i < p.length; i++) {
            p[i] = i;
        }
        super.decompose(qrt);
    }

    
Perform Householder reflection for a minor A(minor, minor) of A.
Params:
minor – minor index
qrt – transposed matrix/** Perform Householder reflection for a minor A(minor, minor) of A.
     * @param minor minor index
     * @param qrt transposed matrix
     */
    @Override
    protected void performHouseholderReflection(int minor, double[][] qrt) {

        double l2NormSquaredMax = 0;
        // Find the unreduced column with the greatest L2-Norm
        int l2NormSquaredMaxIndex = minor;
        for (int i = minor; i < qrt.length; i++) {
            double l2NormSquared = 0;
            for (int j = 0; j < qrt[i].length; j++) {
                l2NormSquared += qrt[i][j] * qrt[i][j];
            }
            if (l2NormSquared > l2NormSquaredMax) {
                l2NormSquaredMax = l2NormSquared;
                l2NormSquaredMaxIndex = i;
            }
        }
        // swap the current column with that with the greated L2-Norm and record in p
        if (l2NormSquaredMaxIndex != minor) {
            double[] tmp1 = qrt[minor];
            qrt[minor] = qrt[l2NormSquaredMaxIndex];
            qrt[l2NormSquaredMaxIndex] = tmp1;
            int tmp2 = p[minor];
            p[minor] = p[l2NormSquaredMaxIndex];
            p[l2NormSquaredMaxIndex] = tmp2;
        }

        super.performHouseholderReflection(minor, qrt);

    }


    
Returns the pivot matrix, P, used in the QR Decomposition of matrix A such that AP = QR.
If no pivoting is used in this decomposition then P is equal to the identity matrix.
Returns:
a permutation matrix./**
     * Returns the pivot matrix, P, used in the QR Decomposition of matrix A such that AP = QR.
     *
     * If no pivoting is used in this decomposition then P is equal to the identity matrix.
     *
     * @return a permutation matrix.
     */
    public RealMatrix getP() {
        if (cachedP == null) {
            int n = p.length;
            cachedP = MatrixUtils.createRealMatrix(n,n);
            for (int i = 0; i < n; i++) {
                cachedP.setEntry(p[i], i, 1);
            }
        }
        return cachedP ;
    }

    
Return the effective numerical matrix rank.
The effective numerical rank is the number of non-negligible
singular values.
This implementation looks at Frobenius norms of the sequence of
bottom right submatrices.  When a large fall in norm is seen,
the rank is returned. The drop is computed as:
  (thisNorm/lastNorm) * rNorm < dropThreshold


where thisNorm is the Frobenius norm of the current submatrix,
lastNorm is the Frobenius norm of the previous submatrix,
rNorm is is the Frobenius norm of the complete matrix

Params:
dropThreshold – threshold triggering rank computation
Returns:
effective numerical matrix rank/**
     * Return the effective numerical matrix rank.
     * <p>The effective numerical rank is the number of non-negligible
     * singular values.</p>
     * <p>This implementation looks at Frobenius norms of the sequence of
     * bottom right submatrices.  When a large fall in norm is seen,
     * the rank is returned. The drop is computed as:</p>
     * <pre>
     *   (thisNorm/lastNorm) * rNorm < dropThreshold
     * </pre>
     * <p>
     * where thisNorm is the Frobenius norm of the current submatrix,
     * lastNorm is the Frobenius norm of the previous submatrix,
     * rNorm is is the Frobenius norm of the complete matrix
     * </p>
     *
     * @param dropThreshold threshold triggering rank computation
     * @return effective numerical matrix rank
     */
    public int getRank(final double dropThreshold) {
        RealMatrix r    = getR();
        int rows        = r.getRowDimension();
        int columns     = r.getColumnDimension();
        int rank        = 1;
        double lastNorm = r.getFrobeniusNorm();
        double rNorm    = lastNorm;
        while (rank < FastMath.min(rows, columns)) {
            double thisNorm = r.getSubMatrix(rank, rows - 1, rank, columns - 1).getFrobeniusNorm();
            if (thisNorm == 0 || (thisNorm / lastNorm) * rNorm < dropThreshold) {
                break;
            }
            lastNorm = thisNorm;
            rank++;
        }
        return rank;
    }

    
Get a solver for finding the A × X = B solution in least square sense.
 Least Square sense means a solver can be computed for an overdetermined system, (i.e. a system with more equations than unknowns, which corresponds to a tall A matrix with more rows than columns). In any case, if the matrix is singular within the tolerance set at construction, an error will be triggered when the solve method will be called. 
Returns:
a solver/**
     * Get a solver for finding the A &times; X = B solution in least square sense.
     * <p>
     * Least Square sense means a solver can be computed for an overdetermined system,
     * (i.e. a system with more equations than unknowns, which corresponds to a tall A
     * matrix with more rows than columns). In any case, if the matrix is singular
     * within the tolerance set at {@link RRQRDecomposition#RRQRDecomposition(RealMatrix,
     * double) construction}, an error will be triggered when
     * the {@link DecompositionSolver#solve(RealVector) solve} method will be called.
     * </p>
     * @return a solver
     */
    @Override
    public DecompositionSolver getSolver() {
        return new Solver(super.getSolver(), this.getP());
    }

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

        
Upper level solver. /** Upper level solver. */
        private final DecompositionSolver upper;

        
A permutation matrix for the pivots used in the QR decomposition /** A permutation matrix for the pivots used in the QR decomposition */
        private RealMatrix p;

        
Build a solver from decomposed matrix.
Params:
upper – upper level solver.
p – permutation matrix/**
         * Build a solver from decomposed matrix.
         *
         * @param upper upper level solver.
         * @param p permutation matrix
         */
        private Solver(final DecompositionSolver upper, final RealMatrix p) {
            this.upper = upper;
            this.p     = p;
        }

        
{@inheritDoc} /** {@inheritDoc} */
        public boolean isNonSingular() {
            return upper.isNonSingular();
        }

        
{@inheritDoc} /** {@inheritDoc} */
        public RealVector solve(RealVector b) {
            return p.operate(upper.solve(b));
        }

        
{@inheritDoc} /** {@inheritDoc} */
        public RealMatrix solve(RealMatrix b) {
            return p.multiply(upper.solve(b));
        }

        
{@inheritDoc}
Throws:
SingularMatrixException – if the decomposed matrix is singular./**
         * {@inheritDoc}
         * @throws SingularMatrixException if the decomposed matrix is singular.
         */
        public RealMatrix getInverse() {
            return solve(MatrixUtils.createRealIdentityMatrix(p.getRowDimension()));
        }
    }
}