//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller
 *  @version 1.6
 *  @date    Tue Apr  2 18:46:44 EDT 2013
 *  @see     LICENSE (MIT style license file).
 *
 *  @title   Line Search Optimizer
 *
 *  @see web.stanford.edu/class/cme304/docs/newton-type-methods.pdf
 */

package scalation.minima

import scala.math.{abs, sqrt}
import scala.util.control.Breaks.{breakable, break}

import scalation.calculus.Differential.{Ⅾ, Δ}
import scalation.math.FunctionS2S

import FunctionSelector._

//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `NewtonRaphson` class is used to find roots (zeros) for a one-dimensional
 *  (scalar) function 'g'.  Depending on the FunctionSelector, it can find zeros
 *  for derivatives or finite differences, which may indicate optima for function 'g'.
 *  Also, for optimization, may pass the derivative of the function, since finding
 *  zeros for the derivative corresponds to finding optima for the function.
 *  @param g     the function to find roots/optima of
 *  @param root  find the root for function, derivative or finite difference
 *  @param eta   the learning rate
 */
class NewtonRaphson (g: FunctionS2S, root: FunctionSelector = Function, eta: Double = 0.5)
{
    private val DEBUG    = true                         // debug flag
    private val EPSILON  = 1E-9                         // a number close to zero
    private val MAX_ITER = 200                          // the maximum number of iterations
    private val f        = root match {                 // find roots for
                           case Function   => g         // function (zeros)
                           case Derivative => Ⅾ (g) _   // derivative (optima)
                           case _          => Δ (g) _   // finite difference (optima, noise)
                           } // match

    //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Solve for/find a root close to the starting point/guess 'x0'.  
     *  This version numerically approximates the derivative.
     *  @param x0  the starting point/guess
     */
    def solve (x0: Double): Double =
    {
        var x    = x0                                         // current point
        var f_x  = -1.0                                       // function value at x
        var df_x = -1.0                                       // derivative value at x

        breakable { for (it <- 0 until MAX_ITER) {
            f_x = f(x)
            if (abs (f_x) < EPSILON) {print (s"it = $it, "); break }  // close to zero => quit
            df_x = Ⅾ (f)(x)               
            if (DEBUG) println (s"it = $it: f($x) = $f_x, df_x = $df_x")
            if (abs (df_x) < EPSILON) df_x = EPSILON          // derivative is too small
        }} // for

        printf ("solution x = %10.5f, f = %10.5f\n", x, f(x))
        x
    } // solve

    //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Solve for/find a root close to the starting point/guess 'x0'.
     *  This version passes in a function for the derivative.
     *  @param x0  the starting point/guess
     *  @param df  the derivative of the function
     */
    def solve (x0: Double, df: FunctionS2S): Double =
    {
        var x    = x0                                         // current point
        var f_x  = -1.0                                       // function value at x
        var df_x = -1.0                                       // derivative value at x

        breakable { for (it <- 0 until MAX_ITER) {
            f_x = f(x)
            if (abs (f_x) < EPSILON) {print (s"it = $it, "); break }  // close to zero => quit
            df_x = df(x)
            if (DEBUG) println (s"it = $it: f($x) = $f_x, df_x = $df_x")
            if (abs (df_x) < EPSILON) df_x = EPSILON          // derivative is too small
            x -= eta * f_x / df_x                             // subtract the ratio
        }} // for

        printf ("solution x = %10.5f, f = %10.5f\n", x, f(x))
        x
    } // solve

} // NewtonRaphson class


//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `NewtonRaphsonTest` object is used to test the `NewtonRaphson` class.
 *  This test numerically approximates the derivative.
 *  > runMain scalation.minima.NewtonRaphsonTest
 */
object NewtonRaphsonTest extends App
{
//  def f (x: Double): Double = 2 * (x*x - 3)
    def f (x: Double): Double = (x - 4) * (x - 4)

    val nr = new NewtonRaphson (f)
    
    nr.solve (0.0)

} // NewtonRaphsonTest object


//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `NewtonRaphsonTest2` object is used to test the `NewtonRaphson` class.
 *  This test passes in a function for the derivative.
 *  > runMain scalation.minima.NewtonRaphsonTest2
 */
object NewtonRaphsonTest2 extends App
{
//  def f (x: Double): Double = 2 * (x*x - 3.0)
    def f (x: Double): Double = (x - 4) * (x - 4)
    def df (x: Double): Double = 2 * (x - 4)

    val nr = new NewtonRaphson (f)
   
    nr.solve (0.0, df)

} // NewtonRaphsonTest object