yat/utility/MultiMinimizer.h

Code
Comments
Other
Rev Date Author Line
4172 12 May 22 peter 1 #ifndef theplu_yat_utility_multi_minimizer_
4172 12 May 22 peter 2 #define theplu_yat_utility_multi_minimizer_
4172 12 May 22 peter 3
4172 12 May 22 peter 4 // $Id$
4172 12 May 22 peter 5
4172 12 May 22 peter 6 /*
4172 12 May 22 peter 7   Copyright (C) 2022 Peter Johansson
4172 12 May 22 peter 8
4172 12 May 22 peter 9   This file is part of the yat library, https://dev.thep.lu.se/yat
4172 12 May 22 peter 10
4172 12 May 22 peter 11   The yat library is free software; you can redistribute it and/or
4172 12 May 22 peter 12   modify it under the terms of the GNU General Public License as
4172 12 May 22 peter 13   published by the Free Software Foundation; either version 3 of the
4172 12 May 22 peter 14   License, or (at your option) any later version.
4172 12 May 22 peter 15
4172 12 May 22 peter 16   The yat library is distributed in the hope that it will be useful,
4172 12 May 22 peter 17   but WITHOUT ANY WARRANTY; without even the implied warranty of
4172 12 May 22 peter 18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4172 12 May 22 peter 19   General Public License for more details.
4172 12 May 22 peter 20
4172 12 May 22 peter 21   You should have received a copy of the GNU General Public License
4172 12 May 22 peter 22   along with yat. If not, see <https://www.gnu.org/licenses/>.
4172 12 May 22 peter 23 */
4172 12 May 22 peter 24
4252 18 Nov 22 peter 25 #include "Exception.h"
4252 18 Nov 22 peter 26 #include "multivariable/f.h"
4172 12 May 22 peter 27
4172 12 May 22 peter 28 #include "Vector.h"
4172 12 May 22 peter 29 #include "yat_assert.h"
4172 12 May 22 peter 30
4172 12 May 22 peter 31 #include <gsl/gsl_multimin.h>
4172 12 May 22 peter 32
4172 12 May 22 peter 33 #include <memory>
4172 12 May 22 peter 34
4172 12 May 22 peter 35 namespace theplu {
4172 12 May 22 peter 36 namespace yat {
4172 12 May 22 peter 37 namespace utility {
4172 12 May 22 peter 38
4172 12 May 22 peter 39   /**
4172 12 May 22 peter 40      \brief Wrapper class around gsl_multimin_fminimizer in GSL.
4172 12 May 22 peter 41
4172 12 May 22 peter 42      Class is the abstract base class for several minimisation
4172 12 May 22 peter 43      algorithms not using the derivative.
4172 12 May 22 peter 44
4172 12 May 22 peter 45      \since New in yat 0.20
4172 12 May 22 peter 46    */
4172 12 May 22 peter 47   class MultiMinimizer
4172 12 May 22 peter 48   {
4172 12 May 22 peter 49   public:
4172 12 May 22 peter 50     /// \brief Copy is not allowed
4172 12 May 22 peter 51     MultiMinimizer(const MultiMinimizer&) = delete;
4172 12 May 22 peter 52
4172 12 May 22 peter 53     /// \brief Assignment is not allowed
4172 12 May 22 peter 54     MultiMinimizer& operator=(const MultiMinimizer&) = delete;
4172 12 May 22 peter 55
4172 12 May 22 peter 56     /**
4172 12 May 22 peter 57        \brief Number of epochs (iterations) used in last minimisation.
4172 12 May 22 peter 58      */
4172 12 May 22 peter 59     unsigned int epochs(void) const;
4172 12 May 22 peter 60
4172 12 May 22 peter 61     /**
4172 12 May 22 peter 62        \brief The size of the first step.
4172 12 May 22 peter 63
4172 12 May 22 peter 64        \see parameter \c step_size in gsl_multimin_fminimizer_set
4172 12 May 22 peter 65      */
4338 14 Apr 23 peter 66     const Vector& step(void) const;
4172 12 May 22 peter 67
4172 12 May 22 peter 68     /**
4172 12 May 22 peter 69        \brief Modify the initial step size.
4172 12 May 22 peter 70
4172 12 May 22 peter 71        The size of the vector cannot be modified but must equal the
4172 12 May 22 peter 72        size set in constructor.
4172 12 May 22 peter 73      */
4338 14 Apr 23 peter 74     Vector& step(void);
4172 12 May 22 peter 75
4252 18 Nov 22 peter 76     /// Abstract class defining the interface for functions
4252 18 Nov 22 peter 77     /// determining when the minimisation stops.
4252 18 Nov 22 peter 78     class Stopper
4252 18 Nov 22 peter 79     {
4252 18 Nov 22 peter 80     public:
4252 18 Nov 22 peter 81       /// return true when search should stop
4252 18 Nov 22 peter 82       virtual bool operator()(const gsl_multimin_fminimizer*)=0;
4252 18 Nov 22 peter 83     };
4252 18 Nov 22 peter 84
4252 18 Nov 22 peter 85
4172 12 May 22 peter 86     /**
4252 18 Nov 22 peter 87        wrapper around
4252 18 Nov 22 peter 88        <a href="https://www.gnu.org/software/gsl/doc/html/multimin.html#stopping-criteria">gsl_multimin_test_size</a>
4252 18 Nov 22 peter 89     */
4252 18 Nov 22 peter 90     class Size : public Stopper
4252 18 Nov 22 peter 91     {
4252 18 Nov 22 peter 92     public:
4252 18 Nov 22 peter 93       /// \param epsabs tolerance
4252 18 Nov 22 peter 94       explicit Size(double epsabs);
4252 18 Nov 22 peter 95
4252 18 Nov 22 peter 96       /**
4252 18 Nov 22 peter 97          \return \c true if <a
4252 18 Nov 22 peter 98          href="https://www.gnu.org/software/gsl/doc/html/multimin.html#stopping-criteria">gsl_multimin_test_size</a>(size,
4252 18 Nov 22 peter 99          epsabs) does not return \c GSL_CONTINUE, where \c size is
4252 18 Nov 22 peter 100          defined by \c gsl_multimin_fminimizer_size and \c epsabs is
4252 18 Nov 22 peter 101          defined in constructor.
4252 18 Nov 22 peter 102        */
4252 18 Nov 22 peter 103       bool operator()(const gsl_multimin_fminimizer*);
4252 18 Nov 22 peter 104     private:
4252 18 Nov 22 peter 105       double epsabs_;
4252 18 Nov 22 peter 106     };
4252 18 Nov 22 peter 107
4252 18 Nov 22 peter 108     /**
4172 12 May 22 peter 109        Function finds an \c x that minimizes the function defined by
4172 12 May 22 peter 110        \c func. It calls gsl_multimin_fminimizer_iterate until either
4252 18 Nov 22 peter 111        \c GSL_ENOPROG is returned or the \c stopper returns true.
4172 12 May 22 peter 112
4172 12 May 22 peter 113        Type Requirements:
4172 12 May 22 peter 114        - \c FUNC must have an operator
4172 12 May 22 peter 115        double operator()(const VectorBase& x)
4172 12 May 22 peter 116      */
4172 12 May 22 peter 117     template<class FUNC>
4338 14 Apr 23 peter 118     void operator()(VectorMutable& x, FUNC& func,
4252 18 Nov 22 peter 119                     Stopper&& stopper);
4172 12 May 22 peter 120
4172 12 May 22 peter 121     /**
4172 12 May 22 peter 122        Same as
4338 14 Apr 23 peter 123        operator()(VectorMutable&, FUNC& func, Stopper&&);
4172 12 May 22 peter 124        but at maximum run \c max_epochs epochs (iterations).
4172 12 May 22 peter 125      */
4172 12 May 22 peter 126     template<class FUNC>
4338 14 Apr 23 peter 127     void operator()(VectorMutable&, FUNC& func,
4252 18 Nov 22 peter 128                     Stopper&& stopper, unsigned int max_epochs);
4172 12 May 22 peter 129   protected:
4172 12 May 22 peter 130     /**
4172 12 May 22 peter 131        \brief Constructor
4172 12 May 22 peter 132
4172 12 May 22 peter 133        \param t defines type of GSL minimizer
4172 12 May 22 peter 134        \param size Number of dimensions in the input space.
4172 12 May 22 peter 135      */
4172 12 May 22 peter 136     MultiMinimizer(const gsl_multimin_fminimizer_type* t, size_t size);
4172 12 May 22 peter 137   private:
4172 12 May 22 peter 138     const gsl_multimin_fminimizer_type* type_;
4172 12 May 22 peter 139     size_t size_;
4172 12 May 22 peter 140     struct GslFree
4172 12 May 22 peter 141     {
4172 12 May 22 peter 142       void operator()(gsl_multimin_fminimizer*) const;
4172 12 May 22 peter 143     };
4172 12 May 22 peter 144     std::unique_ptr<gsl_multimin_fminimizer, GslFree> solver_;
4338 14 Apr 23 peter 145     Vector step_size_;
4172 12 May 22 peter 146     unsigned int epochs_;
4172 12 May 22 peter 147   };
4172 12 May 22 peter 148
4172 12 May 22 peter 149
4172 12 May 22 peter 150   // template implementation
4172 12 May 22 peter 151
4172 12 May 22 peter 152   template<class FUNC>
4338 14 Apr 23 peter 153   void MultiMinimizer::operator()(VectorMutable& x,
4252 18 Nov 22 peter 154                                   FUNC& func,
4252 18 Nov 22 peter 155                                   MultiMinimizer::Stopper&& stopper)
4172 12 May 22 peter 156   {
4172 12 May 22 peter 157     unsigned int max_epochs = std::numeric_limits<unsigned int>::max();
4252 18 Nov 22 peter 158     (*this)(x, func, std::move(stopper), max_epochs);
4172 12 May 22 peter 159   }
4172 12 May 22 peter 160
4172 12 May 22 peter 161
4172 12 May 22 peter 162   template<class FUNC>
4338 14 Apr 23 peter 163   void MultiMinimizer::operator()(VectorMutable& x,
4252 18 Nov 22 peter 164                                   FUNC& func,
4252 18 Nov 22 peter 165                                   MultiMinimizer::Stopper&& stopper,
4172 12 May 22 peter 166                                   unsigned int max_epochs)
4172 12 May 22 peter 167   {
4172 12 May 22 peter 168     YAT_ASSERT(size_ == x.size());
4172 12 May 22 peter 169     gsl_multimin_function gsl_func;
4172 12 May 22 peter 170     gsl_func.n = x.size();
4252 18 Nov 22 peter 171     gsl_func.f = multivariable::f<FUNC>;
4172 12 May 22 peter 172     gsl_func.params = &func;
4172 12 May 22 peter 173
4172 12 May 22 peter 174     gsl_multimin_fminimizer_set(solver_.get(), &gsl_func, x.gsl_vector_p(),
4172 12 May 22 peter 175                                 step_size_.gsl_vector_p());
4172 12 May 22 peter 176
4172 12 May 22 peter 177     int status = 0;
4172 12 May 22 peter 178     for (epochs_=0; epochs_<max_epochs; ++epochs_) {
4172 12 May 22 peter 179       status = gsl_multimin_fminimizer_iterate(solver_.get());
4172 12 May 22 peter 180       if (status) {
4172 12 May 22 peter 181         if (status == GSL_ENOPROG)
4172 12 May 22 peter 182           break;
4252 18 Nov 22 peter 183         throw GSL_error("MultiMinimizer", status);
4172 12 May 22 peter 184       }
4252 18 Nov 22 peter 185       if (stopper(solver_.get()))
4172 12 May 22 peter 186         break;
4252 18 Nov 22 peter 187       //double size = gsl_multimin_fminimizer_size(solver_.get());
4252 18 Nov 22 peter 188       //if (gsl_multimin_test_size(size, epsabs) != GSL_CONTINUE)
4252 18 Nov 22 peter 189       //break;
4172 12 May 22 peter 190     }
4172 12 May 22 peter 191
4172 12 May 22 peter 192     // copy result to passed x
4338 14 Apr 23 peter 193     VectorConstView view(solver_->x);
4172 12 May 22 peter 194     x = view;
4172 12 May 22 peter 195   }
4172 12 May 22 peter 196
4172 12 May 22 peter 197 }}}
4172 12 May 22 peter 198 #endif