//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller
 *  @version 1.6
 *  @date    Tue May 11 16:25:40 EDT 2021
 *  @see     LICENSE (MIT style license file).
 *
 *  @title   Model: Quadratic Spline
 */

package scalation.analytics
package forecaster

import scalation.linalgebra.{MatriD, MatrixD, VectoD, VectorD}
import scalation.plot.Plot
import scalation.util.banner

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `QuadSpline` class fits quadratic splines to time-series data that are equally
 *  spaced in time.  A sliding window consisting of three data points is perfectly fit
 *  to a quadratic curve.
 *  <p>
 *      y_t = a + bt + ct^2
 *  <p>
 *  Note, slope matching and smoothness issues are ignored.
 *  @see wordsandbuttons.online/quadratic_splines_are_useful_too.html
 *  Any time point from t = 3 to the end of time series may be forecasted.
 *  @param y       the time-series
 *  @param hparam  the hyper-parameters
 */
class QuadSpline (y: VectoD, hp: HyperParameter = null)
      extends ForecasterVec (y, 1) with NoFeatureSelectionF
{
    private val DEBUG = true                                          // debug flag
    private val shift = 2                                             // shift center from y_t or y_{t-2} 

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

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Based on three points y_{t-1}, y_t, t_{t+1}, determine values for the
     *  coefficients 'a', 'b' and 'c'.
     *  @param t  the center time point
     */
    def splineFit (t: Int): (Double, Double, Double) =
    {
        val c = 0.5 * (y(t+1) - 2*y(t) + y(t-1))
        val b = 0.5 * (y(t+1) - y(t-1) - 4*c*t)
        val a = y(t) - b*t - c*t*t
        (a, b, c)
    } // splineFit   

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Evaluate the spline function at time point 't', given the coefficients 'a', 'b' and 'c'.
     *  @param t  the time
     *  @param a  the constant term
     *  @param b  the linear term coefficient
     *  @param c  the quadratic term coefficient
     */
    def spline (t: Double, a: Double, b: Double, c: Double): Double =
    {
        a + b*t + c*t*t
    } // spline

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Forecast a one-step ahead value for 'y_t' based on the quadratic curve fit to
     *  the previous three vales: y_{t-3}, y_{t-2}, t_{t-1}.
     *  @param t  the time at which to forecast y
     */
    def forecast1 (t: Int): Double =
    {
        if (t <= shift) y(t)
        else {
            val (a, b, c) = splineFit (t-shift)
            spline (t, a, b, c)
        } // if
    } // forecast1
    
    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Train/fit a `QuadSpline` model to the times-series data in vector 'y_'.
     *  Note: for `QuadSpline` 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): QuadSpline = { super.train (null, y_); this }
    
    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the parameter vector (its null).
     */ 
    def parameter: VectoD = null
    
    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return a vector that is the predictions of a quad spline model, by making forecasts
     *  for all values from time 3 to the end of the time-series.  Note, y_0, y_1 and y_2
     *  can't have forecasts, since they would need a value for y_{-1}.
     */ 
    override def predictAll (): VectoD =
    {   
        val yf = new VectorD (m)                                   // forecasts for all z
        for (t <- 0 to shift)      yf(t) = y(t)                    // copy actual value
        for (t <- shift+1 until m) yf(t) = forecast1(t)            // enter forecasted value
        yf                                                         // return the vector of predicted values
    } // predictAll

    def predictAllz (): VectoD = predictAll () - stats.mu

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** 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) {
            yf(0, k) = y(0)                                              // copy first actual value
            for (t <- 1 until m) {                                       // forecast the rest
                 yf(t, k) = forecast1 (t)                                // FIX - implement for other h beyond 1
            } // for
            if (DEBUG) println (s"forecastAll: yf.col ($k) = ${yf.col (k)}")
        } // for
        yf                                                               // return matrix of forecasted values
    } // forecastAll

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Test the curve produced by the multiple splines.
     */
    def testCurve ()
    {
        println (s"y = $y")
        for (i <- 1 until y.dim - 1) {
            val (a, b, c) = splineFit(i)
            for (j <- -2 to 2) {
                val t   = i + 0.5 * j
                val f_t = spline (t, a, b, c)
                print (s"spline($t) = $f_t \t")
            } // for
            println
        } // for
    } // testCurve

} // QuadSpline class


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `QuadSpline` companion object provides factory methods for the `QuadSpline` class.
 */
object QuadSpline
{
    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Create a `QuadSpline` object.
     *  @param y       the response vector (time series data)
     *  @param hparam  the hyper-parameters
     */
    def apply (y: VectoD, hparam: HyperParameter = null): QuadSpline =
    {
        new QuadSpline (y, hparam)
    } // apply

} // RandomWalk object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `QuadSplineTest` object is used to test the `QuadSpline` class.
 *  Forecasting Fibonacci numbers.
 *  > runMain scalation.analytics.forecaster.QuadSplineTest
 */
object QuadSplineTest extends App
{
    val y = VectorD (1, 2, 3, 5, 8, 13, 21, 34, 55, 89)

    banner ("RandomWalk Model")
    val rw = new RandomWalk (y)
    rw.train (null, y).eval ()
    println (rw.report)
    val yp = rw.predictAll ()
    new Plot (null, y, yp, "RandomWalk: y vs. yp", lines = true)    

    banner ("QuadSpline Model")
    val qs = new QuadSpline (y)
    qs.train (null, y).eval ()
    println (qs.report)
    qs.testCurve ()
    val yf = qs.predictAll ()
    new Plot (null, y, yf, "QuadSpline: y vs. yf", lines = true)    

} // QuadSplineTest object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `QuadSplineTest2` object is used to test the `QuadSpline` class.
 *  > runMain scalation.analytics.forecaster.QuadSplineTest2
 */
object QuadSplineTest2 extends App
{
    // TBD

} // QuadSplineTest object


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

    banner ("RandomWalk Model")
    val rw = new RandomWalk (y)
    rw.train (null, y).eval ()
    println (rw.report)
    val yp = rw.predictAll ()
    new Plot (null, y, yp, "RandomWalk: y vs. yp", lines = true)    

    banner ("QuadSpline Model")
    val qs = new QuadSpline (y)
    qs.train (null, y).eval ()
    println (qs.report)
    val yf = qs.predictAll ()
    new Plot (null, y, yf, "QuadSpline: y vs. yf", lines = true)    

    val mix = (yp + yf) * 0.5
    new Plot (null, y, mix, "Mix: y vs. mix", lines = true)    

} // QuadSplineTest3 object