/*
    * Copyright (C) 2012 Alberto Irurueta Carro (alberto@irurueta.com)
    *
    * Licensed 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 com.irurueta.algebra;
  
  /**
   * This class allows decomposition of matrices using Cholesky decomposition,
   * which consists on retrieving a lower or upper triangular matrix so that input
   * matrix can be decomposed as: A = L * L' = R' * R, where A is provided input
   * matrix, L is a lower triangular matrix and R is an upper triangular matrix.
   * Note: Cholesky decomposition can only be correctly computed on positive
   * definite matrices.
   */
  public class CholeskyDecomposer extends Decomposer {
  
      /**
       * Internal storage of Cholesky decomposition for provided input matrix.
       */
      private Matrix r;
  
      /**
       * Boolean indicating whether provided input matrix is symmetric and
       * positive definite.
       */
      private boolean spd;
  
      /**
       * Constructor of this class.
       */
      public CholeskyDecomposer() {
          super();
          r = null;
          spd = false;
      }
  
      /**
       * Constructor of this class.
       *
       * @param inputMatrix Reference to input matrix to be decomposed.
       */
      public CholeskyDecomposer(final Matrix inputMatrix) {
          super(inputMatrix);
          r = null;
          spd = false;
      }
  
      /**
       * Returns decomposer type corresponding to Cholesky decomposition.
       *
       * @return Decomposer type.
       */
      @Override
      public DecomposerType getDecomposerType() {
          return DecomposerType.CHOLESKY_DECOMPOSITION;
      }
  
      /**
       * Sets reference to input matrix to be decomposed.
       *
       * @param inputMatrix Reference to input matrix to be decomposed.
       * @throws LockedException Exception thrown if attempting to call this
       *                         method while this instance remains locked.
       */
      @Override
      public void setInputMatrix(final Matrix inputMatrix) throws LockedException {
          super.setInputMatrix(inputMatrix);
          r = null;
          spd = false;
      }
  
  
      /**
       * Returns boolean indicating whether decomposition has been computed and
       * results can be retrieved.
       * Attempting to retrieve decomposition results when not available, will
       * probably raise a NotAvailableException.
       *
       * @return Boolean indicating whether decomposition has been computed and
       * results can be retrieved.
       */
      @Override
      public boolean isDecompositionAvailable() {
          return r != null;
      }
  
      /**
       * This method computes Cholesky matrix decomposition, which consists on
      * factoring provided input matrix whenever it is square, symmetric and
      * positive definite into a lower triangulator factor such that it follows
      * next expression: A = L * L'
      * where A is input matrix and L is lower triangular factor (L' is its
      * transposed).
      * Cholesky decomposition can also be computed using Right Cholesky
      * decomposition, in which case A = R' * R, where R is an upper triangular
      * factor equal to L'.
      * Both factors L and R will be accessible once Cholesky decomposition has
      * been computed.
      * Note: During execution of this method, this instance will remain locked,
      * and hence attempting to set some parameters might raise a LockedException
      * Note: After execution of this method, Cholesky decomposition will be
      * available and operations such as retrieving L matrix factor or solving
      * systems of linear equations will be able to be done. Attempting to call
      * any of such operations before calling this method will raise a
      * NotAvailableException because they require computation of Cholesky
      * decomposition first.
      *
      * @throws NotReadyException   Exception thrown if attempting to call this
      *                             method when this instance is not ready (i.e. no input matrix has been
      *                             provided).
      * @throws LockedException     Exception thrown if this decomposer is already
      *                             locked before calling this method. Notice that this method will actually
      *                             lock this instance while it is being executed.
      * @throws DecomposerException Exception thrown if for any reason
      *                             decomposition fails while being executed, like when convergence of
      *                             results cannot be obtained, etc.
      */
     @Override
     public void decompose() throws NotReadyException, LockedException, DecomposerException {
 
         if (isLocked()) {
             throw new LockedException();
         }
 
         if (!isReady()) {
             throw new NotReadyException();
         }
 
         final var rows = inputMatrix.getRows();
         final var columns = inputMatrix.getColumns();
 
         if (rows != columns) {
             throw new DecomposerException();
         }
 
         locked = true;
 
         final Matrix localR;
         try {
             localR = new Matrix(columns, columns);
         } catch (WrongSizeException e) {
             throw new DecomposerException(e);
         }
 
         var localSpd = true;
         double d;
         double s;
 
         // Main loop
         for (var j = 0; j < columns; j++) {
             d = 0.0;
             for (var k = 0; k < j; k++) {
                 s = inputMatrix.getElementAt(k, j);
                 for (var i = 0; i < k; i++) {
                     s = s - localR.getElementAt(i, k) * localR.getElementAt(i, j);
                 }
                 s /= localR.getElementAt(k, k);
                 localR.setElementAt(k, j, s);
                 d += s * s;
                 localSpd &= inputMatrix.getElementAt(k, j) == inputMatrix.getElementAt(j, k);
             }
             d = inputMatrix.getElementAt(j, j) - d;
             localSpd &= d > 0.0;
             // sqrt of max(d, 0.0)
             localR.setElementAt(j, j, Math.sqrt(Math.max(d, 0.0)));
             for (var k = j + 1; k < columns; k++) {
                 localR.setElementAt(k, j, 0.0);
             }
         }
 
         this.spd = localSpd;
         this.r = localR;
 
         locked = false;
     }
 
     /**
      * Returns Cholesky matrix factor corresponding to a Lower triangular matrix
      * following this expression: A = L * L'. Where A is provided input matrix
      * that has been decomposed and L is the left lower triangular matrix
      * factor.
      *
      * @return Returns Cholesky Lower triangular matrix
      * @throws NotAvailableException Exception thrown if attempting to call this
      *                               method before actually computing Cholesky decomposition. To avoid this
      *                               exception call decompose() method first.
      * @see #decompose()
      */
     public Matrix getL() throws NotAvailableException {
         if (!isDecompositionAvailable()) {
             throw new NotAvailableException();
         }
 
         return r.transposeAndReturnNew();
     }
 
     /**
      * Returns Cholesky matrix factor corresponding to an upper triangular
      * matrix following this expression: A = R' * R. Where A is provided input
      * matrix that has been decomposed and R is the right upper triangular
      * matrix factor.
      *
      * @return Returns Cholesky upper triangular matrix
      * @throws NotAvailableException Exception thrown if attempting to call this
      *                               method before actually computing Cholesky decomposition. To avoid this
      *                               exception call decompose() method first.
      * @see #decompose()
      */
     public Matrix getR() throws NotAvailableException {
         if (!isDecompositionAvailable()) {
             throw new NotAvailableException();
         }
 
         return r;
     }
 
     /**
      * Returns boolean indicating whether provided input matrix is
      * Symmetric Positive Definite or not.
      * Notice that if returned value is false, then Cholesky decomposition
      * should be ignored, as Cholesky decomposition can only be computed on
      * symmetric positive definite matrices.
      *
      * @return Boolean indicating whether provided input matrix is symmetric
      * positive definite or not.
      * @throws NotAvailableException Exception thrown if attempting to call this
      *                               method before computing Cholesky decomposition. To avoid this exception
      *                               call decompose() method first.
      * @see #decompose()
      */
     public boolean isSPD() throws NotAvailableException {
         if (!isDecompositionAvailable()) {
             throw new NotAvailableException();
         }
 
         return spd;
     }
 
     /**
      * Solves a linear system of equations of the following form: A * X = B.
      * Where A is the input matrix provided for Cholesky decomposition, X is the
      * solution to the system of equations, and B is the parameters
      * vector/matrix.
      * Note: This method can be reused for different b vectors/matrices without
      * having to recompute Cholesky decomposition on the same input matrix.
      * Note: Provided b matrix must have the same number of rows as provided
      * input matrix A, otherwise an IllegalArgumentException will be raised
      * Note: Provided input matrix A must be square, otherwise a
      * WrongSizeException will be raised.
      * Note: If provided input matrix A is not symmetric positive definite, a
      * NonSymmetricPositiveDefiniteMatrixException will be thrown.
      * Note: In order to be able to execute this method, a Cholesky
      * decomposition must be available, otherwise a NotAvailableException will
      * be raised. In order to avoid this exception call decompose() method first
      * Note: result matrix contains solution of linear system of equations. It
      * will be resized if provided matrix does not have proper size
      *
      * @param b      Parameters of linear system of equations
      * @param result instance where solution X will be stored.
      * @throws com.irurueta.algebra.NotAvailableException                       if decomposition has
      *                                                                          not yet been computed.
      * @throws com.irurueta.algebra.WrongSizeException                          if the number of rows of
      *                                                                          b matrix is not equal to the number of rows of input matrix provided to
      *                                                                          Cholesky decomposer.
      * @throws com.irurueta.algebra.NonSymmetricPositiveDefiniteMatrixException if input matrix provided to Cholesky decomposer is not positive definite.
      */
     public void solve(final Matrix b, final Matrix result) throws NotAvailableException, WrongSizeException,
             NonSymmetricPositiveDefiniteMatrixException {
 
         if (!isDecompositionAvailable()) {
             throw new NotAvailableException();
         }
 
         final var rows = inputMatrix.getRows();
         final var columns = inputMatrix.getColumns();
         final var rowsB = b.getRows();
         final var colsB = b.getColumns();
 
         if (rowsB != rows) {
             throw new WrongSizeException();
         }
 
         if (!isSPD()) {
             throw new NonSymmetricPositiveDefiniteMatrixException();
         }
 
         // resize result matrix if needed
         if (result.getRows() != rowsB || result.getColumns() != colsB) {
             result.resize(rowsB, colsB);
         }
 
         // Copy b into result matrix
         result.copyFrom(b);
 
         final var l = getL();
 
         // Solve L * Y = B
         for (var k = 0; k < columns; k++) {
             for (var j = 0; j < colsB; j++) {
                 for (var i = 0; i < k; i++) {
                     result.setElementAt(k, j, result.getElementAt(k, j)
                             - result.getElementAt(i, j) * l.getElementAt(k, i));
                 }
                 result.setElementAt(k, j, result.getElementAt(k, j) / l.getElementAt(k, k));
             }
         }
 
         // Solve L' * X = Y
         int k2;
         for (var k = columns - 1; k >= 0; k--) {
             k2 = k;
             for (var j = 0; j < colsB; j++) {
                 for (var i = k2 + 1; i < columns; i++) {
                     result.setElementAt(k2, j, result.getElementAt(k2, j)
                             - result.getElementAt(i, j) * l.getElementAt(i, k2));
                 }
                 result.setElementAt(k2, j, result.getElementAt(k2, j) / l.getElementAt(k2, k2));
             }
         }
     }
 
     /**
      * Solves a linear system of equations of the following form: A * X = B.
      * Where A is the input matrix provided for Cholesky decomposition, X is the
      * solution to the system of equations, and B is the parameters
      * vector/matrix.
      * Note: This method can be reused for different b vectors/matrices without
      * having to recompute Cholesky decomposition on the same input matrix.
      * Note: Provided b matrix must have the same number of rows as provided
      * input matrix A, otherwise an IllegalArgumentException will be raised
      * Note: Provided input matrix A must be square, otherwise a
      * WrongSizeException will be raised.
      * Note: If provided input matrix A is not symmetric positive definite, a
      * NonSymmetricPositiveDefiniteMatrixException will be thrown.
      * Note: In order to be able to execute this method, a Cholesky
      * decomposition must be available, otherwise a NotAvailableException will
      * be raised. In order to avoid this exception call decompose() method first
      *
      * @param b Parameters of linear system of equations
      * @return a new matrix containing solution X.
      * @throws com.irurueta.algebra.NotAvailableException                       if decomposition has
      *                                                                          not yet been computed.
      * @throws com.irurueta.algebra.WrongSizeException                          if the number of rows of
      *                                                                          b matrix is not equal to the number of rows of input matrix provided to
      *                                                                          Cholesky decomposer.
      * @throws com.irurueta.algebra.NonSymmetricPositiveDefiniteMatrixException if input matrix provided to Cholesky decomposer is not positive definite.
      */
     public Matrix solve(final Matrix b) throws NotAvailableException, WrongSizeException,
             NonSymmetricPositiveDefiniteMatrixException {
 
         final var columns = inputMatrix.getColumns();
         final var colsB = b.getColumns();
         final var out = new Matrix(columns, colsB);
         solve(b, out);
         return out;
     }
 }