yat/utility/RootFinderDerivative.h

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

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