text: import Matrix from '../matrix';
  import WrapperMatrix2D from '../wrap/WrapperMatrix2D';
  
  import { hypotenuse } from './util';
  
  /**
   * @class QrDecomposition
   * @link https://github.com/lutzroeder/Mapack/blob/master/Source/QrDecomposition.cs
   * @param {Matrix} value
   */
  export default class QrDecomposition {
    constructor(value) {
      value = WrapperMatrix2D.checkMatrix(value);
  
      var qr = value.clone();
      var m = value.rows;
      var n = value.columns;
      var rdiag = new Array(n);
      var i, j, k, s;
  
      for (k = 0; k < n; k++) {
        var nrm = 0;
        for (i = k; i < m; i++) {
          nrm = hypotenuse(nrm, qr.get(i, k));
        }
        if (nrm !== 0) {
          if (qr.get(k, k) < 0) {
            nrm = -nrm;
          }
          for (i = k; i < m; i++) {
            qr.set(i, k, qr.get(i, k) / nrm);
          }
          qr.set(k, k, qr.get(k, k) + 1);
          for (j = k + 1; j < n; j++) {
            s = 0;
            for (i = k; i < m; i++) {
              s += qr.get(i, k) * qr.get(i, j);
            }
            s = -s / qr.get(k, k);
            for (i = k; i < m; i++) {
              qr.set(i, j, qr.get(i, j) + s * qr.get(i, k));
            }
          }
        }
        rdiag[k] = -nrm;
      }
  
      this.QR = qr;
      this.Rdiag = rdiag;
    }
  
    /**
     * Solve a problem of least square (Ax=b) by using the QR decomposition. Useful when A is rectangular, but not working when A is singular.
     * Example : We search to approximate x, with A matrix shape m*n, x vector size n, b vector size m (m > n). We will use :
     * var qr = QrDecomposition(A);
     * var x = qr.solve(b);
     * @param {Matrix} value - Matrix 1D which is the vector b (in the equation Ax = b)
     * @return {Matrix} - The vector x
     */
    solve(value) {
      value = Matrix.checkMatrix(value);
  
      var qr = this.QR;
      var m = qr.rows;
  
      if (value.rows !== m) {
        throw new Error('Matrix row dimensions must agree');
      }
      if (!this.isFullRank()) {
        throw new Error('Matrix is rank deficient');
      }
  
      var count = value.columns;
      var X = value.clone();
      var n = qr.columns;
      var i, j, k, s;
  
      for (k = 0; k < n; k++) {
        for (j = 0; j < count; j++) {
          s = 0;
          for (i = k; i < m; i++) {
            s += qr.get(i, k) * X.get(i, j);
          }
          s = -s / qr.get(k, k);
          for (i = k; i < m; i++) {
            X.set(i, j, X.get(i, j) + s * qr.get(i, k));
          }
        }
      }
      for (k = n - 1; k >= 0; k--) {
        for (j = 0; j < count; j++) {
          X.set(k, j, X.get(k, j) / this.Rdiag[k]);
        }
        for (i = 0; i < k; i++) {
          for (j = 0; j < count; j++) {
            X.set(i, j, X.get(i, j) - X.get(k, j) * qr.get(i, k));
          }
        }
      }
  
      return X.subMatrix(0, n - 1, 0, count - 1);
    }
  
    /**
     *
     * @return {boolean}
     */
    isFullRank() {
      var columns = this.QR.columns;
      for (var i = 0; i < columns; i++) {
        if (this.Rdiag[i] === 0) {
          return false;
        }
      }
      return true;
    }
  
    /**
     *
     * @return {Matrix}
     */
    get upperTriangularMatrix() {
      var qr = this.QR;
      var n = qr.columns;
      var X = new Matrix(n, n);
      var i, j;
      for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
          if (i < j) {
            X.set(i, j, qr.get(i, j));
          } else if (i === j) {
            X.set(i, j, this.Rdiag[i]);
          } else {
            X.set(i, j, 0);
          }
        }
      }
      return X;
    }
  
    /**
     *
     * @return {Matrix}
     */
    get orthogonalMatrix() {
      var qr = this.QR;
      var rows = qr.rows;
      var columns = qr.columns;
      var X = new Matrix(rows, columns);
      var i, j, k, s;
  
      for (k = columns - 1; k >= 0; k--) {
        for (i = 0; i < rows; i++) {
          X.set(i, k, 0);
        }
        X.set(k, k, 1);
        for (j = k; j < columns; j++) {
          if (qr.get(k, k) !== 0) {
            s = 0;
            for (i = k; i < rows; i++) {
              s += qr.get(i, k) * X.get(i, j);
            }
  
            s = -s / qr[k][k];
  
            for (i = k; i < rows; i++) {
              X.set(i, j, X.get(i, j) + s * qr.get(i, k));
            }
          }
        }
      }
      return X;
    }
  }
  