//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller
 *  @version 1.6
 *  @date    Wed Feb 20 17:39:57 EST 2013
 *  @see     LICENSE (MIT style license file).
 *
 *  @title   Model: Non-Linear Regression (non-linear in the parameters)
 */

package scalation.analytics

import scala.collection.mutable.Set
import scala.math.log

import scalation.linalgebra.{FunctionV_2V, MatriD, MatrixD, VectoD, VectorD}
import scalation.minima.QuasiNewton
import scalation.plot.Plot
import scalation.stat.Statistic

import MatrixTransform._

//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `NonLinRegression` class supports non-linear regression.  In this case,
 *  'x' can be multi-dimensional '[1, x1, ... xk]' and the function 'f' is non-linear
 *  in the parameters 'b'.  Fit the parameter vector 'b' in the regression equation
 *  <p>
 *      y  =  f(x, b) + e
 *  <p>
 *  where 'e' represents the residuals (the part not explained by the model).
 *  Use Least-Squares (minimizing the residuals) to fit the parameter vector
 *  'b' by using Non-linear Programming to minimize Sum of Squares Error 'SSE'.
 *  @see www.bsos.umd.edu/socy/alan/stats/socy602_handouts/kut86916_ch13.pdf
 *  @param x       the data/input matrix augmented with a first column of ones
 *  @param y       the response/output vector
 *  @param f       the non-linear function f(x, b) to fit
 *  @param b_init  the initial guess for the parameter vector b
 *  @param fname_  the feature/variable names
 *  @param hparam  the hyper-parameters (currently has none)
 */
class NonLinRegression (x: MatriD, y: VectoD, f: FunctionP2S,
                        b_init: VectorD, fname_ : Strings = null,
                        hparam: HyperParameter = null)
      extends PredictorMat (x, y, fname_, hparam) with NoFeatureSelectionMat
{
    if (y != null && x.dim1 != y.dim) flaw ("constructor", "dimensions of x and y are incompatible")

    private val DEBUG  = false                                    // debug flag

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Function to compute the Sum of Squares Error 'SSE' for given values for
     *  the parameter vector 'b'.
     *  @param b  the parameter vector
     */
    def sseF (b: VectoD): Double =
    {
        val yp = VectorD (for (i <- x.range1) yield f(x(i), b))   // create vector yp of predicted responses
        e = y - yp                                                // residual/error vector
        e dot e                                                   // residual/error sum of squares
    } // sseF

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Train the predictor by fitting the parameter vector (b-vector) in the
     *  non-linear regression equation for the response vector 'y_'.
     *  <p>
     *      y = f(x, b)
     *  <p>
     *  using the least squares method.
     *  Caveat:  Optimizer may converge to an unsatisfactory local optima.
     *           If the regression can be linearized, use linear regression for
     *           starting solution.
     *  @param x_  the training/full data/input matrix
     *  @param y_  the training/full response/output vector
     */
    def train (x_ : MatriD = x, y_ : VectoD = y): NonLinRegression =
    {
        // FIX - currently only works for x_ = x and y_ = y
        val bfgs = new QuasiNewton (sseF)                      // minimize sse using NLP
        b        = bfgs.solve (b_init)                         // estimate for b from optimizer

        this
    } // train

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Predict the value of 'y = f(z)' by evaluating the formula 'y = f(z, b)',
     *  i.e., '(b0, b1) dot (1.0, z1)'.
     *  @param z  the new vector to predict
     */
    override def predict (z: VectoD): Double = f(z, b)

} // NonLinRegression class


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `NonLinRegression` companion object provides factory methods for buidling perceptrons.
 *  Note, 'rescale' is defined in `ModelFactory` in Model.scala.
 */
object NonLinRegression extends ModelFactory
{
    private val DEBUG = false                                          // debug flag
    val drp = (null, null)                                             // default remaining parameters

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Create a `NonLinRegression` with automatic rescaling from a combined data matrix.
     *  @param xy      the combined data/input and response/output matrix
     *  @param f       the non-linear function f(x, b) to fit
     *  @param b_init  the initial guess for the parameter vector b
     *  @param fname   the feature/variable names
     *  @param hparam  the hyper-parameters (currently has none)
     */
    def apply (xy: MatriD, f: FunctionP2S, b_init: VectorD,
               fname: Strings = null, hparam: HyperParameter = null): NonLinRegression =
    {
        var itran: FunctionV_2V = null                                 // inverse transform -> original scale
        val (x, y) = pullResponse (xy)                                 // assumes the last column is the response

/*                                                                     // FIX - function needs bounds
        val x_s = if (rescale) rescaleX (x, f0)
                  else x
        val y_s = if (f0.bounds != null) { val y_i = rescaleY (y, f0); itran = y_i._2; y_i._1 }
                  else y

*/
        val (x_s, y_s) = (x, y)
        if (DEBUG) println (s" scaled: x = $x_s \n scaled y = $y_s")
        new NonLinRegression (x_s, y_s, f, b_init, fname, hparam)
    } // apply

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Create a `NonLinRegression` with automatic rescaling from a data matrix and response vector.
     *  @param x       the data/input matrix
     *  @param y       the response/output vector
     *  @param f       the non-linear function f(x, b) to fit
     *  @param b_init  the initial guess for the parameter vector b
     *  @param fname   the feature/variable names
     *  @param hparam  the hyper-parameters (currently has none)
     */
    def apply (x: MatriD, y: VectoD, f: FunctionP2S, b_init: VectorD,
               fname: Strings, hparam: HyperParameter): NonLinRegression =
    {
        var itran: FunctionV_2V = null                                 // inverse transform -> original scale

/*                                                                     // FIX - function needs bounds
        val x_s = if (rescale) rescaleX (x, f0)
                  else x
        val y_s = if (f0.bounds != null) { val y_i = rescaleY (y, f0); itran = y_i._2; y_i._1 }
                  else y
*/
        val (x_s, y_s) = (x, y)
        if (DEBUG) println (s" scaled: x = $x_s \n scaled y = $y_s")
        new NonLinRegression (x_s, y_s, f, b_init, fname, hparam)
    } // apply

} // NonLinRegression object

import NonLinRegression.drp

//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `NonLinRegressionTest` object tests the `NonLinRegression` class:
 *  y  =  f(x; b)  =  b0 + exp (b1 * x0).
 *  @see www.bsos.umd.edu/socy/alan/stats/socy602_handouts/kut86916_ch13.pdf
 *  Answers:  sse = 49.45929986243339
 *            fit = (VectorD (58.606566327280426, -0.03958645286504356), 0.9874574894685292)
 *            predict (VectorD (50.0)) = 8.09724678182599
 *  FIX: check this example
 */
object NonLinRegressionTest extends App
{
    import scala.math.exp

    // 5 data points: constant term, x1 coordinate, x2 coordinate
    val x = new MatrixD ((15, 1), 2.0, 5.0, 7.0, 10.0, 14.0, 19.0, 26.0, 31.0, 34.0, 38.0, 45.0, 52.0, 53.0, 60.0, 65.0)
    val y = VectorD (54.0, 50.0, 45.0, 37.0, 35.0, 25.0, 20.0, 16.0, 18.0, 13.0, 8.0, 11.0, 8.0, 4.0, 6.0)

    println ("x = " + x)
    println ("y = " + y)

    def f (x: VectoD, b: VectoD): Double = b(0) * exp (b(1) * x(0))   // non-linear regression function
    val b_init = VectorD (4.04, .038)                                   // initial guess for parameter vector b

    val rg = new NonLinRegression (x, y, f, b_init, drp._1, drp._2)
    rg.train ().eval ()
    println ("fit = " + rg.fit)

    val z  = VectorD (1); z(0) = 50.0             // predict y for one point
    val yp = rg.predict (z)
    println ("predict (" + z + ") = " + yp)

} // NonRegressionTest object