Justine Devin / gammalib

Commit 12e667bcc3447b810a2afc8ba288d858d5d82b12

Authored by knodlseder
1 parent f3cdb72c

initial version

Inline Side-by-side
Showing 4 changed files with 624 additions and 0 deletions
include/GFunction.hpp 0 → 100644
  View file @12e667b
  1 +/***************************************************************************
  2 + * GFunction.hpp - Single parameter function abstract base class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2011 by Jurgen Knodlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software; you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation; either version 2 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + ***************************************************************************/
  13 +/**
  14 + * @file GFunction.hpp
  15 + * @brief GFunction abstract virtual base class interface definition.
  16 + * @author J. Knodlseder
  17 + */
  18 +
  19 +#ifndef GFUNCTION_HPP
  20 +#define GFUNCTION_HPP
  21 +
  22 +/* __ Includes ___________________________________________________________ */
  23 +
  24 +
  25 +/***********************************************************************//**
  26 + * @class GFunction
  27 + *
  28 + * @brief GFunction class interface defintion.
  29 + *
  30 + * This class implements the abstract interface for a one parameter function.
  31 + * This function is for example used for integration or numerical computation
  32 + * of derivatives. This class has no members. The only pure virtual method
  33 + * that needs to be implemented by the derived class is the eval() method
  34 + * that provides function evaluation at a given value x, e.g. y=eval(x).
  35 + ***************************************************************************/
  36 +class GFunction {
  37 +
  38 +public:
  39 +
  40 + // Constructors and destructors
  41 + GFunction(void);
  42 + GFunction(const GFunction& func);
  43 + virtual ~GFunction(void);
  44 +
  45 + // Operators
  46 + GFunction& operator= (const GFunction& func);
  47 +
  48 + // Methods
  49 + virtual double eval(double x) = 0;
  50 +
  51 +protected:
  52 + // Protected methods
  53 + void init_members(void);
  54 + void copy_members(const GFunction& func);
  55 + void free_members(void);
  56 +};
  57 +
  58 +#endif /* GFUNCTION_HPP */
... ...
src/numerics/GDerivative.cpp 0 → 100644
  View file @12e667b
  1 +/***************************************************************************
  2 + * GDerivative.cpp - Derivative class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2011 by Jurgen Knodlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software; you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation; either version 2 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + ***************************************************************************/
  13 +/**
  14 + * @file GDerivative.cpp
  15 + * @brief GDerivative class implementation.
  16 + * @author J. Knodlseder
  17 + */
  18 +
  19 +/* __ Includes ___________________________________________________________ */
  20 +#include <cstdlib> // Definition of NULL
  21 +#include <float.h> // for DBL_MAX
  22 +#include "GDerivative.hpp"
  23 +#include "GMatrix.hpp"
  24 +
  25 +/* __ Method name definitions ____________________________________________ */
  26 +#define G_RIDDER "GDerivative::ridder(double&, double&, double&)"
  27 +
  28 +/* __ Macros _____________________________________________________________ */
  29 +
  30 +/* __ Coding definitions _________________________________________________ */
  31 +#define G_VALUE_SIGNAL_BAD_TOLERANCE 1 //!< Signal bad tolerance
  32 +
  33 +/* __ Debug definitions __________________________________________________ */
  34 +#define G_VALUE_DEBUG 0 //!< Debug value method
  35 +
  36 +
  37 +/*==========================================================================
  38 + = =
  39 + = Constructors/destructors =
  40 + = =
  41 + ==========================================================================*/
  42 +
  43 +/***********************************************************************//**
  44 + * @brief Void constructor
  45 + ***************************************************************************/
  46 +GDerivative::GDerivative(void)
  47 +{
  48 + // Initialise private members
  49 + init_members();
  50 +
  51 + // Return
  52 + return;
  53 +}
  54 +
  55 +
  56 +/***********************************************************************//**
  57 + * @brief Function constructor
  58 + *
  59 + * @param[in] func Function pointer.
  60 + *
  61 + * Construct class instance from single parameter function. Note that the
  62 + * function is not copied but just a pointer to it is stored. The function
  63 + * should thus not be destroyed before computation of derivatives is
  64 + * finished.
  65 + ***************************************************************************/
  66 +GDerivative::GDerivative(GFunction* func)
  67 +{
  68 + // Initialise private members
  69 + init_members();
  70 +
  71 + // Set function
  72 + this->function(func);
  73 +
  74 + // Return
  75 + return;
  76 +}
  77 +
  78 +
  79 +/***********************************************************************//**
  80 + * @brief Copy constructor
  81 + *
  82 + * @param[in] dx Derivative.
  83 + ***************************************************************************/
  84 +GDerivative::GDerivative(const GDerivative& dx)
  85 +{
  86 + // Initialise members
  87 + init_members();
  88 +
  89 + // Copy members
  90 + copy_members(dx);
  91 +
  92 + // Return
  93 + return;
  94 +}
  95 +
  96 +
  97 +/***********************************************************************//**
  98 + * @brief Destructor
  99 + ***************************************************************************/
  100 +GDerivative::~GDerivative(void)
  101 +{
  102 + // Free members
  103 + free_members();
  104 +
  105 + // Return
  106 + return;
  107 +}
  108 +
  109 +
  110 +/*==========================================================================
  111 + = =
  112 + = Operators =
  113 + = =
  114 + ==========================================================================*/
  115 +
  116 +/***********************************************************************//**
  117 + * @brief Assignment operator
  118 + *
  119 + * @param[in] dx Derivative.
  120 + ***************************************************************************/
  121 +GDerivative& GDerivative::operator= (const GDerivative& dx)
  122 +{
  123 + // Execute only if object is not identical
  124 + if (this != &dx) {
  125 +
  126 + // Free members
  127 + free_members();
  128 +
  129 + // Initialise members
  130 + init_members();
  131 +
  132 + // Copy members
  133 + copy_members(dx);
  134 +
  135 + } // endif: object was not identical
  136 +
  137 + // Return
  138 + return *this;
  139 +}
  140 +
  141 +
  142 +/*==========================================================================
  143 + = =
  144 + = Public methods =
  145 + = =
  146 + ==========================================================================*/
  147 +
  148 +/***********************************************************************//**
  149 + * @brief Returns derivative
  150 + *
  151 + * @param[in] x Function value.
  152 + *
  153 + * Computes derivative using the Ridders' method. The method starts with a
  154 + * step size h=0.02*|x| (with h >= 1e-6) and iteratively decreases the
  155 + * step size until the error becomes less than 1e-6. A maximum number of
  156 + * 5 iterations is done.
  157 + ***************************************************************************/
  158 +double GDerivative::value(const double& x)
  159 +{
  160 + // Set constants
  161 + const int max_iter = 5;
  162 + const double min_h = 1.0e-6;
  163 + const double max_err = 1.0e-6;
  164 +
  165 + // Initialise result
  166 + double result = 0.0;
  167 +
  168 + // Set initial step size
  169 + double h = 0.02*fabs(x);
  170 + if (h < min_h)
  171 + h = min_h;
  172 +
  173 + // Initialise tolerance
  174 + double err = 1.0;
  175 +
  176 + // Loop over Ridder's method until we have an acceptable tolerance
  177 + int iter = 0;
  178 + for (; iter < max_iter; ++iter) {
  179 +
  180 + // Compute derivative using Ridder's method
  181 + result = ridder(x, h, err);
  182 +
  183 + // Debug option: Show actual results
  184 + #if G_VALUE_DEBUG
  185 + std::cout << "GDerivative::value(";
  186 + std::cout << "iter=" << iter;
  187 + std::cout << ", x=" << x;
  188 + std::cout << ", dy/dx=" << result;
  189 + std::cout << ", h=" << h;
  190 + std::cout << ", err=" << err;
  191 + std::cout << ")" << std::endl;
  192 + #endif
  193 +
  194 + // If uncertainty is below tolerance then exit now
  195 + if (err < max_err)
  196 + break;
  197 +
  198 + // ... otherwise reduce the step size
  199 + h /= 5.0;
  200 +
  201 + } // endfor: tolerance matching loop
  202 +
  203 + // Compile option: signal if we exceed the tolerance
  204 + #if G_VALUE_SIGNAL_BAD_TOLERANCE
  205 + if (err >= max_err) {
  206 + std::cout << "WARNING: GDerivative::value(";
  207 + std::cout << "iter=" << iter;
  208 + std::cout << ", x=" << x;
  209 + std::cout << ", dy/dx=" << result;
  210 + std::cout << ", h=" << h;
  211 + std::cout << ", err=" << err;
  212 + std::cout << "): error exceeds tolerance of ";
  213 + std::cout << max_err << std::endl;
  214 + }
  215 + #endif
  216 +
  217 +
  218 + // Return derivative
  219 + return result;
  220 +}
  221 +
  222 +
  223 +/***********************************************************************//**
  224 + * @brief Returns derivative by Ridders' method
  225 + *
  226 + * @param[in] x Function value.
  227 + * @param[in] h Estimated initial step size.
  228 + * @param[out] err Estimated error in the derivative.
  229 + *
  230 + * @exception GException::invalid_argument
  231 + * Step size of 0 has been specified.
  232 + *
  233 + * Returns the derivative of a function at point x by Ridders's method of
  234 + * polynomial extrapolation. The value h is input as an estimated initial
  235 + * stepsize; it needs not to be small, but rather should be an increment
  236 + * in x over which the function changes substantially.
  237 + *
  238 + * This method has been adopted from Numerical Recipes, 3rd edition, page
  239 + * 321f.
  240 + ***************************************************************************/
  241 +double GDerivative::ridder(const double& x, const double& h, double& err)
  242 +{
  243 + // Set constants
  244 + const int ntab = 10; // Maximum table size
  245 + const double con = 1.4; // Stepsize decreases by con at each iteration
  246 + const double con2 = con * con;
  247 + const double big = DBL_MAX;
  248 + const double safe = 2.0; // Return when error is safe worse than best so far
  249 +
  250 + // Initialise result
  251 + double dx = 0.0;
  252 + err = big;
  253 +
  254 + // Continue only if we have a valid function
  255 + if (m_func != NULL) {
  256 +
  257 + // Allocate Neville table
  258 + GMatrix a(ntab, ntab);
  259 +
  260 + // Throw exception if step size is zero
  261 + if (h <= 0)
  262 + throw GException::invalid_argument(G_RIDDER, "Step size must be positive.");
  263 +
  264 + // Initialise adaptive step size
  265 + double hh = h;
  266 +
  267 + // Perform initial function evaluation
  268 + a(0,0) = (m_func->eval(x+hh) - m_func->eval(x-hh))/(2.0*hh);
  269 +
  270 + // Perform successive computations
  271 + for (int i = 1; i < ntab; ++i) {
  272 +
  273 + // Reduce step size
  274 + hh /= con;
  275 +
  276 + // Try new smaller step size
  277 + a(0,i) = (m_func->eval(x+hh) - m_func->eval(x-hh))/(2.0*hh);
  278 +
  279 + // Compute extrapolations of various orders, requiring no new
  280 + // function evaluations
  281 + double fac = con2;
  282 + for (int j = 1; j <= i; ++j) {
  283 +
  284 + // Compute extrapolation
  285 + a(j,i) = ( a(j-1,i) * fac - a(j-1,i-1) ) / (fac - 1.0);
  286 + fac *= con2;
  287 +
  288 + // If the error is reduced then save the improved answer
  289 + double err1 = fabs(a(j,i) - a(j-1,i));
  290 + double err2 = fabs(a(j,i) - a(j-1,i-1));
  291 + double errt = (err1 > err2) ? err1 : err2;
  292 + if (errt <= err) {
  293 + err = errt;
  294 + dx = a(j,i);
  295 + }
  296 +
  297 + } // endfor: looped over tableau
  298 +
  299 + // If higher order is worse by a significant factor, the quit early
  300 + if (fabs( a(i,i) - a(i-1,i-1) ) >= safe*err)
  301 + break;
  302 +
  303 + } // endfor: looped over computations
  304 +
  305 + } // endif: we had a valid function
  306 +
  307 + // Return derivative
  308 + return dx;
  309 +}
  310 +
  311 +
  312 +/*==========================================================================
  313 + = =
  314 + = Private methods =
  315 + = =
  316 + ==========================================================================*/
  317 +
  318 +/***********************************************************************//**
  319 + * @brief Initialise class members
  320 + ***************************************************************************/
  321 +void GDerivative::init_members(void)
  322 +{
  323 + // Initialise members
  324 + m_func = NULL;
  325 +
  326 + // Return
  327 + return;
  328 +}
  329 +
  330 +
  331 +/***********************************************************************//**
  332 + * @brief Copy class members
  333 + *
  334 + * @param[in] dx Derivative.
  335 + ***************************************************************************/
  336 +void GDerivative::copy_members(const GDerivative& dx)
  337 +{
  338 + // Copy attributes
  339 + m_func = dx.m_func;
  340 +
  341 + // Return
  342 + return;
  343 +}
  344 +
  345 +
  346 +/***********************************************************************//**
  347 + * @brief Delete class members
  348 + ***************************************************************************/
  349 +void GDerivative::free_members(void)
  350 +{
  351 + // Return
  352 + return;
  353 +}
  354 +
  355 +
  356 +/*==========================================================================
  357 + = =
  358 + = Friends =
  359 + = =
  360 + ==========================================================================*/
... ...
src/numerics/GException_numerics.cpp 0 → 100644
  View file @12e667b
  1 +/***************************************************************************
  2 + * GException_numerics.cpp - Numerics exceptions *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2011 by Jurgen Knodlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software; you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation; either version 2 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + ***************************************************************************/
  13 +/**
  14 + * @file GException_numerics.cpp
  15 + * @brief Numerics exception class implementations
  16 + * @author J. Knodlseder
  17 + */
  18 +
  19 +/* __ Includes ___________________________________________________________ */
  20 +#ifdef HAVE_CONFIG_H
  21 +#include <config.h>
  22 +#endif
  23 +#include "GException.hpp"
  24 +#include "GTools.hpp"
  25 +
  26 +
  27 +/***********************************************************************//**
  28 + * @brief <Description>
  29 + *
  30 + * @param[in] origin Name of method that has thrown the exception.
  31 + * @param[in] message Optional error message.
  32 + ***************************************************************************/
  33 +/*
  34 +GException::exception::exception(std::string origin, std::string message)
  35 +{
  36 + // Set origin of exception
  37 + m_origin = origin;
  38 +
  39 + // Set exception message
  40 + m_message = "Exception. "+ message;
  41 +
  42 + // Return
  43 + return;
  44 +}
  45 +*/
0 46 \ No newline at end of file
... ...
src/numerics/GFunction.cpp 0 → 100644
  View file @12e667b
  1 +/***************************************************************************
  2 + * GFunction.cpp - Single parameter function abstract base class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2011 by Jurgen Knodlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software; you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation; either version 2 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + ***************************************************************************/
  13 +/**
  14 + * @file GFunction.cpp
  15 + * @brief GFunction abstract virtual base class implementation.
  16 + * @author J. Knodlseder
  17 + */
  18 +
  19 +/* __ Includes ___________________________________________________________ */
  20 +#include "GFunction.hpp"
  21 +
  22 +/* __ Method name definitions ____________________________________________ */
  23 +
  24 +/* __ Macros _____________________________________________________________ */
  25 +
  26 +/* __ Coding definitions _________________________________________________ */
  27 +
  28 +/* __ Debug definitions __________________________________________________ */
  29 +
  30 +
  31 +/*==========================================================================
  32 + = =
  33 + = Constructors/destructors =
  34 + = =
  35 + ==========================================================================*/
  36 +
  37 +/***********************************************************************//**
  38 + * @brief Void constructor
  39 + ***************************************************************************/
  40 +GFunction::GFunction(void)
  41 +{
  42 + // Initialise members
  43 + init_members();
  44 +
  45 + // Return
  46 + return;
  47 +}
  48 +
  49 +
  50 +/***********************************************************************//**
  51 + * @brief Copy constructor
  52 + *
  53 + * @param[in] func Function.
  54 + ***************************************************************************/
  55 +GFunction::GFunction(const GFunction& func)
  56 +{
  57 + // Initialise members
  58 + init_members();
  59 +
  60 + // Copy members
  61 + copy_members(func);
  62 +
  63 + // Return
  64 + return;
  65 +}
  66 +
  67 +
  68 +/***********************************************************************//**
  69 + * @brief Destructor
  70 + ***************************************************************************/
  71 +GFunction::~GFunction(void)
  72 +{
  73 + // Free members
  74 + free_members();
  75 +
  76 + // Return
  77 + return;
  78 +}
  79 +
  80 +
  81 +/*==========================================================================
  82 + = =
  83 + = Operators =
  84 + = =
  85 + ==========================================================================*/
  86 +
  87 +/***********************************************************************//**
  88 + * @brief Assignment operator
  89 + *
  90 + * @param[in] func Function.
  91 + ***************************************************************************/
  92 +GFunction& GFunction::operator= (const GFunction& func)
  93 +{
  94 + // Execute only if object is not identical
  95 + if (this != &func) {
  96 +
  97 + // Free members
  98 + free_members();
  99 +
  100 + // Initialise members
  101 + init_members();
  102 +
  103 + // Copy members
  104 + copy_members(func);
  105 +
  106 + } // endif: object was not identical
  107 +
  108 + // Return
  109 + return *this;
  110 +}
  111 +
  112 +
  113 +/*==========================================================================
  114 + = =
  115 + = Public methods =
  116 + = =
  117 + ==========================================================================*/
  118 +
  119 +/*==========================================================================
  120 + = =
  121 + = Private methods =
  122 + = =
  123 + ==========================================================================*/
  124 +
  125 +/***********************************************************************//**
  126 + * @brief Initialise class members
  127 + ***************************************************************************/
  128 +void GFunction::init_members(void)
  129 +{
  130 + // Return
  131 + return;
  132 +}
  133 +
  134 +
  135 +/***********************************************************************//**
  136 + * @brief Copy class members
  137 + *
  138 + * @param[in] func Function.
  139 + ***************************************************************************/
  140 +void GFunction::copy_members(const GFunction& func)
  141 +{
  142 + // Return
  143 + return;
  144 +}
  145 +
  146 +
  147 +/***********************************************************************//**
  148 + * @brief Delete class members
  149 + ***************************************************************************/
  150 +void GFunction::free_members(void)
  151 +{
  152 + // Return
  153 + return;
  154 +}
  155 +
  156 +
  157 +/*==========================================================================
  158 + = =
  159 + = Friends =
  160 + = =
  161 + ==========================================================================*/
... ...