yat/utility/Minimizer.h

Code
Comments
Other
Rev Date Author Line
4222 10 Sep 22 peter 1 #ifndef theplu_yat_utility_minimizer_
4222 10 Sep 22 peter 2 #define theplu_yat_utility_minimizer_
4222 10 Sep 22 peter 3
4222 10 Sep 22 peter 4 // $Id$
4222 10 Sep 22 peter 5
4222 10 Sep 22 peter 6 /*
4222 10 Sep 22 peter 7   Copyright (C) 2022 Peter Johansson
4222 10 Sep 22 peter 8
4222 10 Sep 22 peter 9   This file is part of the yat library, https://dev.thep.lu.se/yat
4222 10 Sep 22 peter 10
4222 10 Sep 22 peter 11   The yat library is free software; you can redistribute it and/or
4222 10 Sep 22 peter 12   modify it under the terms of the GNU General Public License as
4222 10 Sep 22 peter 13   published by the Free Software Foundation; either version 3 of the
4222 10 Sep 22 peter 14   License, or (at your option) any later version.
4222 10 Sep 22 peter 15
4222 10 Sep 22 peter 16   The yat library is distributed in the hope that it will be useful,
4222 10 Sep 22 peter 17   but WITHOUT ANY WARRANTY; without even the implied warranty of
4222 10 Sep 22 peter 18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4222 10 Sep 22 peter 19   General Public License for more details.
4222 10 Sep 22 peter 20
4222 10 Sep 22 peter 21   You should have received a copy of the GNU General Public License
4222 10 Sep 22 peter 22   along with yat. If not, see <https://www.gnu.org/licenses/>.
4222 10 Sep 22 peter 23 */
4222 10 Sep 22 peter 24
4222 10 Sep 22 peter 25 #include "Exception.h"
4222 10 Sep 22 peter 26 #include "yat_assert.h"
4237 17 Sep 22 peter 27 #include "univariable/f.h"
4222 10 Sep 22 peter 28
4222 10 Sep 22 peter 29 #include <gsl/gsl_min.h>
4222 10 Sep 22 peter 30 #include <gsl/gsl_errno.h>
4222 10 Sep 22 peter 31
4222 10 Sep 22 peter 32 #include <limits>
4222 10 Sep 22 peter 33 #include <memory>
4222 10 Sep 22 peter 34 #include <sstream>
4222 10 Sep 22 peter 35
4222 10 Sep 22 peter 36 namespace theplu {
4222 10 Sep 22 peter 37 namespace yat {
4222 10 Sep 22 peter 38 namespace utility {
4222 10 Sep 22 peter 39
4222 10 Sep 22 peter 40   /**
4222 10 Sep 22 peter 41      \brief Wrapper class around gsl_min_fminimizer in GSL.
4222 10 Sep 22 peter 42
4222 10 Sep 22 peter 43      Class is the abstract base class for several minimisation
4222 10 Sep 22 peter 44      algorithms.
4222 10 Sep 22 peter 45
4222 10 Sep 22 peter 46      \since New in yat 0.20
4222 10 Sep 22 peter 47    */
4222 10 Sep 22 peter 48   class Minimizer
4222 10 Sep 22 peter 49   {
4222 10 Sep 22 peter 50   public:
4222 10 Sep 22 peter 51     /// \brief Copy is not allowed
4222 10 Sep 22 peter 52     Minimizer(const Minimizer&) = delete;
4222 10 Sep 22 peter 53
4222 10 Sep 22 peter 54     /// \brief Assignment is not allowed
4222 10 Sep 22 peter 55     Minimizer& operator=(const Minimizer&) = delete;
4222 10 Sep 22 peter 56
4222 10 Sep 22 peter 57     /**
4222 10 Sep 22 peter 58        \brief Number of epochs (iterations) used in last minimisation.
4222 10 Sep 22 peter 59      */
4222 10 Sep 22 peter 60     unsigned int epochs(void) const;
4222 10 Sep 22 peter 61
4247 24 Sep 22 peter 62     /// Abstract class defining the interface for functions defining
4247 24 Sep 22 peter 63     /// the stop criteria in the root finding.
4247 24 Sep 22 peter 64     class Stopper
4247 24 Sep 22 peter 65     {
4247 24 Sep 22 peter 66     public:
4247 24 Sep 22 peter 67       /// return true is search should stop
4247 24 Sep 22 peter 68       virtual bool operator()(const gsl_min_fminimizer*)=0;
4247 24 Sep 22 peter 69     };
4222 10 Sep 22 peter 70
4247 24 Sep 22 peter 71     /// wrapper around
4247 24 Sep 22 peter 72     /// <a href="https://www.gnu.org/software/gsl/doc/html/min.html#stopping-parameters">gsl_min_test_interval</a>
4247 24 Sep 22 peter 73     class Interval : public Stopper
4247 24 Sep 22 peter 74     {
4247 24 Sep 22 peter 75     public:
4247 24 Sep 22 peter 76       /**
4247 24 Sep 22 peter 77          \param epsabs absolute tolerance
4247 24 Sep 22 peter 78          \param epsrel relative tolerance
4247 24 Sep 22 peter 79       */
4247 24 Sep 22 peter 80       Interval(double epsabs, double epsrel);
4222 10 Sep 22 peter 81
4247 24 Sep 22 peter 82       /**
4247 24 Sep 22 peter 83          \return true if current search interval satisfies
4222 10 Sep 22 peter 84
4247 24 Sep 22 peter 85          \f$ upper - lower < epsabs + epsrel *
4247 24 Sep 22 peter 86          \inf_{lower \le x \le upper} (|x|) \f$
4247 24 Sep 22 peter 87        */
4247 24 Sep 22 peter 88       bool operator()(const gsl_min_fminimizer*);
4247 24 Sep 22 peter 89     private:
4247 24 Sep 22 peter 90       double epsabs_;
4247 24 Sep 22 peter 91       double epsrel_;
4247 24 Sep 22 peter 92     };
4222 10 Sep 22 peter 93
4222 10 Sep 22 peter 94     /**
4222 10 Sep 22 peter 95        Function finds an \c x that minimizes the function defined by
4247 24 Sep 22 peter 96        \c func. It calls gsl_min_fminimizer_iterate until there is no
4247 24 Sep 22 peter 97        progress as defined by \c stopper.
4222 10 Sep 22 peter 98
4222 10 Sep 22 peter 99        Type Requirements:
4222 10 Sep 22 peter 100        - \c FUNC must have an operator double operator()(double x)
4222 10 Sep 22 peter 101      */
4222 10 Sep 22 peter 102     template<class FUNC>
4247 24 Sep 22 peter 103     double operator()(FUNC& func, double guess, double min, double max,
4247 24 Sep 22 peter 104                       Stopper&& stopper);
4222 10 Sep 22 peter 105
4222 10 Sep 22 peter 106     /**
4222 10 Sep 22 peter 107        Same as operator()(FUNC& func, double guess, double min, double
4247 24 Sep 22 peter 108        max, Stopper&& stopper); but run at maximum \c max_epochs
4247 24 Sep 22 peter 109        epochs (iterations). It is the responsibility of the caller to
4247 24 Sep 22 peter 110        check if maximum number of epochs has been reached, for
4247 24 Sep 22 peter 111        example, with the epochs(void function.
4222 10 Sep 22 peter 112      */
4222 10 Sep 22 peter 113     template<class FUNC>
4222 10 Sep 22 peter 114     double operator()(FUNC& func, double guess, double min, double max,
4247 24 Sep 22 peter 115                       Stopper&& stopper, unsigned int max_epochs);
4222 10 Sep 22 peter 116
4222 10 Sep 22 peter 117   protected:
4222 10 Sep 22 peter 118     /**
4222 10 Sep 22 peter 119        \brief Constructor
4222 10 Sep 22 peter 120
4222 10 Sep 22 peter 121        \param t defines type of GSL minimizer
4222 10 Sep 22 peter 122      */
4222 10 Sep 22 peter 123     Minimizer(const gsl_min_fminimizer_type* t);
4222 10 Sep 22 peter 124   private:
4222 10 Sep 22 peter 125     const gsl_min_fminimizer_type* type_;
4222 10 Sep 22 peter 126     struct GslFree
4222 10 Sep 22 peter 127     {
4222 10 Sep 22 peter 128       void operator()(gsl_min_fminimizer*) const;
4222 10 Sep 22 peter 129     };
4222 10 Sep 22 peter 130     std::unique_ptr<gsl_min_fminimizer, GslFree> solver_;
4222 10 Sep 22 peter 131     unsigned int epochs_;
4222 10 Sep 22 peter 132     double epsabs_;
4222 10 Sep 22 peter 133     double epsrel_;
4222 10 Sep 22 peter 134   };
4222 10 Sep 22 peter 135
4222 10 Sep 22 peter 136
4222 10 Sep 22 peter 137   // template implementation
4222 10 Sep 22 peter 138
4222 10 Sep 22 peter 139   template<class FUNC>
4247 24 Sep 22 peter 140   double Minimizer::operator()(FUNC& func, double guess, double min,double max,
4247 24 Sep 22 peter 141                                Minimizer::Stopper&& stopper)
4222 10 Sep 22 peter 142   {
4222 10 Sep 22 peter 143     unsigned int max_epochs = std::numeric_limits<unsigned int>::max();
4247 24 Sep 22 peter 144     return (*this)(func, guess, min, max, std::move(stopper), max_epochs);
4222 10 Sep 22 peter 145   }
4222 10 Sep 22 peter 146
4222 10 Sep 22 peter 147
4222 10 Sep 22 peter 148   template<class FUNC>
4222 10 Sep 22 peter 149   double Minimizer::operator()(FUNC& func, double guess, double min,
4247 24 Sep 22 peter 150                                double max, Minimizer::Stopper&& stopper,
4247 24 Sep 22 peter 151                                unsigned int max_epochs)
4222 10 Sep 22 peter 152   {
4222 10 Sep 22 peter 153     gsl_function gsl_func;
4237 17 Sep 22 peter 154     gsl_func.function = univariable::f<FUNC>;
4222 10 Sep 22 peter 155     gsl_func.params = &func;
4222 10 Sep 22 peter 156     gsl_min_fminimizer_set(solver_.get(), &gsl_func, guess, min, max);
4222 10 Sep 22 peter 157
4222 10 Sep 22 peter 158     int status = GSL_CONTINUE;
4247 24 Sep 22 peter 159     for (epochs_=0; epochs_<max_epochs && !stopper(solver_.get()); ++epochs_) {
4222 10 Sep 22 peter 160       status = gsl_min_fminimizer_iterate(solver_.get());
4222 10 Sep 22 peter 161       if (status) {
4237 17 Sep 22 peter 162         // GSL returns GSL_FAILURE when the function cannot be improved
4222 10 Sep 22 peter 163         if (status == GSL_FAILURE)
4222 10 Sep 22 peter 164           break;
4222 10 Sep 22 peter 165         if (status == GSL_EBADFUNC) {
4222 10 Sep 22 peter 166           std::ostringstream ss;
4222 10 Sep 22 peter 167           ss << "Minimizer: at x = "
4222 10 Sep 22 peter 168              << gsl_min_fminimizer_x_minimum(solver_.get());
4222 10 Sep 22 peter 169           throw GSL_error(ss.str(), status);
4222 10 Sep 22 peter 170         }
4222 10 Sep 22 peter 171       }
4222 10 Sep 22 peter 172     }
4222 10 Sep 22 peter 173
4222 10 Sep 22 peter 174     return gsl_min_fminimizer_x_minimum(solver_.get());
4222 10 Sep 22 peter 175   }
4222 10 Sep 22 peter 176
4222 10 Sep 22 peter 177 }}}
4222 10 Sep 22 peter 178 #endif

Generated on Sat Oct 14 02:24:54 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)