//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller, Hao Peng
 *  @version 1.6
 *  @date    Thu Jun 13 13:13:26 EDT 2019
 *  @see     LICENSE (MIT style license file).
 *
 *  @title   Model: Simple Exponential Smoothing
 *
 *  @see https://otexts.com/fpp2/ses.html
 */

package scalation.analytics
package forecaster

import scala.collection.mutable.Set

import scalation.linalgebra.{MatriD, MatrixD, VectoD, VectorD}
import scalation.minima.L_BFGS_B
import scalation.plot.Plot
import scalation.random._
import scalation.stat.Statistic
import scalation.util.banner

import Fit._
import RollingValidation.trSize

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `SimpleExpSmoothing` class provide very basic time series analysis using
 *  Simple Exponential Smoothing models.  The forecasted value is the weighted
 *  average the latest value and the previous smoothed value.  The smoothing parameter
 *  'α in [0, 1]' causes the contributions of older values to decay exponentially.
 *  @see Smoothing Equation in section 7.1.
 *  <p>
 *      s_t = α y_t-1 + (1 - α) s_t-1                          // smoothing equation
 *      yf_t = s_t                                             // forecast equation
 *  <p>
 *  where vector 's' is the smoothed version of vector 'y'.
 *  @param y       the response vector (orginal time series data)
 *  @param hparam  the hyper-parameters
 */
class SimpleExpSmoothing (y: VectoD, hparam: HyperParameter = SimpleExpSmoothing.hp)
      extends ForecasterVec (y, n = 1, hparam) with NoFeatureSelection
{
    private val DEBUG     = false                                           // debug flag
    private val TOL       = 1E-4                                            // tolerance
    private val lo        = VectorD.fill (1)(0.0)                           // lower bound on α for optimizer
    private val up        = VectorD.fill (1)(1.5)                           // upper bound on α for optimizer

    private var α         = hparam ("α")                                    // default value for the smoothing parameter
    private var s: VectoD = null                                            // vector of smoothed/leveled values (state)
    private var opt       = true                                            // whehtherf to optimize the smoothing parameter

    modelName = "SimpleExpSmoothing"                                        // name of the modeling technique

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Reset the smoothing parameter α.
     *  @param a  the smoothing parameter
     */
    def reset (a: Double) { α = a }

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Toggle the 'opt' flag that indicates whether optimization should be used to set α.
     */
    def toggleOpt () { opt = ! opt }

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Smooth the time-series data 'y', returning the leveled/smoothed data 'l'.
     *  May be viewed as unoptimized training.
     *  @see Smoothing Equation in section 7.1.
     *  <p>
     *      s_t = α y_t-1 + (1 - α) s_t-1                                   // smoothing equation
     *  <p>
     *  @param a   the smoothing parameter (decay rate for older values)
     *  @param y_  the response/output vector (training/full)
     */
    def smooth (a: Double = α, y_ : VectoD = y): VectoD =
    {
        super.train (null, y_)
        s = new VectorD (y_.dim)
        s(0) = y(0)
        for (t <- 1 until y_.dim) s(t) = a * y_(t-1) + (1 - a) * s(t-1)
        s
    } // smooth

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Train the `SimpleExpSmoothing` model on the time-series data, by finding the value
     *  for the smoothing parameter 'α' that minimizes the sum of squared errors 'sse'.
     *  @param x_null  the data/input matrix (ignored)
     *  @param y_      the response/output vector (training/full)
     */
    override def train (x_null: MatriD, y_ : VectoD): SimpleExpSmoothing =
    {
        //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
        /*  The objective function to be minimized (sum of squared errors) by `L_BFGS_B`.
         *  @param x  the input vector 'VectorD (α)' to be optimized
         */
        def f_obj (x: VectoD): Double = (y_ - smooth (x(0), y_)).normSq     // only one parameter

        super.train (null, y_)
        if (opt) {
            val solver = new L_BFGS_B (f_obj, l = lo, u = up)               // Quasi-Newton optimizer
            val x = solver.solve (VectorD (α), toler = TOL)                 // optimize value for α
            α = x(0)                                                        // pull α from vector result
        } // if
        s = smooth (α)                                                      // vector of smoothed/predicted values, with optimized α
        if (DEBUG) println (s"train: α = $α")
        this
    } // train

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the parameter vector.
     */
    def parameter: VectoD = VectorD (α)

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the vector of predicted values on the training data.
     *  Must call 'smooth' or 'train' first.
     */
    override def predictAll (): VectoD = s

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the vector of centered predicted values on the training data.
     *  Must call 'smooth' or 'train' first.
     */
    def predictAllz (): VectoD = s - stats.mu

    def forecastAll (h: Int): MatriD = ???
    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Forecast values for all 'm' time points and all horizons (1 through 'h'-steps ahead).
     *  Record these in the 'yf' matrix, where
     *  <p>
     *      yf(t, k) = k-steps ahead forecast for y_t
     *  <p>
     *  Note, 'yf.col(0)' is set to 'y' (the actual time-series values).
     *  @param h  the maximum forecasting horizon, number of steps ahead to produce forecasts
     */
    def forecastAll (h: Int): MatriD =
    {
        yf = new MatrixD (m, h+1)                                        // forecasts for all time points t & horizons to h
        yf.setCol (0, y)                                                 // first column is actual values, horizon 0
        for (k <- 1 to h) {
            val c = min (k, p)                                           // cut point from actual to forecasted values
            for (t <- 0 until c) yf(t, k) = y(t)                         // copy first c actual values
            for (t <- c until m) {
                 var sum = δ
                 for (j <- 0 until p) sum += φ(j) * yf(max (0, t-1-j), max (0, k-1-j))
                 yf(t, k) = sum
            } // for
        } // for
        yf                                                               // return matrix of forecasted values
    } // forecastAll

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Produce a vector of size 'h', of 1 through 'h'-steps ahead forecasts for the model.
     *  @see Forecasting Equation in section 7.1.
     *  @param t  the time point from which to make forecasts
     *  @param h  the forecasting horizon, number of steps ahead to produce forecasts
     */
    override def forecast (t: Int, h: Int = 1): VectoD =
    {
        val yf = new VectorD (h+1)
        val sf = new VectorD (h+1)
        
        yf(0) = if (t == 0) y(0) else s(t-1)
        sf(0) = s(t)
        for (k <- 1 to h) {
            yf(k) = sf(k-1)
            sf(k) = α * yf(k-1) + (1 - α) * sf(k-1)
        } // for
        yf.slice (1)
    } // forecast

} // SimpleExpSmoothing class
 

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `SimpleExpSmoothing` companion object provides factory methods for the
 *  `SimpleExpSmoothing` class.
 */
object SimpleExpSmoothing
{
    /** Base hyper-parameter specification for `SimpleExpSmoothing`
     */
    val hp = new HyperParameter;
    hp += ("α", 0.5, 0.5)                                        // default value for the smoothing parameter

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Create a `SimpleExpSmoothing` object.
     *  @param y       the response vector (time series data)
     *  @param hparam  the hyper-parameters
     */
    def apply (y: VectoD, hparam: HyperParameter = hp): SimpleExpSmoothing =
    {
        new SimpleExpSmoothing (y, hparam)
    } // apply

} // SimpleExpSmoothing object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `SimpleExpSmoothingTest` object is used to test the `SimpleExpSmoothing` class.
 *  Test customized smoothing (call 'smooth') versus optimized smoothing (call 'train').
 *  > runMain scalation.analytics.forecaster.SimpleExpSmoothingTest
 */
object SimpleExpSmoothingTest extends App
{
    val m = 50
    val r = Random ()
    val y = VectorD (for (i <- 0 until m) yield i + 10.0 * r.gen)
//  val r2 = RandomVecTrend (m, noise = Uniform (0.0, 10.0))
//  val y2 = r2.gen
//  new Plot (null, y, y2, "RandomVecTrend")
    val h = 3
    println (s"y = $y")

    val ses = new SimpleExpSmoothing (y)                         // smooth time series data: y vs. t

    banner ("Customized Simple Exponential Smoothing")
    ses.smooth (0.5)                                             // use customized parameters, don't train
    ses.eval () 
    println (ses.report)
    val yf = ses.predictAll ()                                   // y forecasted
    println (s"yf = $yf")
    println (s"ses.forecast ($h) = ${ses.forecast (h)}")
    new Plot (null, y, yf, "SES custom: Plot of y, yf vs. t", true)

    banner ("Optimized Simple Exponential Smoothing")
    ses.train (null, y)                                          // train to use optimal α
    ses.eval ()
    println (ses.report)
    val yf2 = ses.predictAll ()                                  // y forecasted
    println (s"yf2 = $yf2")
    println (s"forecast ($h) = ${ses.forecast (h)}")
    new Plot (null, y, yf2, "SES optimize: Plot of y, yf2 vs. t", true)

} // SimpleExpSmoothingTest object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `SimpleExpSmoothingTest2` object is used to test the `SimpleExpSmoothing` class.
 *  Test rolling validation.
 *  > runMain scalation.analytics.forecaster.SimpleExpSmoothingTest2
 */
object SimpleExpSmoothingTest2 extends App
{
    val m = 50
    val r = Random ()
    val y = VectorD (for (i <- 0 until m) yield i + 10.0 * r.gen)
    val h = 3
    println (s"y = $y")

    val ses = new SimpleExpSmoothing (y)                         // smooth time series data: y vs. t

    banner ("Optimized Simple Exponential Smoothing")
    ses.train (null, y)                                          // train to use optimal α
    ses.eval ()
    println (ses.report)
    val yf = ses.predictAll ()
    println (s"yf = $yf")
    println (s"forecast ($h) = ${ses.forecast (h)}")
    new Plot (null, y, yf, "SES optimize: Plot of y, yf2 vs. t", true)

    val yfm = ses.forecastAll (4, 1)
    for (h <- 1 to 4) {                                          // h-steps ahead  forecast
        banner (s"forecastAll h = $h")
        new Plot (null, y, yfm(h), s"SES: Plot y and yfm(${h})", lines = true)

        banner (s"Rolling Validation h = $h")
        val stats = SimpleRollingValidation.crossValidate2 (ses, kt_ = 1, h = h)
//      val stats = RollingValidation.crossValidate2 (ses, kt_ = 5, h = h)
        Fit.showQofStatTable (stats)
    } // for

} // SimpleExpSmoothingTest object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `SimpleExpSmoothingTest3` object is used to test the `SimpleExpSmoothing` class.
 *  Forecasting lake levels.
 *  @see cran.r-project.org/web/packages/fpp/fpp.pdf
 *  > runMain scalation.analytics.forecaster.SimpleExpSmoothingTest3
 */
object SimpleExpSmoothingTest3 extends App
{   
    import ForecasterVec.{t, y}
    
    val ses = new SimpleExpSmoothing (y)                         // time series model SimpleExpSmoothing
    val ar  = new AR (y)                                         // time series model AR(1)
    
    // Build AR(1) and SimpleExpSmoothing models for the time series data
    
    banner (s"Build AR(1) model")
    ar.train (null, y)
    ar.eval ()
    val yf1 = ar.predictAll ()
    println (ar.report)
    new Plot (t, y, yf1, s"Plot of y, AR(1) vs. t", true)

    banner (s"Build SimpleExpSmoothing model")
    ses.train (null, y)
    ses.eval ()
    val yf2 = ses.predictAll ()
    val yf3 = ses.forecastAll (1, 1)(1)
    println (s"yf3 = $yf3")
    assert (yf2 == yf3)
    println (ses.report)
    new Plot (t, y, yf2, s"Plot of y, SimpleExpSmoothing vs. t", true)
    new Plot (t, y, yf3, s"Plot of y, SimpleExpSmoothing3 vs. t", true)

    for (h <- 1 to 4) {                                          // h-steps ahead  forecast
        banner (s"Rolling Validation h = $h")
        val stats = SimpleRollingValidation.crossValidate2 (ses, kt_ = 1, h = h)
//      val stats = RollingValidation.crossValidate2 (ses, kt_ = 5, h = h)
        Fit.showQofStatTable (stats)
    } // for

} // SimpleExpSmoothingTest3 object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `SimpleExpSmoothingTest4` object is used to test the `SimpleExpSmoothing` class.
 *  Forecasting lake levels for several values of the smoothing parameter α.
 *  @see cran.r-project.org/web/packages/fpp/fpp.pdf
 *  > runMain scalation.analytics.forecaster.SimpleExpSmoothingTest4
 */
object SimpleExpSmoothingTest4 extends App
{  
    import ForecasterVec.{t, y}

    val ses = new SimpleExpSmoothing (y)                         // time series model SimpleExpSmoothing
    ses.toggleOpt ()                                             // switch auto optimization off

    for (i <- 0 to 5) {
        val a = i.toDouble / 5.0
        banner (s"Build SimpleExpSmoothing model with α = $a")
        ses.reset (a)
        ses.train (null, y)
        ses.eval ()
        val yf = ses.predictAll ()
        println (ses.report)
        new Plot (t, y, yf, s"Plot of y, SimpleExpSmoothing (a = $a) vs. t", true)
    } // for

} // SimpleExpSmoothingTest4 object