//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller
 *  @version 1.1
 *  @date    Thu Feb 14 14:06:00 EST 2013
 *  @see     LICENSE (MIT style license file).
 *
 *  @see www2.cs.cas.cz/mweb/download/publi/Ples2006.pdf
 *  @see www.math.iit.edu/~fass/477577_Chapter_12.pdf
 *  @see Handbook of Linear Algrbra, Chapter 45
 *  @see cs.fit.edu/~dmitra/SciComp/11Spr/SVD-Presentation-Updated2.ppt
 */

// U N D E R   D E V E L O P M E N T  - to be merged with SVDecomp.scala

package scalation.linalgebra

import math.abs

import scalation.linalgebra.Givens.{givens, givensColUpdate, givensRowUpdate}
import scalation.linalgebra.Householder.house
import scalation.linalgebra.MatrixD.{eye, outer}
import scalation.math.Basic.oneIf
import scalation.util.Error

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `SVD` class is used to compute the Singlar Value Decomposition (SVD) of
 *  matrix 'a' using the Golub-Kahan-Reinsch Algorithm.
 *  Decompose matrix 'a' into the product of three matrices:
 *  <p>
 *      a = u * b * v.t
 *  <p>
 *  where 'u' is a matrix of orthonormal eigenvectors of 'a * a.t',
 *        'v' is a matrix of orthonormal eigenvectors of 'a.t * a' and
 *        'b' is a diagonal matrix of singular values (sqrt of eigenvalues of a.t * a).
 *  @param a  the m-by-n matrix to decompose (requires m >= n)
 */
class SVD (a: MatrixD)
      extends Error
{
    private val EPSILON = 1E-12
    private val m = a.dim1         // number of rows
    private val n = a.dim2         // number of columns

    if (n > m) flaw ("constructor", "SVD requires m >= n")

    private val u = eye (m)        // initialize left unitary matrix to m-by-m identity matrix
    private val v = eye (n)        // initialize right unitary matrix to n-by-n identity matrix

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Use the Householder Bidiagonalization Algorithm to compute unitary
     *  matrices u and v such that u.t * a * v = b where matrix b is bidiagonal.
     *  The b matrix will only have non-zero elements on its main diagonal and
     *  its super diagonal (the diagonal above the main).
     *  This implementation computes b in-place in matrix a.
     *  @see Matrix Computations:  Algorithm 5.4.2 Householder Bidiagonalization
     */
    def bidiagonalize (): Tuple3 [MatrixD, MatrixD, MatrixD] =
    {
        for (j <- 0 until n) {
            val (vu, bu) = house (a(j until m, j))                   // Householder (vector vu, scalar bu)
            val hu = eye (m-j) - outer (vu, vu * bu)                 // Householder matrix hu
            a(j until m, j until n) = hu * a(j until m, j until n)   // zero column j below main diagonal
//          for (i <- j+1 until m) a(i, j) = vu(i-j)                 // save key part of Householder vector vu
            u *= hu.diag (j)                                         // multiply u by Householder matrix hu

            if (j < n-2) {
                val (vv, bv) = house (a(j, j+1 until n))                     // Householder (vector vv, scalar bv)
                val hv = eye (n-j-1) - outer (vv, vv * bv)                   // Householder matrix hv
                a(j until m, j+1 until n) = a(j until m, j+1 until n) * hv   // zero row j past super diagonal
//              for (k <- j+2 until n) a(j, k) = vv(k-j-1)                   // save key part of Householder vector vv
                v *= hv.diag (j+1)                                           // multiply v by Householder matrix hv
            } // if
        } // for

        for (i <- 0 until m; j <- 0 until n if j != i && j != i+1) a(i, j) = 0.0  // clean matrix off diagonals

        (u, a, v)                           // return (u, b as a, v)
    } // bidiagonalize

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Use the Golub-Kahan Bidiagonalization Algorithm to compute unitary matrices
     *  u and v such that u.t * a * v = b where matrix b is bidiagonal (nonzero
     *  elements on the main diagonal and the diagonal above it).  Solve a * v = u * b
     *  by computing column vectors u_k and v_k and diagonals c and d.  Use transposes
     *  of u and v since row access is more efficient than column access.
     *  Caveat: assumes bidiagonals elements are non-negative and need to add
     *  re-orthogonalization steps.
     *  @see web.eecs.utk.edu/~dongarra/etemplates/node198.html
     */
    def bidiagonalization2 (): Tuple2 [VectorD, VectorD] =
    {
        var c  = new VectorD (n-1)                   // above main diagonal (beta)
        var d  = new VectorD (n)                     // main diagonal (alpha)
        var ut = new MatrixD (n, m)                  // transpose of u 
        var vt = new MatrixD (n, m); vt(0, 0) = 1.0  // transpose of v
        val at = a.t                                 // transpose of matrix a

        for (k <- 0 until n) {
            ut(k)  = a * vt(k); if (k > 0) ut(k) -= ut(k-1) * c(k-1)
            d(k)   = ut(k).norm
            ut(k) /= d(k)

            val k1 = k+1
            if (k1 < n) {
                vt(k1)  = at * ut(k) - vt(k) * d(k)
                c(k)    = vt(k1).norm
                vt(k1) /= c(k)
            } // if
        } // for

        println ("ut = " + ut)
        println ("vt = " + vt)

        (c, d)
    } // bidiagonalization2

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Take one step in converting the bidiagoanl matrix 'b' to a diagonal matrix.
     *  @see Matrix Computation: Algorithm 8.6.1 Golub-Kahan Step.
     *  @param b  the bidiagonal matrix to diagonalize
     */
    def diagonalizeStep (b: MatrixD)
    {
        import SVD.trailing
        import Eigen2.eigenvalues

        val tt = trailing (b)                                 // trailing 2-by-2 submatrix of b.t * b
        val td = tt(1, 1)                                     // last diagonal element in b.t * b
        val l  = eigenvalues (tt)                             // the eigenvalues of the submatrix
        val mu = if (abs (td - l(0)) <= abs (td - l(1))) l(0) else l(1)   // pick closest eigenvalue
        var y  = b(0, 0) * b(0, 0) - mu
        var z  = b(0, 0) * b(0, 1)

        for (k <- 0 until n-1) {

            // Givens rotation 1: k, k+1, theta1 (c1, s1)
            val (c1, s1) = givens (y, z)
            givensColUpdate (b, k, k+1, c1, s1)
            y = b(k, k); z = b(k+1, k)

            // Givens rotation 2: k, k+1, theta2 (c2, s2)
            val (c2, s2) = givens (y, z)
            givensRowUpdate (b, k, k+1, c2, s2)
            if (k < n-2) { y = b(k, k+1); z = b(k, k+2) }

        } // for
    } // diagonalizeStep

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the index of the first diagonal entry in 'b' from 'j' until 'k'
     *  that is zero; otherwise -1 (not found).
     *  @param b  the matrix whose diagoanl is being examined
     *  @param j  strart the search here
     *  @param k  end the search here
     */
    def findZero (b: MatrixD, j: Int, k: Int): Int =
    {
        for (i <- j until k if b(i, i) == 0.0) return i
        -1
    } // findZero

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Compute the Singular Value Decomposition (SVD) of matrix 'a'.
     *  @see Matrix Computation: Algorithm 8.6.2 SVD Algorithm.
     */
    def decompose ()
    {
        var (u, b, v) = bidiagonalize ()
        var q = 0
        var p = 0
 
        while (q < n) {
            for (i <- 0 until n-1) {
                if (b(i, i+1) < EPSILON * (abs (b(i, i)) + abs (b(i+1, i+1)))) b(i, i+1) = 0.0
            } // for
            q = 1; p = 0    // FIX

            if (q < n) {
                val k = findZero (b, p, n-q)
                if (k >= 0) {
                    b(k, k+1) = 0.0
                } else {
                    diagonalizeStep (b)
                    b = u.diag (p, q+m-n).t * b * v.diag (p, q)
                } // if
            } // if
        } // while

    } // decompose

} // SVD class


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The SVD companion object.
 */
object SVD
{
    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Compute the trailing 2-by-2 submatrix of 'b.t * b' without multiplying
     *  the full matrices.
     *  @param b  the given bidiagonal
     */
    def trailing (b: MatrixD): MatrixD =
    {
        val n3 = b.dim1 - 1
        val n2 = n3 - 1
        val n1 = n2 - 1
        val b12 = b(n1, n2)
        val b22 = b(n2, n2)
        val b23 = b(n2, n3)
        val b33 = b(n3, n3)
        new MatrixD ((2, 2), b12*b12 + b22*b22,  b22*b23,
                             b22*b23,  b23*b23 + b33*b33) 
    } // trailing

} // SVD object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** This object is used to test the SVD class.
 *  bidiagonalization answer = [ (21.8, -.613), (12.8, 2.24, 0.) ]
 */
object SVDTest extends App
{
    val a = new MatrixD ((4, 3), 1.0,  2.0,  3.0,
                                 4.0,  5.0,  6.0,
                                 7.0,  8.0,  9.0,
                                10.0, 11.0, 12.0)
    println ("a = " + a)
    val svd = new SVD (a)
    val (u, b, v) = svd.bidiagonalize ()
    println ("u     = " + u)
    println ("b     = " + b)    
    println ("v     = " + v)
    println ("ubv.t = " + u * b * v.t)    // should equal the original a

} // SVDTest object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `SVDTest2` object is used to test the SVD companion object.
 */
object SVDTest2 extends App
{
    import SVD.trailing

    val b = new MatrixD ((4, 4), 1.0, 5.0, 0.0, 0.0,    // the bidiagonal matrix
                                 0.0, 2.0, 6.0, 0.0,
                                 0.0, 0.0, 3.0, 7.0,
                                 0.0, 0.0, 0.0, 4.0)
    val n = b.dim1
    println ("b = " + b)
    println ("trailing b.t * b = " + trailing (b))
    println ("check: " + (b.t * b)(n-2 to n, n-2 to n))

} // SVDTest2 object