/*!
 * \file
 * \brief Newton Search optimization algorithms - header file
 * \author Tony Ottosson
 *
 * -------------------------------------------------------------------------
 *
 * IT++ - C++ library of mathematical, signal processing, speech processing,
 *        and communications classes and functions
 *
 * Copyright (C) 1995-2008  (see AUTHORS file for a list of contributors)
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 *
 * -------------------------------------------------------------------------
 */

#ifndef NEWTON_SEARCH_H
#define NEWTON_SEARCH_H

#include <itpp/base/vec.h>
#include <itpp/base/array.h>
#include <limits>


namespace itpp {

  /*!
    \brief Numerical optimization routines
    \addtogroup optimization
  */
  //@{


  //! Newton Search method
  enum Newton_Search_Method {BFGS};

  /*!
    \brief Newton Search

    Newton or Quasi-Newton optimization method that try to minimize the objective function \f$f(\mathbf{x})\f$
    given an initial guess \f$\mathbf{x}\f$.

    The search is stopped when either criterion 1:
    \f[
    \left\| \mathbf{f}'(\mathbf{x})\right\|_{\infty} \leq \varepsilon_1
    \f]
    or criterion 2:
    \f[
    \left\| d\mathbf{x}\right\|_{2} \leq \varepsilon_2 (\varepsilon_2 + \| \mathbf{x} \|_{2} )
    \f]
    is fulfilled. Another possibility is that the search is stopped when the number of function evaluations
    exceeds a threshold (100 per default).

    The default update rule for the inverse of the Hessian matrix is the BFGS algorithm with
    \f$\varepsilon_1 = 10^{-4}\f$ an \f$\varepsilon_2 = 10^{-8}\f$.

  */
  class Newton_Search {
  public:
    //! Default constructor
    Newton_Search();
    //! Destructor
    ~Newton_Search() {};

    //! Set function pointer
    void set_function(double(*function)(const vec&));
    //! Set gradient function pointer
    void set_gradient(vec(*gradient)(const vec&));
    //! Set both function and gradient function pointers
    void set_functions(double(*function)(const vec&), vec(*gradient)(const vec&)) { set_function(function); set_gradient(gradient); }

    //! Set start point \c x for search and approx inverse Hessian at \c x
    void set_start_point(const vec &x, const mat &D);

    //! Set start point \c x for search
    void set_start_point(const vec &x);

    //! Get solution, function value and gradient at solution point
    vec get_solution();

    //! Do the line search
    bool search();
    //! Do the line search and return solution
    bool search(vec &xn);
    //! Set starting point, do the Newton search, and return the solution
    bool search(const vec &x0, vec &xn);

    //! Set stop criterion values
    void set_stop_values(double epsilon_1, double epsilon_2);
    //! Return stop value rho
    double get_epsilon_1() { return stop_epsilon_1; }
    //! Return stop value beta
    double get_epsilon_2() { return stop_epsilon_2; }

    //! Set max number of function evaluations
    void set_max_evaluations(int value);
    //! Return max number of function evaluations
    int get_max_evaluations() { return max_evaluations; }

    //! Set max stepsize
    void set_initial_stepsize(double value);
    //! Return max number of iterations
    double get_initial_stepsize() { return initial_stepsize; }

    //! Set Line search method
    void set_method(const Newton_Search_Method &method);

    //! get function value at solution point
    double get_function_value();
    //! get value of stop criterion 1 at solution point
    double get_stop_1();
    //! get value of stop criterion 2 at solution point
    double get_stop_2();
    //! get number of iterations used to reach solution
    int get_no_iterations();
    //! get number of function evaluations used to reach solution
    int get_no_function_evaluations();

    //! enable trace mode
    void enable_trace() { trace = true; }
    //! disable trace
    void disable_trace() { trace = false; }

    /*! get trace outputs
      \c xvalues are the solutions of every iteration
      \c Fvalues are the function values
      \c ngvalues are the norm(gradient,inf) values
      \c dvalues are the delta values
    */
    void get_trace(Array<vec> & xvalues, vec &Fvalues, vec &ngvalues, vec &dvalues);

  private:
    int n; // dimension of problem, size(x)
    double (*f)(const vec&); // function to minimize
    vec (*df_dx)(const vec&); // df/dx, gradient of f

    // start variables
    vec x_start;
    mat D_start;

    // solution variables
    vec x_end;

    // trace variables
    Array<vec> x_values;
    vec F_values, ng_values, Delta_values;

    Newton_Search_Method method;

    // Parameters
    double initial_stepsize; // opts(1)
    double stop_epsilon_1; // opts(2)
    double stop_epsilon_2; // opt(3)
    int max_evaluations; // opts(4)

    // output parameters
    int no_feval; // number of function evaluations
    int no_iter; // number of iterations
    double F, ng, nh; // function value, stop_1, stop_2 values at solution point

    bool init, finished, trace;
  };



  //! Line Search method
  enum Line_Search_Method {Soft, Exact};

  /*!
    \brief Line Search

    The line search try to minimize the objective function \f$f(\mathbf{x})\f$
    along the direction \f$\mathbf{h}\f$ from the current position \f$\mathbf{x}\f$.

    Hence we look at
    \f[
    \varphi(\alpha) = f(\mathbf{x} + \alpha \mathbf{h})
    \f]
    and try to find an \f$\alpha_s\f$ that minimizes \f$f\f$.

    Two variants are used. Either the soft line search (default) or the exact line
    search.

    The soft line search stops when a point in the acceptable region is found, i.e.
    \f[
    \phi(\alpha_s) \leq \varphi(0) + \alpha_s \rho \varphi'(0)
    \f]
    and
    \f[
    \varphi'(\alpha_s) \geq \beta \varphi'(0),\: \rho < \beta
    \f]
    Default vales are \f$\rho = 10^{-3}\f$ and \f$\beta = 0.99\f$.

    The exact line search
    \f[
    \| \varphi(\alpha_s)\|  \leq \rho \| \varphi'(0) \|
    \f]
    and
    \f[
    b-a \leq \beta b,
    \f]
    where \f$\left[a,b\right]\f$ is the current interval for \f$\alpha_s\f$.
    Default vales are \f$\rho = 10^{-3}\f$ and \f$\beta = 10^{-3}\f$.

    The exact line search can at least in theory give the exact resutl, but it may require
    many extra function evaluations compared to soft line search.
  */
  class Line_Search {
  public:
    //! Default constructor
    Line_Search();
    //! Destructor
    ~Line_Search() {};

    //! Set function pointer
    void set_function(double(*function)(const vec&));
    //! Set gradient function pointer
    void set_gradient(vec(*gradient)(const vec&));
    //! Set both function and gradient function pointers
    void set_functions(double(*function)(const vec&), vec(*gradient)(const vec&)) { set_function(function); set_gradient(gradient); }

    //! Set start point for search
    void set_start_point(const vec &x, double F, const vec &g, const vec &h);

    //! Get solution, function value and gradient at solution point
    void get_solution(vec &xn, double &Fn, vec &gn);

    //! Do the line search
    bool search();
    //! Do the line search and return solution
    bool search(vec &xn, double &Fn, vec &gn);
    //! Set starting point, do the line search, and return the solution
    bool search(const vec &x, double F, const vec &g, const vec &h, vec &xn,
		double &Fn, vec &gn);


    //! return alpha at solution point, xn = x + alpha h
    double get_alpha();
    //! return the slope ratio at solution poin, xn
    double get_slope_ratio();
    //! return number of function evaluations used in search
    int get_no_function_evaluations();


    //! Set stop criterion values
    void set_stop_values(double rho, double beta);
    //! Return stop value rho
    double get_rho() { return stop_rho; }
    //! Return stop value beta
    double get_beta() { return stop_beta; }

    //! Set max number of iterations
    void set_max_iterations(int value);
    //! Return max number of iterations
    int get_max_iterations() { return max_iterations; }

    //! Set max stepsize
    void set_max_stepsize(double value);
    //! Return max number of iterations
    double get_max_stepsize() { return max_stepsize; }

    //! Set Line search method
    void set_method(const Line_Search_Method &method);

    //! enable trace mode
    void enable_trace() { trace = true; }
    //! disable trace
    void disable_trace() { trace = false; }

    /*! get trace outputs
      \c alphavalues are the solutions of every iteration
      \c Fvalues are the function values
      \c dFvalues
    */
    void get_trace(vec &alphavalues, vec &Fvalues, vec &dFvalues);

  private:
    int n; // dimension of problem, size(x)
    double (*f)(const vec&); // function to minimize
    vec (*df_dx)(const vec&); // df/dx, gradient of f

    // start variables
    vec x_start, g_start, h_start;
    double F_start;

    // solution variables
    vec x_end, g_end;
    double F_end;

    // trace variables
    vec alpha_values, F_values, dF_values;

    bool init; // true if functions and starting points are set
    bool finished; // true if functions and starting points are set
    bool trace; // true if trace is enabled

    // Parameters
    Line_Search_Method method;
    double stop_rho; // opts(2)
    double stop_beta; // opts(3)
    int max_iterations; // opts(4)
    double max_stepsize; // opts(5)

    // output parameters
    double alpha; // end value of alpha, info(1)
    double slope_ratio; // slope ratio at xn, info(2)
    int no_feval; // info(3)
  };

  /*!
    \brief Unconstrained minimization

    Unconstrained minimization using a Newton or Quasi-Newton optimization method
    that try to minimize the objective function \f$f(\mathbf{x})\f$ given an initial guess \f$\mathbf{x}\f$.

    The function and the gradient need to be known and supplied.

    The default algorithm is a Quasi-Newton search using BFGS updates of the inverse Hessian matrix.
  */
  vec fminunc(double(*function)(const vec&), vec(*gradient)(const vec&), const vec &x0);

  //@}

} // namespace itpp

#endif // #ifndef NEWTON_SEARCH_H