source: trunk/yat/utility/MultiMinimizer.h @ 4172

#ifndef theplu_yat_utility_multi_minimizer_
#define theplu_yat_utility_multi_minimizer_

// $Id: MultiMinimizer.h 4172 2022-05-12 02:15:31Z peter $

/*
  Copyright (C) 2022 Peter Johansson

  This file is part of the yat library, https://dev.thep.lu.se/yat

  The yat library 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.

  The yat library 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 yat. If not, see <https://www.gnu.org/licenses/>.
*/

#include "multiminimizer/f.h"

#include "Vector.h"
#include "yat_assert.h"

#include <gsl/gsl_multimin.h>

#include <memory>

namespace theplu {
namespace yat {
namespace utility {

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

     Class is the abstract base class for several minimisation
     algorithms not using the derivative.

     \since New in yat 0.20
   */
  class MultiMinimizer
  {
  public:
    /// \brief Copy is not allowed
    MultiMinimizer(const MultiMinimizer&) = delete;

    /// \brief Assignment is not allowed
    MultiMinimizer& operator=(const MultiMinimizer&) = delete;

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

    /**
       \brief The size of the first step.

       \see parameter \c step_size in gsl_multimin_fminimizer_set
     */
    const yat::utility::Vector& step(void) const;

    /**
       \brief Modify the initial step size.

       The size of the vector cannot be modified but must equal the
       size set in constructor.
     */
    yat::utility::Vector& step(void);

    /**
       Function finds an \c x that minimizes the function defined by
       \c func. It calls gsl_multimin_fminimizer_iterate until either
       \c GSL_ENOPROG is returned or the size is smaller than \c
       epsabs, as tested by gsl_multimin_test_size.

       Type Requirements:
       - \c FUNC must have an operator
       double operator()(const VectorBase& x)
     */
    template<class FUNC>
    void operator()(yat::utility::VectorMutable& x, FUNC& func,
                    double epsabs);

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

       \param t defines type of GSL minimizer
       \param size Number of dimensions in the input space.
     */
    MultiMinimizer(const gsl_multimin_fminimizer_type* t, size_t size);
  private:
    const gsl_multimin_fminimizer_type* type_;
    size_t size_;
    struct GslFree
    {
      void operator()(gsl_multimin_fminimizer*) const;
    };
    std::unique_ptr<gsl_multimin_fminimizer, GslFree> solver_;
    yat::utility::Vector step_size_;
    unsigned int epochs_;
  };


  // template implementation

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


  template<class FUNC>
  void MultiMinimizer::operator()(yat::utility::VectorMutable& x,
                                  FUNC& func, double epsabs,
                                  unsigned int max_epochs)
  {
    YAT_ASSERT(size_ == x.size());
    gsl_multimin_function gsl_func;
    gsl_func.n = x.size();
    gsl_func.f = multiminimizer::f<FUNC>;
    gsl_func.params = &func;

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

    int status = 0;
    for (epochs_=0; epochs_<max_epochs; ++epochs_) {
      status = gsl_multimin_fminimizer_iterate(solver_.get());
      if (status) {
        if (status == GSL_ENOPROG)
          break;
        throw yat::utility::GSL_error("MultiMinimizer", status);
      }
      double size = gsl_multimin_fminimizer_size(solver_.get());
      if (gsl_multimin_test_size(size, epsabs) != GSL_CONTINUE)
        break;
    }

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

}}}
#endif