//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller, Naman Fatehpuria
 *  @version 1.6
 *  @date    Mon Jul 29 14:09:25 EDT 2013
 *  @see     LICENSE (MIT style license file).
 *
 *  @title   Data Reduction: Non-negative Matrix Factorization (NMF)
 */

package scalation.analytics

import scala.math.{ceil, min}

import scalation.linalgebra.{MatriD, MatrixD}
import scalation.util.banner

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `NMFactorization` class factors a matrix 'x' into two non negative matrices 
 *  'w' and 'h' such that 'x = wh' approximately.
 *  @see  en.wikipedia.org/wiki/Non-negative_matrix_factorization
 *  @param x      the matrix to be factored 
 *  @param loops  the number of iterations before checking the termination condition  
 *  @param r      factor the m-by-n matrix 'x' in to two factors: an m-by-r and r-by-n matrix  
 */ 
class NMFactorization (x: MatriD, loops: Int = 10, var r: Int = 0)
      extends Reducer
{
    private val DEBUG = true                            // debug flag 
    private val m = x.dim1                              // number of rows in matrix x
    private val n = x.dim2                              // number of columns in matrix x
    private val EPSILON = 1.0E-4                        // number close to zero

    if (r == 0) r = ceil (min (m, n).toDouble / 2).toInt + 1

    private val w  = new MatrixD (m, r)                 // left factor (m-by-r matrix)
    private val h  = new MatrixD (r, n)                 // right factor (r-by-n matrix)
    w.set (1.0); h.set (1.0)                            // initialize all matrix elements to 1

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Reduce the original data matrix by producing a lower dimensionality matrix
     *  that maintains most of the descriptive power of the original matrix.
     */
    def reduce (): MatriD =
    {
        throw new UnsupportedOperationException ("NMFactorization.reduce: should use 'factorReduce', not 'reduce'")
    } // reduce

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Factor and reduce the original matrix 'x' into left 'w' and right 'h' matrices
     *  by iteratively calling the 'update' method until the difference between
     *  'x' and its reduced approximation 'xr = w * h' is sufficiently small.
     */
    override def factorReduce (): (MatrixD, MatrixD) = 
    {
        var xr: MatriD = null                           // holds product of the factors (xr = w*h)

        println ("Original Matrix to Factor x = " + x)
        var converged = false                           // convergence flag
        var i         = 0                               // number of update steps performed
 
        while (! converged) {
            for (k <- 0 until loops) update ()          // perform several update steps
    
            xr = w * h                                  // compute new product of factors
            val dist = (xr - x).norm1                   // compute norm of the difference
            if (DEBUG) {
                println ("Result for iteration " + i + ":"); i += loops
                println ("w  = " + w + "\nh  = " + h)
                println ("xr = " + xr)
                println ("dist = " + dist)
            } // if

            if (dist < EPSILON) converged = true            // quit when distance is small enough
        } // while

        println ("Factors w = " + w + "and     h = " + h)   // print factors
        println ("Product of Factors xr = " + xr)           // print approximation
        (w, h)                                              // return factors
    } // factorReduce

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Perform a multiplicative update step to create better estimates for
     *  factor matrices w and h.
     */
    private def update (): Unit =
    {
        val x_ht = x * h.t
        val h_ht = h * h.t
        for (i <- 0 until m; a <- 0 until r) {    
            w(i, a) *= x_ht(i, a) / (w(i) dot h_ht.col(a))
        } // for

        val wt_x = w.t * x
        val wt_w = w.t * w
        for (b <- 0 until r; j <- 0 until n) {
            h(b, j) *= wt_x(b, j) / (wt_w(b) dot h.col(j))
        } // for              
    } // update

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Approximately recover the original matrix.  The new matrix will have
     *  the same dimensionality, but will have some lose of information.
     */
    def recover (): MatriD = w * h

} // NMFactorization class


//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `NMFactorizationTest` object to test `NMFactorizationTest` class.
 *  > runMain scalation.analytics.NMFactorizationTest
 */
object NMFactorizationTest extends App
{
    val x = new MatrixD ((4, 6), 1.0, 2.0, 3.0,  4.0,  5.0,  6.0,
                                 4.0, 5.0, 6.0,  7.0,  8.0,  6.0,
                                 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
                                 7.0, 8.0, 9.0, 10.0, 11.0, 12.0)

    val nmf = new NMFactorization (x)
    val (w, h) = nmf.factorReduce ()

    banner ("NMF Results: Factor 'x'")
    println ("x = " + x)
    banner ("into")
    println ("w = " + w)
    banner ("and")
    println ("h = " + w)
    banner ("giving reduced rank approximation")
    println ("xr = " + nmf.recover ())

} // NMFactorizationTest object