<commit_before>/*
 Copyright 2016 @author claydonkey (<NAME>)

 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 breeze.linalg

import org.netlib.util.intW
import com.github.fommil.netlib.LAPACK.{ getInstance => lapack }
import breeze.generic.UFunc
import breeze.linalg.hessenberg._
import breeze.linalg.householder._
import breeze.linalg.jacobi._
import breeze.math._
import breeze.numerics._
import scala.util.control.Breaks._
import breeze.linalg.Helper._
import DenseMatrix.canMapValues

/* returns tuple
      (y, z, Q, wR, wI, iLO2, iHI2) where
        y is  upper quasi-triangular matrix T from the Schur decomposition
      z  contains Q*Z. where
*/
object revertSchur extends UFunc {

  implicit object DMC_IMPL_DMC extends Impl2[DenseMatrix[Complex], DenseMatrix[Complex], DenseMatrix[Complex]] {
    def apply(M: DenseMatrix[Complex], T: DenseMatrix[Complex]): DenseMatrix[Complex] = {
      return T * upperTriangular(M) *  conj(T.t)
    }
  }
}

object schur extends UFunc {

  /* Schur decomposition  -computed by first reducing the
 * matrix to Hessenberg form using the class
 * HessenbergDecomposition. The Hessenberg matrix is then reduced
 * to triangular form by performing QR iterations with a single
 * shift. The cost of computing the Schur decomposition depends
 * on the number of iterations; as a rough guide, it may be taken
 * on the number of iterations; as a rough guide, it may be taken
 * to be \f$25n^3\f$ complex flops, or \f$10n^3\f$ complex flops
 * if \a computeU is false.*/

  /*
     *  Hessenberg decomposition of given matrix.
     *
     * The Hessenberg decomposition is computed by bringing the columns of the
     * matrix successively in the required form using Householder reflections
     * (see, e.g., Algorithm 7.4.2 in Golub \& Van Loan, <i>%Matrix
     * Computations</i>).
     */

  //return signature (Q,T,Householder tau, Householder  H)
  implicit object DMD_IMPL_SD extends Impl[DenseMatrix[Double], (DenseMatrix[Complex], DenseMatrix[Complex], DenseVector[Complex], DenseMatrix[Complex])] {
    def apply(M: DenseMatrix[Double]): (DenseMatrix[Complex], DenseMatrix[Complex], DenseVector[Complex], DenseMatrix[Complex]) = {

      new Schur[Double](M) {

        override val hess = hessenberg(M)
        val (s, z, y2, wR, wI, sLO, sHI) = getSchur(M) //LAPACK
        lazy override val T = DenseMatrix.tabulate[Complex](M.cols, M.rows)((i, j) => Complex(s(i, j), 0.0))
        lazy override val Q = hess._1 * DenseMatrix.tabulate[Complex](M.cols, M.rows)((i, j) => Complex(z(i, j), 0.0))
      }.decompose()
    }
  }

  implicit object DMD_IMPL_SI extends Impl[DenseMatrix[Int], (DenseMatrix[Complex], DenseMatrix[Complex], DenseVector[Complex], DenseMatrix[Complex])] {
    def apply(M: DenseMatrix[Int]): (DenseMatrix[Complex], DenseMatrix[Complex], DenseVector[Complex], DenseMatrix[Complex]) = {

      val MD = M.mapValues(_.toDouble)
      new Schur[Double](MD) {

        override val hess = hessenberg(MD)
        val (s, z, y2, wR, wI, sLO, sHI) = getSchur(MD) //LAPACK
        lazy override val T = DenseMatrix.tabulate[Complex](M.cols, M.rows)((i, j) => Complex(s(i, j), 0.0))
        lazy override val Q = hess._1 * DenseMatrix.tabulate[Complex](M.cols, M.rows)((i, j) => Complex(z(i, j), 0.0))
      }.decompose()
    }
  }

  implicit object DMCD_IMPL_SC extends Impl[DenseMatrix[Complex], (DenseMatrix[Complex], DenseMatrix[Complex], DenseVector[Complex], DenseMatrix[Complex])] {
    def apply(M: DenseMatrix[Complex]): (DenseMatrix[Complex], DenseMatrix[Complex], DenseVector[Complex], DenseMatrix[Complex]) = {

      new Schur[Complex](M) {

        override val hess = hessenberg(M)
        val m_maxIterationsPerRow = 30
        val maxIters = m_maxIterationsPerRow * hess._1.rows

        reduceToTriangularForm()
        /**
         * If matT(i+1,i) is neglegible in floating point arithmetic
         * compared to matT(i,i) and matT(j,j), then set it to zero and
         * return true, else return false.
         */
    private    def subdiagonalEntryIsNeglegible(i: Int) =
          {
            val d = norm1(T(i, i)) + norm1(T(i + 1, i + 1))
            val sd = norm1(T(i + 1, i))
            if (isMuchSmallerThan(sd, d)) {
              T(i + 1, i) = Complex(0.0, 0.0)
              true
            } else
              false
          }

        /** Compute the shift in the current QR iteration. */
      private  def computeShift(iu: Int, iter: Int) = {
          if (iter == 10 || iter == 20) {
            // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
            abs(T(iu, iu - 1).real) + abs(T(iu - 1, iu - 2).real)
          }
          // compute the shift as one of the eigenvalues of t, the 2x2
          // diagonal block on the bottom of the active submatrix

          var t = T((iu - 1) to iu, (iu - 1) to iu) // Complex(NormT, 0)
          val normt = Complex(sum(t.mapValues(abs(_))), 0)

          t = t / normt
          val b = t(0, 1) * t(1, 0)
          val c = t(0, 0) - t(1, 1)
          val disc = breeze.numerics.pow((c * c + 4.0 * b), 0.5)
          val det = (t(0, 0) * t(1, 1)) - b
          val trace = t(0, 0) + t(1, 1)
          var eival1 = (trace + disc) / 2.0
          var eival2 = (trace - disc) / 2.0

          if (norm1(eival1) > norm1(eival2))
            eival2 = det / eival1
          else
            eival1 = det / eival2
          // choose the eigenvalue closest to the bottom entry of the diagonal
          if (norm1(eival1 - t(1, 1)) < norm1(eival2 - t(1, 1)))
            normt * eival1
          else
            normt * eival2

        }
        /**A horrbly disFunctional implementation using breaks, iterations  and mutations BLEUURRGGH will change... **/
        // The matrix matT is divided in three parts.
        // Rows 0,...,il-1 are decoupled from the rest because matT(il,il-1) is zero.
        // Rows il,...,iu is the part we are working on (the active submatrix).
        // Rows iu+1,...,end are already brought in triangular form.
   private     def reduceToTriangularForm() = {

          var matnum = 0
          var matnum2 = 0
          var newrot = false
          var iu = T.cols - 1
          val maxIters = m_maxIterationsPerRow * T.rows
          var il = 0
          var iter = 0 // number of iterations we are working on the (iu,iu) element
          var totalIter = 0 // number of iterations for whole matrix

          breakable {
            while (true) {
              // find iu, the bottom row of the active submatrix
              breakable {
                while (iu > 0) {
                  if (subdiagonalEntryIsNeglegible(iu - 1)) {
                  }
                  if (!subdiagonalEntryIsNeglegible(iu - 1)) break
                  iter = 0
                  iu = iu - 1
                }

              }
              // if iu is zero then we are done; the whole matrix is triangularized

              if (iu == 0) break

              // if we spent too many iterations, we give up
              iter = iter + 1
              totalIter = totalIter + 1

              if (totalIter > maxIters) break

              // find il, the top row of the active submatrix
              il = iu - 1
              while (il > 0 && !subdiagonalEntryIsNeglegible(il - 1)) {
                il = il - 1
              }
              /* perform the QR step using Givens rotations. The first rotation
	   *creates a bulge; the (il+2,il) element becomes nonzero. This
	   *bulge is chased down to the bottom of the active submatrix.
	   */
              val givrot1 = makeGivens(T(il, il) - computeShift(iu, iter), T(il + 1, il))
              jacobi(T(il to il + 1, ::)).rotateL(givrot1)
              jacobi(T(0 to (min(il + 2, iu)), il to il + 1)) rotateR (givrot1)
              jacobi(Q(::, il to il + 1)) rotateR (givrot1)

              val idx: Int = 0

              for (idx <- ((il + 1) to iu - 1)) {
                val givrot2 = makeGivens(T(idx, idx - 1), T(idx + 1, idx - 1))
                T(idx, idx - 1) = givrot2.rot
                T(idx + 1, idx - 1) = Complex(0.0, 0.0)
                jacobi(T(idx to idx + 1, ::)) rotateL (givrot2)
                jacobi(T(0 to (min(idx + 2, iu)), idx to idx + 1)) rotateR (givrot2)
                jacobi(Q(::, idx to idx + 1)) rotateR (givrot2)
              }
            }
          }
        }
      }.decompose()
    }
  }
  abstract class Schur[T](val M: DenseMatrix[T]) {
    if (M.rows != M.cols)
      throw new MatrixNotSquareException

    val hess: (DenseMatrix[Complex], DenseMatrix[Complex], Householder)

    lazy val Q = hess._1
    lazy val T = hess._2


    def decompose() = (T, Q, hess._3.tau, hess._3.matrixH)
  }
  //tau =  the householder coeffs
  object getSchur extends Impl[DenseMatrix[Double], (DenseMatrix[Double], DenseMatrix[Double], DenseMatrix[Double], Array[Double], Array[Double], Int, Int)] {
    def apply(X: DenseMatrix[Double]): (DenseMatrix[Double], DenseMatrix[Double], DenseMatrix[Double], Array[Double], Array[Double], Int, Int) = {<commit_msg>error in Schur :: not eq to (0 to cols -1)<commit_after>/*
 Copyright 2016 @author claydonkey (<NAME>)

 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 breeze.linalg

import org.netlib.util.intW
import com.github.fommil.netlib.LAPACK.{ getInstance => lapack }
import breeze.generic.UFunc
import breeze.linalg.hessenberg._
import breeze.linalg.householder._
import breeze.linalg.jacobi._
import breeze.math._
import breeze.numerics._
import scala.util.control.Breaks._
import breeze.linalg.Helper._
import DenseMatrix.canMapValues

/* returns tuple
      (y, z, Q, wR, wI, iLO2, iHI2) where
        y is  upper quasi-triangular matrix T from the Schur decomposition
      z  contains Q*Z. where
*/
object revertSchur extends UFunc {

  implicit object DMC_IMPL_DMC extends Impl2[DenseMatrix[Complex], DenseMatrix[Complex], DenseMatrix[Complex]] {
    def apply(M: DenseMatrix[Complex], T: DenseMatrix[Complex]): DenseMatrix[Complex] = {
      return T * upperTriangular(M) *  conj(T.t)
    }
  }
}

object schur extends UFunc {

  /* Schur decomposition  -computed by first reducing the
 * matrix to Hessenberg form using the class
 * HessenbergDecomposition. The Hessenberg matrix is then reduced
 * to triangular form by performing QR iterations with a single
 * shift. The cost of computing the Schur decomposition depends
 * on the number of iterations; as a rough guide, it may be taken
 * on the number of iterations; as a rough guide, it may be taken
 * to be \f$25n^3\f$ complex flops, or \f$10n^3\f$ complex flops
 * if \a computeU is false.*/

  /*
     *  Hessenberg decomposition of given matrix.
     *
     * The Hessenberg decomposition is computed by bringing the columns of the
     * matrix successively in the required form using Householder reflections
     * (see, e.g., Algorithm 7.4.2 in Golub \& Van Loan, <i>%Matrix
     * Computations</i>).
     */

  //return signature (Q,T,Householder tau, Householder  H)
  implicit object DMD_IMPL_SD extends Impl[DenseMatrix[Double], (DenseMatrix[Complex], DenseMatrix[Complex], DenseVector[Complex], DenseMatrix[Complex])] {
    def apply(M: DenseMatrix[Double]): (DenseMatrix[Complex], DenseMatrix[Complex], DenseVector[Complex], DenseMatrix[Complex]) = {

      new Schur[Double](M) {

        override val hess = hessenberg(M)
        val (s, z, y2, wR, wI, sLO, sHI) = getSchur(M) //LAPACK
        lazy override val T = DenseMatrix.tabulate[Complex](M.cols, M.rows)((i, j) => Complex(s(i, j), 0.0))
        lazy override val Q = hess._1 * DenseMatrix.tabulate[Complex](M.cols, M.rows)((i, j) => Complex(z(i, j), 0.0))
      }.decompose()
    }
  }

  implicit object DMD_IMPL_SI extends Impl[DenseMatrix[Int], (DenseMatrix[Complex], DenseMatrix[Complex], DenseVector[Complex], DenseMatrix[Complex])] {
    def apply(M: DenseMatrix[Int]): (DenseMatrix[Complex], DenseMatrix[Complex], DenseVector[Complex], DenseMatrix[Complex]) = {

      val MD = M.mapValues(_.toDouble)
      new Schur[Double](MD) {

        override val hess = hessenberg(MD)
        val (s, z, y2, wR, wI, sLO, sHI) = getSchur(MD) //LAPACK
        lazy override val T = DenseMatrix.tabulate[Complex](M.cols, M.rows)((i, j) => Complex(s(i, j), 0.0))
        lazy override val Q = hess._1 * DenseMatrix.tabulate[Complex](M.cols, M.rows)((i, j) => Complex(z(i, j), 0.0))
      }.decompose()
    }
  }

  implicit object DMCD_IMPL_SC extends Impl[DenseMatrix[Complex], (DenseMatrix[Complex], DenseMatrix[Complex], DenseVector[Complex], DenseMatrix[Complex])] {
    def apply(M: DenseMatrix[Complex]): (DenseMatrix[Complex], DenseMatrix[Complex], DenseVector[Complex], DenseMatrix[Complex]) = {

      new Schur[Complex](M) {

        override val hess = hessenberg(M)
        val m_maxIterationsPerRow = 30
        val maxIters = m_maxIterationsPerRow * hess._1.rows

        reduceToTriangularForm()
        /**
         * If matT(i+1,i) is neglegible in floating point arithmetic
         * compared to matT(i,i) and matT(j,j), then set it to zero and
         * return true, else return false.
         */
    private    def subdiagonalEntryIsNeglegible(i: Int) =
          {
            val d = norm1(T(i, i)) + norm1(T(i + 1, i + 1))
            val sd = norm1(T(i + 1, i))
            if (isMuchSmallerThan(sd, d)) {
              T(i + 1, i) = Complex(0.0, 0.0)
              true
            } else
              false
          }

        /** Compute the shift in the current QR iteration. */
      private  def computeShift(iu: Int, iter: Int) = {
          if (iter == 10 || iter == 20) {
            // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
            abs(T(iu, iu - 1).real) + abs(T(iu - 1, iu - 2).real)
          }
          // compute the shift as one of the eigenvalues of t, the 2x2
          // diagonal block on the bottom of the active submatrix

          var t = T((iu - 1) to iu, (iu - 1) to iu) // Complex(NormT, 0)
          val normt = Complex(sum(t.mapValues(abs(_))), 0)

          t = t / normt
          val b = t(0, 1) * t(1, 0)
          val c = t(0, 0) - t(1, 1)
          val disc = breeze.numerics.pow((c * c + 4.0 * b), 0.5)
          val det = (t(0, 0) * t(1, 1)) - b
          val trace = t(0, 0) + t(1, 1)
          var eival1 = (trace + disc) / 2.0
          var eival2 = (trace - disc) / 2.0

          if (norm1(eival1) > norm1(eival2))
            eival2 = det / eival1
          else
            eival1 = det / eival2
          // choose the eigenvalue closest to the bottom entry of the diagonal
          if (norm1(eival1 - t(1, 1)) < norm1(eival2 - t(1, 1)))
            normt * eival1
          else
            normt * eival2

        }
        /**A horrbly disFunctional implementation using breaks, iterations  and mutations BLEUURRGGH will change... **/
        // The matrix matT is divided in three parts.
        // Rows 0,...,il-1 are decoupled from the rest because matT(il,il-1) is zero.
        // Rows il,...,iu is the part we are working on (the active submatrix).
        // Rows iu+1,...,end are already brought in triangular form.
   private     def reduceToTriangularForm() = {

          var matnum = 0
          var matnum2 = 0
          var newrot = false
          var iu = T.cols - 1
          val maxIters = m_maxIterationsPerRow * T.rows
          var il = 0
          var iter = 0 // number of iterations we are working on the (iu,iu) element
          var totalIter = 0 // number of iterations for whole matrix

          breakable {
            while (true) {
              // find iu, the bottom row of the active submatrix
              breakable {
                while (iu > 0) {
                  if (subdiagonalEntryIsNeglegible(iu - 1)) {
                  }
                  if (!subdiagonalEntryIsNeglegible(iu - 1)) break
                  iter = 0
                  iu = iu - 1
                }

              }
              // if iu is zero then we are done; the whole matrix is triangularized

              if (iu == 0) break

              // if we spent too many iterations, we give up
              iter = iter + 1
              totalIter = totalIter + 1

              if (totalIter > maxIters) break

              // find il, the top row of the active submatrix
              il = iu - 1
              while (il > 0 && !subdiagonalEntryIsNeglegible(il - 1)) {
                il = il - 1
              }
              /* perform the QR step using Givens rotations. The first rotation
	   *creates a bulge; the (il+2,il) element becomes nonzero. This
	   *bulge is chased down to the bottom of the active submatrix.
	   */
              val givrot1 = makeGivens(T(il, il) - computeShift(iu, iter), T(il + 1, il))
              jacobi(T(il to il + 1, ::)).rotateL(givrot1)
              jacobi(T(0 to (min(il + 2, iu)), il to il + 1)) rotateR (givrot1)
              jacobi(Q(::, il to il + 1)) rotateR (givrot1)

              val idx: Int = 0

              for (idx <- ((il + 1) to iu - 1)) {
                val givrot2 = makeGivens(T(idx, idx - 1), T(idx + 1, idx - 1))
                T(idx, idx - 1) = givrot2.rot
                T(idx + 1, idx - 1) = Complex(0.0, 0.0)
                jacobi(T(idx to idx + 1, idx to T.cols - 1)) rotateL (givrot2)
                jacobi(T(0 to (min(idx + 2, iu)), idx to idx + 1)) rotateR (givrot2)
                jacobi(Q(::, idx to idx + 1)) rotateR (givrot2)
              }
            }
          }
        }
      }.decompose()
    }
  }
  abstract class Schur[T](val M: DenseMatrix[T]) {
    if (M.rows != M.cols)
      throw new MatrixNotSquareException

    val hess: (DenseMatrix[Complex], DenseMatrix[Complex], Householder)

    lazy val Q = hess._1
    lazy val T = hess._2


    def decompose() = (T, Q, hess._3.tau, hess._3.matrixH)
  }
  //tau =  the householder coeffs
  object getSchur extends Impl[DenseMatrix[Double], (DenseMatrix[Double], DenseMatrix[Double], DenseMatrix[Double], Array[Double], Array[Double], Int, Int)] {
    def apply(X: DenseMatrix[Double]): (DenseMatrix[Double], DenseMatrix[Double], DenseMatrix[Double], Array[Double], Array[Double], Int, Int) = {
{'fraction_non_alphanumeric': 0.11493832322809952, 'fraction_numerical': 0.017405394104118753, 'mean_word_length': 2.229743416610398, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6537260', 'n_tokens_mistral': 6391, 'n_tokens_neox': 6004, 'n_words': 1955}