root/Math/GSLMultiRootFinder.h

0001 // @(#)root/mathmore:$Id$
0002 // Author: L. Moneta  03/2011
0003
0004  /**********************************************************************
0005   *                                                                    *
0006   * Copyright (c) 2004 ROOT Foundation,  CERN/PH-SFT                   *
0007   *                                                                    *
0008   * This library is free software; you can redistribute it and/or      *
0009   * modify it under the terms of the GNU General Public License        *
0010   * as published by the Free Software Foundation; either version 2     *
0011   * of the License, or (at your option) any later version.             *
0012   *                                                                    *
0013   * This library is distributed in the hope that it will be useful,    *
0014   * but WITHOUT ANY WARRANTY; without even the implied warranty of     *
0015   * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU   *
0016   * General Public License for more details.                           *
0017   *                                                                    *
0018   * You should have received a copy of the GNU General Public License  *
0019   * along with this library (see file COPYING); if not, write          *
0020   * to the Free Software Foundation, Inc., 59 Temple Place, Suite      *
0021   * 330, Boston, MA 02111-1307 USA, or contact the author.             *
0022   *                                                                    *
0023   **********************************************************************/
0024
0025 // Header file for class GSLMultiRootFinder
0026 //
0027
0028 #ifndef ROOT_Math_GSLMultiRootFinder
0029 #define ROOT_Math_GSLMultiRootFinder
0030
0031
0032
0033 #include "Math/IFunction.h"
0034
0035 #include "Math/WrappedFunction.h"
0036
0037 #include <vector>
0038 #include <utility>
0039 #include <iostream>
0040
0041 namespace ROOT {
0042 namespace Math {
0043
0044
0045    class GSLMultiRootBaseSolver;
0046
0047    /** @defgroup MultiRoot  Multidimensional ROOT finding
0048        Classes for finding the roots of a multi-dimensional system.
0049        @ingroup NumAlgo
0050    */
0051
0052   /**
0053      Class for  Multidimensional root finding algorithms bassed on GSL. This class is used to solve a
0054      non-linear system of equations:
0055
0056      f1(x1,....xn) = 0
0057      f2(x1,....xn) = 0
0058      ..................
0059      fn(x1,....xn) = 0
0060
0061      See the GSL <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html"> online manual</A> for
0062      information on the GSL MultiRoot finding algorithms
0063
0064      The available GSL algorithms require the derivatives of the supplied functions or not (they are
0065      computed internally by GSL). In the first case the user needs to provide a list of multidimensional functions implementing the
0066      gradient interface (ROOT::Math::IMultiGradFunction) while in the second case it is enough to supply a list of
0067      functions impelmenting the ROOT::Math::IMultiGenFunction interface.
0068      The available algorithms requiring derivatives (see also the GSL
0069      <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Algorithms-using-Derivatives.html">documentation</A> )
0070      are the followings:
0071      <ul>
0072          <li><tt>ROOT::Math::GSLMultiRootFinder::kHybridSJ</tt>  with name <i>"HybridSJ"</i>: modified Powell's hybrid
0073      method as implemented in HYBRJ in MINPACK
0074          <li><tt>ROOT::Math::GSLMultiRootFinder::kHybridJ</tt>  with name <i>"HybridJ"</i>: unscaled version of the
0075      previous algorithm</li>
0076          <li><tt>ROOT::Math::GSLMultiRootFinder::kNewton</tt>  with name <i>"Newton"</i>: Newton method </li>
0077          <li><tt>ROOT::Math::GSLMultiRootFinder::kGNewton</tt>  with name <i>"GNewton"</i>: modified Newton method </li>
0078      </ul>
0079      The algorithms without derivatives (see also the GSL
0080      <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Algorithms-without-Derivatives.html">documentation</A> )
0081      are the followings:
0082      <ul>
0083          <li><tt>ROOT::Math::GSLMultiRootFinder::kHybridS</tt>  with name <i>"HybridS"</i>: same as HybridSJ but using
0084      finate difference approximation for the derivatives</li>
0085          <li><tt>ROOT::Math::GSLMultiRootFinder::kHybrid</tt>  with name <i>"Hybrid"</i>: unscaled version of the
0086      previous algorithm</li>
0087          <li><tt>ROOT::Math::GSLMultiRootFinder::kDNewton</tt>  with name <i>"DNewton"</i>: discrete Newton algorithm </li>
0088          <li><tt>ROOT::Math::GSLMultiRootFinder::kBroyden</tt>  with name <i>"Broyden"</i>: Broyden algorithm </li>
0089      </ul>
0090
0091      @ingroup MultiRoot
0092   */
0093
0094
0095  class GSLMultiRootFinder {
0096
0097  public:
0098
0099    /**
0100       enumeration specifying the types of GSL multi root finders
0101       requiring the derivatives
0102
0103    */
0104     enum EDerivType {
0105        kHybridSJ,
0106        kHybridJ,
0107        kNewton,
0108        kGNewton
0109     };
0110     /**
0111        enumeration specifying the types of GSL multi root finders
0112        which do not require the derivatives
0113
0114     */
0115     enum EType {
0116        kHybridS,
0117        kHybrid,
0118        kDNewton,
0119        kBroyden
0120     };
0121
0122
0123
0124     /// create a multi-root finder based on an algorithm not requiring function derivative
0125     GSLMultiRootFinder(EType type);
0126
0127     /// create a multi-root finder based on an algorithm requiring function derivative
0128     GSLMultiRootFinder(EDerivType type);
0129
0130     /*
0131       create a multi-root finder using a string.
0132       The names are those defined in the GSL manuals
0133       after having remived the GSL prefix (gsl_multiroot_fsolver).
0134       Default algorithm  is "hybrids" (without derivative).
0135     */
0136     GSLMultiRootFinder(const char * name = nullptr);
0137
0138     /// destructor
0139     virtual ~GSLMultiRootFinder();
0140
0141  private:
0142     // usually copying is non trivial, so we make this unaccessible
0143     GSLMultiRootFinder(const GSLMultiRootFinder &);
0144     GSLMultiRootFinder & operator = (const GSLMultiRootFinder &);
0145
0146  public:
0147
0148     /// set the type for an algorithm without derivatives
0149     void SetType(EType type) {
0150        fType = type; fUseDerivAlgo = false;
0151     }
0152
0153     /// set the type of algorithm using derivatives
0154     void SetType(EDerivType type) {
0155        fType = type; fUseDerivAlgo = true;
0156     }
0157
0158     /// set the type using a string
0159     void SetType(const char * name);
0160
0161     /*
0162        add the list of functions f1(x1,..xn),...fn(x1,...xn). The list must contain pointers of
0163        ROOT::Math::IMultiGenFunctions. The method requires the
0164        the begin and end of the list iterator.
0165        The list can be any stl container or a simple array of  ROOT::Math::IMultiGenFunctions* or
0166        whatever implementing an iterator.
0167        If using a derivative type algorithm the function pointers must implement the
0168        ROOT::Math::IMultiGradFunction interface
0169     */
0170     template<class FuncIterator>
0171     bool SetFunctionList( FuncIterator begin, FuncIterator end) {
0172        bool ret = true;
0173        for (FuncIterator itr = begin; itr != end; ++itr) {
0174           const ROOT::Math::IMultiGenFunction * f = *itr;
0175           // Using bitwise operator &= require the operand to be a bool
0176           // to have the intended effect here.
0177           ret &= (AddFunction( *f) != 0);
0178        }
0179        return ret;
0180     }
0181
0182     /*
0183       add (set) a single function fi(x1,...xn) which is part of the system of
0184        specifying the begin and end of the iterator.
0185        If using a derivative type algorithm the function must implement the
0186        ROOT::Math::IMultiGradFunction interface
0187        Return the current number of function in the list and 0 if failed to add the function
0188      */
0189     int AddFunction( const ROOT::Math::IMultiGenFunction & func);
0190
0191     /// same method as before but using any function implementing
0192     /// the operator(), so can be wrapped in a IMultiGenFunction interface
0193     template <class Function>
0194     int AddFunction( Function & f, int ndim) {
0195        // no need to care about lifetime of wfunc. It will be cloned inside AddFunction
0196        WrappedMultiFunction<Function &> wfunc(f, ndim);
0197        return AddFunction(wfunc);
0198     }
0199
0200     /**
0201        return the number of sunctions set in the class.
0202        The number must be equal to the dimension of the functions
0203      */
0204     unsigned  int Dim() const { return fFunctions.size(); }
0205
0206     /// clear list of functions
0207     void Clear();
0208
0209     /// return the root X values solving the system
0210     const double * X() const;
0211
0212     /// return the function values f(X) solving the system
0213     /// i.e. they must be close to zero at the solution
0214     const double * FVal() const;
0215
0216     /// return the last step size
0217     const double * Dx() const;
0218
0219
0220     /**
0221        Find the root starting from the point X;
0222        Use the number of iteration and tolerance if given otherwise use
0223        default parameter values which can be defined by
0224        the static method SetDefault...
0225     */
0226     bool Solve(const double * x,  int maxIter = 0, double absTol = 0, double relTol = 0);
0227
0228     /// Return number of iterations
0229     int Iterations() const {
0230        return fIter;
0231     }
0232
0233     /// Return the status of last root finding
0234     int Status() const { return fStatus; }
0235
0236     /// Return the algorithm name used for solving
0237     /// Note the name is available only after having called solved
0238     /// Otherwise an empyty string is returned
0239     const char * Name() const;
0240
0241     /*
0242        set print level
0243        level = 0  quiet (no messages print)
0244              = 1  print only the result
0245              = 3  max debug. Print result at each iteration
0246     */
0247     void SetPrintLevel(int level) { fPrintLevel = level; }
0248
0249     /// return the print level
0250     int PrintLevel() const { return fPrintLevel; }
0251
0252
0253     //-- static methods to set configurations
0254
0255     /// set tolerance (absolute and relative)
0256     /// relative tolerance is only use to verify the convergence
0257     /// do it is a minor parameter
0258     static void SetDefaultTolerance(double abstol, double reltol = 0 );
0259
0260     /// set maximum number of iterations
0261     static void SetDefaultMaxIterations(int maxiter);
0262
0263     /// print iteration state
0264     void PrintState(std::ostream & os = std::cout);
0265
0266
0267  protected:
0268
0269     // return type given a name
0270     std::pair<bool,int> GetType(const char * name);
0271     // clear list of functions
0272     void ClearFunctions();
0273
0274
0275  private:
0276
0277     int fIter;           // current number of iterations
0278     int fStatus;         // current status
0279     int fPrintLevel;     // print level
0280
0281     // int fMaxIter;        // max number of iterations
0282     // double fAbsTolerance;  // absolute tolerance
0283     // double fRelTolerance;  // relative tolerance
0284     int fType;            // type of algorithm
0285     bool fUseDerivAlgo; // algorithm using derivative
0286
0287     GSLMultiRootBaseSolver * fSolver;
0288     std::vector<ROOT::Math::IMultiGenFunction *> fFunctions;   //! transient Vector of the functions
0289
0290
0291  };
0292
0293    // use typedef for most sensible name
0294    typedef GSLMultiRootFinder MultiRootFinder;
0295
0296 } // namespace Math
0297 } // namespace ROOT
0298
0299
0300 #endif /* ROOT_Math_GSLMultiRootFinder */