//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller, Hao Peng
 *  @version 1.6
 *  @date    Sun May 28 13:26:35 EDT 2017
 *  @see     LICENSE (MIT style license file).
 *
 *  @title   Model: Exponential Smoothing
 *
 *  @see www.itl.nist.gov/div898/handbook/pmc/section4/pmc432.htm
 *  @see en.wikipedia.org/wiki/Exponential_smoothing#Triple_exponential_smoothing
 */

package scalation.analytics
package forecaster

// U N D E R   D E V E L O P M E N T

import scala.collection.mutable.Set

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

import RollingValidation.trSize

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `ExpSmoothing` class provide very basic time series analysis capabilities of
 *  Exponential Smoothing models.  `ExpSmoothing` models are often used for forecasting.
 *  Given time series data stored in vector 'y', its next value 'y_t = y(t)'
 *  may be predicted based on prior/smoothed values of 'y':
 *  <p>
 *      y_t = s_t-1 + α (s_t-1 - s_t-2)
 *  <p>
 *  where vector 's' is the smoothed version of vector 'y' and 'α in [0, 1]' is the
 *  smoothing parameter. Trend and seasonality can be factored into the model with
 *  two additional smoothing parameters 'β' and 'γ', respectively.
 *  @param y_      the input vector (orginal time series data)
 *  @param hparam  the hyper-parameters
 */
class ExpSmoothing (y_ : VectoD, hparam: HyperParameter = ExpSmoothing.hp)
      extends ForecasterVec (y_, n = 3, hparam) with NoFeatureSelectionF
{
    private val DEBUG      = false                           // debug flag
    private val ll         = hparam ("ll").toInt             // seasonal period
    private val multi      = hparam ("multi").toInt          // multiplicative/additive seasonality
    private val vSteps     = hparam ("vSteps").toInt         // number of validation steps
    private var α          = hparam ("α")                    // default value for the smoothing parameter
    private var β          = hparam ("β")                    // default value for the trend smoothing parameter
    private var γ          = hparam ("γ")                    // default value for the seasonal change smoothing parameter

    private var y: VectoD  = null                            // the current time series data
    private var n          = 0                               // size of the input vector
    private var mu_        = 0.0                             // the sample mean for current data (mu for original)
    private var sig2_      = 0.0                             // the sample variance for current data
    private var s: VectoD  = null                            // holder for the smoothed data (constant part)
    private var b: VectoD  = null                            // holder for the smoothed data (trend part)
    private var c: VectoD  = null                            // holder for the smoothed data (seasonal part)
    private var nn         = 0                               // number of complete cycles/seasons in data
    private var aa: VectoD = null                            // holder for the average value of 'y' in each cycle/season of data
    private var yp: VectoD = null                            // vector of predicted values

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the model name including current its hyper-parameter.
     */
    override def modelName: String = "ExpSmoothing"

    setTS (y_)                                               // initialize the above parameters

    if (DEBUG) {
        println ("n    = " + n)                              // number of data points
        println ("mu_  = " + mu_)                            // mean
        println ("sig2_ = " + sig2_)                           // variance
    } // for

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Set/change the internal time series.  May be used to set the time series
     *  to a different time window (typically future when new data become available)
     *  in order to produce newer forecast (typically with the new data) without re-training
     *  the model for parameters (use existing parameters from previous training).
     *  @param y_  the new time series
     */
    def setTS (y_ : VectoD)
    {
        y        = y_
        n        = y.dim                                     // size of the input vector
        mu_      = y.mean                                    // the sample mean
        sig2_    = y.variance                                // the sample variance
        s        = new VectorD (n)                           // holder for the smoothed data (constant part)
        b        = new VectorD (n)                           // holder for the smoothed data (trend part)
        c        = new VectorD (n)                           // holder for the smoothed data (seasonal part)
        nn       = n / ll                                    // number of complete cycles/seasons in data
        aa       = new VectorD (nn)                          // holder for the average value of 'y' in each cycle/season of data
        e        = new VectorD (n)                           // initialize the error vector
        yp       = new VectorD (n)                           // initialize vector of predicted values

        if (n < 2*ll) flaw ("ExpSmoothing", "Usually a minimal of 2 full seasons are needed to initialize seasonal parameters")
        init ()
    } // setTS

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Compute initial values of 's', 'b', 'aa' and 'c'.
     */
    def init ()
    {
        s(0) = y(0)

        var sum = 0.0
        for (i <- 0 until ll) sum += (y(ll+i) - y(i)) / ll
        b(0) = sum / ll

        for (j <- 0 until nn) {
            sum = 0.0
            for (i <- 0 until ll) sum += y(ll*j+i)
            aa(j) = sum / ll
        } // for

        for (i <- 0 until ll) {
            sum = 0.0
            for (j <- 0 until nn) sum += y(ll*j+i) / aa(j)
            c(i) = sum / nn
        } // for
    } // init

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Smooth the times series data.
     *  @param input  the input vector of 'α', 'β' and 'γ' to be optimized
     */
    def smooth (input: VectoD = VectorD (α, β, γ)): VectoD =
    {
        α = input(0)
        β = input(1)
        γ = input(2)

        val α_1 = 1.0 - α
        val β_1 = 1.0 - β
        val γ_1 = 1.0 - γ

        if (multi == 1) {
            for (t <- 1 until n) {
                if (t - ll >= 0) {
                    s(t) = α * y(t) / c(t-ll) + α_1 * (s(t-1) + b(t-1))
                    b(t) = β * (s(t) - s(t-1)) + β_1 * b(t-1)
                    c(t) = γ * y(t) / s(t) + γ_1 * c(t-ll)
                } else {
                    s(t) = α * y(t) + α_1 * (s(t-1) + b(t-1))
                    b(t) = β * (s(t) - s(t-1)) + β_1 * b(t-1)
                } // if
            } // for
        } else {
            for (t <- 1 until n) {
                if (t - ll >= 0) {
                    s(t) = α * (y(t) - c(t-ll)) + α_1 * (s(t-1) + b(t-1))
                    b(t) = β * (s(t) - s(t-1)) + β_1 * b(t-1)
                    c(t) = γ * (y(t) - s(t)) + γ_1 * c(t-ll)
                } else {
                    s(t) = α * y(t) + α_1 * (s(t-1) + b(t-1))
                    b(t) = β * (s(t) - s(t-1)) + β_1 * b(t-1)
                } // if
            } // for
        } // if

        for (i <- vSteps until n) {
            yp(i) = forecast (vSteps, i - vSteps).last
            e(i)  = y(i) - yp(i)
        } // for
        s
    } // smooth

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** The objective function to be minimized (sum of squared errors).
     *  @param input  the input vector of 'α', 'β' and 'γ' to be optimized
     */
    def f_obj (input: VectoD): Double = {smooth (input); e.normSq }

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Train the `ExpSmoothing` model to times series data, by finding the value for
     *  the smoothing parameter 'α', 'β' and 'γ' 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): ExpSmoothing =
    {
        super.train (null, y_)
        val l = VectorD.fill(3)(0.0)
        val u = VectorD.fill(3)(1.0)
        val solver = new L_BFGS_B (f_obj, l=l, u=u)
        val result = solver.solve (VectorD (α, β, γ), 1E-4)     // optimized value for α, β, γ
        smooth (result)
        if (DEBUG) println(s"(α, β, γ) = ${(α, β, γ)}")
        this
    } // train

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

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the vector of predicted values on the training data.
     */
    override def predictAll (): VectoD = yp

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the vector of centered predicted values on the training data.
     */
    override def predictAllz (): VectoD = yp - stats.mu

    def forecastAll (h: Int): MatriD = ???

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Produce a vector of size 'h', of 1 through 'h'-steps ahead forecasts for the model.
     *  @see ams.sunysb.edu/~zhu/ams586/Forecasting.pdf
     *  @param t  the time point from which to make forecasts (defaults to the last time point)
     *  @param h  the forecasting horizon, number of steps ahead to produce forecasts
     */
    override def forecast (t: Int = y.dim-1, h: Int = 1): VectoD =
    {
        val yf = new VectorD (h)
        for (i <- 1 to h) {
            if (multi == 1) {
                yf(i-1) = if (t-ll >= -1) (s(t) + i * b(t)) * c(t-ll+1+(i-1)%ll) else s(t) + i * b(t)
            } else {
                yf(i-1) = if (t-ll >= -1) s(t) + i * b(t) + c(t-ll+1+(i-1)%ll) else s(t) + i * b(t)
            } // if
        } // for
        yf
    } // forecast

} // ExpSmoothing class
 

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `ExpSmoothing` companion object provides factory methods for the `ExpSmoothing` class.
 */
object ExpSmoothing
{
    /** Base hyper-parameter specification for `ExpSmoothing`
     */
    val hp = new HyperParameter
    hp += ("ll", 1, 1)                        // seasonal period
    hp += ("multi", 0, 0)                     // whether to use multiplicative seasonality or not
                                              // if false, use additive seasonality
    hp += ("vSteps", 1, 1)                    // number of steps ahead within-sample forecast sse to minimize
    hp += ("α", 0.3, 0.3)                     // default value for the smoothing parameter
    hp += ("β", 0.1, 0.1)                     // default value for the trend smoothing parameter
    hp += ("γ", 0.1, 0.1)                     // default value for the seasonal change smoothing parameter

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

} // ExpSmoothing object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `ExpSmoothingTest` object is used to test the `ExpSmoothing` class.
 *  > runMain scalation.analytics.forecaster.ExpSmoothingTest
 */
object ExpSmoothingTest extends App
{
    val n = 100
    val r = Random ()
    val t = VectorD.range (0, n)
    val y = VectorD (for (i <- 0 until n) yield t(i) + 10.0 * r.gen)
    val h = 3
    val ts = new ExpSmoothing (y)                                 // smooth time series data: y vs. t

    banner ("Customized Exponential Smoothing")
    ts.smooth ()                                                  // use customized parameters
    ts.eval ()
    val s = ts.predictAll
    println (s"fit = ${ts.fit}")
    println (s"predict (s) = ${ts.forecast (h)}")
    new Plot (t, y, s, "Plot of y, s vs. t")

    banner ("Optimized Exponential Smoothing")
    ts.train (null, y).eval ()                                    // use optimal α
    val s2 = ts.predictAll
    println (s"fit = ${ts.fit}")
    println (s"predict (s2) = ${ts.forecast (h)}")
    new Plot (t, y, s2, "Plot of y, s2 vs. t")

} // ExpSmoothingTest object