//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller
 *  @version 1.6
 *  @date    Sat Jun 15 13:19:05 EDT 2019
 *  @see     LICENSE (MIT style license file).
 *
 *  @title   Model: Extreme-Learning Machines with Quadratic Cross Regression
 *
 *  @see     www.ntu.edu.sg/home/egbhuang/pdf/ELM-Unified-Learning.pdf
 *  @see     www.sciencedirect.com/science/article/pii/S0893608014002214
 */

package scalation.analytics

import scala.collection.mutable.Set

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

import ActivationFun._
import Fit._
import Initializer._
import MatrixTransform._
import PredictorMat2._

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `QuadXELM_3L1` class supports multi-output, 3-layer (input, hidden and output)
 *  Extreme-Learning Machines.  It can be used for both classification and prediction,
 *  depending on the activation functions used.  Given several input vectors and output
 *  vectors (training data), fit the parameters 'a' and 'b' connecting the layers,
 *  so that for a new input vector 'v', the net can predict the output value, i.e.,
 *  <p>
 *      yp = forms (f1 (b * f0 (a * v)))
 *  <p>
 *  where 'f0' and 'f1' are the activation functions and the parameter 'a' and 'b'
 *  are the parameters between input-hidden and hidden-output layers.
 *  Unlike `NeuralNet_2L` which adds input 'x0 = 1' to account for the intercept/bias,
 *  `QuadXELM_3L1` explicitly adds bias.
 *  'forms' expands a vector to include quadratic forms with cross terms.
 *  Note: only uses 'f_id' implicitly, use `ELM_3L` for other options for 'f1'.
 *  @param x       the m-by-n input matrix (training data consisting of m input vectors)
 *  @param y       the m output vector (training data consisting of m output scalars)
 *  @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 itran   the inverse transformation function returns responses to original scale
 */
class QuadXELM_3L1 (x: MatriD, y: VectoD,
              private var nz: Int = -1,
              fname_ : Strings = null, hparam: HyperParameter = null,
              f0: AFF = f_tanh,
              val itran: FunctionV_2V = null)
      extends PredictorMat (x, y, fname_, hparam)
{
    private val DEBUG = false                                              // debug flag

    if (nz < 1) nz = 3 * n + 1                                             // default number of nodes for hidden layer
    val df_m = compute_df_m (nz)                                           // degrees of freedom for model (first output only)
    resetDF (df_m, x.dim1 - df_m)                                          // degrees of freedom for (model, error)
 
    private val s = 8                                                      // random number stream to use (0 - 999)

    private var a = new NetParam (weightMat3 (n, nz, s),
                                  weightVec3 (nz, s))                      // parameters (weights & biases) in to hid (fixed)

    println (s"Create an QuadXELM_3L1 with $n input, $nz hidden and 1 output nodes: df_m = $df_m")

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Compute the degrees of freedom for the model (based on 'n, nz_, ny = 1').
     *  Rough extimate based on total number of parameters - 1.
     *  @param nz_  the number of nodes in the hidden layer
     */
    def compute_df_m (nz_ : Int): Int = nz_

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the parameters 'b'.  Since the 'a' weights are fixed, only return 'b'.
     */
    def parameters: VectoD = b

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given training data 'x_' and 'y_', with parameters 'a' fixed, fit parameters 'b'.
     *  Use matrix factorization in `QuadXRegression` to find optimal values for the
     *  parameters/weights 'b'.
     *  @param x_  the training/full data/input matrix
     *  @param y_  the training/full response/output vector
     */
    def train (x_ : MatriD = x, y_ : VectoD = y): QuadXELM_3L1 =
    {
        val z = f0.fM (a * x_)                                             // Z  = f(XA)
        if (DEBUG) println (s"train: x_ = $x_, \n z = $z, \n y_ = $y_")
        val rg = new QuadXRegression (z, y_)                               // QuadXRegression
        rg.train ()
        b = rg.parameter
        this
    } // train

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Build a sub-model that is restricted to the given columns of the data matrix.
     *  @param x_cols  the columns that the new model is restricted to
     */
    def buildModel (x_cols: MatriD): QuadXELM_3L1 =
    {
        new QuadXELM_3L1 (x_cols, y, -1, null, hparam, f0, itran)
    } // buildModel

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given a new input vector 'v', predict the output/response vector 'f(v)'.
     *  @param v  the new input vector
     */
    override def predict (v: VectoD): Double = b dot QuadXRegression.forms (f0.fV (a dot v), v.dim, QuadXRegression.numTerms (v.dim))
//  override def predict (v: VectoD): Double = b dot f0.fV (a dot v)

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given an input matrix 'v', predict the output/response matrix 'f(v)'.
     *  @param v  the input matrix
     */
    override def predict (v: MatriD = x): VectoD = QuadXRegression.allForms (f0.fM (a * v)) * b
//  override def predict (v: MatriD = x): VectoD = f0.fM (a * v) * b

} // QuadXELM_3L1 class


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `QuadXELM_3L1` companion object provides factory functions for buidling three-layer
 *  (one hidden layer) extreme learning machines.
 *  Note, 'rescale' is defined in `ModelFactory` in Model.scala.
 */
object QuadXELM_3L1 extends ModelFactory
{
    private val DEBUG = false                                          // debug flag
    val drp = (-1, null, null, f_tanh)                                 // default remaining parameters

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Create an `QuadXELM_3L1` for a combined data-response matrix.
     *  @param xy         the combined input/data and output/response matrix
     *  @param nz         the number of nodes in hidden layer
     *  @param fname      the feature/variable names
     *  @param hparam     the hyper-parameters
     *  @param f0         the activation function family for layers 1->2 (input to hidden)
     */
    def apply (xy: MatriD, nz: Int = -1,
               fname: Strings = null, hparam: HyperParameter = null,
               f0: AFF = f_tanh): QuadXELM_3L1 =
    {
        var itran: FunctionV_2V = null                                 // inverse transform -> orginal 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 = y

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

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Create an `QuadXELM_3L1` for a data matrix and response vector.
     *  @param x       the input/data matrix
     *  @param y       the output/response vector
     *  @param nz      the number of nodes in hidden layer
     *  @param fname   the feature/variable names
     *  @param hparam  the hyper-parameters
     *  @param f0      the activation function family for layers 1->2 (input to hidden)
     */
    def apply (x: MatriD, y: VectoD, nz: Int,
               fname: Strings, hparam: HyperParameter,
               f0: AFF): QuadXELM_3L1 =
    {
        var itran: FunctionV_2V = null                                 // inverse transform -> orginal scale

        val x_s = if (rescale) rescaleX (x, f0)
                  else x
        val y_s = y

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

} // QuadXELM_3L1 object


//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `QuadXELM_3L1Test` object tests the multi-collinearity method in the
 *  `QuadXELM_3L1` class using the following regression equation on the Blood
 *  Pressure dataset.  It also applies forward selection and backward elimination.
 *  <p>
 *      y = b dot x = b_0 + b_1*x_1 + b_2*x_2 + b_3*x_3 + b_4 * x_4
 *  <p>
 *  @see online.stat.psu.edu/online/development/stat501/12multicollinearity/05multico_vif.html
 *  @see online.stat.psu.edu/online/development/stat501/data/bloodpress.txt
 *  > runMain scalation.analytics.QuadXELM_3L1Test
 */
object QuadXELM_3L1Test extends App                                       // FIX - fails
{
    import ExampleBPressure._
    val x = ox                                                       // use ox for intercept

    println ("model: y = b_0 + b_1*x1 + b_2*x_ + b3*x3 + b4*x42")
    println ("x = " + x)
    println ("y = " + y)

    banner ("Parameter Estimation and Quality of Fit")
    val elm = new QuadXELM_3L1 (x, y, -1, ofname)
    println (elm.analyze ().report)
//  println (elm.summary)

    banner ("Collinearity Diagnostics")
    println ("corr (x) = " + corr (x))                               // correlations of column vectors in x

    banner ("Multi-collinearity Diagnostics")
    println ("vif      = " + elm.vif ())                             // test multi-colinearity (VIF)

    banner ("Forward Selection Test")
    elm.forwardSelAll ()

/*
    banner ("Backward Elimination Test")
    val bcols = Set (0) ++ Array.range (1, x.dim2)                   // start with all x_j in model
    for (l <- 1 until x.dim2  - 1) {
        val (x_j, b_j, fit_j) = elm.backwardElim (bcols)             // eliminate least predictive variable
        println (s"backward model: remove x_j = $x_j with b = $b_j \n fit = $fit_j")
        bcols -= x_j
    } // for
*/

} // QuadXELM_3L1Test object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `QuadXELM_3L1Test2` object trains a extreme learning machine on the `ExampleBasketBall` dataset.
 *  > runMain scalation.analytics.QuadXELM_3L1Test2
 */
object QuadXELM_3L1Test2 extends App
{
    import ExampleBasketBall._
    banner ("QuadXELM_3L1 vs. Regession - ExampleBasketBall")

    println ("x = " + x)
    println ("y = " + y)

    banner ("Regression")
    val rg = Regression (oxy)
    println (rg.analyze ().report)

    banner ("prediction")                                            // not currently rescaling
    val yq = rg.predict ()                                           // scaled predicted output values for all x
    println ("target output:    y  = " + y)
    println ("predicted output: yq = " + yq)
    println ("error:            e  = " + (y - yq))

    banner ("QuadXELM_3L1 with scaled y values")

    val elm = QuadXELM_3L1 (xy)                                      // factory function automatically rescales
//  val elm = new QuadXELM_3L1 (x, y))                               // constructor does not automatically rescale
    println (elm.analyze ().report)

    banner ("prediction")
    val yp = elm.predict ()                                          // scaled predicted output values for all x
    println ("target output:    y  = " + y)
    println ("predicted output: yp = " + yp)
    println ("error:            e  = " + (y - yp))

} // QuadXELM_3L1Test2 object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `QuadXELM_3L1Test3` object trains a extreme learning machine on the `ExampleAutoMPG` dataset.
 *  > runMain scalation.analytics.QuadXELM_3L1Test3
 */
object QuadXELM_3L1Test3 extends App
{
    import ExampleAutoMPG._
    banner ("QuadXELM_3L1 vs. Regession - ExampleAutoMPG")

    banner ("Regression")
    val rg = Regression (oxy)
    println (rg.analyze ().report)

/*
    banner ("prediction")                                            // not currently rescaling
    val yq = rg.predict ()                                           // scaled predicted output values for all x
    println ("target output:    y  = " + y)
    println ("predicted output: yq = " + yq)
    println ("error:            e  = " + (y - yq))
*/

    banner ("QuadXELM_3L1 with scaled y values")

//  val elm = new QuadXELM_3L1 (x, y))                               // constructor does not automatically rescale
    val elm = QuadXELM_3L1 (xy)                                      // factory function automatically rescales
    println (elm.analyze ().report)

    banner ("prediction")
    val yp = elm.predict ()                                          // scaled predicted output values for all x
    println ("target output:    y  = " + y)
    println ("predicted output: yp = " + yp)
    println ("error:            e  = " + (y - yp))

} // QuadXELM_3L1Test3 object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `QuadXELM_3L1Test4` object trains a extreme learning machine on the `ExampleAutoMPG` dataset.
 *  It test cross-validation.
 *  > runMain scalation.analytics.QuadXELM_3L1Test4
 */
object QuadXELM_3L1Test4 extends App
{
    import ExampleAutoMPG._
    banner ("QuadXELM_3L1 cross-validation - ExampleAutoMPG")

    banner ("QuadXELM_3L1 with scaled y values")

    var elm: QuadXELM_3L1 = null
    for (nz <- 11 to 31 by 2) {
        elm = QuadXELM_3L1 (xy, nz = nz)                             // factory function automatically rescales
//      elm = new QuadXELM_3L1 (x, Seq (y))                          // constructor does not automatically rescale
        println (elm.analyze ().report)
    } // for

//  banner ("cross-validation")
//  elm.crossVal ()

} // QuadXELM_3L1Test4 object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `QuadXELM_3L1Test5` object trains a extreme learning machine on the `ExampleAutoMPG` dataset.
 *  This tests forward feature/variable selection.
 *  FIX (1) missing intercept/bias, (2) R^2 cv too high
 *  > runMain scalation.analytics.QuadXELM_3L1Test5
 */
object QuadXELM_3L1Test5 extends App
{
    import ExampleAutoMPG._
    val n   = x.dim2                                                 // number of parameters/variables
    val nz  = 17                                                     // number of nodes in hidden layer
    banner ("QuadXELM_3L1 feature selection - ExampleAutoMPG")

    banner ("QuadXELM_3L1 with scaled y values")
    val elm = QuadXELM_3L1 (xy, nz = nz)                             // factory function automatically rescales
//  val elm = new QuadXELM_3L1 (x, y)                                // constructor does not automatically rescale
    println (elm.analyze ().report)
    val ft = elm.fit                                                 // quality of fit for first output

    banner ("Forward Selection Test")
    val (cols, rSq) = elm.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")

} // QuadXELM_3L1Test5 object