//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller
 *  @version 1.1
 *  @date    Tue Sep 17 17:12:07 EDT 2013
 *  @see     LICENSE (MIT style license file).
 */

package scalation.linalgebra

import math.{abs, hypot}

import scalation.linalgebra.MatrixD.eye

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `Givens` objects has methods for determinng values 'c = cos(theta)' and
 *  's = sin(theta) for Givens rotation matrices as well as methods for applying
 *  Givens rotations.
 */
object Givens
{
    private val EPSILON = 1E-12            // a value close to zero
    private val DEBUG   = false            // debug flag

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Create the values for a Givens 2-by-2 rotation matrix.  Given scalars
     *  'y' and 'z', efficiently compute 'c = cos(theta)' and 's = sin(theta)'
     *  that can be used to form the rotation matrix.
     *  @see Algorithm 5.1.3 in Matrix Computation.
     *  @param y  the first scalar
     *  @param z  the second scalar
     */
    def givens (y: Double, z: Double): Tuple2 [Double, Double] =
    {
        val r = hypot (y, z)    // radius: sqrt (y~^2, z~^2) avoiding over/under-flow

        if (z == 0.0) (1.0, 0.0)
        else (y / r, -z / r)
/*
        else if (abs (z) > abs (y)) (z / r, -y / r)
        else (y / r, -z / r)
        if (zz > yy) {
            r = -y/z; s = 1.0 / sqrt (1.0 + r*r); c = s * r
        } else if (zz > EPSILON) {
            r = -z/y; c = 1.0 / sqrt (1.0 + r*r); s = c * r
        } // if
        (c, s)
//      (s, c)              // FIX - worked for SVD.bak2 when flipped?
*/
    } // givens

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Efficiently perform a Givens row update: 'a = g (i, k, theta).t * a'.
     *  The update just affects two row.
     *  @param a  the matrix to update
     *  @param i  the first row
     *  @param k  the second row
     *  @param c  the cos(theta)
     *  @param s  the sin(theta)
     */
    def givensRowUpdate (a: MatrixD, i: Int, k: Int, c: Double, s: Double)
    {
        for (j <- 0 until a.dim2) {
            val r1 = a(i, j)
            val r2 = a(k, j)
            a(i, j) = c * r1 - s * r2
            a(k, j) = s * r1 + c * r2
        } // for
    } // givensRowUpdate

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Efficiently perform a Givens column update: a = a * g (i, k, theta).
     *  The update just affects two columns.
     *  @param a  the matrix to update
     *  @param j  the first column
     *  @param k  the second column
     *  @param c  the cos(theta)
     *  @param s  the sin(theta)
     */
    def givensColUpdate (a: MatrixD, j: Int, k: Int, c: Double, s: Double)
    {
        for (i <- 0 until a.dim1) {
            val r1 = a(i, j)
            val r2 = a(i, k)
            a(i, j) = c * r1 - s * r2
            a(k, j) = s * r1 + c * r2
        } // for
    } // givensColUpdate

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return a Givens rotation matrix with angle 'theta = atan(s/c)'.  The
     *  2-by-2 rotation is embedded in an identity matrix of dimension 'n'
     *  @param i  the first diagonal position i, i
     *  @param k  the second diagonal position k, k
     *  @param c  the cosine of theta
     *  @param s  the sine of theta
     *  @param n  the dimension of the resulting rotation matrix
     */
    def givensRo (i: Int, k: Int, c: Double, s: Double, n: Int): MatrixD =
    {
        val b = eye (n)
        b(i, i) =  c; b(i, k) =  s
        b(k, i) = -s; b(k, k) =  c
        if (DEBUG) println ("givensRo: b = " + b)
        b
    }  // givensRo

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return a transposed Givens rotation matrix with angle 'theta = atan(s/c)'.
     *  The 2-by-2 rotation is embedded in an identity matrix of dimension 'n'
     *  @param i  the first diagonal position i, i
     *  @param k  the second diagonal position k, k
     *  @param c  the cosine of theta
     *  @param s  the sine of theta
     *  @param n  the dimension of the resulting rotation matrix
     */
    def givensRoT (i: Int, k: Int, c: Double, s: Double, n: Int): MatrixD =
    {
        val b = eye (n)
        b(i, i) =  c; b(i, k) = -s
        b(k, i) =  s; b(k, k) =  c
        if (DEBUG) println ("givensRoT: b = " + b)
        b
    }  // givensRoT

} // Givens object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `GivensTest` object tests the `Givens` object by using two rotations to
 *  turn matrix 'a' into an upper triangular matrix.
 *  @see http://en.wikipedia.org/wiki/Givens_rotation
 */
object GivensTest extends App
{
    import Givens._

    var a = new MatrixD ((3, 3), 6.0, 5.0, 0.0,
                                 5.0, 1.0, 4.0,
                                 0.0, 4.0, 3.0)
    println ("a = " + a)

    val (c1, s1) = givens (a(0, 0), a(1, 0))    // rotate (1, 0) 5 -> 0 
    val g1 = givensRoT (0, 1, c1, s1, 3)
    a = g1 * a
    println ("(c1, s1) = " + (c1, s1))
    println ("rotation 1: a = " + a)

    val (c2, s2) = givens (a(1, 1), a(2, 1))    // rotate (2, 1) 4 -> 0 
    val g2 = givensRoT (1, 2, c2, s2, 3)
    a = g2 * a
    println ("(c2, s2) = " + (c2, s2))
    println ("rotation 2: a = " + a)

} // GivensTest object