yat/utility/RootFinder.h

Code
Comments
Other
Rev Date Author Line
4237 17 Sep 22 peter 1 #ifndef theplu_yat_utility_root_finder_
4237 17 Sep 22 peter 2 #define theplu_yat_utility_root_finder_
4237 17 Sep 22 peter 3
4237 17 Sep 22 peter 4 // $Id$
4237 17 Sep 22 peter 5
4237 17 Sep 22 peter 6 /*
4237 17 Sep 22 peter 7   Copyright (C) 2022 Peter Johansson
4237 17 Sep 22 peter 8
4237 17 Sep 22 peter 9   This file is part of the yat library, https://dev.thep.lu.se/yat
4237 17 Sep 22 peter 10
4237 17 Sep 22 peter 11   The yat library is free software; you can redistribute it and/or
4237 17 Sep 22 peter 12   modify it under the terms of the GNU General Public License as
4237 17 Sep 22 peter 13   published by the Free Software Foundation; either version 3 of the
4237 17 Sep 22 peter 14   License, or (at your option) any later version.
4237 17 Sep 22 peter 15
4237 17 Sep 22 peter 16   The yat library is distributed in the hope that it will be useful,
4237 17 Sep 22 peter 17   but WITHOUT ANY WARRANTY; without even the implied warranty of
4237 17 Sep 22 peter 18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4237 17 Sep 22 peter 19   General Public License for more details.
4237 17 Sep 22 peter 20
4237 17 Sep 22 peter 21   You should have received a copy of the GNU General Public License
4237 17 Sep 22 peter 22   along with yat. If not, see <https://www.gnu.org/licenses/>.
4237 17 Sep 22 peter 23 */
4237 17 Sep 22 peter 24
4237 17 Sep 22 peter 25 #include "univariable/f.h"
4237 17 Sep 22 peter 26 #include "yat_assert.h"
4237 17 Sep 22 peter 27
4237 17 Sep 22 peter 28 #include <gsl/gsl_errno.h>
4237 17 Sep 22 peter 29 #include <gsl/gsl_roots.h>
4237 17 Sep 22 peter 30
4243 18 Sep 22 peter 31 #include <limits>
4237 17 Sep 22 peter 32 #include <memory>
4237 17 Sep 22 peter 33
4237 17 Sep 22 peter 34 namespace theplu {
4237 17 Sep 22 peter 35 namespace yat {
4237 17 Sep 22 peter 36 namespace utility {
4237 17 Sep 22 peter 37
4237 17 Sep 22 peter 38   /**
4237 17 Sep 22 peter 39      Wrapper class around
4237 17 Sep 22 peter 40      <a href="https://www.gnu.org/software/gsl/doc/html/roots.html">
4237 17 Sep 22 peter 41      GSL One-Dimensional Root-Finding</a>
4237 17 Sep 22 peter 42
4237 17 Sep 22 peter 43      \since New in yat 0.20
4237 17 Sep 22 peter 44    */
4237 17 Sep 22 peter 45   class RootFinder
4237 17 Sep 22 peter 46   {
4237 17 Sep 22 peter 47   public:
4237 17 Sep 22 peter 48     /// \brief Copy is not allowed
4237 17 Sep 22 peter 49     RootFinder(const RootFinder&) = delete;
4237 17 Sep 22 peter 50
4237 17 Sep 22 peter 51     /// \brief Assignment is not allowed
4237 17 Sep 22 peter 52     RootFinder& operator=(const RootFinder&) = delete;
4237 17 Sep 22 peter 53
4237 17 Sep 22 peter 54     /**
4237 17 Sep 22 peter 55        \brief Number of epochs (iterations) used in last minimisation.
4237 17 Sep 22 peter 56      */
4237 17 Sep 22 peter 57     unsigned int epochs(void) const;
4237 17 Sep 22 peter 58
4237 17 Sep 22 peter 59     /**
4237 17 Sep 22 peter 60        \return name of the type of root finder
4237 17 Sep 22 peter 61      */
4237 17 Sep 22 peter 62     const char* name(void) const;
4237 17 Sep 22 peter 63
4237 17 Sep 22 peter 64     /// Abstract class defining the interface for functions defining
4237 17 Sep 22 peter 65     /// the stop criteria in the root finding.
4237 17 Sep 22 peter 66     class Stopper
4237 17 Sep 22 peter 67     {
4237 17 Sep 22 peter 68     public:
4237 17 Sep 22 peter 69       /// return true is search should stop
4237 17 Sep 22 peter 70       virtual bool operator()(const gsl_root_fsolver*)=0;
4237 17 Sep 22 peter 71     };
4237 17 Sep 22 peter 72
4237 17 Sep 22 peter 73     /// wrapper around
4237 17 Sep 22 peter 74     /// <a href="https://www.gnu.org/software/gsl/doc/html/roots.html#search-stopping-parameters">gsl_root_test_interval</a>
4237 17 Sep 22 peter 75     class Interval : public Stopper
4237 17 Sep 22 peter 76     {
4237 17 Sep 22 peter 77     public:
4237 17 Sep 22 peter 78       /**
4237 17 Sep 22 peter 79          \param epsabs absolute tolerance
4237 17 Sep 22 peter 80          \param epsrel relative tolerance
4237 17 Sep 22 peter 81       */
4237 17 Sep 22 peter 82       Interval(double epsabs, double epsrel);
4237 17 Sep 22 peter 83
4237 17 Sep 22 peter 84       /**
4237 17 Sep 22 peter 85          \return true if current search interval satisfies
4237 17 Sep 22 peter 86
4237 17 Sep 22 peter 87          \f$ upper - lower < epsabs + epsrel *
4237 17 Sep 22 peter 88          \inf_{lower \le x \le upper} (|x|) \f$
4237 17 Sep 22 peter 89        */
4237 17 Sep 22 peter 90       bool operator()(const gsl_root_fsolver*);
4237 17 Sep 22 peter 91     private:
4237 17 Sep 22 peter 92       double epsabs_;
4237 17 Sep 22 peter 93       double epsrel_;
4237 17 Sep 22 peter 94     };
4237 17 Sep 22 peter 95
4237 17 Sep 22 peter 96     /// wrapper around
4237 17 Sep 22 peter 97     /// <a href="https://www.gnu.org/software/gsl/doc/html/roots.html#search-stopping-parameters">gsl_root_test_delta</a>
4237 17 Sep 22 peter 98     class Delta : public Stopper
4237 17 Sep 22 peter 99     {
4237 17 Sep 22 peter 100     public:
4237 17 Sep 22 peter 101       /**
4237 17 Sep 22 peter 102          \param epsabs absolute tolerance
4237 17 Sep 22 peter 103          \param epsrel relative tolerance
4237 17 Sep 22 peter 104       */
4237 17 Sep 22 peter 105       Delta(double epsabs, double epsrel);
4237 17 Sep 22 peter 106
4237 17 Sep 22 peter 107       /**
4237 17 Sep 22 peter 108          \return true if search sequence x_0, x_1 satisfies
4237 17 Sep 22 peter 109          \f$ |x_1 - x_0| < epsabs + epsrel |x_1| \f$
4237 17 Sep 22 peter 110        */
4237 17 Sep 22 peter 111       bool operator()(const gsl_root_fsolver*);
4237 17 Sep 22 peter 112     private:
4237 17 Sep 22 peter 113       double epsabs_;
4237 17 Sep 22 peter 114       double epsrel_;
4237 17 Sep 22 peter 115       std::unique_ptr<double> prev_;
4237 17 Sep 22 peter 116     };
4237 17 Sep 22 peter 117
4237 17 Sep 22 peter 118     /// wrapper around
4237 17 Sep 22 peter 119     /// <a href="https://www.gnu.org/software/gsl/doc/html/roots.html#search-stopping-parameters">gsl_root_test_residual</a>
4237 17 Sep 22 peter 120     template<class FUNC>
4237 17 Sep 22 peter 121     class Residual : public Stopper
4237 17 Sep 22 peter 122     {
4237 17 Sep 22 peter 123     public:
4237 17 Sep 22 peter 124       /**
4237 17 Sep 22 peter 125          \param func same functor as used in RootFinder::operator().
4237 17 Sep 22 peter 126          \param epsabs absolute tolerance
4237 17 Sep 22 peter 127       */
4237 17 Sep 22 peter 128       Residual(const FUNC& func, double epsabs)
4237 17 Sep 22 peter 129         : func_(func), epsabs_(epsabs)
4237 17 Sep 22 peter 130       {}
4237 17 Sep 22 peter 131
4237 17 Sep 22 peter 132       /**
4237 17 Sep 22 peter 133          \return true if current residual \a f satisfies
4237 17 Sep 22 peter 134          \f$ |f| < epsabs \f$
4237 17 Sep 22 peter 135        */
4237 17 Sep 22 peter 136       bool operator()(const gsl_root_fsolver*);
4237 17 Sep 22 peter 137     private:
4237 17 Sep 22 peter 138       const FUNC& func_;
4237 17 Sep 22 peter 139       double epsabs_;
4237 17 Sep 22 peter 140     };
4237 17 Sep 22 peter 141
4237 17 Sep 22 peter 142     /**
4237 17 Sep 22 peter 143        Find a value \a x within range [\c lower, \c upper] such that
4237 17 Sep 22 peter 144        func(x) = 0. Function keeps iterating until \c stopper returns
4237 17 Sep 22 peter 145        \c true.
4237 17 Sep 22 peter 146      */
4237 17 Sep 22 peter 147     template<class FUNC>
4237 17 Sep 22 peter 148     double operator()(FUNC& func, double lower, double upper,Stopper&& stopper);
4237 17 Sep 22 peter 149
4237 17 Sep 22 peter 150     /**
4237 17 Sep 22 peter 151        Same as operator()(4) but do maximum \c max_epochs iterations.
4237 17 Sep 22 peter 152      */
4237 17 Sep 22 peter 153     template<class FUNC>
4237 17 Sep 22 peter 154     double operator()(FUNC& func, double lower, double upper,
4237 17 Sep 22 peter 155                       Stopper&& stopper, unsigned int max_epochs);
4237 17 Sep 22 peter 156
4237 17 Sep 22 peter 157   protected:
4237 17 Sep 22 peter 158     /**
4237 17 Sep 22 peter 159        \brief Constructor
4237 17 Sep 22 peter 160
4237 17 Sep 22 peter 161        \param t defines type of GSL minimizer
4237 17 Sep 22 peter 162      */
4237 17 Sep 22 peter 163     RootFinder(const gsl_root_fsolver_type* t);
4237 17 Sep 22 peter 164   private:
4237 17 Sep 22 peter 165     const gsl_root_fsolver_type* type_;
4237 17 Sep 22 peter 166     struct GslFree
4237 17 Sep 22 peter 167     {
4237 17 Sep 22 peter 168       void operator()(gsl_root_fsolver*) const;
4237 17 Sep 22 peter 169     };
4237 17 Sep 22 peter 170     std::unique_ptr<gsl_root_fsolver, GslFree> solver_;
4237 17 Sep 22 peter 171     unsigned int epochs_;
4237 17 Sep 22 peter 172   };
4237 17 Sep 22 peter 173
4237 17 Sep 22 peter 174
4237 17 Sep 22 peter 175   // template implementation
4237 17 Sep 22 peter 176   template<class FUNC>
4237 17 Sep 22 peter 177   double RootFinder::operator()(FUNC& func, double lower, double upper,
4237 17 Sep 22 peter 178                                 RootFinder::Stopper&& stopper)
4237 17 Sep 22 peter 179   {
4237 17 Sep 22 peter 180     unsigned int max_epochs = std::numeric_limits<unsigned int>::max();
4237 17 Sep 22 peter 181     return (*this)(func, lower, upper, std::move(stopper), max_epochs);
4237 17 Sep 22 peter 182   }
4237 17 Sep 22 peter 183
4237 17 Sep 22 peter 184
4237 17 Sep 22 peter 185   template<class FUNC>
4237 17 Sep 22 peter 186   double
4237 17 Sep 22 peter 187   RootFinder::operator()(FUNC& func, double lower, double upper,
4237 17 Sep 22 peter 188                          RootFinder::Stopper&& stopper,
4237 17 Sep 22 peter 189                          unsigned int max_epochs)
4237 17 Sep 22 peter 190   {
4237 17 Sep 22 peter 191     gsl_function gsl_func;
4237 17 Sep 22 peter 192     gsl_func.function = univariable::f<FUNC>;
4237 17 Sep 22 peter 193     gsl_func.params = &func;
4237 17 Sep 22 peter 194     gsl_root_fsolver_set(solver_.get(), &gsl_func, lower, upper);
4237 17 Sep 22 peter 195
4237 17 Sep 22 peter 196     for (epochs_=0; epochs_<max_epochs && !stopper(solver_.get()); ++epochs_) {
4237 17 Sep 22 peter 197       int status = gsl_root_fsolver_iterate(solver_.get());
4237 17 Sep 22 peter 198       if (status) {
4237 17 Sep 22 peter 199         std::ostringstream ss;
4237 17 Sep 22 peter 200         ss << "RootFinder: at x = "
4237 17 Sep 22 peter 201            << gsl_root_fsolver_root(solver_.get());
4237 17 Sep 22 peter 202         throw GSL_error(ss.str(), status);
4237 17 Sep 22 peter 203
4237 17 Sep 22 peter 204       }
4237 17 Sep 22 peter 205     }
4237 17 Sep 22 peter 206
4237 17 Sep 22 peter 207     return gsl_root_fsolver_root(solver_.get());
4237 17 Sep 22 peter 208   }
4237 17 Sep 22 peter 209
4237 17 Sep 22 peter 210
4237 17 Sep 22 peter 211   template<class FUNC>
4237 17 Sep 22 peter 212   bool RootFinder::Residual<FUNC>::operator()(const gsl_root_fsolver* s)
4237 17 Sep 22 peter 213   {
4237 17 Sep 22 peter 214     double x = gsl_root_fsolver_root(s);
4237 17 Sep 22 peter 215     double f = func_(x);
4237 17 Sep 22 peter 216     return gsl_root_test_residual(f, epsabs_) == GSL_SUCCESS;
4237 17 Sep 22 peter 217   }
4237 17 Sep 22 peter 218
4237 17 Sep 22 peter 219 }}}
4237 17 Sep 22 peter 220 #endif

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