//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller
 *  @version 1.5
 *  @date    Mon Sep  9 13:30:41 EDT 2013
 *  @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, VectorD}
import scalation.util.banner

import ActivationFun._
import MatrixTransform._
import Optimizer._

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `Perceptron` class supports single-output, 2-layer (input and output)
 *  Neural-Networks.  Although perceptrons are typically used for classification,
 *  this class is used for prediction.  Given several input vectors and output
 *  values (training data), fit the weights/parameters 'b' connecting the layers,
 *  so that for a new input vector 'z', the net can predict the output value, i.e.,
 *  <p>
 *      z = f (b dot z)
 *  <p>
 *  The parameter vector 'b' (w) gives the weights between input and output layers.
 *  Note, 'b0' is treated as the bias, so 'x0' must be 1.0.
 *  @param x       the input m-by-n matrix (training data consisting of m input vectors)
 *  @param y       the output m-vector (training data consisting of m output values)
 *  @param fname_  the feature/variable names
 *  @param hparam  the hyper-parameters
 *  @param f1      the activation function family for layers 1->2 (input to output)
 */
class Perceptron (x: MatriD, y: VectoD, fname_ : Strings = null,
                  hparam: HyperParameter = Optimizer.hp,
                  f1: AFF = f_sigmoid)
      extends PredictorMat (x, y, fname_, hparam)
{
    private val DEBUG     = false                                    // debug flag
    private var eta       = hparam ("eta")                           // the learning/convergence rate (requires adjustment)
    private var bSize     = hparam ("bSize").toInt                   // the batch size
    private val maxEpochs = hparam ("maxEpochs").toInt               // the maximum number of training epcochs/iterations
    private val n         = x.dim2                                   // dimensionality of the input
    private val _1        = VectorD.one (m)                          // vector of all ones

    if (y.dim != m) flaw ("constructor", "dimensions of x and y are incompatible")

    println (s"Create a Perceptron with $n input nodes and 1 output node")

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Set the initial parameter/weight vector 'b' manually before training.
     *  This is mainly for testing purposes.
     *  @param w0  the initial weights for b
     */
    def setWeights (w0: VectoD)
    {
        b = w0
        if (DEBUG) {
            eval ()
            println ("b      = " + parameter)
            println ("fitMap = " + fitMap)
        } // if
    } // setWeights

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Reset the learning rate 'eta'.
     *  @param eta  the learning rate
     */
    def reset (eta_ : Double) { eta = eta_ }

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given training data 'x' and 'yy', fit the parameter/weight vector 'b'.
     *  Minimize the error in the prediction by adjusting the weight vector 'b'.
     *  The error 'e' is simply the difference between the target value 'yy' and the
     *  predicted value 'yp'.  Minimize the dot product of error with itself using
     *  gradient-descent (move in the opposite direction of the gradient).
     *  Iterate over several epochs (no batching).
     *  @param yy  the output vector
     */
    def train0 (yy: VectoD = y): Perceptron =
    {
        println (s"train0: eta = $eta")
        if (b == null) b = weightVec (n)                         // initialize parameters/weights
        var sse0 = Double.MaxValue

        for (epoch <- 1 until maxEpochs) {                       // epoch-th learning phase
            val yp  = f1.fV (x * b)                              // vector of predicted outputs Yp = f (X b)
            val e   = yp - yy                                    // -e where e is the error vector
//          val dy  = yp * (_1 - yp) * e                         // delta y (for sigmoid only)
            val dy  = f1.dV (yp) * e                             // delta y
            val bup = x.t * dy * eta                             // change in parameters/weights
            b      -= bup                                        // update the parameters/weights

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

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given training data 'x' and 'yy', fit the parameter/weight vector 'b'.
     *  Minimize the error in the prediction by adjusting the weight vector 'b'.
     *  Iterate over several epochs, where each epoch divides the training set into
     *  'nbat' batches.  Each batch is used to update the weights.
     *  @param yy  the output vector
     */
    def train (yy: VectoD = y): Perceptron =
    {
        if (yy.dim < 2 * bSize) flaw ("train", "not enough data for batching - use 'train0'")
        println (s"train: eta = $eta")
        if (b == null) b = weightVec (n)                         // initialize parameters/weights
        val result = optimize (x, y, b, eta, bSize, maxEpochs, f1)
        println (s"result = (sse, ending_epoch_) = $result")
        this
    } // train

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given training data 'x' and 'yy', fit the parameter/weight vector 'b'.
     *  Minimize the error in the prediction by adjusting the weight vector 'b'.
     *  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.
     *  @param yy  the output vector
     */
    override def train2 (yy: VectoD = y): Perceptron =
    {
        if (yy.dim < 2 * bSize) flaw ("train2", "not enough data for batching - use 'train0'")
        val etaI = (0.25 * eta, 4.0 * eta)                       // quarter to four times eta
        println (s"train2: etaI = $etaI")
        if (b == null) b = weightVec (n)                         // initialize parameters/weights
        val result = optimizeI (x, y, b, etaI, bSize, maxEpochs, f1)
        println (s"result = (sse, ending_epoch_) = $result")
        this
    } // train2

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Evaluate the quality of fit for the current parameter/weight vector 'b'.
     */
    override def eval ()
    { 
        val yp = f1.fV (x * b)                                   // predicted output yp
        e = y - yp                                               // difference between actual and predicted output
        diagnose (e)                                             // compute diagonostics
    } // eval

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given a new input vector 'z', predict the output/response value 'f(z)'.
     *  @param z  the new input vector
     */
    override def predict (z: VectoD): Double = f1.f (b dot z)

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given a new input matrix 'z', predict the output/response value 'f(z)'.
     *  @param z  the new input matrix
     */
    override def predict (z: MatriD = x): VectoD = f1.fV (z * b)

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** 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: VectoD) => new Perceptron (x, y), k, rando)
    } // crossVal

} // Perceptron class


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `Perceptron` companion object provides factory methods for buidling perceptrons.
 */
object Perceptron
{
    import PredictorMat.pullResponse

    private val DEBUG = false                                    // debug flag

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Create a `Perceptron` for a combined data matrix.
     *  @param xy       the combined input and output matrix
     *  @param fname    the feature/variable names
     *  @param hparam   the hyper-parameters
     *  @param f1       the activation function family for layers 1->2 (input to output)
     *  @param rescale  the whether to rescale output y to match activation function
     */
    def apply (xy: MatriD, fname: Strings = null,
               hparam: HyperParameter = Optimizer.hp,
               f1: AFF = f_sigmoid, rescale: Boolean = true): Perceptron =
    {
        val xy_s = if (rescale && f1.bounds != null) {
                       val (min_xy, max_xy) = (min (xy), max (xy))
                       scale (xy, min_xy, max_xy, f1.bounds)
                   } else {
                       xy
                   } // if
        val (x, y) = pullResponse (xy_s)
        setCol2One (x)
        if (DEBUG) println ("scaled: x = " + x + "\nscaled y = " + y)
        new Perceptron (x, y, fname, hparam, f1)
    } // apply

} // Perceptron object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `PerceptronTest` object is used to test the `Perceptron` 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.PerceptronTest
 */
object PerceptronTest 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 = VectorD (0.5, 0.30, 0.6)                 // training data - output vector

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

    val nn = new Perceptron (x, y)                   // create a Perceptron

    banner ("PerceptronTest: Set the parameter vector b manually")

    val b = VectorD (0.0, 0.5, 0.5)                  // set weight vector b manually
    nn.setWeights (b)
    println ("b      = " + nn.parameter)
    println ("fitMap = " + nn.fitMap)

    var yp = nn.predict ()                           // predicted output value
    println ("target output:    y   = " + y)
    println ("predicted output: yp  = " + yp)

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

        nn.reset (eta)
        nn.train0 ().eval ()                         // fit the weights using training data
        println ("b      = " + nn.parameter)
        println ("fitMap = " + nn.fitMap)

        yp = nn.predict ()                           // predicted output value
        println ("target output:    y   = " + y)
        println ("predicted output: yp  = " + yp)
    } // for

    banner ("PerceptronTest: Compare with Linear Regression")

    val rg = new Regression (x, y)                  // create a Regression model
    rg.train ().eval ()
    println ("b      = " + rg.parameter)
    println ("fitMap = " + rg.fitMap)

    yp = rg.predict ()                              // predicted output value
    println ("target output:    y   = " + y)
    println ("predicted output: yp  = " + yp)

} // PerceptronTest object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `PerceptronTest2` object trains a perceptron on a small dataset of
 *  temperatures from counties in Texas where the variables/factors to consider
 *  are Latitude (x1), Elevation (x2) and Longitude (x3).  The regression equation
 *  is the following:
 *  <p>
 *      y  =  sigmoid (w dot x)  =  sigmoid (w0 + w1*x1 + w2*x2 + w3*x3)
 *  <p>
 *  This test case illustrates how to transform the columns of the matrix
 *  so that the 'sigmoid' activation function can work effectively.
 *  Since the dataset is very small, should use 'train0' which does no batching.
 *  > runMain scalation.analytics.PerceptronTest2
 */
object PerceptronTest2 extends App
{
    // 16 data points:        Constant      x1      x2       x3     y
    //                                     Lat    Elev     Long  Temp        County
    val xy = new  MatrixD ((16, 5), 1.0, 29.767,   41.0,  95.367, 56.0,    // Harris
                                    1.0, 32.850,  440.0,  96.850, 48.0,    // Dallas
                                    1.0, 26.933,   25.0,  97.800, 60.0,    // Kennedy
                                    1.0, 31.950, 2851.0, 102.183, 46.0,    // Midland
                                    1.0, 34.800, 3840.0, 102.467, 38.0,    // Deaf Smith
                                    1.0, 33.450, 1461.0,  99.633, 46.0,    // Knox
                                    1.0, 28.700,  815.0, 100.483, 53.0,    // Maverick
                                    1.0, 32.450, 2380.0, 100.533, 46.0,    // Nolan
                                    1.0, 31.800, 3918.0, 106.400, 44.0,    // El Paso
                                    1.0, 34.850, 2040.0, 100.217, 41.0,    // Collington
                                    1.0, 30.867, 3000.0, 102.900, 47.0,    // Pecos
                                    1.0, 36.350, 3693.0, 102.083, 36.0,    // Sherman
                                    1.0, 30.300,  597.0,  97.700, 52.0,    // Travis
                                    1.0, 26.900,  315.0,  99.283, 60.0,    // Zapata
                                    1.0, 28.450,  459.0,  99.217, 56.0,    // Lasalle
                                    1.0, 25.900,   19.0,  97.433, 62.0)    // Cameron

    val (x, y) = PredictorMat.pullResponse (xy)
    println ("x = " + x)

    banner ("Perceptron with scaled y values")
    val nn  = Perceptron (xy)                        // factory function automatically rescales
//  val nn  = new Perceptron (x, y)                  // constructor does not automatically rescale
    val eta = 0.5                                    // try several values
    nn.reset (eta)
    nn.train0 ().eval ()                             // fit the weights using training data

    println ("b      = " + nn.parameter)
    println ("fitMap = " + nn.fitMap)

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

    banner ("unscaled prediction")
    val (ymin, ymax) = (y.min (), y.max ())
    val ypu = unscaleV (yp, ymin, ymax, (0, 1))      // unscaled predicted output values for all x
    println ("target output:   y   = " + y)
    println ("unscaled output: ypu = " + ypu)

} // PerceptronTest2 object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `PerceptronTest3` object trains a perceptron on a small dataset with variables
 *  x1 and x2.  The regression equation is the following:
 *  <p>
 *      y  =  sigmoid (b dot x)  =  sigmoid (b0 + b1*x1 + b2*x2)
 *  <p>
 *  Does not call the 'train' method; improvements steps for sigmoid are explicitly in code below.
 *  > runMain scalation.analytics.PerceptronTest3
 */
object PerceptronTest3 extends App
{
    import PredictorMat.pullResponse

    // 9 data points:        Constant    x1    x2     y
    val xy = new MatrixD ((9, 4), 1.0,  0.0,  0.0,  0.5,
                                  1.0,  0.0,  0.5,  0.3,
                                  1.0,  0.0,  1.0,  0.2,

                                  1.0,  0.5,  0.0,  0.8,
                                  1.0,  0.5,  0.5,  0.5,
                                  1.0,  0.5,  1.0,  0.3,

                                  1.0,  1.0,  0.0,  1.0,
                                  1.0,  1.0,  0.5,  0.8,
                                  1.0,  1.0,  1.0,  0.5)

    private val _1 = VectorD.one (xy.dim1)           // vector of all ones

    println (s"xy = $xy")
    val (x, y) = pullResponse (xy)                   // input matrix, output vector
    val sst = Fit.sstf (y)                           // sum of squares total

    val b = VectorD (0.1, 0.2, 0.1)                  // initial weights/parameters

    val eta = 2.0                                    // try several values
    val hp = Optimizer.hp.updateReturn ("eta", eta)
    val nn = Perceptron (xy, null, hp)               // create a perceptron
    banner ("initialize")

    nn.setWeights (b)
    var yp  = nn.predict ()
    var e   = y - yp
    var sse = e dot e
    println (s"b   = $b")
    println (s"y   = $y")
    println (s"yp  = $yp")
    println (s"e   = $e")
    println (s"sse = $sse")
    println (s"R^2 = ${1 - sse/sst}")
 
    for (epoch <- 1 to 10) {
    banner (s"improvement step $epoch")
        val bup = x.t * (e * yp * (_1 - yp)) * eta   // adjust the parameters/weights (for sigmoid)
        b += bup
        nn.setWeights (b)
        yp  = nn.predict ()
        e   = y - yp
        sse = e dot e
        println (s"bup = $bup")
        println (s"b   = $b")
        println (s"y   = $y")
        println (s"yp  = $yp")
        println (s"e   = $e")
        println (s"sse = $sse")
        println (s"R^2 = ${1 - sse/sst}")
    } // for
 
} // PerceptronTest3 object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `PerceptronTest4` object trains a perceptron on a small dataset with variables
 *  x1 and x2.  The regression equation is the following:
 *  <p>
 *      y  =  sigmoid (b dot x)  =  sigmoid (b0 + b1*x1 + b2*x2)
 *  <p>
 *  This version calls the 'train0' method.
 *  > runMain scalation.analytics.PerceptronTest4
 */
object PerceptronTest4 extends App
{
    // 16 data points:        Constant    x1    x2     y
    val xy = new  MatrixD ((9, 4), 1.0,  0.0,  0.0,  0.5,
                                   1.0,  0.0,  0.5,  0.3,
                                   1.0,  0.0,  1.0,  0.2,

                                   1.0,  0.5,  0.0,  0.8,
                                   1.0,  0.5,  0.5,  0.5,
                                   1.0,  0.5,  1.0,  0.3,

                                   1.0,  1.0,  0.0,  1.0,
                                   1.0,  1.0,  0.5,  0.8,
                                   1.0,  1.0,  1.0,  0.5)

    println (s"xy = $xy")

    val eta = 2.0                                    // try several values
    val hp = Optimizer.hp.updateReturn ("eta", eta)
    val nn = Perceptron (xy, null, hp)               // create a perceptron

    nn.train0 ().eval ()                             // fit the weights using training data

    println ("b      = " + nn.parameter)
    println ("fitMap = " + nn.fitMap)

} // PerceptronTest4 object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `PerceptronTest5` object trains a perceptron on a small dataset with variables
 *  > runMain scalation.analytics.PerceptronTest5
 */
object PerceptronTest5 extends App
{
    val x = new MatrixD ((6, 2), 1, 1, 1, 2, 1, 3, 1, 4, 1, 5, 1, 6)
    val y = VectorD (0, 0, 1, 0, 1, 1)

    val eta = 0.3                                    // try several values
    val hp = Optimizer.hp.updateReturn ("eta", eta)
    val nn = new Perceptron (x, y, null, hp)         // create a perceptron

    nn.setWeights (VectorD (-3.1, 0.9))              // start search from this point

    nn.train0 ().eval ()                             // fit the weights using training data

    println ("b      = " + nn.parameter)
    println ("fitMap = " + nn.fitMap)

    val yp = nn.predict ()                           // predicted output value
    val ypr = yp.map (math.round (_))
    println ("target output:    y   = " + y)
    println ("predicted output: ypr = " + ypr)
    println ("predicted output: yp  = " + yp)

} // PerceptronTest5 object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `PerceptronTest6` object trains a perceptron on a small dataset with 4 variables.
 *  > runMain scalation.analytics.PerceptronTest6
 */
object PerceptronTest6 extends App
{
    val xy = new MatrixD ((9, 4), 1, 0, 0, 0,
                                  1, 0, 1, 0,
                                  1, 0, 2, 0,
                                  1, 1, 0, 0,
                                  1, 1, 1, 1,
                                  1, 1, 2, 1,
                                  1, 2, 0, 0,
                                  1, 2, 1, 1,
                                  1, 2, 2, 1)

    val eta = 0.1                                    // try several values
    val hp = Optimizer.hp.updateReturn ("eta", eta)
    val nn = Perceptron (xy, null, hp)               // create a perceptron

//  nn.setWeights (VectorD (-3.1, 0.9))              // start search from this point

    nn.train0 ().eval ()                             // fit the weights using training data

    println ("b      = " + nn.parameter)
    println ("fitMap = " + nn.fitMap)

    val (x, y) = PredictorMat.pullResponse (xy)
    val yp = nn.predict ()                           // predicted output value
    val ypr = yp.map (math.round (_))
    println ("target output:    y   = " + y)
    println ("predicted output: ypr = " + ypr)
    println ("predicted output: yp  = " + yp)

} // PerceptronTest6 object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `PerceptronTest7` object trains a perceptron on the `ExampleBasketBall dataset.
 *  > runMain scalation.analytics.PerceptronTest7
 */
object PerceptronTest7 extends App
{
    import ExampleBasketBall._

    println ("ox = " + ox)
    println ("y  = " + y)

    banner ("Regression")
    val rg = Regression (oxy)
    rg.train ().eval ()

    println ("b      = " + rg.parameter)
    println ("fitMap = " + rg.fitMap)

    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 ("Perceptron with scaled y values")
//  val nn  = Perceptron (oxy, f1 = f_id)            // factory function automatically rescales
    val nn  = new Perceptron (ox, y, f1 = f_id)      // constructor does not automatically rescale
    val eta = 0.000001                               // try several values
    nn.reset (eta)
    nn.train ().eval ()                              // fit the weights using training data

    println ("b      = " + nn.parameter)
    println ("fitMap = " + nn.fitMap)

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

    banner ("unscaled prediction")
    val (ymin, ymax) = (y.min (), y.max ())
    val ypu = unscaleV (yp, ymin, ymax, (0, 1))      // unscaled predicted output values for all x
    println ("target output:   y   = " + y)
    println ("unscaled output: ypu = " + ypu)

} // PerceptronTest7 object