//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller
 *  @version 1.6
 *  @date    Sat Jun 13 01:27:00 EST 2017
 *  @see     LICENSE (MIT style license file).
 *
 *  @title   Model: Simple Moving Average with Derivatives (experimental)
 */

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

package scalation.analytics
package forecaster

import scala.collection.mutable.Set
import scala.math.{max, min}

import scalation.linalgebra.{MatriD, MatrixD, VectoD, VectorD}
import scalation.plot.Plot
import scalation.random.{Normal, Random}
import scalation.stat.{Statistic, vectorD2StatVector}
import scalation.util.banner

import Fit._
import ForecasterVec.testInSamp
import RollingValidation.trSize

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `MovingAverageD` class provides basic time series analysis capabilities.
 *  For a 'MovingAverageD' model with the time series data stored in vector 'y', the
 *  next value 'y_t = y(t)' may be predicted based on the mean of prior values
 *  of 'y' and its noise:
 *  <p>
 *      y_t = mean (y_t-1, ..., y_t-q) + e_t
 *  <p>
 *  where 'e' is the noise vector and 'q' is the number of prior values used to
 *  compute the mean.
 *  @param y       the response vector (time series data)
 *  @param hparam  the hyper-parameters
 */
class MovingAverageD (y: VectoD, hparam: HyperParameter = MovingAverage.hp)
      extends ForecasterVec (y, 1, hparam) with NoFeatureSelectionF
{
    private   val DEBUG = true                                           // debug flag
    private   val q     = hparam("q").toInt                              // take mean of last q values

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the model name including its current hyper-parameter, e.g., MovingAverageD(3)
     */
    override def modelName: String = s"MovingAverageD($q)"

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Train/fit an `MovingAverageD` model to the times series data in vector 'y_'.
     *  Note: for `MovingAverageD` there are no parameters to train.
     *  @param x_null  the data/input matrix (ignored)
     *  @param y_      the response/output vector (currently only works for y)
     */
    override def train (x_null: MatriD, y_ : VectoD): MovingAverageD =
    {
        super.train (null, y_)
        this
    } // train

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the parameter vector (its null).
     */
    def parameter: VectoD = null

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return a vector that is the predictions (zero-centered) of a moving average
     *  model, with side-effect of updating the error.
     */
    def predictAllz (): VectoD =
    {
        val zp   = new VectorD (m)
        val sumq = new SumQueue (q)

        for (i <- z.range) { sumq += z(i); zp(i) = sumq.mean }
        zp
    } // predictAllz

    def forecastAll (h: Int): MatriD = forecastAll (0, h)

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Produce a vector of size 'h', of 1 through 'h'-steps ahead forecasts for the model.
     *  @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)
        val sumq  = new SumQueue (q)
        val sumdq = new SumQueue (q)
        val sum2q = new SumQueue (q)
        for (i <- max (0, t - q - 1) until t) sumq  += y(i)
        for (i <- max (1, t - q - 1) until t) sumdq += (y(i) - y(i-1)) / 2.0
        for (i <- max (2, t - q - 1) until t) sum2q += (y(i) - 2 * y(i-1) + y(t-2)) / 4.0

//      for (i <- yf.range) { yf(i) = sumq.mean; sumq += yf(i) }
        for (i <- yf.range) { yf(i) = sumq.mean + sumdq.mean + sum2q.mean;
                              sumq  += yf(i);
                              sumdq += (yf(i) - yf(max (0, i-1))) / 2.0;
                              sum2q += (yf(i) - 2 * yf(max(0, i-1)) + yf(max (0, i-2))) / 4.0  }
        yf
    } // forecast

} // MovingAverageD class


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `MovingAverageD` companion object provides factory methods for the `MovingAverageD` class.
 */
object MovingAverageD
{
    /** Base hyper-parameter specification for `MovingAverageD`
     */
    val hp = new HyperParameter
    hp += ("q", 3, 3)                                                    // number of prior values for mean

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

} // MovingAverageD object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `MovingAverageDTest` object is used to test the `MovingAverageD` class.
 *  > runMain scalation.analytics.forecaster.MovingAverageDTest
 */
object MovingAverageDTest extends App
{
    val ran = Random ()
    val n = 100
    val t = VectorD.range (0, n)
    val y = VectorD (for (i <- 0 until n) yield t(i) + 10.0 * ran.gen)

    banner ("Build AR(1) Model")
    val ar = new AR (y)                                                  // time series model
    ar.train (null, y).eval ()                                           // train for AR(1) model
    println (ar.report)
    new Plot (t, y, ar.predictAll (), "Plot of y, AR(1) vs. t", true)

    banner ("Build MovingAverageD Model")
    val rw = new MovingAverageD (y)                                       // time series model
    rw.train (null, y).eval ()                                           // train for MovingAverage model
    println (rw.report)
    new Plot (t, y, rw.predictAll (), "Plot of y, MovingAverageD vs. t", true)

} // MovingAverageDTest object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `MovingAverageDTest2` object is used to test the `MovingAverageD` class.
 *  > runMain scalation.analytics.forecaster.MovingAverageDTest2
 */
object MovingAverageDTest2 extends App
{
    val noise = Normal (0.0, 1.0)
    val n = 30
    val t = VectorD.range (0, n)
    val y = VectorD (for (i <- 0 until n) yield i + noise.gen)

    banner ("Build AR(1) Model")
    val ar = new AR (y)                                                  // time series model
    ar.train (null, y).eval ()                                           // train for AR(1) model
    println (ar.report)
    new Plot (t, y, ar.predictAll (), "Plot of y, AR(1) vs. t", true)

    banner ("Build MovingAverageD Model")
    val rw = new MovingAverageD (y)                                       // time series model
    rw.train (null, y).eval ()                                           // train for MovingAverage model
    println (rw.report)
    new Plot (t, y, rw.predictAll (), s"Plot of y, MovingAverageD vs. t", true)

    banner ("Make Forecasts")
    val steps = 10                                                       // number of steps for the forecasts
    val rw_f = rw.forecast (steps)
    println (s"$steps-step ahead forecasts using MovingAverageD model = $rw_f")
    val tf = VectorD.range (n, n + steps)
    new Plot (tf, rw_f, null, s"Plot MovingAverageD forecasts vs. t", true)

} // MovingAverageDTest2 object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `MovingAverageDTest3` object is used to test the `MovingAverageD` class.
 *  Forecasting lake levels.
 *  @see cran.r-project.org/web/packages/fpp/fpp.pdf
 *  > runMain scalation.analytics.forecaster.MovingAverageDTest3
 */
object MovingAverageDTest3 extends App
{
    import ForecasterVec.y

    banner ("Build AR(1) Model")
    val ar = new AR (y)                                                  // time series model
    ar.train (null, y).eval ()                                           // train for AR(1) model
    testInSamp (y, 1, ar, 1)
    println (ar.summary)

    banner ("Build MovingAverageD Model")
    val ma = new MovingAverageD (y)                                      // time series model
    ma.train (null, y).eval ()                                           // train for MovingAverage model
    testInSamp (y, 0, ma, 1)
//  println (ma.summary)

    banner ("Make Forecasts")
    val steps = 2                                                        // number of steps for the forecasts
    val yf = ma.forecast (steps)
    println (s"$steps-step ahead forecasts using MovingAverageD model = $yf")

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

} // MovingAverageDTest3 object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `MovingAverageDTest4` object is used to test the `MovingAverageD` class.
 *  > runMain scalation.analytics.forecaster.MovingAverageDTest4
 */
object MovingAverageDTest4 extends App
{
    val path = BASE_DIR + "travelTime.csv"
    val data = MatrixD (path)

    val (t, y) = (data.col(0), data.col(1))

    banner ("Build AR(1) Model")
    val ar = new AR (y)                                                  // time series model
    ar.train (null, y).eval ()                                           // train for AR(1) model
    println (ar.report)
    new Plot (t, y, ar.predictAll (), "Plot of y, AR(1) vs. t", true)

    banner (s"Build MovingAverageD model")
    val rw = new MovingAverageD (y)                                       // time series model
    rw.train (null, y).eval ()                                           // train for MovingAverage model
    println (rw.report)
    new Plot (t, y, rw.predictAll (), s"Plot of y, MovingAverageD vs. t", true)

    banner ("Make Forecasts")
    val steps = 1                                                        // number of steps for the forecasts
    val rw_f = rw.forecast (steps)
    println (s"$steps-step ahead forecasts using MovingAverageD model = $rw_f")

} // MovingAverageDTest4 object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `MovingAverageDTest5` object is used to test the `MovingAverageD` class.
 *  > runMain scalation.analytics.forecaster.MovingAverageDTest5
 */
object MovingAverageDTest5 extends App
{
    val sig2  = 10000.0
    val noise = Normal (0.0, sig2)
    val n = 50
    val t = VectorD.range (0, n)
    val y = VectorD (for (i <- 0 until n) yield 40 * (i-1) - (i-2) * (i-2) + noise.gen)

    banner ("Build AR(1) Model")
    val ar = new AR (y)                                                  // time series model
    ar.train (null, y).eval ()                                           // train for AR(1) model
    println (ar.report)
    new Plot (t, y, ar.predictAll (), "Plot of y, AR(1) vs. t", true)

    banner ("Build MovingAverageD model")
    val rw = new MovingAverageD (y)                                       // time series model
    rw.train (null, y).eval ()                                           // train for MovingAverage model
    println (rw.report)
    new Plot (t, y, rw.predictAll (), s"Plot of y, MovingAverageD vs. t", true)

    banner ("Make Forecasts")
    val steps = 2                                                        // number of steps for the forecasts
    val rw_f = rw.forecast (steps)
    println (s"$steps-step ahead forecasts using MovingAverageD model = $rw_f")

} // MovingAverageDTest5 object