//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller, Hao Peng
 *  @version 1.5
 *  @date    Fri Mar 16 15:13:38 EDT 2018
 *  @see     LICENSE (MIT style license file).
 *
 *  @see     hebb.mit.edu/courses/9.641/2002/lectures/lecture03.pdf
 */

package scalation.analytics

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

import ActivationFun._
import Optimizer._

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `NeuralNet_3L` class supports multi-output, 3-layer (input, hidden and output)
 *  Neural-Networks.  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 weights/parameters 'aa' and 'bb' connecting the layers,
 *  so that for a new input vector 'v', the net can predict the output value, i.e.,
 *  <p>
 *      yp = f2 (bb * f1V (aa * v))
 *  <p>
 *  where 'f1' and 'f2' are the activation functions and the parameter matrices
 *  'aa' and 'bb' gives the  weights between input-hidden and hidden-output layers.
 *  Note, if 'a0' is to be treated as bias/intercept, 'x0' must be 1.0.
 *  @param x          the m-by-nx input matrix (training data consisting of m input vectors)
 *  @param y          the m-by-ny output matrix (training data consisting of m output vectors)
 *  @param nz         the number of nodes in hidden layer
 *  @param eta_       the learning/convergence rate
 *  @param bSize      the mini-batch size
 *  @param maxEpochs  the maximum number of training epochs/iterations
 *  @param f1         the activation function family for layers 1->2 (input to hidden)
 *  @param f2         the activation function family for layers 2->3 (hidden to output)
 */
class NeuralNet_3L (x: MatriD, y: MatriD,
                    private var nz: Int = -1,
                    eta_ : Double  = hp ("eta"),
                    bSize: Int     = hp ("bSize").toInt,
                    maxEpochs: Int = hp ("maxEpochs").toInt,
                    f1: AFF = f_sigmoid,
                    f2: AFF = f_sigmoid)
      extends NeuralNet (x, y, eta_)                                  // sets eta in parent class
{
    private val DEBUG      = false                                    // debug flag
    private var aa: MatriD = null                                     // weights/parameters (in to hid) - new MatrixD (nx, nz)
    private var bb: MatriD = null                                     // weights/parameters (hid to out) - new MatrixD (nz, ny)

    // Guidelines for setting the number of nodes in hidden layer, e.g.,
    // 2 nx + 1, nx + 1, (nx + ny) / 2, sqrt (nx ny)
    if (nz < 1) nz = nx + 1

    println (s"Create a NeuralNet_3L with $nx input, $nz hidden and $ny output nodes")

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the weight matrices 'aa' and 'bb'.
     */
    def weights: Array [MatriD] = Array (aa, bb)

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given training data 'x' and 'y', fit the parameter/weight matrices 'aa' and 'bb'.
     *  Iterate over several epochs (no batching).
     */
    def train0 (): NeuralNet_3L =
    {
        println (s"train0: eta = $eta")
        if (aa == null) aa = weightMat (nx, nz)                         // initialize parameter/weight matrix aa
        if (bb == null) bb = weightMat (nz, ny)                         // initialize parameter/weight matrix bb
        var sse0 = Double.MaxValue                                      // hold prior value of sse

        for (epoch <- 1 to maxEpochs) {                                 // iterate over each epoch
            var z  = f1.fM (x * aa); /* z.setCol (0, _1) */             // Z  = f (X A) (and opt. set hidden bias)
            var yp = f2.fM (z * bb)                                     // Yp = f (Z B)
            var ee = yp - y                                             // -E where E is the error matrix
            val dy = f2.dM (yp) ** ee                                   // delta y
            val dz = f1.dM (z) ** (dy * bb.t)                           // delta z
            val (aup, bup) = (x.t * dz * eta, z.t * dy * eta)           // changes to weights
            aa -= aup; bb -= bup                                        // update weight matrices

            val sse = sseF (y, f2.fM (f1.fM (x * aa) * bb))
            if (DEBUG) println (s"train0: weights for $epoch th epoch: sse = $sse")
            if (sse > sse0) return this                                 // return early if moving up
            sse0 = sse                                                  // save prior sse
            this
        } // for
        this
    } // train0

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given training data 'x' and 'y', fit the parameter/weight matrices 'aa' and 'bb'.
     *  Iterate over several epochs, where each epoch divides the training set into
     *  'nbat' batches.  Each batch is used to update the weights. 
     */
    def train (): NeuralNet_3L =
    {
        if (aa == null) aa = weightMat (nx, nz)                                  // initialize parameter/weight matrix aa
        if (bb == null) bb = weightMat (nz, ny)                                  // initialize parameter/weight matrix bb
        val epochs = optimize3 (x, y, aa, bb, eta, bSize, maxEpochs, f1, f2)     // optimize parameter/weight matrices aa, bb
        println (s"ending epoch = $epochs")
        this
    } // train

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given training data 'x' and 'y', fit the parameter/weight matrices 'aa' and 'bb'.
     *  Iterate over several epochs, where each epoch divides the training set into
     *  'nbat' batches.  Each batch is used to update the weights. 
     *  This version preforms an interval search for the best 'eta' value. 
     */
    def train2 (): NeuralNet_3L =
    {
        if (aa == null) aa = weightMat (nx, nz)                                  // initialize parameter/weight matrix aa
        if (bb == null) bb = weightMat (nz, ny)                                  // initialize parameter/weight matrix bb
        val etaI = (0.25 * eta, 4.0 * eta)                                       // quarter to four times eta
        val epochs = optimize3I (x, y, aa, bb, etaI, bSize, maxEpochs, f1, f2)   // optimize parameter/weight matrices aa, bb
        println (s"ending epoch = $epochs")
        this
    } // train2

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given a new input vector 'v', predict the output/response vector 'f(v)'.
     *  @param v  the new input vector
     */
    def predictV (v: VectoD): VectoD = f2.fV (bb dot f1.fV (aa dot v))

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given an input matrix 'x', predict the output/response matrix 'f(x)'.
     *  @param x  the input matrix
     */
    def predict (x: MatriD): MatriD = f2.fM (f1.fM (x * aa) * bb)

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Perform 'k'-fold cross-validation.
     *  @param k      the number of folds
     *  @param rando  whether to use randomized cross-validation
     */
    def crossVal (k: Int = 10, rando: Boolean = true)
    {
        crossValidate ((x: MatriD, y: MatriD) => new NeuralNet_3L (x, y), k, rando) 
    } // crossVal

} // NeuralNet_3L class


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `NeuralNet_3LTest` object is used to test the `NeuralNet_3L` class.
 *  @see www4.rgu.ac.uk/files/chapter3%20-%20bp.pdf
 *  > runMain scalation.analytics.NeuralNet_3LTest
 */
object NeuralNet_3LTest extends App
{
    val s = 0                                          // random number stream to use
    val x = new MatrixD ((3, 3), 1.0, 0.35, 0.9,       // training data - input matrix (m vectors)
                                 1.0, 0.20, 0.7,
                                 1.0, 0.40, 0.95)
    val y = new MatrixD ((3, 2), 0.5, 0.4,             // training data - output matrix (m vectors)
                                 0.3, 0.3,
                                 0.6, 0.5)

    println ("input  matrix x = " + x)
    println ("output matrix y = " + y)

    val nn = new NeuralNet_3L (x, y, 3, bSize = 1)      // create a NeuralNet_3L - FIX

//  for (eta <- 0.5 to 10.0 by 0.5) {
    for (i <- 1 to 20) {
        val eta = i * 0.5
        banner (s"NeuralNet_3LTest: Fit the parameter matrix bb using optimization with learning rate $eta")

        nn.reset (eta)
        nn.train ().eval ()                             // fit the weights using training data
        println ("(aa, bb) = " + nn.weights.deep)
        nn.fitMap ()

//      yp = nn.predict (x)                             // predicted output values
//      println ("target output:    y   = " + y)
//      println ("predicted output: yp  = " + yp)
        println ("yp0 = " + nn.predict (x(0)))          // predicted output values for row 0
    } // for

    banner ("NeuralNet_3LTest: Compare with Linear Regression - first column of y")

    val y0  = y.col(0)                                  // use first column of matrix y
    val rg0 = new Regression (x, y0)                    // create a Regression model
    rg0.train ().eval ()
    println ("b      = " + rg0.parameter)
    println ("fitMap = " + rg0.fitMap)

    val y0p = rg0.predict (x)                           // predicted output value
    println ("target output:    y0  = " + y0)
    println ("predicted output: y0p = " + y0p)

    banner ("NeuralNet_3LTest: Compare with Linear Regression - second column of y")

    val y1 = y.col(1)                                   // use second column of matrix y
    val rg1 = new Regression (x, y1)                    // create a Regression model
    rg1.train ().eval ()
    println ("b      = " + rg1.parameter)
    println ("fitMap = " + rg1.fitMap)

    val y1p = rg1.predict (x)                           // predicted output value
    println ("target output:    y1  = " + y1)
    println ("predicted output: y1p = " + y1p)

} // NeuralNet_3LTest object