//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller
 *  @version 1.6
 *  @date    Sat Jun 13 13:37:06 EDT 2020
 *  @see     LICENSE (MIT style license file).
 *
 *  @title   Model: Feed-Forward Neural Network for Time Series
 */

package scalation.analytics
package forecaster

import scala.collection.mutable.{ArrayBuffer, Set}

import scalation.linalgebra.{FunctionV_2V, MatriD, MatrixD, VectoD, VectorD}
import scalation.linalgebra.VectorD.one
import scalation.math.{double_exp, noDouble}
import scalation.plot.PlotM
import scalation.stat.Statistic
import scalation.util.banner

import ActivationFun._
//import ImputeBackward.impute
import ImputeMean.impute
import MatrixTransform._
import PredictorMat2.{rescaleX, rescaleY}

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `NeuralNet_3L1_4TS` class uses a neural network to fit the lagged data. 
 *  Lag columns ranging from 'lag1' (inclusive) to 'lag2' (inclusive) are added before
 *  delegating the problem to the `NeuralNet_3L1` class.  A constant term for intercept
 *  can be added (@see 'allForms' method) but must not include intercept (column of ones)
 *  in initial data matrix.
 *  @param x_      the m-by-nx input/data matrix (training data consisting of m input vectors)
 *                     the initial input/data matrix (before lag term expansion)
 *  @param y       the output/response m vector (training data consisting of m output values)
 *  @param nz      the number of nodes in hidden layer (-1 => use default formula)
 *  @param fname_  the feature/variable names (if null, use x_j's)
 *  @param hparam  the hyper-parameters for the model/network
 *  @param f0      the activation function family for layers 1->2 (input to hidden)
 *  @param f1      the activation function family for layers 2->3 (hidden to output)
 *  @param itran   the inverse transformation function returns responses to original scale
 */
class NeuralNet_3L1_4TS (x_ : MatriD, y: VectoD, nz: Int = -1,
                         fname_ : Strings = null,
                         hparam: HyperParameter = NeuralNet_3L1_4TS.hp,
                         f0: AFF = f_tanh, f1: AFF = f_id,
                         itran: FunctionV_2V = null)
      extends NeuralNet_3L1 (Regression4TS.allForms (x_, hparam ("lag1").toInt, hparam ("lag2").toInt, true),
                             y, nz, fname_, hparam, f0, f1, itran)
      with  ForecasterMat
{
    private val DEBUG = false                                          // debug flag

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the model name including its current hyper-parameter.
     */
    override def modelName: String = s"NeuralNet_3L1_4TS($nx)"

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Compute the error (difference between actual and forecasted) and useful
     *  diagnostics for the test dataset.
     *  @param ym  the mean of the actual response/output vector (full/training)
     *  @param yy  the actual response/output vector (full/testing)
     *  @param yf  the forecasted response/output vector (full/testing)
     */
    override def eval (ym: Double, yy: VectoD, yf: VectoD): NeuralNet_3L1_4TS =
    {
        if (DEBUG) {
            println (s"eval: ym     = $ym")
            println (s"eval: yy.dim = ${yy.dim}")
            println (s"eval: yf.dim = ${yf.dim}")
        } // if
        val e = yy - yf                                                // compute residual/error vector e, ee(0) ?
        fit.diagnose (e, yy, yf, null, ym)                             // compute diagnostics
        this
    } // eval

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Produce a forecast for 'h' steps ahead into the future.
     *  Note: the forecasts for time 't = 0, ... , h-1' will be duplicates.  
     *  @param xe  the relevant expanded data matrix
     *  @param t   the time for the forecast
     *  @param h   the forecasting horizon, number of steps ahead to produce forecast
     */
    def forecast (xe: MatriD, t: Int, h: Int = 1): Double = 
    {
        val tt = math.max (0, t+1-h)                                   // time @ row h-1 back in matrix
        predict (xe(tt))
    } // forecast

} // NeuralNet_3L1_4TS class


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `NeuralNet_3L1_4TS` companion object provides factory functions and functions
 *  for creating functional forms.
 */
object NeuralNet_3L1_4TS extends ModelFactory
{
    /** Base hyper-parameter specification for `NeuralNet_3L1_4TS`
     */
    val hp = Optimizer.hp ++ new HyperParameter
    hp += ("lag1",   1, 1)                                        // first lag included (inclusive)
    hp += ("lag2",   2, 2)                                        // last lag included (inclusive)

    private val DEBUG = false                                     // debug flag
    
    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Create a `NeuralNet_3L1_4TS` object from a combined data-response matrix.
     *  @param xy      the m-by-(nx+1) combined input/data matrix and output/response vector 
     *  @param nz      the number of nodes in hidden layer (-1 => use default formula)
     *  @param fname   the feature/variable names (if null, use x_j's)
     *  @param hparam  the hyper-parameters for the model/network
     *  @param f0      the activation function family for layers 1->2 (input to hidden)
     *  @param f1      the activation function family for layers 2->3 (hidden to output)
     */
    def apply (xy: MatriD, nz: Int = -1, fname: Strings = null, hparam: HyperParameter = hp,
               f0: AFF = f_tanh, f1: AFF = f_id): NeuralNet_3L1_4TS =
    {
        var itran: FunctionV_2V = null                                 // inverse transform -> original scale
        val (x, y) = pullResponse (xy)                                 // assumes the last column is the response

        val x_s = if (rescale) rescaleX (x, f0)
                  else x
        val y_s = if (f1.bounds != null) { val y_i = rescaleY (y, f1); itran = y_i._2; y_i._1 }
                  else y

        if (DEBUG) println (s" scaled: x = $x_s \n scaled y = $y_s")
        new NeuralNet_3L1_4TS (x_s, y_s, nz, fname, hparam, f0, f1, itran)
    } // apply 

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Create a `NeuralNet_3L1_4TS` object from a data matrix and a response vector.
     *  This factory function provides data rescaling.
     *  @see `ModelFactory`
     *  @param x       the initial input/data matrix (before expansion)
     *  @param y       the output/response m-vector
     *  @param nz      the number of nodes in hidden layer (-1 => use default formula)
     *  @param fname   the feature/variable names (if null, use x_j's)
     *  @param hparam  the hyper-parameters for the model/network
     *  @param f0      the activation function family for layers 1->2 (input to hidden)
     *  @param f1      the activation function family for layers 2->3 (hidden to output)
     */
    def apply (x: MatriD, y: VectoD, nz: Int, fname: Strings, hparam: HyperParameter,
               f0: AFF, f1: AFF): NeuralNet_3L1_4TS =
    {
        val hp2 = if (hparam == null) hp else hparam
        var itran: FunctionV_2V = null                                 // inverse transform -> original scale

        val x_s = if (rescale) rescaleX (x, f0)
                  else x
        val y_s = if (f1.bounds != null) { val y_i = rescaleY (y, f1); itran = y_i._2; y_i._1 }
                  else y

        if (DEBUG) println (s" scaled: x = $x_s \n scaled y = $y_s")
        new NeuralNet_3L1_4TS (x_s, y_s, nz, fname, hp2, f0, f1, itran)
    } // apply

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** The number of terms include current value and lag one value.
     *  when there are no cross-terms.
     *  @param k  number of features/predictor variables (not counting intercept)
     */
    override def numTerms (k: Int) =  2 * k + 1                 // FIX

} // NeuralNet_3L1_4TS object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `NeuralNet_3L1_4TSTest` object is used to test the `NeuralNet_3L1_4TS` class.
 *  > runMain scalation.analytics.forecaster.NeuralNet_3L1_4TSTest
 */
object NeuralNet_3L1_4TSTest extends App
{
    import ExampleBPressure.{x01 => x, y}

    val rg4 = new NeuralNet_3L1_4TS (x, y)
    rg4.analyze ()
    val nTerms = NeuralNet_3L1_4TS.numTerms (2)
    println (s"x = ${rg4.getX}")
    println (s"y = $y")

    println (s"nTerms = $nTerms")
    println (rg4.report)

    banner ("Forward Selection Test")
    rg4.forwardSelAll (cross = false)

    banner ("Backward Elimination Test")
    rg4.backwardElimAll (cross = false)

} // NeuralNet_3L1_4TSTest object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `NeuralNet_3L1_4TSTest2` object is used to test the `NeuralNet_3L1_4TS` class.
 *  The 'x' matrix in one dimensional.
 *  > runMain scalation.analytics.forecaster.NeuralNet_3L1_4TSTest2
 */
object NeuralNet_3L1_4TSTest2 extends App
{
    import scalation.random.Normal
    import scalation.plot.Plot

    val (m, n) = (400, 1)
    val noise = new Normal (0, 10 * m * m)
    val x = new MatrixD (m, n)
    val y = new VectorD (m)
    val t = VectorD.range (0, m)

    for (i <- x.range1) { 
        x(i, 0) = i + 1
        y(i) = i*i + i + noise.gen
    } // for

    banner ("Regression")
    val ox = VectorD.one (y.dim) +^: x
    val rg = new Regression (ox, y)
    rg.analyze ()
    println (rg.report)
    val yp = rg.predict ()
    val e  = rg.residual

    banner ("NeuralNet_3L1_4TS")
    val hp2 = NeuralNet_3L1_4TS.hp.updateReturn (("lag1", 1.0), ("lag2", 2.0))
    val rg4 = new NeuralNet_3L1_4TS (x, y, hparam = hp2)
    rg4.analyze ()
    println (rg4.report)
    val yp4 = rg4.predict ()
    val e4  = rg4.residual

    val x0 = x.col(0)
    new Plot (x0, y, null, "y vs x")
    new Plot (t, y, yp, "y and yp vs t")
    new Plot (t, y, yp4, "y and yp4 vs t")
    new Plot (x0, e, null, "e vs x")
    new Plot (x0, e4, null, "e4 vs x")

} // NeuralNet_3L1_4TSTest2 object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `NeuralNet_3L1_4TSTest3` object is used to test the `NeuralNet_3L1_4TS` class.
 *  The 'x' matrix in two dimensional.
 *  > runMain scalation.analytics.forecaster.NeuralNet_3L1_4TSTest3
 */
object NeuralNet_3L1_4TSTest3 extends App
{
    import scalation.random.Normal
    import scalation.stat.StatVector.corr

    val s      = 20
    val grid   = 1 to s
    val (m, n) = (s*s, 2)
    val noise  = new Normal (0, 10 * s * s)
    val x = new MatrixD (m, n)
    val y = new VectorD (m)

    var k = 0
    for (i <- grid; j <- grid) {
        x(k) = VectorD (i, j)
        y(k) = x(k, 0)~^2 + 2 * x(k, 1) +  noise.gen
        k += 1
    } // for

    banner ("Regression")
    val ox = VectorD.one (y.dim) +^: x
    val rg = new Regression (ox, y)
    rg.analyze ()
    println (rg.report)

    banner ("NeuralNet_3L1_4TS")
    val rg4 = new NeuralNet_3L1_4TS (x, y)
    rg4.analyze ()
    println (rg4.report)

    banner ("Multi-collinearity Check")
    val x4 = rg4.getX
    println (corr (x4.asInstanceOf [MatrixD]))
    println (s"vif = ${rg4.vif ()}")

} // NeuralNet_3L1_4TSTest3 object


//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `NeuralNet_3L1_4TSTest4` object tests the `NeuralNet_3L1_4TS` class using the AutoMPG
 *  dataset.  It illustrates using the `Relation` class for reading the data
 *  from a .csv file "auto-mpg.csv".  Assumes no missing values.
 *  It also combines feature selection with cross-validation and plots
 *  R^2, R^2 Bar and R^2 cv vs. the instance index.
 *  > runMain scalation.analytics.forecaster.NeuralNet_3L1_4TSTest4
 */
object NeuralNet_3L1_4TSTest4 extends App
{
    import scalation.columnar_db.Relation

    banner ("auto_mpg relation")
    val auto_tab = Relation (BASE_DIR + "auto-mpg.csv", "auto_mpg", null, -1)
    auto_tab.show ()

    banner ("auto_mpg (x, y) dataset")
    val (x, y) = auto_tab.toMatriDD (ArrayBuffer.range (1, 6), 0)
    println (s"x = $x")
    println (s"y = $y")

    banner ("auto_mpg regression")
    val rg4 = new NeuralNet_3L1_4TS (x, y)
    rg4.analyze ()
    val n  = x.dim2                                                  // number of variables
    val nt = NeuralNet_3L1_4TS.numTerms (n)                           // number of terms
    println (s"n = $n, nt = $nt")
    println (rg4.report)

    banner ("Forward Selection Test")
    val (cols, rSq) = rg4.forwardSelAll ()                           // R^2, R^2 Bar, R^2 cv
    val k = cols.size
    val t = VectorD.range (1, k)                                     // instance index
    new PlotM (t, rSq, lines = true)
    println (s"rSq = $rSq")
    
/*
    banner ("Forward Selection Test")
    val qx  = rg4.getX                                               // get matrix with all columns
    val rg  = new Regression (qx, y)                                 // regression on all columns
    val rSq = new MatrixD (qx.dim2-1, 3)                             // R^2, R^2 Bar,  R^2 cv
    
    val fcols = Set (0)                                              // start with x_0 in model
    for (l <- 1 until nt) {
        val (x_j, b_j, fit_j) = rg.forwardSel (fcols)                // add most predictive variable
        println (s"forward model: add x_j = $x_j with b = $b_j \n fit = $fit_j")
        if (x_j == -1) {
            println (s"the 'forwardSel' could not find a variable to add: x_j = $x_j")
        } else { 
            fcols += x_j                                             // add variable x_j
            val x_cols = qx.selectCols (fcols.toArray)               // qx projected onto cols_j columns
            rSq(l-1)   = Fit.qofVector (fit_j, rg.crossVal (x_cols))   // use main model, 'rg'
//          val rg_j   = new Regression (x_cols, y)                  // new model, regress with x_j added
//          rSq(l-1)   = Fit.qofVector (fit_j, rg_j.crossVal ())     // use new model, rg_j
        } // if
    } // for

    val k = fcols.size
    for (l <- k until nt) rSq(l-1) = rSq(l-2)
    println (s"rSq = $rSq")
    val t = VectorD.range (1, k)                                     // instance index
    new PlotM (t, rSq.t, lines = true)
*/

} // NeuralNet_3L1_4TSTest4 object


//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `NeuralNet_3L1_4TSTest5` object tests `NeuralNet_3L1_4TS` class using the following
 *  regression equation.
 *  <p>
 *      y = b dot x = b_0 + b_1*x1 + b_2*x_2.
 *  <p>
 *  > runMain scalation.analytics.forecaster.NeuralNet_3L1_4TSTest5
 */
object NeuralNet_3L1_4TSTest5 extends App
{
    // 4 data points:             x1 x2  y
    val xy = new MatrixD ((9, 3), 2, 1, 0.4,               // 4-by-3 matrix
                                  3, 1, 0.5,
                                  4, 1, 0.6,
                                  2, 2, 1.0,
                                  3, 2, 1.1,
                                  4, 2, 1.2,
                                  2, 3, 2.0,
                                  3, 3, 2.1,
                                  4, 3, 2.2)

    println ("model: y = b0 + b1*x1 b2*x1^2 + b3*x2 + b4*x2^2")
    println (s"xy = $xy")

    val oxy = VectorD.one (xy.dim1) +^: xy
    val xy_ = oxy.selectCols (Array (0, 2, 3))
    println (s"xy_ = $xy_")

    banner ("SimpleRegression")
    val srg  = SimpleRegression (xy_)
    srg.analyze ()
    println (srg.report)
    println (s"predict = ${srg.predict ()}")

    banner ("Regression")
    val rg  = Regression (oxy)
    rg.analyze ()
    println (rg.report)
    println (s"predict = ${rg.predict ()}")

    banner ("NeuralNet_3L1_4TS")
    val rg4 = NeuralNet_3L1_4TS (xy)
    rg4.analyze ()
    println (rg4.report)
    println (s"predict = ${rg4.predict ()}")

} // NeuralNet_3L1_4TSTest5 object