//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller
 *  @version 1.6
 *  @date    Mon Feb  2 18:18:15 EST 2015
 *  @see     LICENSE (MIT style license file).
 *
 *  @title   Model: Trigonometric Regression
 */

package scalation.analytics

import scala.collection.mutable.{ListMap, Map, Set}
import scala.math.{cos, Pi, sin}

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

import Fit.showQofStatTable
import MatrixTransform._
import RegTechnique._

//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `TrigRegression` class supports trigonometric regression.  In this case,
 *  't' is expanded to '[1, sin (wt), cos (wt), sin (2wt), cos (2wt), ...]'.
 *  Fit the parameter vector 'b' in the regression equation
 *  <p>
 *      y  =  b dot x + e  =  b_0 + b_1 sin (wt)  + b_2 cos (wt)  +
 *                                  b_3 sin (2wt) + b_4 cos (2wt) + ... + e
 *  <p>
 *  where 'e' represents the residuals (the part not explained by the model).
 *  Use Least-Squares (minimizing the residuals) to solve for the parameter vector 'b'
 *  using the Normal Equations:
 *  <p>
 *      x.t * x * b  =  x.t * y
 *      b  =  fac.solve (.)
 *  <p>
 *  @see link.springer.com/article/10.1023%2FA%3A1022436007242#page-1
 *  @param t          the initial data/input m-by-1 matrix: t_i expands to x_i
 *  @param y          the response/ouput vector
 *  @param ord        the order (k), maximum multiplier in the trig function (kwt)
 *  @param fname_     the feature/variable names
 *  @param hparam     the hyper-parameters
 *  @param technique  the technique used to solve for b in x.t*x*b = x.t*y
 */
class TrigRegression (t: MatriD, y: VectoD, ord: Int, fname_ : Strings = null, hparam: HyperParameter = null,
                      technique: RegTechnique = QR)
      extends Regression (TrigRegression.allForms (t, ord), y, fname_, hparam, technique)
      with ExpandableForms
{
    private val w  = (2.0 * Pi) / (t.max() - t.min())            // base displacement angle in radians
    private val n0 = 1                                           // number of terms/columns originally
    private val nt = TrigRegression.numTerms (2 * ord)           // number of terms/columns after expansion

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Expand the vector 'z' into a vector of that includes additional terms,
     *  i.e., add polynomial terms.
     *  @param z  the un-expanded vector
     */
    def expand (z: VectoD): VectoD = TrigRegression.forms (z, n0, nt, w)

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given the scalar 'z', expand it and predict the response value.
     *  @param z  the un-expanded scalar
     */
    def predict (z: Double): Double = predict_ex (VectorD (z))

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given the vector 'z', expand it and predict the response value.
     *  @param z  the un-expanded vector
     */
    def predict_ex (z: VectoD): Double = predict (expand (z))

} // TrigRegression class


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `TrigRegression` companion object provides factory functions and functions
 *  for creating functional forms.
 */
object TrigRegression extends ModelFactory
{
    val drp = (null, null, QR)                                   // default remaining parameters
    private val DEBUG = false                                    // debug flag

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Create a `TrigRegression` object from a combined data-response matrix.
     *  @param xy         the initial combined data-response matrix (before polynomial term expansion)
     *  @param ord        the order (k) of the polynomial (max degree)
     *  @param fname_     the feature/variable names
     *  @param hparam     the hyper-parameters
     *  @param technique  the technique used to solve for b in x.t*x*b = x.t*y
     */
    def apply (xy: MatriD, ord: Int, fname: Strings = null, hparam: HyperParameter = null,
               technique: RegTechnique = Cholesky): TrigRegression =
    {
        val n = xy.dim2
        if (n < 2) {
            flaw ("apply", s"dim2 = $n of the 'xy' matrix must be at least 2")
            null
        } else {
            val (x, y) = pullResponse (xy)
            new TrigRegression (x, y, ord, fname, hparam, technique)
        } // if
    } // apply

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Create a `TrigRegression` object from a combined data-response matrix.
     *  @param t          the initial data/input vector: t_i expands to x_i = [1, t_i, t_i^2, ... t_i^k]
     *  @param y          the response/ouput vector
     *  @param ord        the order (k) of the polynomial (max degree)
     *  @param fname_     the feature/variable names
     *  @param hparam     the hyper-parameters
     *  @param technique  the technique used to solve for b in x.t*x*b = x.t*y
     */
    def apply (t: VectoD, y: VectoD, ord: Int, fname: Strings, hparam: HyperParameter,
               technique: RegTechnique): TrigRegression =
    {
        new TrigRegression (MatrixD (Seq (t)), y, ord, fname, hparam, technique)
    } // apply

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Create a `TrigRegression` object from a data matrix and a response vector.
     *  This factory function provides data rescaling.
     *  @see `ModelFactory`
     *  @param x          the initial data/input matrix (before polynomial term expansion)
     *  @param y          the response/output m-vector
     *  @param ord        the order (k) of the polynomial (max degree)
     *  @param fname      the feature/variable names (use null for default)
     *  @param hparam     the hyper-parameters (use null for default)
     *  @param technique  the technique used to solve for b in x.t*x*b = x.t*y
     */
    def apply (x: MatriD, y: VectoD, ord: Int, fname: Strings, hparam: HyperParameter,
               technique: RegTechnique): TrigRegression =
    {
        val n = x.dim2
        if (n < 1) {
            flaw ("apply", s"dim2 = $n of the 'x' matrix must be at least 1")
            null
        } else if (rescale) {                                    // normalize the x matrix
            val (mu_x, sig_x) = (x.mean, stddev (x))
            val xn = normalize (x, (mu_x, sig_x))
            new TrigRegression (xn, y, ord, fname, hparam, technique)
        } else {
            new TrigRegression (x, y, ord, fname, hparam, technique)
       } // if
    } // apply

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given a 1-vector/point 'v', compute the values for all of its trigonmetric
     *  forms/terms, returning them as a vector.
     *  '[1, sin (wt), cos (wt), sin (2wt), cos (2wt), ...]'.
     *  @param v   the vector/point (i-th row of t) for creating forms/terms
     *  @param k   number of features/predictor variables (not counting intercept) = 1
     *  @param nt  the number of terms
     *  @param w   the base displacement angle in radians
     */
    def forms (v: VectoD, k: Int, nt: Int, w: Double): VectoD =
    {
        val wt = w * v(0)
        val u  = new VectorD (nt)
        u(0)   = 1.0
        for (j <- 1 to nt / 2) {
            u(2*j-1) = sin (j * wt)
            u(2*j)   = cos (j * wt)
        } // for
        u
    } // forms

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Create all forms/terms for each row/point placing them in a new matrix.
     *  @param x    the original un-expanded input/data matrix
     *  @param ord  the order (max degree) of the polynomial
     */
    def allForms (x: MatriD, ord: Int): MatriD =
    {
        val t  = x.col (0)                                           // first and only column of x
        val w  = (2.0 * Pi) / (t.max() - t.min())                    // base displacement angle in radians
        val k  = 1
        val nt = numTerms (2 * ord)
        println (s"allForms: create expanded data matrix with nt = $nt columns from k = $k columns")
        val xe = new MatrixD (x.dim1, nt)
        for (i <- x.range1) xe(i) = forms (x(i), k, nt, w)           // vector with values for all forms/terms
        if (DEBUG) println (s"expanded data matrix xe = $xe")
        xe                                                           // expanded matrix
    } // allForms

} // TrigRegression object

import TrigRegression.drp

//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `TrigRegressionTest` object tests `TrigRegression` class using the following
 *  regression equation.
 *  <p>
 *      y  =  b dot x  =  b_0 + b_1 sin wt + b_2 cos wt + ... b_2k-1 sin kwt + b_2k cos kwt + e
 *  <p>
 *  The data is generated from a noisy cubic function.
 *  > runMain scalation.analytics.TrigRegressionTest
 */
object TrigRegressionTest extends App
{
    import scalation.random.Normal

    val noise = Normal (0.0, 10000.0)
    val t     = VectorD.range (0, 100)
    val y     = new VectorD (t.dim)
    for (i <- 0 until 100) { val x = (i - 40)/2.0; y(i) = 1000.0 + x + x*x + x*x*x + noise.gen }

    println ("t = " + t)
    println ("y = " + y)

    for (harmonics <- 2 to 16 by 2) {
        val trg = TrigRegression (t, y, harmonics, drp._1, drp._2, drp._3)
        trg.train ().eval ()
        println (trg.report)

        val z   = 10.5                                  // predict y for one point
        val yp1 = trg.predict (z)
        println (s"predict ($z) = $yp1")

        banner ("TrigRegressionTest: test predictions for harmonics = $harmonics")
        val yp = t.map (trg.predict (_))
//      println (s" y = $y \n yp = $yp")
        new Plot (t, y, yp, s"TrigRegression: harmonics = $harmonics")

        if (harmonics == 16) {
           val stats = trg.crossValidate ()
           showQofStatTable (stats)
        } // if
    } // for

} // TrigRegressionTest object


//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `TrigRegressionTest2` object tests `TrigRegression` class using the following
 *  regression equation.
 *  <p>
 *      y  =  b dot x  =  b_0 + b_1 sin wt + b_2 cos wt + ... b_2k-1 sin kwt + b_2k cos kwt + e
 *  <p>
 *  The data is generated from periodic noisy cubic functions.
 *  > runMain scalation.analytics.TrigRegressionTest2
 */
object TrigRegressionTest2 extends App
{
    import scalation.random.Normal

    val noise = Normal (0.0, 10.0)
    val t     = VectorD.range (0, 200)
    val y     = new VectorD (t.dim)
    for (i <- 0 until 5) {
        for (j <- 0 until 20) { val x = j - 4;  y(40*i+j) = 100.0 + x + x*x + x*x*x + noise.gen }
        for (j <- 0 until 20) { val x = 16 - j; y(40*i+20+j) = 100.0 + x + x*x + x*x*x + noise.gen }
    } // for

    println ("t = " + t)
    println ("y = " + y)

    for (harmonics <- 2 to 16 by 2) {
        val trg = TrigRegression (t, y, harmonics, drp._1, drp._2, drp._3)
        trg.train ().eval ()
        println (trg.report)

        val z   = 10.5                                  // predict y for one point
        val yp1 = trg.predict (z)
        println (s"predict ($z) = $yp1")

        banner ("TrigRegressionTest: test predictions for harmonics = $harmonics")
        val yp = t.map (trg.predict (_))
//      println (s" y = $y \n yp = $yp")
        new Plot (t, y, yp, s"TrigRegression: harmonics = $harmonics")

        if (harmonics == 16) {
           val stats = trg.crossValidate ()
           showQofStatTable (stats)
        } // if
    } // for

} // TrigRegressionTest2 object