source: trunk/yat/utility/MultiMinimizerDerivative.h @ 4252

#ifndef theplu_yat_utility_multi_minimizer_derivative
#define theplu_yat_utility_multi_minimizer_derivative

// $Id: MultiMinimizerDerivative.h 4252 2022-11-18 02:54:04Z peter $
//
// Copyright (C) 2022 Peter Johansson
//
// 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 3 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, see <https://www.gnu.org/licenses/>.

#include "multivariable/df.h"
#include "multivariable/f.h"

#include <yat/utility/Vector.h>
#include <yat/utility/VectorMutable.h>
#include <yat/utility/version.h>

#include <gsl/gsl_multimin.h>

#include <cassert>
#include <memory>

namespace theplu {
namespace yat {
namespace utility {

  /**
     \brief Wrapper class around gsl_multimin_fdfminimizer in GSL.

     Class is the abstract base class for several minimisation
     algorithms using the derivative.
  */
  class MultiMinimizerDerivative
  {
  public:
    /// Copy is not allowed
    MultiMinimizerDerivative(const MultiMinimizerDerivative&) = delete;

    /// Assignment not allowed
    MultiMinimizerDerivative&
    operator=(const MultiMinimizerDerivative&) = delete;

    /**
       \brief Number of epochs (iterations) used in last minimisation.
     */
    unsigned int epochs(void) const;

    /**
       The size of the first trial step is given by \c step_size

       \see gsl_multimin_fdfminimizer_set
     */
    double step_size(void) const;

    /**
       The size of the first trial step is given by \c step_size

       \see gsl_multimin_fdfminimizer_set
     */
    void step_size(double ss);

    /**
       The accuracy of the line minimization is specified by \c
       tolerance

       \see <a href="https://www.gnu.org/software/gsl/doc/html/multimin.html#initializing-the-multidimensional-minimizer">gsl_multimin_fdfminimizer_set</a>
     */
    double tolerance(void) const;

    /**
       The accuracy of the line minimization is specified by \c
       tolerance

       \see <a href="https://www.gnu.org/software/gsl/doc/html/multimin.html#initializing-the-multidimensional-minimizer">gsl_multimin_fdfminimizer_set</a>
     */
    void tolerance(double tol);

    /// Abstract class defining the interface for functions defining
    /// the stop criteria in the minimization
    class Stopper
    {
    public:
      /// return true is search should stop
      virtual bool operator()(const gsl_multimin_fdfminimizer*)=0;
    };


    /// wrapper around
    /// <a href="https://www.gnu.org/software/gsl/doc/html/multimin.html#stopping-criteria">gsl_multimin_test_gradient</a>
    class Gradient : public Stopper
    {
    public:
      /// \param epsabs tolerance
      explicit Gradient(double epsabs);

      /**
         \return \c true if <a
         href="https://www.gnu.org/software/gsl/doc/html/multimin.html#stopping-criteria">gsl_multimin_test_gradient</a>(gradient,
         epsabs) does not return \c GSL_CONTINUE, where \c gradient is
         defined by \c gsl_multimin_fdfminimizer_gradient and \c
         epsabs is defined in constructor.
       */
      bool operator()(const gsl_multimin_fdfminimizer*);
    private:
      double epsabs_;
    };


    /**
       Function finds an \c x that minimizes the function defined by
       \c func. It calls gsl_multimin_fdfminimizer_iterate until
       either \c GSL_ENOPROG is returned or the norm of the gradient
       is smaller than \c epsabs, as tested by
       gsl_multimin_test_gradient.

       Type Requirements:
       - \c FUNC must have an operator defining the function
       double operator()(const VectorBase& x)
       - \c FUNC must have an operator defining the gradient
       void operator()(const VectorBase& x, VectorMutable&)
    */
    template<class FUNC>
    void operator()(yat::utility::VectorMutable&, FUNC& func,
                    MultiMinimizerDerivative::Stopper&& stopper);

    /**
       Same as
       operator()(yat::utility::VectorMutable&, FUNC& func, double epsabs)
       but iterate at maximum \c max_epochs iterations.
     */
    template<class FUNC>
    void operator()(yat::utility::VectorMutable&, FUNC& func,
                    MultiMinimizerDerivative::Stopper&& stopper,
                    unsigned int max_epochs);
  protected:
    /**
       \brief Constructor

       \param t defines type of GSL minimizer
       \param size Number of dimensions in the input space.
     */
    MultiMinimizerDerivative(const gsl_multimin_fdfminimizer_type* t,
                             size_t size);
  private:
    const gsl_multimin_fdfminimizer_type* type_;
    unsigned int epochs_;
    size_t size_;
    struct GslFree
    {
      void operator()(gsl_multimin_fdfminimizer*) const;
    };
    std::unique_ptr<gsl_multimin_fdfminimizer, GslFree> solver_;
    double step_size_;
    double tol_;
  };


  // template implementation

  template<class FUNC>
  void
  MultiMinimizerDerivative::operator()(yat::utility::VectorMutable& x,
                                       FUNC& func,
                                       MultiMinimizerDerivative::Stopper&& stopper)
  {
    unsigned int max_epochs = std::numeric_limits<unsigned int>::max();
    (*this)(x, func, std::move(stopper), max_epochs);
  }


  template<class FUNC>
  void
  MultiMinimizerDerivative::operator()(yat::utility::VectorMutable& x,
                                       FUNC& func,
                                       MultiMinimizerDerivative::Stopper&& stopper,
                                       unsigned int max_epochs)
  {
    assert(size_ == x.size());
    gsl_multimin_function_fdf gsl_func;
    gsl_func.n = size_;
    gsl_func.f = multivariable::f<FUNC>;
    gsl_func.df = multivariable::df<FUNC>;
    gsl_func.fdf = multivariable::fdf<FUNC>;
    gsl_func.params = &func;

    gsl_multimin_fdfminimizer_set(solver_.get(), &gsl_func, x.gsl_vector_p(),
                                  step_size_, tol_);

    int status = 0;
    for (epochs_=0; epochs_<max_epochs; ++epochs_) {
      status = gsl_multimin_fdfminimizer_iterate(solver_.get());
      if (status) {
        if (status == GSL_ENOPROG)
          break;
        throw yat::utility::GSL_error("MultiMinimizerDerivative", status);
      }

      if (stopper(solver_.get()))
        break;
      //if (gsl_multimin_test_gradient(solver_->gradient,epsabs) != GSL_CONTINUE)
      //break;
    }

    // copy result to passed x
    yat::utility::VectorConstView view(solver_->x);
    x = view;
  }

}}}
#endif