//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller, Michael E. Cotterell
 *  @version 1.6
 *  @date    Sun Jan 11 19:05:20 EST 2015
 *  @see     LICENSE (MIT style license file).
 *
 *  @title   Model: Simple Exponential Regression
 */

package scalation.analytics

import scala.collection.mutable.Set
import scala.math.{exp, log, sqrt}

import scalation.linalgebra.{MatriD, MatrixD, VectoD, VectorD}
import scalation.math.noDouble
import scalation.minima.QuasiNewton
import scalation.stat.Statistic
import scalation.plot.Plot

import Fit._
import Initializer._
import MatrixTransform._

//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `SimpleExpRegression` class supports exponential regression.  In this case,
 *  'x' is multi-dimensional [1, x_1, ... x_k].  Fit the parameter vector 'b' in the
 *  exponential regression equation
 *  <p>
 *      log (mu (x))  =  b dot x  =  b_0 + b_1 * x_1 + ... b_k * x_k
 *  <p>
 *  @see www.stat.uni-muenchen.de/~leiten/Lehre/Material/GLM_0708/chapterGLM.pdf 
 *  @param x       the data/input matrix
 *  @param y       the response/output vector
 *  @param fname_  the feature/variable names
 *  @param hparam  the hyper-parameters (currently none)
 *  @param nonneg  whether to check that responses are nonnegative
 */
class SimpleExpRegression (x: MatriD, y: VectoD,
                     fname_ : Strings = null, hparam: HyperParameter = null,
                     nonneg: Boolean = true)
      extends PredictorMat (x, y, fname_, hparam) with NoFeatureSelectionMat
{
    if (nonneg && ! y.isNonnegative) flaw ("constructor", "response vector y must be nonnegative")

    private val DEBUG     = true                                             // debug flag
    private val cutoff    = 1E-4                                             // cutoff threshold
    private var eta       = 0.00001                                              // the learning/convergence rate (requires adjustment)
    private val maxEpochs = 100                                              // the maximum number of training epcochs/iterations
//  private var eta       = hparam ("eta")                                   // the learning/convergence rate (requires adjustment)
//  private val maxEpochs = hparam ("maxEpochs").toInt                       // the maximum number of training epcochs/iterations

    def f (xi: Double, b: VectoD): Double = b(0) * exp (b(1) * xi)

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Train the predictor by fitting the parameter vector (b-vector) in the
     *  simple exponential regression equation.
     *  @param x_  the training/full data/input matrix
     *  @param y_  the training/full response/output vector
     */
    def train (x_ : MatriD = x, y_ : VectoD = y): SimpleExpRegression =
    {
        val xx = x_.col (0)
        if (b == null) b = weightVec (2)                                     // initial random/guessed parameter vector
        b(0) = 50.0
        b(1) = -0.03

        for (epoch <- 1 until maxEpochs) {                                   // iterate over the learning epochs
            val yp = VectorD (for (i <- xx.range) yield f (xx(i), b) )       // y-predicted vector
            e      = y_ - yp                                                 // error vector
            val d0 = (e dot yp) / b(0)                                       // delta 0
            val d1 = e dot yp * xx                                           // delta 1
            b(0)  += eta * d0                                                // update to first parameter
            b(1)  += eta * d1                                                // update to second parameter
            if (DEBUG) println (s"train: parameters for $epoch th epoch: b = $b \n yp = $yp")
            if (VectorD (d0, d1).norm < cutoff) return this                  // stopping rule, may refine
        } // for
        this
    } // train

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Predict the value of 'y = f(z)' by evaluating the formula 'y = exp (b dot z)',
     *  e.g., 'exp (b_0, b_1, b_2) dot (1, z_1, z_2)'.
     *  @param z  the new vector to predict
     */
    override def predict (z: VectoD): Double = f (z(0), b)

    def predictVec (z: VectoD): VectoD = VectorD (for (i <- z.range) yield f (z(i), b))

} // SimpleExpRegression class


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `SimpleExpRegression` companion object provides factory apply functions and a testing method.
 */
object SimpleExpRegression extends ModelFactory
{
    val drp = (null, null, true)                                   // default remaining parameters

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Create an `SimpleExpRegression` object from a combined data-response matrix.
     *  The last column is assumed to be the response column.
     *  @param xy         the combined data-response matrix (predictors and response)
     *  @param fname      the feature/variable names
     *  @param hparam     the hyper-parameters
     *  @param nonneg  whether to check that responses are nonnegative
     */
    def apply (xy: MatriD, fname: Strings = null, hparam: HyperParameter = null,
               nonneg: Boolean = true): SimpleExpRegression =
    {
        val n = xy.dim2
        if (n < 2) {
            flaw ("apply", s"dim2 = $n of the 'xy' matrix must be at least 2")
            null
        } else {
            val (x, y) = pullResponse (xy)
            new SimpleExpRegression (x, y, fname, hparam, nonneg)
       } // if
    } // apply

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Create an `SimpleExpRegression` object from a data matrix and a response vector.
     *  This factory function provides data rescaling.
     *  @see `ModelFactory`
     *  @param x          the data/input m-by-n matrix
     *                        (augment with a first column of ones to include intercept in model)
     *  @param y          the response/output m-vector
     *  @param fname      the feature/variable names (use null for default)
     *  @param hparam     the hyper-parameters (use null for default)
     *  @param nonneg  whether to check that responses are nonnegative
     */
    def apply (x: MatriD, y: VectoD, fname: Strings, hparam: HyperParameter,
               nonneg: Boolean): SimpleExpRegression =
    {
        val n = x.dim2
        if (n < 1) {
            flaw ("apply", s"dim2 = $n of the 'x' matrix must be at least 1")
            null
        } else if (rescale) {                                   // normalize the x matrix
            val (mu_x, sig_x) = (x.mean, stddev (x))
            val xn = normalize (x, (mu_x, sig_x))
            new SimpleExpRegression (xn, y, fname, hparam, nonneg)
        } else {
            new SimpleExpRegression (x, y, fname, hparam, nonneg)
       } // if
    } // apply

} // SimpleExpRegression object


//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `SimpleExpRegressionTest` object tests `SimpleExpRegression` class using the following
 *  exponential regression problem.
 *  > runMain scalation.analytics.SimpleExpRegressionTest
 */
object SimpleExpRegressionTest extends App
{
    val x = new MatrixD ((5, 3), 1.0, 36.0,  66.0,               // 5-by-3 matrix
                                 1.0, 37.0,  68.0,
                                 1.0, 47.0,  64.0,
                                 1.0, 32.0,  53.0,
                                 1.0,  1.0, 101.0)
    val y = VectorD (745.0, 895.0, 442.0, 440.0, 1598.0)
    val z = VectorD (1.0, 20.0, 80.0)

    println ("x = " + x)

    val erg = new SimpleExpRegression (x, y)
    erg.train ().eval ()
    println (erg.report)

    val yp = erg.predict (x)
    println ("y  = " + y)
    println ("yp = " + yp)

    val t = VectorD.range (0, y.dim)
    new Plot (t, y, yp, "SimpleExpRegressionTest")

    val yp2 = erg.predict (z)
    println (s"predict ($z) = $yp2")

} // SimpleExpRegressionTest object


//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `SimpleExpRegressionTest2` object has a basic test for the `SimpleExpRegression` class.
 *  > runMain scalation.analytics.SimpleExpRegressionTest2
 */
object SimpleExpRegressionTest2 extends App
{
    import scalation.random.{Uniform, Exponential, Random}

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Test `SimpleExpRegression` by simulating 'n'-many observations.
     *  @param n  number of observations
     *  @param k  number of variables
     *  @return   (actual, estimated, r^2)
     */
    def test (n: Int = 10000, k: Int = 5): (Int, Int, VectorD, VectoD, Double) =
    {
        val u = new Uniform (0, 1)     // uniform random
        val e = new Exponential (1)    // exponential error
        val r = new Random ()

        val x = new MatrixD (n, k)     // data matrix
        val y = new VectorD (x.dim1)   // response vector
        val b = new VectorD (k)        // known coefficients

        for (i <- b.range) b(i) = 1 + r.gen * 6

        for (i <- x.range1; j <- x.range2) x(i, j) = u.gen

        for (i <- y.range) y(i) = exp (x(i) dot b) * e.gen

        val erg = new SimpleExpRegression (x, y)
        erg.train ().eval ()

        (n, k, b, erg.parameter, erg.fit(0))
    } // test

    val tests = Array.ofDim [(Int, Int, VectorD, VectoD, Double)] (10)

    for (k <- tests.indices) tests(k) = test (1000, k + 1)

    tests.foreach {
        case (n, k, actual, fit, rSquared) => {
            actual.setFormat ("%.4f, ")
            fit.setFormat ("%.4f, ")
            println ("nobs = %d, regressors = %d, R^2 = %f\nactual = %s\nfir    = %s\n".format(n, k, rSquared, actual, fit))
        } // case
    } // foreach

} // SimpleExpRegressionTest2


//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `SimpleExpRegressionTest3` object has a basic test for the `SimpleExpRegression` class.
 *  @see http://www.cnachtsheim-text.csom.umn.edu/kut86916_ch13.pdf
 *  <p>
 *      y = 58.6065 exp (-.03959 x) + e
 *  <p>
 *  > runMain scalation.analytics.SimpleExpRegressionTest3
 */
object SimpleExpRegressionTest3 extends App
{
    val xy = new MatrixD ((15, 2), 2, 54,
                                   5, 50,
                                   7, 45,
                                  10, 37,
                                  14, 35,
                                  19, 25,
                                  26, 20,
                                  31, 16,
                                  34, 18,
                                  38, 13,
                                  45, 8,
                                  52, 11,
                                  53, 8,
                                  60, 4,
                                  65, 6)

    println ("xy = " + xy)
    val x = xy.col(0)
    val y = xy.col(1)

    val erg = new SimpleExpRegression (MatrixD (Seq (x)), y)
    erg.train ().eval ()
    println (erg.report)

    val yp = erg.predictVec (x)
    println ("y  = " + y)
    println ("yp = " + yp)

    new Plot (null, y, yp, "SimpleExpRegressionTest")

} // SimpleExpRegressionTest3