//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller
 *  @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_2L` class supports multi-output, 2-layer (input 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 'bb' connecting the layers,
 *  so that for a new input vector 'z', the net can predict the output value, i.e.,
 *  <p>
 *      yp_j = f (bb dot z)
 *  <p>
 *  where 'f' is the activation function and the parameter matrix 'bb' gives the
 *  weights between input and output layers.
 *  No batching is used for this algorithm.
 *  Note, 'b0' is treated as the bias, so '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 eta_       the learning/convergence rate
 *  @param bSize      the batch size
 *  @param maxEpochs  the maximum number of training epochs/iterations
 *  @param f1         the activation function family for layers 1->2 (input to output)
 */
class NeuralNet_2L (x: MatriD, y: MatriD,
                    eta_ : Double  = hp ("eta"),
                    bSize: Int     = hp ("bSize").toInt,
                    maxEpochs: Int = hp ("maxEpochs").toInt,
                    f1: AFF = f_sigmoid)
      extends NeuralNet (x, y, eta_)                                     // sets eta in parent class
{
    private val DEBUG   = false                                          // debug flag

    private var bb: MatriD = null                                        // parameters/weights - new MatrixD (nx, ny)
    private var yp: MatriD = null                                        // prediction matrix
    private var ee: MatriD = null                                        // error matrix

    println (s"Create a NeuralNet_2L with $nx input and $ny output nodes")

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the weight matrix 'bb' (array of 1).
     */
    def weights: Array [MatriD] = Array (bb)

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given training data 'x' and 'y', fit the parameter/weight matrix 'bb'.
     *  Minimize the error in the prediction by adjusting the weight matrix 'bb'.
     *  The error 'ee' is simply the difference between the target value 'y' and the
     *  predicted value 'yp'.  Minimize the dot product of error with itself using
     *  gradient-descent. specifically move in the opposite direction of the gradient.
     *  Iterate over several epochs (no batching).
     */
    def train0 (): NeuralNet_2L =
    {
        println (s"train0: eta = $eta")
        if (bb == null) bb = weightMat (nx, ny)                         // initialize parameters/weights
        var sse0 = Double.MaxValue                                      // hold prior value of sse

        for (epoch <- 1 to maxEpochs) {                                 // iterate over each epoch

            val yp  = f1.fM (x * bb)                                    // Yp = f (X B)
            val ee  = yp - y                                            // -E where E is the error matrix
            val dy  = f1.dM (yp) * ee                                   // delta y
            val bup = x.t * dy * eta                                    // change in input-output weights
            bb -= bup                                                   // update weights

            val sse = sseF (y, f1.fM (x * bb))                          // recompute sum of squared errors
            if (DEBUG) println (s"train0: weights for $epoch th epoch: sse = $sse")
            if (sse > sse0) return this                                 // return when sse increases
            sse0 = sse                                                  // save prior sse
        } // for
        this
    } // train0

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given training data 'x' and 'y', fit the parameter/weight matrix 'bb'.
     *  Minimize the error in the prediction by adjusting the weight matrix '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_2L =
    {
        if (bb == null) bb = weightMat (nx, ny)                          // initialize parameter/weight matrix bb
        val epochs = optimize2 (x, y, bb, eta, bSize, maxEpochs, f1)     // optimize parameter/weight matrix bb
        println (s"ending epoch = $epochs")
        this
    } // train

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given training data 'x' and 'y', fit the parameter/weight matrix 'bb'.
     *  Minimize the error in the prediction by adjusting the weight matrix '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_2L =
    {
        if (bb == null) bb = weightMat (nx, ny)                          // initialize parameter/weight matrix bb
        val etaI = (0.25 * eta, 4.0 * eta)                               // quarter to four times eta
        val epochs = optimize2I (x, y, bb, etaI, bSize, maxEpochs, f1)   // optimize parameter/weight matrix bb
        println (s"ending epoch = $epochs")
        this
    } // train2

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

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given an input matrix 'z', predict the output/response matrix 'f(z)'.
     *  @param z  the input matrix
     */
    def predict (z: MatriD): MatriD = f1.fM (z * 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_2L (x, y), k, rando)
    } // crossVal

} // NeuralNet_2L class


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `NeuralNet_2LTest` object is used to test the `NeuralNet_2L` class.  For this
 *  test, training data is used to fit the weights before using them for prediction.
 *  @see www4.rgu.ac.uk/files/chapter3%20-%20bp.pdf
 *  > runMain scalation.analytics.NeuralNet_2LTest
 */
object NeuralNet_2LTest 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_2L (x, y)                  // create a NeuralNet_2L

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

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

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

    banner ("NeuralNet_2LTest: 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_2LTest: 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_2LTest object