[1230] | 1 | // |
---|
| 2 | // ******************************************************************** |
---|
| 3 | // * License and Disclaimer * |
---|
| 4 | // * * |
---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
| 7 | // * conditions of the Geant4 Software License, included in the file * |
---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
| 9 | // * include a list of copyright holders. * |
---|
| 10 | // * * |
---|
| 11 | // * Neither the authors of this software system, nor their employing * |
---|
| 12 | // * institutes,nor the agencies providing financial support for this * |
---|
| 13 | // * work make any representation or warranty, express or implied, * |
---|
| 14 | // * regarding this software system or assume any liability for its * |
---|
| 15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
| 16 | // * for the full disclaimer and the limitation of liability. * |
---|
| 17 | // * * |
---|
| 18 | // * This code implementation is the result of the scientific and * |
---|
| 19 | // * technical work of the GEANT4 collaboration. * |
---|
| 20 | // * By using, copying, modifying or distributing the software (or * |
---|
| 21 | // * any work based on the software) you agree to acknowledge its * |
---|
| 22 | // * use in resulting scientific publications, and indicate your * |
---|
| 23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
| 24 | // ******************************************************************** |
---|
| 25 | // |
---|
| 26 | // |
---|
[807] | 27 | /** |
---|
| 28 | * \file |
---|
| 29 | * \brief Provides the headers for the general c2_function algebra which supports |
---|
| 30 | * fast, flexible operations on piecewise-twice-differentiable functions |
---|
| 31 | * |
---|
| 32 | * \author Created by R. A. Weller and Marcus H. Mendenhall on 7/9/05. |
---|
| 33 | * \author Copyright 2005 __Vanderbilt University__. All rights reserved. |
---|
| 34 | * |
---|
[1230] | 35 | * \version c2_function.hh,v 1.238 2008/05/22 12:45:19 marcus Exp |
---|
| 36 | * \see \ref c2_factory "Factory Functions" for information on constructing things in here |
---|
[807] | 37 | */ |
---|
| 38 | |
---|
[1230] | 39 | #ifndef __has_c2_function_hh |
---|
| 40 | #define __has_c2_function_hh 1 |
---|
[807] | 41 | |
---|
[1230] | 42 | // MSVC does not automatically define numerical constants such as M_PI without this. |
---|
| 43 | // this came from the msdn website, so it should be right... |
---|
| 44 | #ifdef _MSC_VER |
---|
| 45 | #define _USE_MATH_DEFINES |
---|
| 46 | #define c2_isnan _isnan |
---|
| 47 | #define c2_isfinite _finite |
---|
| 48 | #else |
---|
| 49 | #define c2_isnan std::isnan |
---|
| 50 | #define c2_isfinite std::isfinite |
---|
| 51 | #endif |
---|
| 52 | |
---|
[807] | 53 | #include <cmath> |
---|
| 54 | #include <vector> |
---|
[1230] | 55 | #include <utility> |
---|
[807] | 56 | #include <string> |
---|
[1230] | 57 | #include <stdexcept> |
---|
| 58 | #include <typeinfo> |
---|
| 59 | #include <sstream> |
---|
[807] | 60 | |
---|
| 61 | /// \brief the exception class for c2_function operations. |
---|
| 62 | class c2_exception : public std::exception { |
---|
| 63 | public: |
---|
| 64 | /// \brief construct the exception with an error message |
---|
| 65 | /// \param msgcode the message |
---|
| 66 | c2_exception(const char msgcode[]) : info(msgcode) { } |
---|
| 67 | virtual ~c2_exception() throw() { } |
---|
| 68 | /** Returns a C-style character string describing the general cause |
---|
| 69 | * of the current error. */ |
---|
| 70 | virtual const char* what() const throw() { return info.c_str(); } |
---|
| 71 | private: |
---|
| 72 | std::string info; |
---|
| 73 | }; |
---|
| 74 | |
---|
| 75 | // put these forward references here, and with a bogus typename to make swig happy. |
---|
[1230] | 76 | template <typename float_type> class c2_composed_function_p; |
---|
| 77 | template <typename float_type> class c2_sum_p; |
---|
| 78 | template <typename float_type> class c2_diff_p; |
---|
| 79 | template <typename float_type> class c2_product_p; |
---|
| 80 | template <typename float_type> class c2_ratio_p; |
---|
| 81 | template <typename float_type> class c2_piecewise_function_p; |
---|
| 82 | template <typename float_type> class c2_quadratic_p; |
---|
| 83 | template <typename float_type> class c2_ptr; |
---|
| 84 | /** |
---|
| 85 | \defgroup abstract_classes Abstract Classes |
---|
| 86 | \defgroup arithmetic_functions Arithmetic Functions |
---|
| 87 | \defgroup math_functions Mathemetical Functions |
---|
| 88 | \defgroup parametric_functions Parametric Families of Functions |
---|
| 89 | \defgroup interpolators Interpolating Functions |
---|
| 90 | \defgroup containers Functions which are containers for, or functions of, other functions |
---|
| 91 | \defgroup factories Factory classes which reduce silly template typing |
---|
| 92 | \defgroup transforms Classes which provide coordinate system transformations, wih derivatives |
---|
| 93 | */ |
---|
[807] | 94 | |
---|
[1230] | 95 | /// \brief structure used to hold evaluated function data at a point. |
---|
| 96 | /// |
---|
| 97 | /// Contains all the information for the function at one point. |
---|
| 98 | template <typename float_type> class c2_fblock |
---|
| 99 | { |
---|
| 100 | public: |
---|
| 101 | /// \brief the abscissa |
---|
| 102 | float_type x; |
---|
| 103 | /// \brief the value of the function at \a x |
---|
| 104 | float_type y; |
---|
| 105 | /// \brief the derivative at \a x |
---|
| 106 | float_type yp; |
---|
| 107 | /// \brief the second derivative at \a x |
---|
| 108 | float_type ypp; |
---|
| 109 | /// flag, filled in by c2_function::fill_fblock(), indicating the derivative is NaN of Inf |
---|
| 110 | bool ypbad; |
---|
| 111 | /// flag, filled in by c2_function::fill_fblock(), indicating the second derivative is NaN of Inf |
---|
| 112 | bool yppbad; |
---|
| 113 | }; |
---|
| 114 | |
---|
[807] | 115 | /** |
---|
| 116 | \brief the parent class for all c2_functions. |
---|
[1230] | 117 | \ingroup abstract_classes |
---|
[807] | 118 | c2_functions know their value, first, and second derivative at almost every point. |
---|
| 119 | They can be efficiently combined with binary operators, via c2_binary_function, |
---|
[1230] | 120 | composed via c2_composed_function_, |
---|
[807] | 121 | have their roots found via find_root(), |
---|
| 122 | and be adaptively integrated via partial_integrals() or integral(). |
---|
| 123 | They also can carry information with them about how to find 'interesting' points on the function. |
---|
| 124 | This information is set with set_sampling_grid() and extracted with get_sampling_grid(). |
---|
| 125 | |
---|
| 126 | Particularly important subclasses are the interpolating functions classes, |
---|
| 127 | interpolating_function , lin_log_interpolating_function, log_lin_interpolating_function, |
---|
| 128 | log_log_interpolating_function, and arrhenius_interpolating_function, |
---|
| 129 | as well as the template functions |
---|
[1230] | 130 | inverse_integrated_density_function(). |
---|
[807] | 131 | |
---|
[1230] | 132 | For a discussion of memory management, see \ref memory_management |
---|
[807] | 133 | |
---|
| 134 | */ |
---|
| 135 | template <typename float_type=double> class c2_function { |
---|
| 136 | public: |
---|
| 137 | /// \brief get versioning information for the header file |
---|
| 138 | /// \return the CVS Id string |
---|
| 139 | const std::string cvs_header_vers() const { return |
---|
[1230] | 140 | "c2_function.hh,v 1.238 2008/05/22 12:45:19 marcus Exp"; |
---|
[807] | 141 | } |
---|
| 142 | |
---|
| 143 | /// \brief get versioning information for the source file |
---|
| 144 | /// \return the CVS Id string |
---|
| 145 | const std::string cvs_file_vers() const ; |
---|
| 146 | |
---|
| 147 | public: |
---|
| 148 | /// \brief destructor |
---|
[1230] | 149 | virtual ~c2_function() { |
---|
| 150 | if(sampling_grid && !no_overwrite_grid) delete sampling_grid; |
---|
| 151 | if(root_info) delete root_info; |
---|
| 152 | if(owner_count) { |
---|
| 153 | std::ostringstream outstr; |
---|
| 154 | outstr << "attempt to delete an object with non-zero ownership in class "; |
---|
| 155 | outstr << typeid(*this).name() << std::endl; |
---|
| 156 | throw c2_exception(outstr.str().c_str()); |
---|
| 157 | } |
---|
| 158 | } |
---|
[807] | 159 | |
---|
| 160 | /// \brief get the value and derivatives. |
---|
| 161 | /// |
---|
| 162 | /// There is required checking for null pointers on the derivatives, |
---|
| 163 | /// and most implementations should operate faster if derivatives are not needed. |
---|
| 164 | /// \param[in] x the point at which to evaluate the function |
---|
| 165 | /// \param[out] yprime the first derivative (if pointer is non-null) |
---|
| 166 | /// \param[out] yprime2 the second derivative (if pointer is non-null) |
---|
| 167 | /// \return the value of the function |
---|
| 168 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) =0 ; // { return 0; }; |
---|
| 169 | |
---|
| 170 | /// \brief evaluate the function in the classic way, ignoring derivatives. |
---|
| 171 | /// \param x the point at which to evaluate |
---|
| 172 | /// \return the value of the function |
---|
| 173 | inline float_type operator () (float_type x) const throw(c2_exception) |
---|
| 174 | { return value_with_derivatives(x, (float_type *)0, (float_type *)0); } |
---|
| 175 | |
---|
| 176 | /// \brief get the value and derivatives. |
---|
| 177 | /// |
---|
| 178 | /// \param[in] x the point at which to evaluate the function |
---|
| 179 | /// \param[out] yprime the first derivative (if pointer is non-null) |
---|
| 180 | /// \param[out] yprime2 the second derivative (if pointer is non-null) |
---|
| 181 | /// \return the value of the function |
---|
| 182 | inline float_type operator () (float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
| 183 | { return value_with_derivatives(x, yprime, yprime2); } |
---|
| 184 | |
---|
| 185 | /// \brief solve f(x)==value very efficiently, with explicit knowledge of derivatives of the function |
---|
| 186 | /// |
---|
| 187 | /// find_root solves by iterated inverse quadratic extrapolation for a solution to f(x)=y. It |
---|
| 188 | /// includes checks against bad convergence, so it should never be able to fail. Unlike typical |
---|
| 189 | /// secant method or fancier Brent's method finders, this does not depend in any strong wasy on the |
---|
| 190 | /// brackets, unless the finder has to resort to successive approximations to close in on a root. |
---|
| 191 | /// Often, it is possible to make the brackets equal to the domain of the function, if there is |
---|
| 192 | /// any clue as to where the root lies, as given by the parameter \a start. |
---|
| 193 | /// \param lower_bracket the lower bound for the search |
---|
| 194 | /// \param upper_bracket the upper bound for the search. Function sign must be |
---|
| 195 | /// opposite to that at \a lower_bracket |
---|
| 196 | /// \param start starting value for the search |
---|
| 197 | /// \param value the value of the function being sought (solves f(x) = \a value) |
---|
| 198 | /// \param[out] error If pointer is zero, errors raise exception. Otherwise, returns error here. |
---|
| 199 | /// \param[out] final_yprime If pointer is not zero, return derivative of function at root |
---|
| 200 | /// \param[out] final_yprime2 If pointer is not zero, return second derivative of function at root |
---|
| 201 | /// \return the position of the root. |
---|
[1230] | 202 | /// \see \ref rootfinder_subsec "Root finding sample" |
---|
[807] | 203 | float_type find_root(float_type lower_bracket, float_type upper_bracket, float_type start, |
---|
| 204 | float_type value, int *error=0, |
---|
| 205 | float_type *final_yprime=0, float_type *final_yprime2=0 ) const throw(c2_exception) ; // solve f(x)=value |
---|
| 206 | |
---|
| 207 | /// \brief for points in xgrid, adaptively return Integral[f(x),{x,xgrid[i],xgrid[i+1]}] and return in vector, along with sum |
---|
| 208 | /// |
---|
| 209 | /// partial_integrals uses a method with an error O(dx**10) with full information from the derivatives, |
---|
| 210 | /// and falls back to lower order methods if informed of incomplete derivatives. |
---|
| 211 | /// It uses exact midpoint splitting of the intervals for recursion, resulting in no recomputation of the function |
---|
| 212 | /// during recursive descent at previously computed points. |
---|
| 213 | /// \param xgrid points between which to evaluate definite integrals. |
---|
| 214 | /// \param partials if non-NULL, a vector in which to receive the partial integrals. |
---|
| 215 | /// It will automatically be sized apprpropriately, if provided, to contain \a n - 1 elements where \a n is the length of \a xgrid |
---|
| 216 | /// \param abs_tol the absolute error bound for each segment |
---|
| 217 | /// \param rel_tol the fractional error bound for each segment. |
---|
| 218 | /// If the error is smaller than either the relative or absolute tolerance, the integration step is finished. |
---|
| 219 | /// \param derivs number of derivatives to trust, which sets the order of the integrator. The order |
---|
| 220 | /// is 3*\a derivs + 4. \a derivs can be 0, 1, or 2. |
---|
| 221 | /// \param adapt if true, use recursive adaptation, otherwise do simple evaluation on the grid provided |
---|
| 222 | /// with no error checking. |
---|
| 223 | /// \param extrapolate if true, use simple Richardson extrapolation on the final 2 steps to reduce the error. |
---|
[1230] | 224 | /// \return sum of partial integrals, which is the definite integral from the first value in \a xgrid to the last. |
---|
[807] | 225 | float_type partial_integrals(std::vector<float_type> xgrid, std::vector<float_type> *partials = 0, |
---|
[1230] | 226 | float_type abs_tol=1e-12, float_type rel_tol=1e-12, int derivs=2, bool adapt=true, bool extrapolate=true) |
---|
| 227 | const throw(c2_exception); |
---|
[807] | 228 | |
---|
| 229 | /// \brief a fully-automated integrator which uses the information provided by the get_sampling_grid() function |
---|
| 230 | /// to figure out what to do. |
---|
| 231 | /// |
---|
| 232 | /// It returns the integral of the function over the domain requested |
---|
| 233 | /// with error tolerances as specified. It is just a front-end to partial_integrals() |
---|
| 234 | /// |
---|
| 235 | /// \param xmin lower bound of the domain for integration |
---|
| 236 | /// \param xmax upper bound of the domain for integration |
---|
| 237 | /// \param partials if non-NULL, a vector in which to receive the partial integrals. |
---|
[1230] | 238 | /// It will automatically be sized appropriately, if provided, to contain \a n - 1 elements where \a n is the length of \a xgrid |
---|
[807] | 239 | /// \param abs_tol the absolute error bound for each segment |
---|
| 240 | /// \param rel_tol the fractional error bound for each segment. |
---|
| 241 | /// If the error is smaller than either the relative or absolute tolerance, the integration step is finished. |
---|
| 242 | /// \param derivs number of derivatives to trust, which sets the order of the integrator. The order |
---|
| 243 | /// is 3*\a derivs + 4. \a derivs can be 0, 1, or 2. |
---|
| 244 | /// \param adapt if true, use recursive adaptation, otherwise do simple evaluation on the grid provided |
---|
| 245 | /// with no error checking. |
---|
| 246 | /// \param extrapolate if true, use simple Richardson extrapolation on the final 2 steps to reduce the error. |
---|
[1230] | 247 | /// \return sum of partial integrals, which is the definite integral from the first value in \a xgrid to the last. |
---|
[807] | 248 | float_type integral(float_type xmin, float_type xmax, std::vector<float_type> *partials = 0, |
---|
[1230] | 249 | float_type abs_tol=1e-12, float_type rel_tol=1e-12, int derivs=2, bool adapt=true, bool extrapolate=true) |
---|
| 250 | const throw(c2_exception); |
---|
[807] | 251 | |
---|
[1230] | 252 | /// \brief create a c2_piecewise_function_p from c2_connector_function_p segments which |
---|
| 253 | /// is a representation of the parent function to the specified accuracy, but maybe much cheaper to evaluate |
---|
| 254 | /// |
---|
| 255 | /// This method has three modes, depending on the \a derivs flag. |
---|
| 256 | /// |
---|
| 257 | /// If \a derivs is 2, |
---|
| 258 | /// it computes a c2_piecewise_function_p representation of its parent function, which may be a much faster |
---|
| 259 | /// function to use in codes if the parent function is expensive. If \a xvals and \a yvals are non-null, |
---|
| 260 | /// it will also fill them in with the function values at each grid point the adaptive algorithm chooses. |
---|
| 261 | /// |
---|
| 262 | /// If \a derivs is 1, this does not create the connectors, |
---|
| 263 | /// and returns an null pointer, but will fill in the \a xvals and \a yvals |
---|
| 264 | /// vectors with values of the function at points such that the linear interpolation error between the points |
---|
| 265 | /// is bounded by the tolerance values given. Because it uses derivative information from the function to manage the |
---|
| 266 | /// error control, it is almost completely free of issues with missing periods of oscillatory functions, |
---|
| 267 | /// even with no information provided in the sampling grid. |
---|
| 268 | /// This is typically useful for sampling a function for plotting. |
---|
| 269 | /// |
---|
| 270 | /// If \a derivs is 0, this does something very like what it does if \a derivs = 1, but without derivatives. |
---|
| 271 | /// Instead, to compute the intermediate value of the function for error control, it just uses |
---|
| 272 | /// 3-point parabolic interpolation. This is useful amost exclusively for converting a non-c2_function, |
---|
| 273 | /// with no derivatives, but wrapped in a c2_classic_function wrapper, into a table of values to seed an interpolating_function_p. |
---|
| 274 | /// Note, however, that without derivatives, this is very susceptible to missing periods of oscillatory |
---|
| 275 | /// functions, so it is important to set a sampling grid which isn't too much coarser than the typical oscillations. |
---|
| 276 | /// |
---|
| 277 | /// \note the \a sampling_grid of the returned function matches the \a sampling_grid of its parent. |
---|
| 278 | /// \see \ref sample_function_for_plotting "Adaptive Sampling Examples" |
---|
| 279 | /// \param xmin lower bound of the domain for sampling |
---|
| 280 | /// \param xmax upper bound of the domain for sampling |
---|
| 281 | /// \param abs_tol the absolute error bound for each segment |
---|
| 282 | /// \param rel_tol the fractional error bound for each segment. |
---|
| 283 | /// \param derivs if 0 or 1, return a useless function, but fill in the \a xvals and \a yvals vectors (if non-null). |
---|
| 284 | /// Also, if 0 or 1, tolerances refer to linear interpolation, not high-order interpolation. |
---|
| 285 | /// If 2, return a full piecewise collection of c2_connector_function_p segments. See discussion above. |
---|
| 286 | /// \param [in,out] xvals vector of abscissas at which the function was actually sampled (if non-null) |
---|
| 287 | /// \param [in,out] yvals vector of function values corresponding to \a xvals (if non-null) |
---|
| 288 | /// \return a new, sampled representation, if \a derivs is 2. A null pointer if \a derivs is 0 or 1. |
---|
| 289 | c2_piecewise_function_p<float_type> *adaptively_sample(float_type xmin, float_type xmax, |
---|
| 290 | float_type abs_tol=1e-12, float_type rel_tol=1e-12, |
---|
| 291 | int derivs=2, std::vector<float_type> *xvals=0, std::vector<float_type> *yvals=0) const throw(c2_exception); |
---|
| 292 | |
---|
[807] | 293 | /// \brief return the lower bound of the domain for this function as set by set_domain() |
---|
| 294 | inline float_type xmin() const { return fXMin; } |
---|
| 295 | /// \brief return the upper bound of the domain for this function as set by set_domain() |
---|
| 296 | inline float_type xmax() const { return fXMax; } |
---|
| 297 | /// \brief set the domain for this function. |
---|
| 298 | void set_domain(float_type xmin, float_type xmax) { fXMin=xmin; fXMax=xmax; } |
---|
| 299 | |
---|
| 300 | /// \brief this is a counter owned by the function but which can be used to monitor efficiency of algorithms. |
---|
| 301 | /// |
---|
[1230] | 302 | /// It is not maintained automatically in general! The root finder, integrator, and sampler do increment it. |
---|
[807] | 303 | /// \return number of evaluations logged since last reset. |
---|
[1230] | 304 | volatile size_t get_evaluations() const { return evaluations; } |
---|
[807] | 305 | /// \brief reset the counter |
---|
| 306 | void reset_evaluations() const { evaluations=0; } // evaluations are 'invisible' to constant |
---|
| 307 | /// \brief count evaluations |
---|
| 308 | inline void increment_evaluations() const { evaluations++; } |
---|
| 309 | |
---|
| 310 | /// \brief check that a vector is monotonic, throw an exception if not, and return a flag if it is reversed |
---|
| 311 | /// |
---|
| 312 | /// \param data a vector of data points which are expected to be monotonic. |
---|
| 313 | /// \param message an informative string to include in an exception if this throws c2_exception |
---|
| 314 | /// \return true if in decreasing order, false if increasing |
---|
[1230] | 315 | bool check_monotonicity(const std::vector<float_type> &data, const char message[]) const throw(c2_exception); |
---|
[807] | 316 | |
---|
| 317 | /// \brief establish a grid of 'interesting' points on the function. |
---|
| 318 | /// |
---|
| 319 | /// The sampling grid describes a reasonable initial set of points to look at the function. |
---|
| 320 | /// this should generally be set at a scale which is quite coarse, and sufficient for initializing |
---|
| 321 | /// adaptive integration or possibly root bracketing. For sampling a function to build a new interpolating |
---|
| 322 | /// function, one may want to refine this for accuracy. However, interpolating_functions themselves |
---|
| 323 | /// return their original X grid by default, so refining the grid in this case might be a bad idea. |
---|
| 324 | /// \param grid a vector of abscissas. The contents is copied into an internal vector, so the \a grid can be discarded after passingin. |
---|
| 325 | virtual void set_sampling_grid(const std::vector<float_type> &grid) throw(c2_exception); |
---|
| 326 | |
---|
[1230] | 327 | /// \brief get the sampling grid, which may be a null pointer |
---|
| 328 | /// \return pointer to the sampling grid |
---|
| 329 | std::vector<float_type> *get_sampling_grid_pointer() const { return sampling_grid; } |
---|
| 330 | |
---|
[807] | 331 | /// \brief return the grid of 'interesting' points along this function which lie in the region requested |
---|
| 332 | /// |
---|
| 333 | /// if a sampling grid is defined, work from there, otherwise return vector of (xmin, xmax) |
---|
| 334 | /// \param xmin the lower bound for which the function is to be sampled |
---|
| 335 | /// \param xmax the upper bound for which the function is to be sampled |
---|
[1230] | 336 | /// \param [in,out] grid filled vector containing the samplng grid. |
---|
| 337 | virtual void get_sampling_grid(float_type xmin, float_type xmax, std::vector<float_type> &grid) const ; |
---|
[807] | 338 | |
---|
| 339 | /// \brief clean up endpoints on a grid of points |
---|
| 340 | /// \param[in,out] result the sampling grid with excessively closely space endpoints removed. |
---|
| 341 | /// The grid is modified in place. |
---|
| 342 | void preen_sampling_grid(std::vector<float_type> *result) const; |
---|
| 343 | /// \brief refine a grid by splitting each interval into more intervals |
---|
[1230] | 344 | /// \param [in,out] grid the grid to refine in place |
---|
[807] | 345 | /// \param refinement the number of new steps for each old step |
---|
[1230] | 346 | void refine_sampling_grid(std::vector<float_type> &grid, size_t refinement) const; |
---|
[807] | 347 | |
---|
| 348 | /// \brief create a new c2_function from this one which is normalized on the interval |
---|
| 349 | /// \param xmin lower bound of the domain for integration |
---|
| 350 | /// \param xmax upper bound of the domain for integration |
---|
| 351 | /// \param norm the desired integral for the function over the region |
---|
| 352 | /// \return a new c2_function with the desired \a norm. |
---|
[1230] | 353 | c2_function<float_type> &normalized_function(float_type xmin, float_type xmax, float_type norm=1.0) const throw(c2_exception); |
---|
[807] | 354 | /// \brief create a new c2_function from this one which is square-normalized on the interval |
---|
| 355 | /// \param xmin lower bound of the domain for integration |
---|
| 356 | /// \param xmax upper bound of the domain for integration |
---|
| 357 | /// \param norm the desired integral for the function over the region |
---|
| 358 | /// \return a new c2_function with the desired \a norm. |
---|
[1230] | 359 | c2_function<float_type> &square_normalized_function(float_type xmin, float_type xmax, float_type norm=1.0) const throw(c2_exception); |
---|
[807] | 360 | /// \brief create a new c2_function from this one which is square-normalized with the provided \a weight on the interval |
---|
| 361 | /// \param xmin lower bound of the domain for integration |
---|
| 362 | /// \param xmax upper bound of the domain for integration |
---|
| 363 | /// \param weight a c2_function providing the weight |
---|
| 364 | /// \param norm the desired integral for the function over the region |
---|
| 365 | /// \return a new c2_function with the desired \a norm. |
---|
[1230] | 366 | c2_function<float_type> &square_normalized_function( |
---|
| 367 | float_type xmin, float_type xmax, const c2_function<float_type> &weight, float_type norm=1.0) |
---|
| 368 | const throw(c2_exception); |
---|
[807] | 369 | |
---|
[1230] | 370 | /// \brief factory function to create a c2_sum_p from a regular algebraic expression. |
---|
| 371 | /// \param rhs the right-hand term of the sum |
---|
| 372 | /// \return a new c2_function |
---|
| 373 | c2_sum_p<float_type> &operator + (const c2_function<float_type> &rhs) const |
---|
| 374 | { return *new c2_sum_p<float_type>(*this, rhs); } |
---|
| 375 | /// \brief factory function to create a c2_diff_p from a regular algebraic expression. |
---|
| 376 | /// \param rhs the right-hand term of the difference |
---|
| 377 | /// \return a new c2_function |
---|
| 378 | c2_diff_p<float_type> &operator - (const c2_function<float_type> &rhs) const |
---|
| 379 | { return *new c2_diff_p<float_type>(*this, rhs); } |
---|
| 380 | /// \brief factory function to create a c2_product_p from a regular algebraic expression. |
---|
| 381 | /// \param rhs the right-hand term of the product |
---|
| 382 | /// \return a new c2_function |
---|
| 383 | c2_product_p<float_type> &operator * (const c2_function<float_type> &rhs) const |
---|
| 384 | { return *new c2_product_p<float_type>(*this, rhs); } |
---|
| 385 | /// \brief factory function to create a c2_ratio_p from a regular algebraic expression. |
---|
| 386 | /// \param rhs the right-hand term of the ratio (the denominator) |
---|
| 387 | /// \return a new c2_function |
---|
| 388 | c2_ratio_p<float_type> &operator / (const c2_function<float_type> &rhs) const |
---|
| 389 | { return *new c2_ratio_p<float_type>(*this, rhs); } |
---|
| 390 | /// \brief compose this function outside another. |
---|
| 391 | /// \param inner the inner function |
---|
| 392 | /// \return the composed function |
---|
| 393 | /// \anchor compose_operator |
---|
| 394 | c2_composed_function_p<float_type> & operator ()(const c2_function<float_type> &inner) const |
---|
| 395 | { return *new c2_composed_function_p<float_type>((*this), inner); } |
---|
| 396 | |
---|
| 397 | /// \brief Find out where a calculation ran into trouble, if it got a nan. |
---|
| 398 | /// If the most recent computation did not return a nan, this is undefined. |
---|
| 399 | /// \return \a x value of point at which something went wrong, if integrator (or otherwise) returned a nan. |
---|
| 400 | float_type get_trouble_point() const { return bad_x_point; } |
---|
[807] | 401 | |
---|
[1230] | 402 | /// \brief increment our reference count. Destruction is only legal if the count is zero. |
---|
| 403 | void claim_ownership() const { owner_count++; } |
---|
| 404 | /// \brief decrement our reference count. Do not destroy at zero. |
---|
| 405 | /// \return final owner count, to check whether object should disappear. |
---|
| 406 | size_t release_ownership_for_return() const throw(c2_exception) { |
---|
| 407 | if(!owner_count) { |
---|
| 408 | std::ostringstream outstr; |
---|
| 409 | outstr << "attempt to release ownership of an unowned function in class "; |
---|
| 410 | outstr << typeid(*this).name() << std::endl; |
---|
| 411 | throw c2_exception(outstr.str().c_str()); |
---|
| 412 | } |
---|
| 413 | owner_count--; |
---|
| 414 | return owner_count; |
---|
| 415 | } |
---|
| 416 | /// \brief decrement our reference count. If the count reaches zero, destroy ourself. |
---|
| 417 | void release_ownership() const throw(c2_exception) { |
---|
| 418 | if(!release_ownership_for_return()) delete this; |
---|
| 419 | } |
---|
| 420 | /// \brief get the reference count, mostly for debugging |
---|
| 421 | /// \return the count |
---|
| 422 | size_t count_owners() const { return owner_count; } |
---|
[807] | 423 | |
---|
| 424 | protected: |
---|
| 425 | c2_function(const c2_function<float_type> &src) : sampling_grid(0), |
---|
| 426 | no_overwrite_grid(false), |
---|
[1230] | 427 | fXMin(src.fXMin), fXMax(src.fXMax), root_info(0), owner_count(0) |
---|
[807] | 428 | {} // copy constructor only copies domain, and is only for internal use |
---|
| 429 | c2_function() : |
---|
| 430 | sampling_grid(0), no_overwrite_grid(0), |
---|
| 431 | fXMin(-std::numeric_limits<float_type>::max()), |
---|
[1230] | 432 | fXMax(std::numeric_limits<float_type>::max()), root_info(0), owner_count(0) |
---|
[807] | 433 | {} // prevent accidental naked construction (impossible any since this has pure virtual methods) |
---|
| 434 | |
---|
| 435 | // this should only be called very early on, by a constructor, before anyone else |
---|
| 436 | // sets a sampling grid, or it will leak memory |
---|
| 437 | virtual void set_sampling_grid_pointer(std::vector<float_type> &grid) |
---|
| 438 | { |
---|
| 439 | if (sampling_grid && !no_overwrite_grid) delete sampling_grid; // grid was ours, lose it. |
---|
| 440 | sampling_grid=&grid; no_overwrite_grid=1; |
---|
| 441 | } |
---|
| 442 | |
---|
[1230] | 443 | std::vector<float_type> * sampling_grid; |
---|
| 444 | bool no_overwrite_grid; |
---|
| 445 | |
---|
[807] | 446 | float_type fXMin, fXMax; |
---|
[1230] | 447 | mutable size_t evaluations; |
---|
| 448 | /// \brief this point may be used to record where a calculation ran into trouble |
---|
| 449 | mutable float_type bad_x_point; |
---|
| 450 | public: |
---|
| 451 | /// \brief fill in a c2_fblock<float_type>... a shortcut for the integrator & sampler |
---|
| 452 | /// \param [in,out] fb the block to fill in with information |
---|
| 453 | inline void fill_fblock(c2_fblock<float_type> &fb) const throw(c2_exception) |
---|
| 454 | { |
---|
| 455 | fb.y=value_with_derivatives(fb.x, &fb.yp, &fb.ypp); |
---|
| 456 | fb.ypbad=c2_isnan(fb.yp) || !c2_isfinite(fb.yp); |
---|
| 457 | fb.yppbad=c2_isnan(fb.ypp) || !c2_isfinite(fb.ypp); |
---|
| 458 | increment_evaluations(); |
---|
| 459 | } |
---|
| 460 | |
---|
| 461 | private: |
---|
| 462 | /// \brief the data element for the internal recursion stack for the sampler and integrator |
---|
| 463 | struct recur_item { |
---|
| 464 | c2_fblock<float_type> f1; size_t depth; |
---|
| 465 | float_type previous_estimate, abs_tol, step_sum; |
---|
| 466 | bool done; |
---|
| 467 | size_t f0index, f2index; |
---|
| 468 | }; |
---|
[807] | 469 | |
---|
[1230] | 470 | |
---|
| 471 | /// \brief structure used to pass information recursively in integrator. |
---|
[807] | 472 | /// |
---|
| 473 | /// the \a abs_tol is scaled by a factor of two at each division. |
---|
| 474 | /// Everything else is just passed down. |
---|
[1230] | 475 | struct c2_integrate_recur { |
---|
| 476 | c2_fblock<float_type> *f0, *f1; |
---|
| 477 | float_type abs_tol, rel_tol, eps_scale, extrap_coef, extrap2, dx_tolerance, abs_tol_min; |
---|
| 478 | std::vector< recur_item > *rb_stack; |
---|
| 479 | int derivs; |
---|
| 480 | bool adapt, extrapolate, inited; |
---|
[807] | 481 | }; |
---|
[1230] | 482 | |
---|
| 483 | /// \brief structure used to pass information recursively in sampler. |
---|
| 484 | /// |
---|
| 485 | struct c2_sample_recur { |
---|
| 486 | c2_fblock<float_type> *f0, *f1; |
---|
| 487 | float_type abs_tol, rel_tol, dx_tolerance, abs_tol_min; |
---|
| 488 | int derivs; |
---|
| 489 | c2_piecewise_function_p<float_type> *out; |
---|
| 490 | std::vector<float_type> *xvals, *yvals; |
---|
| 491 | std::vector< recur_item > *rb_stack; |
---|
| 492 | bool inited; |
---|
| 493 | }; |
---|
| 494 | |
---|
| 495 | /// \brief structure used to hold root bracketing information |
---|
| 496 | /// |
---|
| 497 | struct c2_root_info { |
---|
| 498 | c2_fblock<float_type> lower, upper; |
---|
| 499 | bool inited; |
---|
| 500 | }; |
---|
[807] | 501 | |
---|
| 502 | /// \brief Carry out the recursive subdivision and integration. |
---|
| 503 | /// |
---|
| 504 | /// This passes information recursively through the \a recur block pointer |
---|
| 505 | /// to allow very efficient recursion. |
---|
| 506 | /// \param rb a pointer to the recur struct. |
---|
[1230] | 507 | float_type integrate_step(struct c2_integrate_recur &rb) const throw(c2_exception); |
---|
[807] | 508 | |
---|
[1230] | 509 | /// \brief Carry out the recursive subdivision for sampling. |
---|
| 510 | /// |
---|
| 511 | /// This passes information recursively through the \a recur block pointer |
---|
| 512 | /// to allow very efficient recursion. |
---|
| 513 | /// \param rb a pointer to the recur struct. |
---|
| 514 | void sample_step(struct c2_sample_recur &rb) const throw(c2_exception); |
---|
| 515 | |
---|
| 516 | /// this carry a memory of the last root bracketing, |
---|
[807] | 517 | /// to avoid the necessity of evaluating the function on the brackets every time |
---|
| 518 | /// if the brackets have not been changed. |
---|
[1230] | 519 | /// it is declared as a pointer, since many c2_functions may never need one allocated |
---|
| 520 | mutable struct c2_root_info *root_info; |
---|
[807] | 521 | |
---|
[1230] | 522 | mutable size_t owner_count; |
---|
[807] | 523 | }; |
---|
| 524 | |
---|
[1230] | 525 | /// \brief a container into which any conventional c-style function can be dropped, |
---|
| 526 | /// to create a degenerate c2_function without derivatives. |
---|
| 527 | /// Mostly useful for sampling into interpolating functions. |
---|
| 528 | /// construct a reference to this with c2_classic_function() |
---|
| 529 | /// \ingroup containers |
---|
| 530 | /// The factory function c2_factory::classic_function() creates *new c2_classic_function_p() |
---|
| 531 | template <typename float_type=double> class c2_classic_function_p : public c2_function<float_type> { |
---|
| 532 | public: |
---|
| 533 | /// \brief construct the container |
---|
| 534 | /// \param c_func a pointer to a conventional c-style function |
---|
| 535 | c2_classic_function_p(const float_type (*c_func)(float_type)) : c2_function<float_type>(), func(c_func) {} |
---|
| 536 | |
---|
| 537 | /// \copydoc c2_function::value_with_derivatives |
---|
| 538 | /// Uses the internal function pointer set by set_function(). |
---|
| 539 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
| 540 | { |
---|
| 541 | if(!func) throw c2_exception("c2_classic_function called with null function"); |
---|
| 542 | if(yprime) *yprime=0; |
---|
| 543 | if(yprime2) *yprime2=0; |
---|
| 544 | return func(x); |
---|
| 545 | } |
---|
| 546 | ~c2_classic_function_p() { } |
---|
| 547 | |
---|
| 548 | protected: |
---|
| 549 | /// \brief pointer to our function |
---|
| 550 | const float_type (*func)(float_type); |
---|
| 551 | }; |
---|
| 552 | |
---|
| 553 | /// \brief create a container for a c2_function which handles the reference counting. |
---|
| 554 | /// \ingroup containers |
---|
| 555 | /// It is useful as a smart container to hold a c2_function and keep the reference count correct. |
---|
| 556 | /// The recommended way for a class to store a c2_function which is handed in from the outside |
---|
| 557 | /// is for it to have a c2_ptr member into which the passed-in function is stored. |
---|
| 558 | /// This way, when the class instance is deleted, it will automatically dereference any function |
---|
| 559 | /// which it was handed. |
---|
| 560 | /// |
---|
| 561 | /// This class contains a copy constructor and operator=, to make it fairly easy to make |
---|
| 562 | /// a std::vector of these objects, and have it work as expected. |
---|
| 563 | template <typename float_type> class c2_const_ptr { |
---|
| 564 | public: |
---|
| 565 | /// \brief construct the container with no function |
---|
| 566 | c2_const_ptr() : func(0) {} |
---|
| 567 | /// \brief construct the container with a pre-defined function |
---|
| 568 | /// \param f the function to store |
---|
| 569 | c2_const_ptr(const c2_function<float_type> &f) : func(0) |
---|
| 570 | { set_function(&f); } |
---|
| 571 | /// \brief copy constructor |
---|
| 572 | /// \param src the container to copy |
---|
| 573 | c2_const_ptr(const c2_const_ptr<float_type> &src) : func(0) |
---|
| 574 | { set_function(src.get_ptr()); } |
---|
| 575 | /// \brief fill the container with a new function, or clear it with a null pointer |
---|
| 576 | /// \param f the function to store, releasing any previously held function |
---|
| 577 | void set_function(const c2_function<float_type> *f) |
---|
| 578 | { |
---|
| 579 | if(func) func->release_ownership(); |
---|
| 580 | func=f; |
---|
| 581 | if(func) func->claim_ownership(); |
---|
| 582 | } |
---|
| 583 | |
---|
| 584 | /// \brief fill the container from another container |
---|
| 585 | /// \param f the container to copy |
---|
| 586 | void operator =(const c2_const_ptr<float_type> &f) |
---|
| 587 | { set_function(f.get_ptr()); } |
---|
| 588 | /// \brief fill the container with a function |
---|
| 589 | /// \param f the function |
---|
| 590 | void operator =(const c2_function<float_type> &f) |
---|
| 591 | { set_function(&f); } |
---|
| 592 | /// \brief release the function without destroying it, so it can be returned from a function |
---|
| 593 | /// |
---|
| 594 | /// This is usually the very last line of a function before the return statement, so that |
---|
| 595 | /// any exceptions that happen during execution of the function will cause proper cleanup. |
---|
| 596 | /// Once the function has been released from its container this way, it is an orhpaned object |
---|
| 597 | /// until the caller claims it, so it could get lost if an exception happens. |
---|
| 598 | void release_for_return() throw(c2_exception) |
---|
| 599 | { |
---|
| 600 | if(func) func->release_ownership_for_return(); |
---|
| 601 | func=0; |
---|
| 602 | } |
---|
| 603 | /// \brief clear the function |
---|
| 604 | /// |
---|
| 605 | /// Any attempt to use this c2_plugin_function_p throws an exception if the saved function is cleared. |
---|
| 606 | void unset_function(void) { set_function(0); } |
---|
| 607 | /// \brief destructor |
---|
| 608 | ~c2_const_ptr() { set_function(0); } |
---|
| 609 | |
---|
| 610 | /// \brief get a reference to our owned function |
---|
| 611 | inline const c2_function<float_type> &get() const throw(c2_exception) |
---|
| 612 | { |
---|
| 613 | if(!func) throw c2_exception("c2_ptr accessed uninitialized"); |
---|
| 614 | return *func; |
---|
| 615 | } |
---|
| 616 | /// \brief get an unchecked pointer to our owned function |
---|
| 617 | inline const c2_function<float_type> *get_ptr() const { return func; } |
---|
| 618 | /// \brief get a checked pointer to our owned function |
---|
| 619 | inline const c2_function<float_type> *operator -> () const |
---|
| 620 | { return &get(); } |
---|
| 621 | /// \brief check if we have a valid function |
---|
| 622 | bool valid() const { return func != 0; } |
---|
| 623 | |
---|
| 624 | /// \brief type coercion operator which lets us use a pointer as if it were a const c2_function |
---|
| 625 | operator const c2_function<float_type>& () const { return this->get(); } |
---|
| 626 | |
---|
| 627 | /// \brief convenience operator to make us look like a function |
---|
| 628 | /// \param x the value at which to evaluate the contained function |
---|
| 629 | /// \return the evaluated function |
---|
| 630 | /// \note If you using this repeatedly, do const c2_function<float_type> &func=ptr; |
---|
| 631 | /// and use func(x). Calling this operator wastes some time, since it checks the validity of the |
---|
| 632 | /// pointer every time. |
---|
| 633 | float_type operator()(float_type x) const throw(c2_exception) { return get()(x); } |
---|
| 634 | /// \brief convenience operator to make us look like a function |
---|
| 635 | /// \param x the value at which to evaluate the contained function |
---|
| 636 | /// \param yprime the derivative |
---|
| 637 | /// \param yprime2 the second derivative |
---|
| 638 | /// \return the evaluated function |
---|
| 639 | /// \note If you using this repeatedly, do const c2_function<float_type> &func=ptr; |
---|
| 640 | /// and use func(x). Calling this operator wastes some time, since it checks the validity of the |
---|
| 641 | /// pointer every time. |
---|
| 642 | float_type operator()(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
| 643 | { return get().value_with_derivatives(x, yprime, yprime2); } |
---|
| 644 | /// \brief factory function to create a c2_sum_p from a regular algebraic expression. |
---|
| 645 | /// \param rhs the right-hand term of the sum |
---|
| 646 | /// \return a new c2_function |
---|
| 647 | c2_sum_p<float_type> &operator + (const c2_function<float_type> &rhs) const throw(c2_exception) |
---|
| 648 | { return *new c2_sum_p<float_type>(get(), rhs); } |
---|
| 649 | /// \brief factory function to create a c2_diff_p from a regular algebraic expression. |
---|
| 650 | /// \param rhs the right-hand term of the difference |
---|
| 651 | /// \return a new c2_function |
---|
| 652 | c2_diff_p<float_type> &operator - (const c2_function<float_type> &rhs) const throw(c2_exception) |
---|
| 653 | { return *new c2_diff_p<float_type>(get(), rhs); } |
---|
| 654 | /// \brief factory function to create a c2_product_p from a regular algebraic expression. |
---|
| 655 | /// \param rhs the right-hand term of the product |
---|
| 656 | /// \return a new c2_function |
---|
| 657 | c2_product_p<float_type> &operator * (const c2_function<float_type> &rhs) const throw(c2_exception) |
---|
| 658 | { return *new c2_product_p<float_type>(get(), rhs); } |
---|
| 659 | /// \brief factory function to create a c2_ratio_p from a regular algebraic expression. |
---|
| 660 | /// \param rhs the right-hand term of the ratio (the denominator) |
---|
| 661 | /// \return a new c2_function |
---|
| 662 | c2_ratio_p<float_type> &operator / (const c2_function<float_type> &rhs) const throw(c2_exception) |
---|
| 663 | { return *new c2_ratio_p<float_type>(get(), rhs); } |
---|
| 664 | /// \brief compose this function outside another. |
---|
| 665 | /// \param inner the inner function |
---|
| 666 | /// \return the composed function |
---|
| 667 | c2_composed_function_p<float_type> & operator ()(const c2_function<float_type> &inner) const throw(c2_exception) |
---|
| 668 | { return *new c2_composed_function_p<float_type>(get(), inner); } |
---|
| 669 | |
---|
| 670 | protected: |
---|
| 671 | const c2_function<float_type> * func; |
---|
| 672 | }; |
---|
| 673 | |
---|
| 674 | /// \brief create a container for a c2_function which handles the reference counting. |
---|
| 675 | /// \ingroup containers |
---|
| 676 | /// |
---|
| 677 | /// \see c2_const_ptr and \ref memory_management "Use of c2_ptr for memory management" |
---|
| 678 | |
---|
| 679 | template <typename float_type> class c2_ptr : public c2_const_ptr<float_type > |
---|
| 680 | { |
---|
| 681 | public: |
---|
| 682 | /// \brief construct the container with no function |
---|
| 683 | c2_ptr() : c2_const_ptr<float_type>() {} |
---|
| 684 | /// \brief construct the container with a pre-defined function |
---|
| 685 | /// \param f the function to store |
---|
| 686 | c2_ptr(c2_function<float_type> &f) : |
---|
| 687 | c2_const_ptr<float_type>() { set_function(&f); } |
---|
| 688 | /// \brief copy constructor |
---|
| 689 | /// \param src the container to copy |
---|
| 690 | c2_ptr(const c2_ptr<float_type> &src) : |
---|
| 691 | c2_const_ptr<float_type>() { set_function(src.get_ptr()); } |
---|
| 692 | /// \brief get a checked pointer to our owned function |
---|
| 693 | inline c2_function<float_type> &get() const throw(c2_exception) |
---|
| 694 | { return *const_cast<c2_function<float_type>*>(&c2_const_ptr<float_type>::get()); } |
---|
| 695 | /// \brief get an unchecked pointer to our owned function |
---|
| 696 | inline c2_function<float_type> *get_ptr() const |
---|
| 697 | { return const_cast<c2_function<float_type>*>(this->func); } |
---|
| 698 | /// \brief get a checked pointer to our owned function |
---|
| 699 | inline c2_function<float_type> *operator -> () const |
---|
| 700 | { return &get(); } |
---|
| 701 | /// \brief fill the container from another container |
---|
| 702 | /// \param f the container to copy |
---|
| 703 | void operator =(const c2_ptr<float_type> &f) |
---|
| 704 | { set_function(f.get_ptr()); } |
---|
| 705 | /// \brief fill the container with a function |
---|
| 706 | /// \param f the function |
---|
| 707 | void operator =(c2_function<float_type> &f) |
---|
| 708 | { set_function(&f); } |
---|
| 709 | private: |
---|
| 710 | /// \brief hidden non-const-safe version of operator= |
---|
| 711 | void operator =(const c2_const_ptr<float_type> &f) { } |
---|
| 712 | /// \brief hidden non-const-safe version of operator= |
---|
| 713 | void operator =(const c2_function<float_type> &f) { } |
---|
| 714 | }; |
---|
| 715 | |
---|
| 716 | /// \brief create a non-generic container for a c2_function which handles the reference counting. |
---|
| 717 | /// \ingroup containers |
---|
| 718 | /// |
---|
| 719 | /// \see c2_const_ptr and \ref memory_management "Use of c2_ptr for memory management" |
---|
| 720 | /// |
---|
| 721 | /// \note Overuse of this class will generate massive bloat. Use c2_ptr and c2_const_ptr if you don't _really_ need specific pointer types. |
---|
| 722 | /// \see \ref memory_management "Use of c2_ptr for memory management" |
---|
| 723 | /* |
---|
| 724 | template <typename float_type, template <typename> class c2_class > class c2_typed_ptr : public c2_const_ptr<float_type> { |
---|
| 725 | public: |
---|
| 726 | /// \brief construct the container with no function |
---|
| 727 | c2_typed_ptr() : c2_ptr<float_type>() {} |
---|
| 728 | /// \brief construct the container with a pre-defined function |
---|
| 729 | /// \param f the function to store |
---|
| 730 | c2_typed_ptr(c2_class<float_type> &f) |
---|
| 731 | : c2_const_ptr<float_type>() { this->set_function(&f); } |
---|
| 732 | /// \brief copy constructor |
---|
| 733 | /// \param src the container to copy |
---|
| 734 | c2_typed_ptr(const c2_typed_ptr<float_type, c2_class> &src) |
---|
| 735 | : c2_const_ptr<float_type>() { this->set_function(src.get_ptr()); } |
---|
| 736 | |
---|
| 737 | /// \brief get a reference to our owned function |
---|
| 738 | inline c2_class<float_type> &get() const throw(c2_exception) |
---|
| 739 | { |
---|
| 740 | return *static_cast<c2_class<float_type> *>(const_cast<c2_function<float_type>*>(&c2_const_ptr<float_type>::get())); |
---|
| 741 | } |
---|
| 742 | /// \brief get a checked pointer to our owned function |
---|
| 743 | inline c2_class<float_type> *operator -> () const |
---|
| 744 | { return &get(); } |
---|
| 745 | /// \brief get an unchecked pointer to our owned function |
---|
| 746 | inline c2_class<float_type> *get_ptr() const |
---|
| 747 | { return static_cast<c2_class<float_type> *>(const_cast<c2_function<float_type>*>(this->func)); } |
---|
| 748 | /// \brief type coercion operator which lets us use a pointer as if it were a c2_function |
---|
| 749 | operator c2_class<float_type>&() const { return get(); } |
---|
| 750 | /// \brief fill the container from another container |
---|
| 751 | /// \param f the container to copy |
---|
| 752 | void operator =(const c2_typed_ptr<float_type, c2_class> &f) |
---|
| 753 | { set_function(f.get_ptr()); } |
---|
| 754 | /// \brief fill the container with a function |
---|
| 755 | /// \param f the function |
---|
| 756 | void operator =(c2_class<float_type> &f) |
---|
| 757 | { set_function(&f); } |
---|
| 758 | private: |
---|
| 759 | /// \brief hidden downcasting version of operator= |
---|
| 760 | void operator =(const c2_const_ptr<float_type> &f) { } |
---|
| 761 | /// \brief hidden downcasting version of operator=. Use an explicit dynamic_cast<c2_class<float_type>&>(f) if you need to try this. |
---|
| 762 | void operator =(const c2_function<float_type> &f) { } |
---|
| 763 | }; |
---|
| 764 | */ |
---|
[807] | 765 | /// \brief a container into which any other c2_function can be dropped, to allow expressions |
---|
| 766 | /// with replacable components. |
---|
[1230] | 767 | /// \ingroup containers |
---|
[807] | 768 | ///It is useful for plugging different InterpolatingFunctions into a c2_function expression. |
---|
| 769 | ///It saves a lot of effort in other places with casting away const declarations. |
---|
| 770 | /// |
---|
| 771 | /// It is also useful as a wrapper for a function if it is necessary to have a copy of a function |
---|
| 772 | /// which has a different domain or sampling grid than the parent function. This can be |
---|
[1230] | 773 | /// be used, for example, to patch badly-behaved functions with c2_piecewise_function_p by |
---|
[807] | 774 | /// taking the parent function, creating two plugins of it with domains on each side of the |
---|
| 775 | /// nasty bit, and then inserting a nice function in the hole. |
---|
[1230] | 776 | /// |
---|
| 777 | /// This can also be used as a fancier c2_ptr which allows direct evaluation |
---|
| 778 | /// instead of having to dereference the container first. |
---|
| 779 | /// |
---|
| 780 | /// The factory function c2_factory::plugin_function() creates *new c2_plugin_function_p() |
---|
| 781 | template <typename float_type=double> class c2_plugin_function_p : |
---|
| 782 | public c2_function<float_type> { |
---|
[807] | 783 | public: |
---|
| 784 | /// \brief construct the container with no function |
---|
[1230] | 785 | c2_plugin_function_p() : c2_function<float_type>(), func() {} |
---|
[807] | 786 | /// \brief construct the container with a pre-defined function |
---|
[1230] | 787 | c2_plugin_function_p(c2_function<float_type> &f) : |
---|
| 788 | c2_function<float_type>(),func(f) { } |
---|
| 789 | /// \brief fill the container with a new function, or clear it with a null pointer |
---|
| 790 | /// and copy our domain |
---|
| 791 | void set_function(c2_function<float_type> *f) |
---|
| 792 | { |
---|
| 793 | func.set_function(f); |
---|
| 794 | if(f) set_domain(f->xmin(), f->xmax()); |
---|
[807] | 795 | } |
---|
| 796 | /// \copydoc c2_function::value_with_derivatives |
---|
| 797 | /// Uses the internal function pointer set by set_function(). |
---|
| 798 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
| 799 | { |
---|
[1230] | 800 | if(!func.valid()) throw c2_exception("c2_plugin_function_p called uninitialized"); |
---|
| 801 | return func->value_with_derivatives(x, yprime, yprime2); |
---|
[807] | 802 | } |
---|
| 803 | /// \brief destructor |
---|
[1230] | 804 | ~c2_plugin_function_p() { } |
---|
[807] | 805 | |
---|
[1230] | 806 | /// \brief clear our function |
---|
| 807 | void unset_function() { func.unset_function(); } |
---|
| 808 | |
---|
| 809 | virtual void get_sampling_grid(float_type xmin, float_type xmax, std::vector<float_type> &grid) const |
---|
| 810 | { |
---|
| 811 | if(!func.valid()) throw c2_exception("c2_plugin_function_p called uninitialized"); |
---|
| 812 | if(this->sampling_grid) c2_function<float_type>::get_sampling_grid(xmin, xmax, grid); |
---|
| 813 | else func->get_sampling_grid(xmin, xmax, grid); |
---|
| 814 | } |
---|
[807] | 815 | protected: |
---|
[1230] | 816 | c2_ptr<float_type> func; |
---|
| 817 | }; |
---|
| 818 | |
---|
| 819 | /// \brief a c2_plugin_function_p which promises not to fiddle with the plugged function. |
---|
| 820 | /// \ingroup containers |
---|
| 821 | /// |
---|
| 822 | /// The factory function c2_factory::const_plugin_function() creates *new c2_const_plugin_function_p() |
---|
| 823 | template <typename float_type=double> class c2_const_plugin_function_p : public c2_plugin_function_p<float_type> { |
---|
| 824 | public: |
---|
| 825 | /// \brief construct the container with no function |
---|
| 826 | c2_const_plugin_function_p() : c2_plugin_function_p<float_type>() {} |
---|
| 827 | /// \brief construct the container with a pre-defined function |
---|
| 828 | c2_const_plugin_function_p(const c2_function<float_type> &f) : |
---|
| 829 | c2_plugin_function_p<float_type>() { set_function(&f); } |
---|
| 830 | /// \brief fill the container with a new function, or clear it with a null pointer |
---|
| 831 | void set_function(const c2_function<float_type> *f) |
---|
| 832 | { c2_plugin_function_p<float_type>::set_function(const_cast<c2_function<float_type>*>(f)); } |
---|
| 833 | /// \brief destructor |
---|
| 834 | ~c2_const_plugin_function_p() { } |
---|
[807] | 835 | |
---|
[1230] | 836 | /// \brief get a const reference to our owned function, for direct access |
---|
| 837 | const c2_function<float_type> &get() const throw(c2_exception) |
---|
| 838 | { return this->func.get(); } |
---|
[807] | 839 | }; |
---|
| 840 | |
---|
| 841 | /// \brief Provides support for c2_function objects which are constructed from two other c2_function |
---|
| 842 | /// objects. |
---|
[1230] | 843 | /// \ingroup abstract_classes |
---|
[807] | 844 | template <typename float_type=double> class c2_binary_function : public c2_function<float_type> { |
---|
| 845 | public: |
---|
| 846 | /// \brief function to manage the binary operation, used by c2_binary_function::value_with_derivatives() |
---|
| 847 | /// |
---|
| 848 | |
---|
| 849 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw (c2_exception) |
---|
| 850 | { |
---|
[1230] | 851 | if(stub) throw c2_exception("attempt to evaluate a c2_binary_function stub"); |
---|
| 852 | return this->combine(*Left.get_ptr(), *Right.get_ptr(), x, yprime, yprime2); |
---|
[807] | 853 | } |
---|
| 854 | |
---|
[1230] | 855 | /// \brief destructor releases ownership of member functions |
---|
[807] | 856 | /// |
---|
[1230] | 857 | virtual ~c2_binary_function() { } |
---|
| 858 | |
---|
[807] | 859 | protected: |
---|
| 860 | /// \brief construct the binary function |
---|
[1230] | 861 | /// \param combiner pointer to the function which actualy knows how to execute the binary |
---|
[807] | 862 | /// \param left the c2_function to be used in the left side of the binary relation |
---|
| 863 | /// \param right the c2_function to be used in the right side of the binary relation |
---|
[1230] | 864 | c2_binary_function( |
---|
| 865 | float_type (*combiner)(const c2_function<float_type> &left, const c2_function<float_type> &right, |
---|
| 866 | float_type x, float_type *yprime, float_type *yprime2), |
---|
| 867 | const c2_function<float_type> &left, const c2_function<float_type> &right) : |
---|
| 868 | c2_function<float_type>(), combine(combiner), Left(left), Right(right), stub(false) |
---|
[807] | 869 | { |
---|
[1230] | 870 | set_domain( |
---|
| 871 | (left.xmin() > right.xmin()) ? left.xmin() : right.xmin(), |
---|
| 872 | (left.xmax() < right.xmax()) ? left.xmax() : right.xmax() |
---|
| 873 | ); |
---|
[807] | 874 | } |
---|
| 875 | |
---|
| 876 | /// \brief construct a 'stub' c2_binary_function, which provides access to the combine() function |
---|
| 877 | /// \note Do not evaluate a 'stub' ever. It is only used so that combine() can be called |
---|
[1230] | 878 | c2_binary_function( |
---|
| 879 | float_type (*combiner)(const c2_function<float_type> &left, const c2_function<float_type> &right, |
---|
| 880 | float_type x, float_type *yprime, float_type *yprime2) |
---|
| 881 | ) : c2_function<float_type>(), combine(combiner), Left(), Right(), stub(true) { } |
---|
[807] | 882 | |
---|
[1230] | 883 | public: |
---|
| 884 | float_type (* const combine)(const c2_function<float_type> &left, const c2_function<float_type> &right, |
---|
| 885 | float_type x, float_type *yprime, float_type *yprime2); |
---|
| 886 | |
---|
| 887 | protected: |
---|
| 888 | const c2_const_ptr<float_type> Left, Right; |
---|
| 889 | /// \brief if true, we don't own any functions, we are just a source of a combining function. |
---|
| 890 | bool stub; |
---|
| 891 | |
---|
[807] | 892 | }; |
---|
| 893 | |
---|
| 894 | /// \brief Create a very lightweight method to return a scalar multiple of another function. |
---|
[1230] | 895 | /// \ingroup containers \ingroup arithmetic_functions \ingroup parametric_functions |
---|
| 896 | /// |
---|
| 897 | /// The factory function c2_factory::scaled_function() creates *new c2_scaled_function_p |
---|
| 898 | template <typename float_type=double> class c2_scaled_function_p : public c2_function<float_type> { |
---|
[807] | 899 | public: |
---|
| 900 | /// \brief construct the function with its scale factor. |
---|
| 901 | /// |
---|
| 902 | /// \param outer the function to be scaled |
---|
| 903 | /// \param scale the multiplicative scale factor |
---|
[1230] | 904 | c2_scaled_function_p(const c2_function<float_type> &outer, float_type scale) : |
---|
| 905 | c2_function<float_type>(), func(outer), yscale(scale) { } |
---|
[807] | 906 | |
---|
| 907 | /// \brief set a new scale factor |
---|
| 908 | /// \param scale the new factor |
---|
| 909 | void reset(float_type scale) { yscale=scale; } |
---|
| 910 | |
---|
| 911 | /// \copydoc c2_function::value_with_derivatives |
---|
| 912 | /// |
---|
| 913 | /// provide our own value_with_derivatives which bypasses the combiner for quicker operation |
---|
| 914 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw (c2_exception) |
---|
| 915 | { |
---|
| 916 | float_type y=this->func->value_with_derivatives(x, yprime, yprime2); |
---|
| 917 | if(yprime) (*yprime)*=yscale; |
---|
| 918 | if(yprime2) (*yprime2)*=yscale; |
---|
| 919 | return y*yscale; |
---|
| 920 | } |
---|
| 921 | |
---|
| 922 | protected: |
---|
[1230] | 923 | c2_scaled_function_p<float_type>() : func() {} // hide default constructor, since its use is almost always an error. |
---|
| 924 | /// \brief the scaling factor for the function |
---|
| 925 | const c2_const_ptr<float_type> func; |
---|
[807] | 926 | float_type yscale; |
---|
| 927 | }; |
---|
| 928 | |
---|
| 929 | /// \brief A container into which any other c2_function can be dropped. |
---|
[1230] | 930 | /// \ingroup containers |
---|
[807] | 931 | /// It allows a function to be pre-evaluated at a point, and used at multiple places in an expression |
---|
| 932 | /// efficiently. If it is re-evaluated at the previous point, it returns the remembered values; |
---|
| 933 | /// otherwise, it re-evauates the function at the new point. |
---|
| 934 | /// |
---|
[1230] | 935 | /// The factory function c2_factory::cached_function() creates *new c2_cached_function_p |
---|
| 936 | template <typename float_type=double> class c2_cached_function_p : public c2_function<float_type> { |
---|
[807] | 937 | public: |
---|
| 938 | /// \brief construct the container |
---|
| 939 | /// |
---|
| 940 | /// \param f the function to be cached |
---|
[1230] | 941 | c2_cached_function_p(const c2_function<float_type> &f) : c2_function<float_type>(), |
---|
| 942 | func(f), init(false) {} |
---|
[807] | 943 | /// \copydoc c2_function::value_with_derivatives |
---|
| 944 | /// |
---|
| 945 | /// Checks to see if the function is being re-evaluated at the previous point, and |
---|
| 946 | /// returns remembered values if so. |
---|
| 947 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
| 948 | { |
---|
| 949 | if(!init || x != x0) { |
---|
| 950 | y=this->func->value_with_derivatives(x, &yp, &ypp); |
---|
| 951 | x0=x; |
---|
| 952 | init=true; |
---|
| 953 | } |
---|
| 954 | if(yprime) *yprime=yp; |
---|
| 955 | if(yprime2) *yprime2=ypp; |
---|
| 956 | return y; |
---|
| 957 | } |
---|
| 958 | |
---|
| 959 | protected: |
---|
[1230] | 960 | c2_cached_function_p() : func() {} // hide default constructor, since its use is almost always an error. |
---|
| 961 | const c2_const_ptr<float_type> func; |
---|
[807] | 962 | mutable bool init; |
---|
| 963 | mutable float_type x0, y, yp, ypp; |
---|
| 964 | |
---|
| 965 | }; |
---|
| 966 | |
---|
| 967 | /// \brief Provides function composition (nesting) |
---|
[1230] | 968 | /// \ingroup arithmetic_functions |
---|
[807] | 969 | /// This allows evaluation of \a f(g(x)) where \a f and \a g are c2_function objects. |
---|
| 970 | /// |
---|
[1230] | 971 | /// This should always be constructed using \ref compose_operator "c2_function::operator()" |
---|
| 972 | template <typename float_type=double> class c2_composed_function_p : public c2_binary_function<float_type> { |
---|
[807] | 973 | public: |
---|
| 974 | |
---|
| 975 | /// \brief construct \a outer( \a inner (x)) |
---|
| 976 | /// \note See c2_binary_function for discussion of ownership. |
---|
| 977 | /// \param outer the outer function |
---|
| 978 | /// \param inner the inner function |
---|
[1230] | 979 | c2_composed_function_p(const c2_function<float_type> &outer, const c2_function<float_type> &inner) : |
---|
| 980 | c2_binary_function<float_type>(combine, outer, inner) { this->set_domain(inner.xmin(), inner.xmax()); } |
---|
[807] | 981 | /// \brief Create a stub just for the combiner to avoid statics. |
---|
[1230] | 982 | c2_composed_function_p() : c2_binary_function<float_type>(combine) {} |
---|
[807] | 983 | |
---|
[1230] | 984 | /// \brief execute math necessary to do composition |
---|
| 985 | static float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, |
---|
| 986 | float_type x, float_type *yprime, float_type *yprime2) throw(c2_exception) |
---|
[807] | 987 | { |
---|
[1230] | 988 | float_type y0, y1; |
---|
[807] | 989 | if(yprime || yprime2) { |
---|
[1230] | 990 | float_type yp0, ypp0, yp1, ypp1; |
---|
| 991 | y0=right.value_with_derivatives(x, &yp0, &ypp0); |
---|
| 992 | y1=left.value_with_derivatives(y0, &yp1, &ypp1); |
---|
| 993 | if(yprime) *yprime=yp1*yp0; |
---|
| 994 | if(yprime2) *yprime2=ypp0*yp1+yp0*yp0*ypp1; |
---|
[807] | 995 | } else { |
---|
[1230] | 996 | y0=right(x); |
---|
| 997 | y1=left(y0); |
---|
[807] | 998 | } |
---|
| 999 | return y1; |
---|
| 1000 | } |
---|
| 1001 | }; |
---|
| 1002 | |
---|
| 1003 | /// \brief create a c2_function which is the sum of two other c2_function objects. |
---|
[1230] | 1004 | /// \ingroup arithmetic_functions |
---|
| 1005 | /// This should always be constructed using c2_function::operator+() |
---|
| 1006 | template <typename float_type=double> class c2_sum_p : public c2_binary_function<float_type> { |
---|
[807] | 1007 | public: |
---|
| 1008 | /// \brief construct \a left + \a right |
---|
[1230] | 1009 | /// \param left the left function |
---|
[807] | 1010 | /// \param right the right function |
---|
[1230] | 1011 | c2_sum_p(const c2_function<float_type> &left, const c2_function<float_type> &right) : c2_binary_function<float_type>(combine, left, right) {} |
---|
[807] | 1012 | /// \brief Create a stub just for the combiner to avoid statics. |
---|
[1230] | 1013 | c2_sum_p() : c2_binary_function<float_type>(combine) {} ; // create a stub just for the combiner to avoid statics |
---|
[807] | 1014 | |
---|
[1230] | 1015 | /// \brief execute math necessary to do addition |
---|
| 1016 | static float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, |
---|
| 1017 | float_type x, float_type *yprime, float_type *yprime2) throw(c2_exception) |
---|
[807] | 1018 | { |
---|
[1230] | 1019 | float_type y0, y1; |
---|
[807] | 1020 | if(yprime || yprime2) { |
---|
[1230] | 1021 | float_type yp0, ypp0, yp1, ypp1; |
---|
| 1022 | y0=left.value_with_derivatives(x, &yp0, &ypp0); |
---|
| 1023 | y1=right.value_with_derivatives(x, &yp1, &ypp1); |
---|
| 1024 | if(yprime) *yprime=yp0+yp1; |
---|
| 1025 | if(yprime2) *yprime2=ypp0+ypp1; |
---|
[807] | 1026 | } else { |
---|
[1230] | 1027 | y0=left(x); |
---|
| 1028 | y1=right(x); |
---|
[807] | 1029 | } |
---|
| 1030 | return y0+y1; |
---|
| 1031 | } |
---|
| 1032 | }; |
---|
| 1033 | |
---|
| 1034 | |
---|
| 1035 | /// \brief create a c2_function which is the difference of two other c2_functions. |
---|
[1230] | 1036 | /// \ingroup arithmetic_functions |
---|
| 1037 | /// This should always be constructed using c2_function::operator-() |
---|
| 1038 | template <typename float_type=double> class c2_diff_p : public c2_binary_function<float_type> { |
---|
[807] | 1039 | public: |
---|
| 1040 | /// \brief construct \a left - \a right |
---|
| 1041 | /// \param left the left function |
---|
| 1042 | /// \param right the right function |
---|
[1230] | 1043 | c2_diff_p(const c2_function<float_type> &left, const c2_function<float_type> &right) : c2_binary_function<float_type>(combine, left, right) {} |
---|
[807] | 1044 | /// \brief Create a stub just for the combiner to avoid statics. |
---|
[1230] | 1045 | c2_diff_p() : c2_binary_function<float_type>(combine) {} ; // create a stub just for the combiner to avoid statics |
---|
[807] | 1046 | |
---|
[1230] | 1047 | /// \brief execute math necessary to do subtraction |
---|
| 1048 | static float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, |
---|
| 1049 | float_type x, float_type *yprime, float_type *yprime2) throw(c2_exception) |
---|
[807] | 1050 | { |
---|
[1230] | 1051 | float_type y0, y1; |
---|
[807] | 1052 | if(yprime || yprime2) { |
---|
[1230] | 1053 | float_type yp0, ypp0, yp1, ypp1; |
---|
| 1054 | y0=left.value_with_derivatives(x, &yp0, &ypp0); |
---|
| 1055 | y1=right.value_with_derivatives(x, &yp1, &ypp1); |
---|
| 1056 | if(yprime) *yprime=yp0-yp1; |
---|
| 1057 | if(yprime2) *yprime2=ypp0-ypp1; |
---|
[807] | 1058 | } else { |
---|
[1230] | 1059 | y0=left(x); |
---|
| 1060 | y1=right(x); |
---|
[807] | 1061 | } |
---|
| 1062 | return y0-y1; |
---|
| 1063 | } |
---|
| 1064 | }; |
---|
| 1065 | |
---|
| 1066 | |
---|
| 1067 | /// \brief create a c2_function which is the product of two other c2_functions. |
---|
[1230] | 1068 | /// \ingroup arithmetic_functions |
---|
| 1069 | /// This should always be constructed using c2_function::operator*() |
---|
| 1070 | template <typename float_type=double> class c2_product_p : public c2_binary_function<float_type> { |
---|
[807] | 1071 | public: |
---|
| 1072 | /// \brief construct \a left * \a right |
---|
| 1073 | /// \param left the left function |
---|
| 1074 | /// \param right the right function |
---|
[1230] | 1075 | c2_product_p(const c2_function<float_type> &left, const c2_function<float_type> &right) : c2_binary_function<float_type>(combine, left, right) {} |
---|
[807] | 1076 | /// \brief Create a stub just for the combiner to avoid statics. |
---|
[1230] | 1077 | c2_product_p() : c2_binary_function<float_type>(combine) {} ; // create a stub just for the combiner to avoid statics |
---|
[807] | 1078 | |
---|
[1230] | 1079 | /// \brief execute math necessary to do multiplication |
---|
| 1080 | static float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, |
---|
| 1081 | float_type x, float_type *yprime, float_type *yprime2) throw(c2_exception) |
---|
[807] | 1082 | { |
---|
[1230] | 1083 | float_type y0, y1; |
---|
[807] | 1084 | if(yprime || yprime2) { |
---|
[1230] | 1085 | float_type yp0, ypp0, yp1, ypp1; |
---|
| 1086 | y0=left.value_with_derivatives(x, &yp0, &ypp0); |
---|
| 1087 | y1=right.value_with_derivatives(x, &yp1, &ypp1); |
---|
| 1088 | if(yprime) *yprime=y1*yp0+y0*yp1; |
---|
| 1089 | if(yprime2) *yprime2=ypp0*y1+2.0*yp0*yp1+ypp1*y0; |
---|
[807] | 1090 | } else { |
---|
[1230] | 1091 | y0=left(x); |
---|
| 1092 | y1=right(x); |
---|
[807] | 1093 | } |
---|
| 1094 | return y0*y1; |
---|
| 1095 | } |
---|
| 1096 | }; |
---|
| 1097 | |
---|
| 1098 | |
---|
| 1099 | /// \brief create a c2_function which is the ratio of two other c2_functions. |
---|
[1230] | 1100 | /// \ingroup arithmetic_functions |
---|
| 1101 | /// This should always be constructed using c2_function::operator/() |
---|
| 1102 | template <typename float_type=double> class c2_ratio_p : public c2_binary_function<float_type> { |
---|
[807] | 1103 | public: |
---|
| 1104 | /// \brief construct \a left / \a right |
---|
| 1105 | /// \param left the left function |
---|
| 1106 | /// \param right the right function |
---|
[1230] | 1107 | c2_ratio_p(const c2_function<float_type> &left, const c2_function<float_type> &right) : c2_binary_function<float_type>(combine, left, right) {} |
---|
[807] | 1108 | /// \brief Create a stub just for the combiner to avoid statics. |
---|
[1230] | 1109 | c2_ratio_p() : c2_binary_function<float_type>(combine) {} ; // create a stub just for the combiner to avoid statics |
---|
[807] | 1110 | |
---|
[1230] | 1111 | /// \brief execute math necessary to do division |
---|
| 1112 | static float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, |
---|
| 1113 | float_type x, float_type *yprime, float_type *yprime2) throw(c2_exception) |
---|
[807] | 1114 | { |
---|
[1230] | 1115 | float_type y0, y1; |
---|
[807] | 1116 | if(yprime || yprime2) { |
---|
[1230] | 1117 | float_type yp0, ypp0, yp1, ypp1; |
---|
| 1118 | y0=left.value_with_derivatives(x, &yp0, &ypp0); |
---|
| 1119 | y1=right.value_with_derivatives(x, &yp1, &ypp1); |
---|
| 1120 | if(yprime) *yprime=(yp0*y1-y0*yp1)/(y1*y1); // first deriv of ratio |
---|
| 1121 | if(yprime2) *yprime2=(y1*y1*ypp0+y0*(2*yp1*yp1-y1*ypp1)-2*y1*yp0*yp1)/(y1*y1*y1); |
---|
[807] | 1122 | } else { |
---|
[1230] | 1123 | y0=left(x); |
---|
| 1124 | y1=right(x); |
---|
[807] | 1125 | } |
---|
| 1126 | return y0/y1; |
---|
| 1127 | } |
---|
| 1128 | |
---|
| 1129 | }; |
---|
| 1130 | |
---|
[1230] | 1131 | /// \brief a c2_function which is constant |
---|
| 1132 | /// \ingroup parametric_functions |
---|
| 1133 | /// |
---|
| 1134 | /// The factory function c2_factory::constant() creates *new c2_constant_p() |
---|
| 1135 | template <typename float_type> class c2_constant_p : public c2_function<float_type> { |
---|
[807] | 1136 | public: |
---|
[1230] | 1137 | c2_constant_p(float_type x) : c2_function<float_type>(), value(x) {} |
---|
[807] | 1138 | void reset(float_type val) { value=val; } |
---|
| 1139 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
| 1140 | { if(yprime) *yprime=0; if(yprime2) *yprime2=0; return value; } |
---|
| 1141 | |
---|
| 1142 | private: |
---|
| 1143 | float_type value; |
---|
| 1144 | }; |
---|
| 1145 | |
---|
[1230] | 1146 | /// \brief a transformation of a coordinate, including an inverse |
---|
| 1147 | /// \ingroup transforms |
---|
| 1148 | template <typename float_type> class c2_transformation { |
---|
| 1149 | public: |
---|
| 1150 | /// \brief initialize all our function pointers |
---|
| 1151 | /// \param transformed true if this function is not the identity |
---|
| 1152 | /// \param xin input X transform |
---|
| 1153 | /// \param xinp input X transform derivative |
---|
| 1154 | /// \param xinpp input X transform second derivative |
---|
| 1155 | /// \param xout output X transform, which MUST be the inverse of \a xin |
---|
| 1156 | c2_transformation(bool transformed, |
---|
| 1157 | float_type (*xin)(float_type), float_type (*xinp)(float_type), float_type (*xinpp)(float_type), float_type (*xout)(float_type) |
---|
| 1158 | ) : |
---|
| 1159 | fTransformed(transformed), fHasStaticTransforms(true), |
---|
| 1160 | pIn(xin), pInPrime(xinp), pInDPrime(xinpp), pOut(xout) { } |
---|
| 1161 | |
---|
| 1162 | /// \brief initialize all our function pointers so that only the (overridden) virtual functions can be called without an error |
---|
| 1163 | /// \param transformed true if this function is nonlinear |
---|
| 1164 | c2_transformation(bool transformed) : |
---|
| 1165 | fTransformed(transformed), fHasStaticTransforms(false), |
---|
| 1166 | pIn(report_error), pInPrime(report_error), pInDPrime(report_error), pOut(report_error) { } |
---|
| 1167 | /// \brief the destructor |
---|
| 1168 | virtual ~c2_transformation() { } |
---|
| 1169 | /// \brief flag to indicate if this transform is not the identity |
---|
| 1170 | const bool fTransformed; |
---|
| 1171 | /// \brief flag to indicate if the static function pointers can be used for efficiency |
---|
| 1172 | const bool fHasStaticTransforms; |
---|
| 1173 | |
---|
| 1174 | /// \note the pointers to functions allow highly optimized access when static functions are available. |
---|
| 1175 | /// They are only used inside value_with_derivatives(), which is assumed to be the most critical routine. |
---|
| 1176 | /// \brief non-virtual pointer to input X transform |
---|
| 1177 | float_type (* const pIn)(float_type); |
---|
| 1178 | /// \brief non-virtual pointer to input X transform derivative |
---|
| 1179 | float_type (* const pInPrime)(float_type); |
---|
| 1180 | /// \brief non-virtual pointer to input X transform second derivative |
---|
| 1181 | float_type (* const pInDPrime)(float_type); |
---|
| 1182 | /// \brief non-virtual pointer to output X transform |
---|
| 1183 | float_type (* const pOut)(float_type); |
---|
| 1184 | |
---|
| 1185 | /// \brief virtual input X transform |
---|
| 1186 | virtual float_type fIn(float_type x) const { return pIn(x); } |
---|
| 1187 | /// \brief virtual input X transform derivative |
---|
| 1188 | virtual float_type fInPrime(float_type x) const { return pInPrime(x); } |
---|
| 1189 | /// \brief virtual input X transform second derivative |
---|
| 1190 | virtual float_type fInDPrime(float_type x) const { return pInDPrime(x); } |
---|
| 1191 | /// \brief virtual output X transform |
---|
| 1192 | virtual float_type fOut(float_type x) const { return pOut(x); } |
---|
| 1193 | |
---|
| 1194 | protected: |
---|
| 1195 | /// \brief utility function for unimplemented conversion |
---|
| 1196 | static float_type report_error(float_type x) { throw c2_exception("use of improperly constructed axis transform"); return x; } |
---|
| 1197 | /// \brief utility function f(x)=x useful in axis transforms |
---|
| 1198 | static float_type ident(float_type x) { return x; } |
---|
| 1199 | /// \brief utility function f(x)=1 useful in axis transforms |
---|
| 1200 | static float_type one(float_type) { return 1; } |
---|
| 1201 | /// \brief utility function f(x)=0 useful in axis transforms |
---|
| 1202 | static float_type zero(float_type) { return 0; } |
---|
| 1203 | /// \brief utility function f(x)=1/x useful in axis transforms |
---|
| 1204 | static float_type recip(float_type x) { return 1.0/x; } |
---|
| 1205 | /// \brief utility function f(x)=-1/x**2 useful in axis transforms |
---|
| 1206 | static float_type recip_prime(float_type x) { return -1/(x*x); } |
---|
| 1207 | /// \brief utility function f(x)=2/x**3 useful in axis transforms |
---|
| 1208 | static float_type recip_prime2(float_type x) { return 2/(x*x*x); } |
---|
| 1209 | |
---|
| 1210 | }; |
---|
| 1211 | |
---|
| 1212 | /// \brief the identity transform |
---|
| 1213 | /// \ingroup transforms |
---|
| 1214 | template <typename float_type> class c2_transformation_linear : public c2_transformation<float_type> { |
---|
| 1215 | public: |
---|
| 1216 | /// \brief constructor |
---|
| 1217 | c2_transformation_linear() : c2_transformation<float_type>(false, this->ident, this->one, this->zero, this->ident) { } |
---|
| 1218 | /// \brief destructor |
---|
| 1219 | ~c2_transformation_linear() { } |
---|
| 1220 | }; |
---|
| 1221 | /// \brief log axis transform |
---|
| 1222 | /// \ingroup transforms |
---|
| 1223 | template <typename float_type> class c2_transformation_log : public c2_transformation<float_type> { |
---|
| 1224 | public: |
---|
| 1225 | /// \brief constructor |
---|
| 1226 | c2_transformation_log() : c2_transformation<float_type>(true, std::log, this->recip, this->recip_prime, std::exp) { } |
---|
| 1227 | /// \brief destructor |
---|
| 1228 | ~c2_transformation_log() { } |
---|
| 1229 | }; |
---|
| 1230 | /// \brief reciprocal axis transform |
---|
| 1231 | /// \ingroup transforms |
---|
| 1232 | template <typename float_type> class c2_transformation_recip : public c2_transformation<float_type> { |
---|
| 1233 | public: |
---|
| 1234 | /// \brief constructor |
---|
| 1235 | c2_transformation_recip() : c2_transformation<float_type>(true, this->recip, this->recip_prime, this->recip_prime2, this->recip) { } |
---|
| 1236 | /// \brief destructor |
---|
| 1237 | ~c2_transformation_recip() { } |
---|
| 1238 | }; |
---|
| 1239 | |
---|
| 1240 | /// \brief a transformation of a function in and out of a coordinate space, using 2 c2_transformations |
---|
| 1241 | /// |
---|
| 1242 | /// This class is a container for two axis transforms, but also provides the critical evaluate() |
---|
| 1243 | /// function which converts a result in internal coordinates (with derivatives) into the external representation |
---|
| 1244 | /// \ingroup transforms |
---|
| 1245 | template <typename float_type> |
---|
| 1246 | class c2_function_transformation { |
---|
| 1247 | public: |
---|
| 1248 | /// \brief construct this from two c2_transformation instances |
---|
| 1249 | /// \param xx the X axis transform |
---|
| 1250 | /// \param yy the Y axis transform |
---|
| 1251 | c2_function_transformation( |
---|
| 1252 | const c2_transformation<float_type> &xx, const c2_transformation<float_type> &yy) : |
---|
| 1253 | isIdentity(!(xx.fTransformed || yy.fTransformed)), X(xx), Y(yy) { } |
---|
| 1254 | /// \brief destructor |
---|
| 1255 | virtual ~c2_function_transformation() { delete &X; delete &Y; } |
---|
| 1256 | /// \brief evaluate the transformation from internal coordinates to external coordinates |
---|
| 1257 | /// \param xraw the value of \a x in external cordinates at which the transform is taking place |
---|
| 1258 | /// \param y the value of the function in internal coordinates |
---|
| 1259 | /// \param yp0 the derivative in internal coordinates |
---|
| 1260 | /// \param ypp0 the second derivative in internal coordinates |
---|
| 1261 | /// \param [out] yprime pointer to the derivative, or NULL, in external coordinates |
---|
| 1262 | /// \param [out] yprime2 pointer to the second derivative, or NULL, in external coordinates |
---|
| 1263 | /// \return the value of the function in external coordinates |
---|
| 1264 | virtual float_type evaluate(float_type xraw, |
---|
| 1265 | float_type y, float_type yp0, float_type ypp0, |
---|
| 1266 | float_type *yprime, float_type *yprime2) const; |
---|
| 1267 | /// \brief flag indicating of the transform is the identity, and can be skipped for efficiency |
---|
| 1268 | const bool isIdentity; |
---|
| 1269 | /// \brief the X axis transform |
---|
| 1270 | const c2_transformation<float_type> &X; |
---|
| 1271 | /// \brief the Y axis transform |
---|
| 1272 | const c2_transformation<float_type> &Y; |
---|
| 1273 | }; |
---|
| 1274 | |
---|
| 1275 | /// \brief a transformation of a function in and out of lin-lin space |
---|
| 1276 | /// |
---|
| 1277 | /// \ingroup transforms |
---|
| 1278 | template <typename float_type> class c2_lin_lin_function_transformation : |
---|
| 1279 | public c2_function_transformation<float_type> { |
---|
| 1280 | public: |
---|
| 1281 | c2_lin_lin_function_transformation() : |
---|
| 1282 | c2_function_transformation<float_type>( |
---|
| 1283 | *new c2_transformation_linear<float_type>, |
---|
| 1284 | *new c2_transformation_linear<float_type> |
---|
| 1285 | ) { } |
---|
| 1286 | virtual ~c2_lin_lin_function_transformation() { } |
---|
| 1287 | }; |
---|
| 1288 | |
---|
| 1289 | /// \brief a transformation of a function in and out of log-log space |
---|
| 1290 | /// |
---|
| 1291 | /// \ingroup transforms |
---|
| 1292 | template <typename float_type> class c2_log_log_function_transformation : |
---|
| 1293 | public c2_function_transformation<float_type> { |
---|
| 1294 | public: |
---|
| 1295 | c2_log_log_function_transformation() : |
---|
| 1296 | c2_function_transformation<float_type>( |
---|
| 1297 | *new c2_transformation_log<float_type>, |
---|
| 1298 | *new c2_transformation_log<float_type> |
---|
| 1299 | ) { } |
---|
| 1300 | virtual ~c2_log_log_function_transformation() { } |
---|
| 1301 | }; |
---|
| 1302 | |
---|
| 1303 | /// \brief a transformation of a function in and out of lin-log space |
---|
| 1304 | /// |
---|
| 1305 | /// \ingroup transforms |
---|
| 1306 | template <typename float_type> class c2_lin_log_function_transformation : |
---|
| 1307 | public c2_function_transformation<float_type> { |
---|
| 1308 | public: |
---|
| 1309 | c2_lin_log_function_transformation() : |
---|
| 1310 | c2_function_transformation<float_type>( |
---|
| 1311 | *new c2_transformation_linear<float_type>, |
---|
| 1312 | *new c2_transformation_log<float_type> |
---|
| 1313 | ) { } |
---|
| 1314 | virtual ~c2_lin_log_function_transformation() { } |
---|
| 1315 | }; |
---|
| 1316 | |
---|
| 1317 | /// \brief a transformation of a function in and out of log-lin space |
---|
| 1318 | /// |
---|
| 1319 | /// \ingroup transforms |
---|
| 1320 | template <typename float_type> class c2_log_lin_function_transformation : |
---|
| 1321 | public c2_function_transformation<float_type> { |
---|
| 1322 | public: |
---|
| 1323 | c2_log_lin_function_transformation() : |
---|
| 1324 | c2_function_transformation<float_type>( |
---|
| 1325 | *new c2_transformation_log<float_type>, |
---|
| 1326 | *new c2_transformation_linear<float_type> |
---|
| 1327 | ) { } |
---|
| 1328 | virtual ~c2_log_lin_function_transformation() { } |
---|
| 1329 | }; |
---|
| 1330 | |
---|
| 1331 | /// \brief a transformation of a function in and out of Arrhenuis (1/x vs. log(y)) space |
---|
| 1332 | /// |
---|
| 1333 | /// \ingroup transforms |
---|
| 1334 | template <typename float_type> class c2_arrhenius_function_transformation : |
---|
| 1335 | public c2_function_transformation<float_type> { |
---|
| 1336 | public: |
---|
| 1337 | c2_arrhenius_function_transformation() : |
---|
| 1338 | c2_function_transformation<float_type>( |
---|
| 1339 | *new c2_transformation_recip<float_type>, |
---|
| 1340 | *new c2_transformation_log<float_type> |
---|
| 1341 | ) { } |
---|
| 1342 | virtual ~c2_arrhenius_function_transformation() { } |
---|
| 1343 | }; |
---|
| 1344 | |
---|
[807] | 1345 | /** |
---|
| 1346 | \brief create a cubic spline interpolation of a set of (x,y) pairs |
---|
[1230] | 1347 | \ingroup interpolators |
---|
[807] | 1348 | This is one of the main reasons for c2_function objects to exist. |
---|
| 1349 | |
---|
| 1350 | It provides support for cubic spline interpolation of data provides from tables of \a x, \a y pairs. |
---|
| 1351 | It supports automatic, transparent linearization of the data before storing in its tables (through |
---|
| 1352 | subclasses such as |
---|
| 1353 | log_lin_interpolating_function, lin_log_interpolating_function, and |
---|
| 1354 | log_log_interpolating_function) to permit very high accuracy representations of data which have a suitable |
---|
| 1355 | structure. It provides utility functions LinearInterpolatingGrid() and LogLogInterpolatingGrid() |
---|
| 1356 | to create grids for mapping other functions onto a arithmetic or geometric grid. |
---|
| 1357 | |
---|
| 1358 | In its simplest form, an untransformed cubic spline of a data set, using natural boundary conditions |
---|
| 1359 | (vanishing second derivative), is created as: \n |
---|
| 1360 | \code |
---|
[1230] | 1361 | c2_ptr<double> c2p; |
---|
| 1362 | c2_factory<double> c2; |
---|
[807] | 1363 | std::vector<double> xvals(10), yvals(10); |
---|
| 1364 | // < fill in xvals and yvals > |
---|
[1230] | 1365 | c2p myfunc=c2.interpolating_function().load(xvals, yvals,true,0,true,0); |
---|
[807] | 1366 | // and it can be evaluated at a point for its value only by: |
---|
| 1367 | double y=myfunc(x); |
---|
| 1368 | // or it can be evaluated with its derivatives by |
---|
| 1369 | double yprime, yprime2; |
---|
| 1370 | double y=myfunc(x,&yprime, &yprime2); |
---|
| 1371 | \endcode |
---|
[1230] | 1372 | |
---|
| 1373 | The factory function c2_factory::interpolating_function() creates *new interpolating_function_p() |
---|
[807] | 1374 | */ |
---|
| 1375 | |
---|
[1230] | 1376 | template <typename float_type=double> class interpolating_function_p : public c2_function<float_type> { |
---|
[807] | 1377 | public: |
---|
[1230] | 1378 | /// \brief an empty linear-linear cubic-spline interpolating_function_p |
---|
[807] | 1379 | /// |
---|
| 1380 | /// lots to say here, but see Numerical Recipes for a discussion of cubic splines. |
---|
[1230] | 1381 | /// |
---|
| 1382 | interpolating_function_p() : c2_function<float_type>(), |
---|
| 1383 | fTransform(*new c2_lin_lin_function_transformation<float_type>) { } |
---|
| 1384 | |
---|
| 1385 | /// \brief an empty cubic-spline interpolating_function_p with a specific transform |
---|
| 1386 | /// |
---|
| 1387 | interpolating_function_p(const c2_function_transformation<float_type> &transform) : c2_function<float_type>(), |
---|
| 1388 | fTransform(transform) { } |
---|
| 1389 | |
---|
| 1390 | /// \brief do the dirty work of constructing the spline from a function. |
---|
[807] | 1391 | /// \param x the list of abscissas. Must be either strictly increasing or strictly decreasing. |
---|
| 1392 | /// Strictly increasing is preferred, as less memory is used since a copy is not required for the sampling grid. |
---|
| 1393 | /// \param f the list of function values. |
---|
| 1394 | /// \param lowerSlopeNatural if true, set y''(first point)=0, otherwise compute it from \a lowerSope |
---|
| 1395 | /// \param lowerSlope derivative of the function at the lower bound, used only if \a lowerSlopeNatural is false |
---|
| 1396 | /// \param upperSlopeNatural if true, set y''(last point)=0, otherwise compute it from \a upperSope |
---|
| 1397 | /// \param upperSlope derivative of the function at the upper bound, used only if \a upperSlopeNatural is false |
---|
[1230] | 1398 | /// \param splined if true (default), use cubic spline, if false, use linear interpolation. |
---|
| 1399 | /// \return the same interpolating function, filled |
---|
| 1400 | interpolating_function_p<float_type> & load(const std::vector<float_type> &x, const std::vector<float_type> &f, |
---|
| 1401 | bool lowerSlopeNatural, float_type lowerSlope, |
---|
| 1402 | bool upperSlopeNatural, float_type upperSlope, bool splined=true |
---|
| 1403 | ) throw(c2_exception); |
---|
| 1404 | |
---|
| 1405 | /// \brief do the dirty work of constructing the spline from a function. |
---|
| 1406 | /// \param data std::vector of std::pairs of x,y. Will be sorted into x increasing order in place. |
---|
| 1407 | /// \param lowerSlopeNatural if true, set y''(first point)=0, otherwise compute it from \a lowerSope |
---|
| 1408 | /// \param lowerSlope derivative of the function at the lower bound, used only if \a lowerSlopeNatural is false |
---|
| 1409 | /// \param upperSlopeNatural if true, set y''(last point)=0, otherwise compute it from \a upperSope |
---|
| 1410 | /// \param upperSlope derivative of the function at the upper bound, used only if \a upperSlopeNatural is false |
---|
| 1411 | /// \param splined if true (default), use cubic spline, if false, use linear interpolation. |
---|
| 1412 | /// \return the same interpolating function, filled |
---|
| 1413 | interpolating_function_p<float_type> & load_pairs( |
---|
| 1414 | std::vector<std::pair<float_type, float_type> > &data, |
---|
| 1415 | bool lowerSlopeNatural, float_type lowerSlope, |
---|
| 1416 | bool upperSlopeNatural, float_type upperSlope, bool splined=true |
---|
| 1417 | ) throw(c2_exception); |
---|
| 1418 | |
---|
| 1419 | /// \brief do the dirty work of constructing the spline from a function. |
---|
| 1420 | /// \param func a function without any requirement of valid derivatives to sample into an interpolating function. |
---|
| 1421 | /// Very probably a c2_classic_function. |
---|
| 1422 | /// \param xmin the lower bound of the region to sample |
---|
| 1423 | /// \param xmax the upper bound of the region to sample |
---|
| 1424 | /// \param abs_tol the maximum absolute error permitted when linearly interpolating the points. |
---|
| 1425 | /// the real error will be much smaller, since this uses cubic splines at the end. |
---|
| 1426 | /// \param rel_tol the maximum relative error permitted when linearly interpolating the points. |
---|
| 1427 | /// the real error will be much smaller, since this uses cubic splines at the end. |
---|
| 1428 | /// \param lowerSlopeNatural if true, set y'(first point) from 3-point parabola, otherwise compute it from \a lowerSope |
---|
| 1429 | /// \param lowerSlope derivative of the function at the lower bound, used only if \a lowerSlopeNatural is false |
---|
| 1430 | /// \param upperSlopeNatural if true, set y'(last point) from 3-point parabola, otherwise compute it from \a upperSope |
---|
| 1431 | /// \param upperSlope derivative of the function at the upper bound, used only if \a upperSlopeNatural is false |
---|
| 1432 | /// \return the same interpolating function, filled |
---|
| 1433 | /// \note If the interpolator being filled has a log vertical axis, put the desired relative error in |
---|
| 1434 | /// \a abs_tol, and 0 in \a rel_tol since the absolute error on the log of a function is the relative error |
---|
| 1435 | /// on the function itself. |
---|
| 1436 | interpolating_function_p<float_type> & sample_function(const c2_function<float_type> &func, |
---|
| 1437 | float_type xmin, float_type xmax, float_type abs_tol, float_type rel_tol, |
---|
| 1438 | bool lowerSlopeNatural, float_type lowerSlope, |
---|
| 1439 | bool upperSlopeNatural, float_type upperSlope |
---|
| 1440 | ) throw(c2_exception); |
---|
[807] | 1441 | |
---|
[1230] | 1442 | |
---|
| 1443 | /// \brief initialize from a grid of points and a c2_function (un-normalized) to an |
---|
| 1444 | /// interpolator which, when evaluated with a uniform random variate on [0,1] returns random numbers |
---|
| 1445 | /// distributed as the input function. |
---|
| 1446 | /// \see \ref random_subsec "Arbitrary random generation" |
---|
| 1447 | /// inverse_integrated_density starts with a probability density std::vector, generates the integral, |
---|
| 1448 | /// and generates an interpolating_function_p of the inverse function which, when evaluated using a uniform random on [0,1] returns values |
---|
| 1449 | /// with a density distribution equal to the input distribution |
---|
| 1450 | /// If the data are passed in reverse order (large X first), the integral is carried out from the big end. |
---|
| 1451 | /// \param bincenters the positions at which to sample the function \a binheights |
---|
| 1452 | /// \param binheights a function which describes the density of the random number distribution to be produced. |
---|
| 1453 | /// \return an initialized interpolator, which |
---|
| 1454 | /// if evaluated randomly with a uniform variate on [0,1] produces numbers |
---|
| 1455 | /// distributed according to \a binheights |
---|
| 1456 | interpolating_function_p<float_type> & load_random_generator_function( |
---|
| 1457 | const std::vector<float_type> &bincenters, const c2_function<float_type> &binheights) |
---|
| 1458 | throw(c2_exception); |
---|
| 1459 | |
---|
| 1460 | /// \brief initialize from a grid of points and an std::vector of probability densities (un-normalized) to an |
---|
| 1461 | /// interpolator which, when evaluated with a uniform random variate on [0,1] returns random numbers |
---|
| 1462 | /// distributed as the input histogram. |
---|
| 1463 | /// \see \ref random_subsec "Arbitrary random generation" |
---|
| 1464 | /// inverse_integrated_density starts with a probability density std::vector, generates the integral, |
---|
| 1465 | /// and generates an interpolating_function_p of the inverse function which, when evaluated using a uniform random on [0,1] returns values |
---|
| 1466 | /// with a density distribution equal to the input distribution |
---|
| 1467 | /// If the data are passed in reverse order (large X first), the integral is carried out from the big end. |
---|
| 1468 | /// \param bins if \a bins .size()==\a binheights .size(), the centers of the bins. \n |
---|
| 1469 | /// if \a bins .size()==\a binheights .size()+1, the edges of the bins |
---|
| 1470 | /// \param binheights a vector which describes the density of the random number distribution to be produced. |
---|
| 1471 | /// Note density... the numbers in the bins are not counts, but counts/unit bin width. |
---|
| 1472 | /// \return an initialized interpolator, which |
---|
| 1473 | /// if evaluated randomly with a uniform variate on [0,1] produces numbers |
---|
| 1474 | /// distributed according to \a binheights |
---|
| 1475 | interpolating_function_p<float_type> & load_random_generator_bins( |
---|
| 1476 | const std::vector<float_type> &bins, const std::vector<float_type> &binheights) |
---|
| 1477 | throw(c2_exception); |
---|
[807] | 1478 | |
---|
[1230] | 1479 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception); |
---|
[807] | 1480 | |
---|
[1230] | 1481 | /// \brief destructor |
---|
| 1482 | virtual ~interpolating_function_p() { delete &fTransform; } |
---|
| 1483 | |
---|
| 1484 | /// \brief create a new, empty interpolating function of this type (virtual constructor) |
---|
| 1485 | virtual interpolating_function_p<float_type> &clone() const throw(c2_exception) |
---|
| 1486 | { return *new interpolating_function_p<float_type>(); } |
---|
| 1487 | |
---|
[807] | 1488 | /// \brief retrieve copies of the x & y tables from which this was built |
---|
| 1489 | /// |
---|
| 1490 | /// This is often useful in the creation of new interpolating functions with transformed data. |
---|
[1230] | 1491 | /// The vectors will have their sizes set correctly on return. |
---|
[807] | 1492 | /// \param [in, out] xvals the abscissas |
---|
| 1493 | /// \param [in, out] yvals the ordinates |
---|
| 1494 | void get_data(std::vector<float_type> &xvals, std::vector<float_type> &yvals) const throw() ; |
---|
| 1495 | |
---|
| 1496 | /// \brief enable extrapolation of the function below the tabulated data. |
---|
| 1497 | /// |
---|
| 1498 | /// This allows the interpolator to be extrapolated outside the bounds of the data, |
---|
| 1499 | /// using whatever derivatives it already had at the lower bound. |
---|
| 1500 | /// \param bound the abscissa to which the function should be extended. |
---|
| 1501 | void set_lower_extrapolation(float_type bound); |
---|
| 1502 | /// \brief enable extrapolation of the function above the tabulated data. |
---|
| 1503 | /// |
---|
| 1504 | /// This allows the interpolator to be extrapolated outside the bounds of the data, |
---|
| 1505 | /// using whatever derivatives it already had at the upper bound. |
---|
| 1506 | /// \param bound the abscissa to which the function should be extended. |
---|
| 1507 | void set_upper_extrapolation(float_type bound); |
---|
| 1508 | |
---|
| 1509 | // these functions correctly combine the interpolating function with another interpolating function |
---|
| 1510 | // preserving the X bounds and mapping functions of the host (left hand) function. |
---|
| 1511 | |
---|
[1230] | 1512 | /// \brief create a new interpolating_function_p which is the \a source |
---|
[807] | 1513 | /// function applied to every point in the interpolating tables |
---|
| 1514 | /// |
---|
| 1515 | /// This carefully manages the derivative of the composed function at the two ends. |
---|
| 1516 | /// \param source the function to apply |
---|
[1230] | 1517 | /// \return a new interpolating_function_p with the same mappings for x and y |
---|
| 1518 | interpolating_function_p <float_type> & unary_operator(const c2_function<float_type> &source) const; |
---|
[807] | 1519 | |
---|
[1230] | 1520 | /// \brief create a new interpolating_function_p which is the parent interpolating_function_p |
---|
[807] | 1521 | /// combined with \a rhs using \a combiner at every point in the interpolating tables |
---|
| 1522 | /// |
---|
| 1523 | /// This carefully manages the derivative of the composed function at the two ends. |
---|
| 1524 | /// \param rhs the function to apply |
---|
| 1525 | /// \param combining_stub a function which defines which binary operation to use. |
---|
[1230] | 1526 | /// \return a new interpolating_function_p with the same mappings for x and y |
---|
| 1527 | interpolating_function_p <float_type> & binary_operator(const c2_function<float_type> &rhs, |
---|
| 1528 | const c2_binary_function<float_type> *combining_stub |
---|
[807] | 1529 | ) const; |
---|
[1230] | 1530 | /// \brief produce a newly resampled interpolating_function_p which is the specified sum. |
---|
[807] | 1531 | /// \param rhs the function to add, pointwise |
---|
[1230] | 1532 | /// \return a new interpolating_function_p |
---|
| 1533 | interpolating_function_p <float_type> & add_pointwise (const c2_function<float_type> &rhs) const { |
---|
| 1534 | return binary_operator(rhs, new c2_sum_p<float_type>()); } |
---|
| 1535 | /// \brief produce a newly resampled interpolating_function_p which is the specified difference. |
---|
[807] | 1536 | /// \param rhs the function to subtract, pointwise |
---|
[1230] | 1537 | /// \return a new interpolating_function_p |
---|
| 1538 | interpolating_function_p <float_type> & subtract_pointwise (const c2_function<float_type> &rhs) const { |
---|
| 1539 | return binary_operator(rhs, new c2_diff_p<float_type>()); } |
---|
| 1540 | /// \brief produce a newly resampled interpolating_function_p which is the specified product. |
---|
[807] | 1541 | /// \param rhs the function to multiply, pointwise |
---|
[1230] | 1542 | /// \return a new interpolating_function_p |
---|
| 1543 | interpolating_function_p <float_type> & multiply_pointwise (const c2_function<float_type> &rhs) const { |
---|
| 1544 | return binary_operator(rhs, new c2_product_p<float_type>()); } |
---|
| 1545 | /// \brief produce a newly resampled interpolating_function_p which is the specified ratio. |
---|
[807] | 1546 | /// \param rhs the function to divide, pointwise |
---|
[1230] | 1547 | /// \return a new interpolating_function_p |
---|
| 1548 | interpolating_function_p <float_type> & divide_pointwise (const c2_function<float_type> &rhs) const { |
---|
| 1549 | return binary_operator(rhs, new c2_ratio_p<float_type>()); } |
---|
| 1550 | /// \brief copy data from another interpolating function. This only makes sense if the source |
---|
| 1551 | /// function has the same transforms as the destination. |
---|
| 1552 | /// \param rhs interpolating_function_p to copy from |
---|
| 1553 | void clone_data(const interpolating_function_p <float_type> &rhs) { |
---|
| 1554 | Xraw=rhs.Xraw; X=rhs.X; F=rhs.F; y2=rhs.y2; |
---|
| 1555 | set_sampling_grid_pointer(Xraw); |
---|
| 1556 | } |
---|
| 1557 | |
---|
| 1558 | const c2_function_transformation<float_type> &fTransform; |
---|
[807] | 1559 | |
---|
[1230] | 1560 | protected: |
---|
| 1561 | /// \brief create the spline coefficients |
---|
| 1562 | void spline( |
---|
[807] | 1563 | bool lowerSlopeNatural, float_type lowerSlope, |
---|
[1230] | 1564 | bool upperSlopeNatural, float_type upperSlope |
---|
| 1565 | ) throw(c2_exception); |
---|
[807] | 1566 | |
---|
[1230] | 1567 | // This is for sorting the data. It must be static if it's going to be a class member. |
---|
| 1568 | static bool comp_pair(std::pair<float_type,float_type> const &i, std::pair<float_type,float_type> const &j) {return i.first<j.first;} |
---|
| 1569 | |
---|
[807] | 1570 | std::vector<float_type> Xraw, X, F, y2; |
---|
[1230] | 1571 | c2_const_ptr<float_type> sampler_function; |
---|
| 1572 | bool xInverted; |
---|
| 1573 | mutable size_t lastKLow; |
---|
[807] | 1574 | }; |
---|
| 1575 | |
---|
[1230] | 1576 | /// \brief A spline with X transformed into log space. |
---|
| 1577 | /// \ingroup interpolators |
---|
[807] | 1578 | /// Most useful for functions looking like y=log(x) or any other function with a huge X dynamic range, |
---|
| 1579 | /// and a slowly varying Y. |
---|
[1230] | 1580 | /// |
---|
| 1581 | /// The factory function c2_factory::log_lin_interpolating_function() creates *new log_lin_interpolating_function_p() |
---|
| 1582 | template <typename float_type=double> class log_lin_interpolating_function_p : public interpolating_function_p <float_type> { |
---|
[807] | 1583 | public: |
---|
[1230] | 1584 | /// \brief an empty log-linear cubic-spline interpolating_function_p |
---|
| 1585 | /// |
---|
| 1586 | log_lin_interpolating_function_p() : interpolating_function_p<float_type>(*new c2_log_lin_function_transformation<float_type>) |
---|
| 1587 | { } |
---|
| 1588 | /// \brief create a new, empty interpolating function of this type (virtual constructor) |
---|
| 1589 | virtual interpolating_function_p<float_type> &clone() const throw(c2_exception) |
---|
| 1590 | { return *new log_lin_interpolating_function_p<float_type>(); } |
---|
[807] | 1591 | }; |
---|
| 1592 | |
---|
| 1593 | |
---|
[1230] | 1594 | /// \brief A spline with Y transformed into log space. |
---|
| 1595 | /// \ingroup interpolators |
---|
| 1596 | /// Most useful for functions looking like y=exp(x) |
---|
[807] | 1597 | /// |
---|
[1230] | 1598 | /// The factory function c2_factory::lin_log_interpolating_function() creates *new lin_log_interpolating_function_p() |
---|
| 1599 | template <typename float_type=double> class lin_log_interpolating_function_p : public interpolating_function_p <float_type> { |
---|
[807] | 1600 | public: |
---|
[1230] | 1601 | /// \brief an empty linear-log cubic-spline interpolating_function_p |
---|
| 1602 | /// |
---|
| 1603 | lin_log_interpolating_function_p() : interpolating_function_p<float_type>(*new c2_lin_log_function_transformation<float_type>) |
---|
| 1604 | { } |
---|
| 1605 | /// \brief create a new, empty interpolating function of this type (virtual constructor) |
---|
| 1606 | virtual interpolating_function_p<float_type> &clone() const throw(c2_exception) |
---|
| 1607 | { return *new lin_log_interpolating_function_p<float_type>(); } |
---|
[807] | 1608 | }; |
---|
| 1609 | |
---|
| 1610 | |
---|
[1230] | 1611 | /// \brief A spline with X and Y transformed into log space. |
---|
| 1612 | /// \ingroup interpolators |
---|
| 1613 | /// Most useful for functions looking like y=x^n or any other function with a huge X and Y dynamic range. |
---|
[807] | 1614 | /// |
---|
[1230] | 1615 | /// The factory function c2_factory::log_log_interpolating_function() creates *new log_log_interpolating_function_p() |
---|
| 1616 | template <typename float_type=double> class log_log_interpolating_function_p : public interpolating_function_p <float_type> { |
---|
| 1617 | public: |
---|
| 1618 | /// \brief an empty log-log cubic-spline interpolating_function_p |
---|
| 1619 | /// |
---|
| 1620 | log_log_interpolating_function_p() : interpolating_function_p<float_type>(*new c2_log_log_function_transformation<float_type>) |
---|
| 1621 | { } |
---|
| 1622 | /// \brief create a new, empty interpolating function of this type (virtual constructor) |
---|
| 1623 | virtual interpolating_function_p<float_type> &clone() const throw(c2_exception) |
---|
| 1624 | { return *new log_log_interpolating_function_p<float_type>(); } |
---|
[807] | 1625 | }; |
---|
| 1626 | |
---|
| 1627 | |
---|
[1230] | 1628 | /// \brief A spline with X in reciprocal space and Y transformed in log space. |
---|
| 1629 | /// \ingroup interpolators |
---|
[807] | 1630 | /// Most useful for thermodynamic types of data where Y is roughly A*exp(-B/x). |
---|
| 1631 | /// Typical examples are reaction rate data, and thermistor calibration data. |
---|
[1230] | 1632 | /// |
---|
| 1633 | /// The factory function c2_factory::arrhenius_interpolating_function() creates *new arrhenius_interpolating_function_p() |
---|
| 1634 | template <typename float_type=double> class arrhenius_interpolating_function_p : public interpolating_function_p <float_type> { |
---|
[807] | 1635 | public: |
---|
[1230] | 1636 | /// \brief an empty arrhenius cubic-spline interpolating_function_p |
---|
| 1637 | /// |
---|
| 1638 | arrhenius_interpolating_function_p() : interpolating_function_p<float_type>(*new c2_arrhenius_function_transformation<float_type>) |
---|
| 1639 | { } |
---|
| 1640 | /// \brief create a new, empty interpolating function of this type (virtual constructor) |
---|
| 1641 | virtual interpolating_function_p<float_type> &clone() const throw(c2_exception) |
---|
| 1642 | { return *new arrhenius_interpolating_function_p<float_type>(); } |
---|
[807] | 1643 | }; |
---|
| 1644 | |
---|
| 1645 | /// \brief compute sin(x) with its derivatives. |
---|
[1230] | 1646 | /// \ingroup math_functions |
---|
[807] | 1647 | /// |
---|
[1230] | 1648 | /// The factory function c2_factory::sin() creates *new c2_sin_p |
---|
| 1649 | template <typename float_type=double> class c2_sin_p : public c2_function<float_type> { |
---|
[807] | 1650 | public: |
---|
[1230] | 1651 | /// \brief constructor. |
---|
| 1652 | c2_sin_p() : c2_function<float_type>() {} |
---|
[807] | 1653 | |
---|
| 1654 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
| 1655 | { float_type q=std::sin(x); if(yprime) *yprime=std::cos(x); if(yprime2) *yprime2=-q; return q; } |
---|
| 1656 | |
---|
| 1657 | /// \brief return a grid dynamically, suitable for use with trig functions with period 2*pi |
---|
| 1658 | /// \param xmin the lower bound for the grid |
---|
| 1659 | /// \param xmax upper bound for the grid |
---|
[1230] | 1660 | /// \param [in, out] grid the sampling grid. |
---|
| 1661 | virtual void get_sampling_grid(float_type xmin, float_type xmax, std::vector<float_type> &grid) const; |
---|
[807] | 1662 | }; |
---|
[1230] | 1663 | |
---|
[807] | 1664 | /// \brief compute cos(x) with its derivatives. |
---|
[1230] | 1665 | /// \ingroup math_functions |
---|
[807] | 1666 | /// |
---|
[1230] | 1667 | /// The factory function c2_factory::cos() creates *new c2_cos_p |
---|
| 1668 | template <typename float_type=double> class c2_cos_p : public c2_sin_p<float_type> { |
---|
[807] | 1669 | public: |
---|
[1230] | 1670 | /// \brief constructor. |
---|
| 1671 | c2_cos_p() : c2_sin_p<float_type>() {} |
---|
[807] | 1672 | |
---|
| 1673 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
| 1674 | { float_type q=std::cos(x); if(yprime) *yprime=-std::sin(x); if(yprime2) *yprime2=-q; return q; } |
---|
| 1675 | }; |
---|
[1230] | 1676 | |
---|
[807] | 1677 | /// \brief compute tan(x) with its derivatives. |
---|
[1230] | 1678 | /// \ingroup math_functions |
---|
[807] | 1679 | /// |
---|
[1230] | 1680 | /// The factory function c2_factory::tan() creates *new c2_tan_p |
---|
| 1681 | template <typename float_type=double> class c2_tan_p : public c2_function<float_type> { |
---|
[807] | 1682 | public: |
---|
[1230] | 1683 | /// \brief constructor. |
---|
| 1684 | c2_tan_p() : c2_function<float_type>() {} |
---|
[807] | 1685 | |
---|
| 1686 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
| 1687 | { |
---|
| 1688 | float_type c=std::cos(x), s=std::sin(x); |
---|
| 1689 | float_type t=s/c; |
---|
| 1690 | float_type yp=1/(c*c); |
---|
| 1691 | if(yprime) *yprime=yp; if(yprime2) *yprime2=2*t*yp; |
---|
| 1692 | return t; |
---|
| 1693 | } |
---|
| 1694 | }; |
---|
[1230] | 1695 | |
---|
[807] | 1696 | /// \brief compute log(x) with its derivatives. |
---|
[1230] | 1697 | /// \ingroup math_functions |
---|
[807] | 1698 | /// |
---|
[1230] | 1699 | /// The factory function c2_factory::log() creates *new c2_log_p |
---|
| 1700 | template <typename float_type=double> class c2_log_p : public c2_function<float_type> { |
---|
[807] | 1701 | public: |
---|
[1230] | 1702 | /// \brief constructor. |
---|
| 1703 | c2_log_p() : c2_function<float_type>() {} |
---|
[807] | 1704 | |
---|
| 1705 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
| 1706 | { if(yprime) *yprime=1.0/x; if(yprime2) *yprime2=-1.0/(x*x); return std::log(x); } |
---|
| 1707 | }; |
---|
[1230] | 1708 | |
---|
[807] | 1709 | /// \brief compute exp(x) with its derivatives. |
---|
[1230] | 1710 | /// \ingroup math_functions |
---|
[807] | 1711 | /// |
---|
[1230] | 1712 | /// The factory function c2_factory::exp() creates *new c2_exp_p |
---|
| 1713 | template <typename float_type=double> class c2_exp_p : public c2_function<float_type> { |
---|
[807] | 1714 | public: |
---|
[1230] | 1715 | /// \brief constructor. |
---|
| 1716 | c2_exp_p() : c2_function<float_type>() {} |
---|
[807] | 1717 | |
---|
| 1718 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
| 1719 | { float_type q=std::exp(x); if(yprime) *yprime=q; if(yprime2) *yprime2=q; return q; } |
---|
| 1720 | }; |
---|
[1230] | 1721 | |
---|
[807] | 1722 | /// \brief compute sqrt(x) with its derivatives. |
---|
[1230] | 1723 | /// \ingroup math_functions |
---|
[807] | 1724 | /// |
---|
[1230] | 1725 | /// The factory function c2_factory::sqrt() creates *new c2_sqrt_p() |
---|
| 1726 | template <typename float_type=double> class c2_sqrt_p : public c2_function<float_type> { |
---|
[807] | 1727 | public: |
---|
[1230] | 1728 | /// \brief constructor. |
---|
| 1729 | c2_sqrt_p() : c2_function<float_type>() {} |
---|
[807] | 1730 | |
---|
| 1731 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
| 1732 | { float_type q=std::sqrt(x); if(yprime) *yprime=0.5/q; if(yprime2) *yprime2=-0.25/(x*q); return q; } |
---|
| 1733 | }; |
---|
[1230] | 1734 | |
---|
[807] | 1735 | /// \brief compute scale/x with its derivatives. |
---|
[1230] | 1736 | /// \ingroup parametric_functions |
---|
[807] | 1737 | /// |
---|
[1230] | 1738 | /// The factory function c2_factory::recip() creates *new c2_recip_p |
---|
| 1739 | template <typename float_type=double> class c2_recip_p : public c2_function<float_type> { |
---|
[807] | 1740 | public: |
---|
[1230] | 1741 | /// \brief constructor. |
---|
| 1742 | c2_recip_p(float_type scale) : c2_function<float_type>(), rscale(scale) {} |
---|
[807] | 1743 | |
---|
| 1744 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
| 1745 | { |
---|
| 1746 | float_type q=1.0/x; |
---|
| 1747 | float_type y=rscale*q; |
---|
| 1748 | if(yprime) *yprime=-y*q; |
---|
| 1749 | if(yprime2) *yprime2=2*y*q*q; |
---|
| 1750 | return y; |
---|
| 1751 | } |
---|
| 1752 | /// \brief reset the scale factor |
---|
| 1753 | /// \param scale the new numerator |
---|
| 1754 | void reset(float_type scale) { rscale=scale; } |
---|
| 1755 | private: |
---|
| 1756 | float_type rscale; |
---|
| 1757 | }; |
---|
[1230] | 1758 | |
---|
[807] | 1759 | /// \brief compute x with its derivatives. |
---|
[1230] | 1760 | /// \ingroup math_functions |
---|
[807] | 1761 | /// |
---|
[1230] | 1762 | /// The factory function c2_factory::identity() creates *new c2_identity_p |
---|
| 1763 | template <typename float_type=double> class c2_identity_p : public c2_function<float_type> { |
---|
[807] | 1764 | public: |
---|
[1230] | 1765 | /// \brief constructor. |
---|
| 1766 | c2_identity_p() : c2_function<float_type>() {} |
---|
[807] | 1767 | |
---|
| 1768 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
| 1769 | { if(yprime) *yprime=1.0; if(yprime2) *yprime2=0; return x; } |
---|
| 1770 | }; |
---|
| 1771 | |
---|
| 1772 | /** |
---|
| 1773 | \brief create a linear mapping of another function |
---|
[1230] | 1774 | \ingroup parametric_functions |
---|
[807] | 1775 | for example, given a c2_function \a f |
---|
| 1776 | \code |
---|
[1230] | 1777 | c2_function<double> &F=c2_linear<double>(1.2, 2.0, 3.0)(f); |
---|
[807] | 1778 | \endcode |
---|
| 1779 | produces a new c2_function F=2.0+3.0*(\a f - 1.2) |
---|
[1230] | 1780 | |
---|
| 1781 | The factory function c2_factory::linear() creates *new c2_linear_p |
---|
[807] | 1782 | */ |
---|
[1230] | 1783 | template <typename float_type=double> class c2_linear_p : public c2_function<float_type> { |
---|
[807] | 1784 | public: |
---|
| 1785 | /// \brief Construct the operator f=y0 + slope * (x-x0) |
---|
| 1786 | /// \param x0 the x offset |
---|
| 1787 | /// \param y0 the y-intercept i.e. f(x0) |
---|
| 1788 | /// \param slope the slope of the mapping |
---|
[1230] | 1789 | c2_linear_p(float_type x0, float_type y0, float_type slope) : |
---|
| 1790 | c2_function<float_type>(), xint(x0), intercept(y0), m(slope) {} |
---|
[807] | 1791 | /// \brief Change the slope and intercepts after construction. |
---|
| 1792 | /// \param x0 the x offset |
---|
| 1793 | /// \param y0 the y-intercept |
---|
| 1794 | /// \param slope the slope of the mapping |
---|
| 1795 | void reset(float_type x0, float_type y0, float_type slope) { xint=x0; intercept=y0; m=slope; } |
---|
| 1796 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
| 1797 | { if(yprime) *yprime=m; if(yprime2) *yprime2=0; return m*(x-xint)+intercept; } |
---|
| 1798 | |
---|
| 1799 | private: |
---|
| 1800 | float_type xint, intercept, m; |
---|
| 1801 | protected: |
---|
[1230] | 1802 | c2_linear_p() {} // do not allow naked construction... it is usually an accident. |
---|
[807] | 1803 | }; |
---|
| 1804 | |
---|
| 1805 | /** |
---|
| 1806 | \brief create a quadratic mapping of another function |
---|
[1230] | 1807 | \ingroup parametric_functions |
---|
[807] | 1808 | for example, given a c2_function \a f |
---|
| 1809 | \code |
---|
[1230] | 1810 | c2_function<double> &F=c2_quadratic<double>(1.2, 2.0, 3.0, 4.0)(f); |
---|
[807] | 1811 | \endcode |
---|
| 1812 | produces a new c2_function F=2.0 + 3.0*(f-1.2) + 4.0*(f-1.2)^2 |
---|
| 1813 | |
---|
| 1814 | note that the parameters are overdetermined, but allows the flexibility of two different representations |
---|
| 1815 | |
---|
[1230] | 1816 | The factory function c2_factory::quadratic() creates *new c2_quadratic_p |
---|
[807] | 1817 | */ |
---|
[1230] | 1818 | template <typename float_type=double> class c2_quadratic_p : public c2_function<float_type> { |
---|
[807] | 1819 | public: |
---|
| 1820 | /// \brief Construct the operator |
---|
| 1821 | /// \param x0 the center around which the powers are computed |
---|
| 1822 | /// \param y0 the value of the function at \a x = \a x0 |
---|
| 1823 | /// \param xcoef the scale on the (\a x - \a x0) term |
---|
| 1824 | /// \param x2coef the scale on the (\a x - \a x0)^2 term |
---|
[1230] | 1825 | c2_quadratic_p(float_type x0, float_type y0, float_type xcoef, float_type x2coef) : |
---|
| 1826 | c2_function<float_type>(), intercept(y0), center(x0), a(x2coef), b(xcoef) {} |
---|
| 1827 | /// \brief Modify the coefficients after construction |
---|
[807] | 1828 | /// \param x0 the new center around which the powers are computed |
---|
| 1829 | /// \param y0 the new value of the function at \a x = \a x0 |
---|
| 1830 | /// \param xcoef the new scale on the (\a x - \a x0) term |
---|
| 1831 | /// \param x2coef the new scale on the (\a x - \a x0)^2 term |
---|
| 1832 | void reset(float_type x0, float_type y0, float_type xcoef, float_type x2coef) { intercept=y0; center=x0; a=x2coef; b=xcoef; } |
---|
| 1833 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
| 1834 | { float_type dx=x-center; if(yprime) *yprime=2*a*dx+b; if(yprime2) *yprime2=2*a; return a*dx*dx+b*dx+intercept; } |
---|
| 1835 | |
---|
| 1836 | private: |
---|
| 1837 | float_type intercept, center, a, b; |
---|
| 1838 | protected: |
---|
[1230] | 1839 | c2_quadratic_p() {} // do not allow naked construction... it is usually an accident. |
---|
[807] | 1840 | }; |
---|
| 1841 | |
---|
| 1842 | /** |
---|
| 1843 | \brief create a power law mapping of another function |
---|
[1230] | 1844 | \ingroup parametric_functions |
---|
[807] | 1845 | for example, given a c2_function \a f |
---|
| 1846 | \code |
---|
[1230] | 1847 | c2_power_law_p<double> PLaw(1.2, 2.5); |
---|
| 1848 | c2_composed_function_p<double> &F=PLaw(f); |
---|
[807] | 1849 | \endcode |
---|
| 1850 | produces a new c2_function F=1.2 * f^2.5 |
---|
| 1851 | |
---|
[1230] | 1852 | The factory function c2_factory::power_law() creates *new c2_power_law_p |
---|
[807] | 1853 | */ |
---|
[1230] | 1854 | template <typename float_type=double> class c2_power_law_p : public c2_function<float_type> { |
---|
[807] | 1855 | public: |
---|
| 1856 | /// \brief Construct the operator |
---|
| 1857 | /// \param scale the multipler |
---|
| 1858 | /// \param power the exponent |
---|
[1230] | 1859 | c2_power_law_p(float_type scale, float_type power) : |
---|
| 1860 | c2_function<float_type>(), a(scale), b(power) {} |
---|
[807] | 1861 | /// \brief Modify the mapping after construction |
---|
| 1862 | /// \param scale the new multipler |
---|
| 1863 | /// \param power the new exponent |
---|
| 1864 | void reset(float_type scale, float_type power) { a=scale; b=power; } |
---|
| 1865 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
| 1866 | { float_type q=a*std::pow(x,b-2); if(yprime) *yprime=b*q*x; if(yprime2) *yprime2=b*(b-1)*q; return q*x*x; } |
---|
| 1867 | |
---|
| 1868 | private: |
---|
| 1869 | float_type a, b; |
---|
| 1870 | protected: |
---|
[1230] | 1871 | c2_power_law_p() {} // do not allow naked construction... it is usually an accident. |
---|
[807] | 1872 | }; |
---|
| 1873 | |
---|
| 1874 | /** |
---|
| 1875 | \brief create the formal inverse function of another function |
---|
[1230] | 1876 | \ingroup containers |
---|
[807] | 1877 | for example, given a c2_function \a f |
---|
| 1878 | \code |
---|
| 1879 | c2_inverse_function<double> inv(f); |
---|
| 1880 | a=f(x); |
---|
| 1881 | x1=inv(a); |
---|
| 1882 | \endcode |
---|
| 1883 | will return x1=x to machine precision. The important part of this |
---|
| 1884 | is that the resulting function is a first-class c2_function, so it knows its |
---|
| 1885 | derivatives, too, unlike the case of a simple root-finding inverse. This means |
---|
| 1886 | it can be integrated (for example) quite efficiently. |
---|
[1230] | 1887 | |
---|
| 1888 | \see \ref combined_inversion_hinting_sampling |
---|
| 1889 | |
---|
| 1890 | The factory function c2_factory::inverse_function() creates *new c2_inverse_function_p |
---|
| 1891 | */ |
---|
| 1892 | template <typename float_type=double> class c2_inverse_function_p : public c2_function<float_type> { |
---|
[807] | 1893 | public: |
---|
| 1894 | /// \brief Construct the operator |
---|
| 1895 | /// \param source the function to be inverted |
---|
[1230] | 1896 | c2_inverse_function_p(const c2_function<float_type> &source); |
---|
[807] | 1897 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception); |
---|
| 1898 | |
---|
| 1899 | /// \brief give the function a hint as to where to look for its inverse |
---|
| 1900 | /// \param hint the likely value of the inverse, which defaults to whatever the evaluation returned. |
---|
| 1901 | void set_start_hint(float_type hint) const { start_hint=hint; } |
---|
| 1902 | |
---|
| 1903 | /// \brief get the starting hint. |
---|
| 1904 | /// |
---|
| 1905 | /// This is virtual so if there is a better way, this can be easily overridden. |
---|
| 1906 | /// It is used in value_with_derivatives() to guess where to start the root finder. |
---|
| 1907 | /// \param x the abscissa for which an estimate is needed |
---|
[1230] | 1908 | virtual float_type get_start_hint(float_type x) const |
---|
| 1909 | { return hinting_function.valid()? hinting_function(x) : start_hint; } |
---|
| 1910 | |
---|
| 1911 | /// \brief set or unset the approximate function used to start the root finder |
---|
| 1912 | /// \anchor set_hinting_function_discussion |
---|
| 1913 | /// A hinting function is mostly useful if the evaluation of this inverse is |
---|
| 1914 | /// going to be carried out in very non-local order, so the root finder has to start over |
---|
| 1915 | /// for each step. If most evaluations are going to be made in fairly localized clusters (scanning |
---|
| 1916 | /// through the function, for example), the default mechanism used (which just remembers the last point) |
---|
| 1917 | /// is almost certainly faster. |
---|
| 1918 | /// |
---|
| 1919 | /// Typically, the hinting function is likely to be set up by creating the inverse function, |
---|
| 1920 | /// and then adaptively sampling an interpolating function from it, and then using the result |
---|
| 1921 | /// to hint it. Another way, if the parent function is already an interpolating function, is just to create |
---|
| 1922 | /// a version of the parent with the x & y coordinates reversed. |
---|
| 1923 | /// |
---|
| 1924 | /// \see \ref combined_inversion_hinting_sampling |
---|
| 1925 | /// |
---|
| 1926 | /// \param hint_func the function that is an approximate inverse of the parent of this inverse_function |
---|
| 1927 | void set_hinting_function(const c2_function<float_type> *hint_func) |
---|
| 1928 | { hinting_function.set_function(hint_func); } |
---|
| 1929 | /// \brief set the hinting function from a pointer. |
---|
| 1930 | /// |
---|
| 1931 | /// See \ref set_hinting_function_discussion "discussion" |
---|
| 1932 | /// \param hint_func the container holding the function |
---|
| 1933 | void set_hinting_function(const c2_const_ptr<float_type> hint_func) |
---|
| 1934 | { hinting_function=hint_func; } |
---|
| 1935 | |
---|
[807] | 1936 | protected: |
---|
[1230] | 1937 | c2_inverse_function_p() {} // do not allow naked construction... it is usually an accident. |
---|
| 1938 | mutable float_type start_hint; |
---|
| 1939 | const c2_const_ptr<float_type> func; |
---|
| 1940 | c2_const_ptr<float_type> hinting_function; |
---|
[807] | 1941 | }; |
---|
| 1942 | |
---|
| 1943 | /** |
---|
| 1944 | \brief |
---|
[1230] | 1945 | An interpolating_function_p which is the cumulative integral of a histogram. |
---|
| 1946 | \ingroup interpolators |
---|
[807] | 1947 | Note than binedges should be one element longer than binheights, since the lower & upper edges are specified. |
---|
| 1948 | Note that this is a malformed spline, since the second derivatives are all zero, so it has less continuity. |
---|
| 1949 | Also, note that the bin edges can be given in backwards order to generate the |
---|
| 1950 | reversed accumulation (starting at the high end) |
---|
| 1951 | */ |
---|
| 1952 | |
---|
[1230] | 1953 | template <typename float_type=double> class accumulated_histogram : public interpolating_function_p <float_type> { |
---|
[807] | 1954 | public: |
---|
| 1955 | /// \brief Construct the integrated histogram |
---|
| 1956 | /// \param binedges the edges of the bins in \a binheights. It should have one more element than \a binheights |
---|
| 1957 | /// \param binheights the number of counts in each bin. |
---|
| 1958 | /// \param normalize if true, normalize integral to 1 |
---|
| 1959 | /// \param inverse_function if true, drop zero channels, and return inverse function for random generation |
---|
| 1960 | /// \param drop_zeros eliminate null bins before integrating, so integral is strictly monotonic. |
---|
| 1961 | accumulated_histogram(const std::vector<float_type>binedges, const std::vector<float_type> binheights, |
---|
| 1962 | bool normalize=false, bool inverse_function=false, bool drop_zeros=true); |
---|
| 1963 | |
---|
| 1964 | }; |
---|
| 1965 | |
---|
| 1966 | /** |
---|
[1230] | 1967 | \anchor inverse_integrated_density_bins |
---|
| 1968 | \brief construct from a grid of points and an std::vector of probability densities (un-normalized) |
---|
| 1969 | \see \ref random_subsec "Arbitrary random generation" |
---|
| 1970 | \ingroup interpolators |
---|
| 1971 | inverse_integrated_density starts with a probability density std::vector, generates the integral, |
---|
| 1972 | and generates an interpolating_function_p of the inverse function which, when evaluated using a uniform random on [0,1] returns values |
---|
| 1973 | with a density distribution equal to the input distribution |
---|
| 1974 | If the data are passed in reverse order (large X first), the integral is carried out from the big end. |
---|
| 1975 | |
---|
| 1976 | \param bins if \a bins .size()==\a binheights .size(), the centers of the bins. \n |
---|
| 1977 | if \a bins .size()==\a binheights .size()+1, the edges of the bins |
---|
| 1978 | \param binheights a vector which describes the density of the random number distribution to be produced. |
---|
| 1979 | Note density... the numbers in the bins are not counts, but counts/unit bin width. |
---|
| 1980 | \return an interpolating_function_p of the type requested in the template which, |
---|
| 1981 | if evaluated randomly with a uniform variate on [0,1] produces numbers |
---|
| 1982 | distributed according to \a binheights |
---|
[807] | 1983 | */ |
---|
| 1984 | |
---|
[1230] | 1985 | template <typename float_type, typename Final> |
---|
| 1986 | interpolating_function_p<float_type> & inverse_integrated_density_bins( |
---|
| 1987 | const std::vector<float_type> &bins, const std::vector<float_type> &binheights) |
---|
| 1988 | throw(c2_exception); |
---|
[807] | 1989 | |
---|
| 1990 | /** |
---|
[1230] | 1991 | \anchor inverse_integrated_density_function |
---|
| 1992 | \brief construct from a grid of points and a c2_function of probability densities (un-normalized) |
---|
| 1993 | \see \ref random_subsec "Arbitrary random generation" |
---|
| 1994 | \ingroup interpolators |
---|
| 1995 | inverse_integrated_density starts with a probability density std::vector, generates the integral, |
---|
| 1996 | and generates an interpolating_function_p of the inverse function which, when evaluated using a uniform random on [0,1] returns values |
---|
| 1997 | with a density distribution equal to the input distribution |
---|
| 1998 | If the data are passed in reverse order (large X first), the integral is carried out from the big end. |
---|
| 1999 | |
---|
| 2000 | \param bincenters the centers of the bins. |
---|
| 2001 | \param binheights a c2_function which describes the density of the random number distribution to be produced. |
---|
| 2002 | \return an interpolating_function_p of the type requested in the template which, |
---|
| 2003 | if evaluated randomly with a uniform variate on [0,1] produces numbers |
---|
| 2004 | distributed according to \a binheights |
---|
| 2005 | */ |
---|
| 2006 | template <typename float_type, typename Final> |
---|
| 2007 | interpolating_function_p<float_type> & inverse_integrated_density_function( |
---|
| 2008 | const std::vector<float_type> &bincenters, const c2_function<float_type> &binheights) |
---|
| 2009 | throw(c2_exception); |
---|
[807] | 2010 | |
---|
| 2011 | /// \brief create a c2_function which smoothly connects two other c2_functions. |
---|
[1230] | 2012 | /// \ingroup parametric_functions |
---|
[807] | 2013 | /// This takes two points and generates a polynomial which matches two c2_function arguments |
---|
| 2014 | /// at those two points, with two derivatives at each point, and an arbitrary value at the center of the |
---|
| 2015 | /// region. It is useful for splicing together functions over rough spots (0/0, for example). |
---|
| 2016 | /// |
---|
| 2017 | /// If \a auto_center is true, the value at the midpoint is computed so that the resulting polynomial is |
---|
| 2018 | /// of order 5. If \a auto_center is false, the value \a y1 is used at the midpoint, resulting in a |
---|
| 2019 | /// polynomial of order 6. |
---|
[1230] | 2020 | /// |
---|
| 2021 | /// This is usually used in conjunction with c2_piecewise_function_p to assemble an apparently seamless |
---|
| 2022 | /// function from a series of segments. |
---|
| 2023 | /// \see \ref piecewise_applications_subsec "Sample Applications" and \ref c2_function::adaptively_sample() "Adaptive sampling" |
---|
| 2024 | /// |
---|
| 2025 | /// The factory function c2_factory::connector_function() creates *new c2_connector_function_p |
---|
| 2026 | template <typename float_type=double> class c2_connector_function_p : public c2_function<float_type> { |
---|
[807] | 2027 | public: |
---|
[1230] | 2028 | /// \brief construct the container from two functions |
---|
[807] | 2029 | /// \param x0 the point at which to match \a f1 and its derivatives |
---|
[1230] | 2030 | /// \param f0 the function on the left side to be connected |
---|
[807] | 2031 | /// \param x2 the point at which to match \a f2 and its derivatives |
---|
[1230] | 2032 | /// \param f2 the function on the right side to be connected |
---|
[807] | 2033 | /// \param auto_center if true, no midpoint value is specified. If false, match the value \a y1 at the midpoint |
---|
| 2034 | /// \param y1 the value to match at the midpoint, if \a auto_center is false |
---|
[1230] | 2035 | /// \return a c2_function with domain (\a x0,\a x2) which smoothly connects \a f0(x0) and \a f2(x2) |
---|
| 2036 | c2_connector_function_p(float_type x0, const c2_function<float_type> &f0, float_type x2, const c2_function<float_type> &f2, |
---|
[807] | 2037 | bool auto_center, float_type y1); |
---|
[1230] | 2038 | /// \brief construct the container from numerical values |
---|
| 2039 | /// \param x0 the position of the left edge |
---|
| 2040 | /// \param y0 the function derivative on the left boundary |
---|
| 2041 | /// \param yp0 the function second derivative on the left boundary |
---|
| 2042 | /// \param ypp0 the function value on the left boundary |
---|
| 2043 | /// \param x2 the position of the right edge |
---|
| 2044 | /// \param y2 the function derivative on the right boundary |
---|
| 2045 | /// \param yp2 the function second derivative on the right boundary |
---|
| 2046 | /// \param ypp2 the function value on the right boundary |
---|
| 2047 | /// \param auto_center if true, no midpoint value is specified. If false, match the value \a y1 at the midpoint |
---|
| 2048 | /// \param y1 the value to match at the midpoint, if \a auto_center is false |
---|
| 2049 | /// \return a c2_function with domain (\a x0,\a x2) which smoothly connects the points described |
---|
| 2050 | /// \anchor c2_connector_raw_init_docs |
---|
| 2051 | c2_connector_function_p( |
---|
| 2052 | float_type x0, float_type y0, float_type yp0, float_type ypp0, |
---|
| 2053 | float_type x2, float_type y2, float_type yp2, float_type ypp2, |
---|
| 2054 | bool auto_center, float_type y1); |
---|
| 2055 | /// \brief construct the container from c2_fblock<float_type> objects |
---|
| 2056 | /// \param fb0 the left edge |
---|
| 2057 | /// \param fb2 the right edge |
---|
| 2058 | /// \param auto_center if true, no midpoint value is specified. If false, match the value \a y1 at the midpoint |
---|
| 2059 | /// \param y1 the value to match at the midpoint, if \a auto_center is false |
---|
| 2060 | /// \return a c2_function with domain (\a fb0.x,\a fb2.x) which smoothly connects \a fb0 and \a fb2 |
---|
| 2061 | c2_connector_function_p( |
---|
| 2062 | const c2_fblock<float_type> &fb0, |
---|
| 2063 | const c2_fblock<float_type> &fb2, |
---|
| 2064 | bool auto_center, float_type y1); |
---|
| 2065 | |
---|
[807] | 2066 | /// \brief destructor |
---|
[1230] | 2067 | virtual ~c2_connector_function_p(); |
---|
[807] | 2068 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw (c2_exception); |
---|
| 2069 | protected: |
---|
[1230] | 2070 | /// \brief fill container numerically |
---|
| 2071 | void init( |
---|
| 2072 | const c2_fblock<float_type> &fb0, |
---|
| 2073 | const c2_fblock<float_type> &fb2, |
---|
| 2074 | bool auto_center, float_type y1); |
---|
| 2075 | |
---|
| 2076 | float_type fhinv, fy1, fa, fb, fc, fd, fe, ff; |
---|
[807] | 2077 | }; |
---|
| 2078 | |
---|
| 2079 | /// \brief create a c2_function which is a piecewise assembly of other c2_functions. |
---|
[1230] | 2080 | /// \ingroup containers |
---|
[807] | 2081 | /// The functions must have increasing, non-overlapping domains. Any empty space |
---|
| 2082 | /// between functions will be filled with a linear interpolation. |
---|
[1230] | 2083 | /// |
---|
| 2084 | /// \note If you want a smooth connection, instead of the default linear interpolation, |
---|
| 2085 | /// create a c2_connector_function_p to bridge the gap. The linear interpolation is intended |
---|
| 2086 | /// to be a barely intelligent bridge, and may never get used by anyone. |
---|
| 2087 | /// |
---|
[807] | 2088 | /// \note The creation of the container results in the creation of an explicit sampling grid. |
---|
| 2089 | /// If this is used with functions with a large domain, or which generate very dense sampling grids, |
---|
| 2090 | /// it could eat a lot of memory. Do not abuse this by using functions which can generate gigantic grids. |
---|
| 2091 | /// |
---|
[1230] | 2092 | /// \see \ref piecewise_applications_subsec "Sample Applications" \n |
---|
| 2093 | /// c2_plugin_function_p page \n |
---|
| 2094 | /// c2_connector_function_p page \n |
---|
| 2095 | /// \ref c2_function::adaptively_sample() "Adaptive sampling" |
---|
| 2096 | /// |
---|
| 2097 | /// The factory function c2_factory::piecewise_function() creates *new c2_piecewise_function_p |
---|
| 2098 | template <typename float_type=double> class c2_piecewise_function_p : public c2_function<float_type> { |
---|
[807] | 2099 | public: |
---|
| 2100 | /// \brief construct the container |
---|
[1230] | 2101 | c2_piecewise_function_p(); |
---|
[807] | 2102 | /// \brief destructor |
---|
[1230] | 2103 | virtual ~c2_piecewise_function_p(); |
---|
[807] | 2104 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw (c2_exception); |
---|
| 2105 | /// \brief append a new function to the sequence |
---|
| 2106 | /// |
---|
| 2107 | /// This takes a c2_function, and appends it onto the end of the piecewise collection. |
---|
| 2108 | /// The domain of the function (which MUST be set) specifies the place it will be used in |
---|
| 2109 | /// the final function. If the domain exactly abuts the domain of the previous function, it |
---|
| 2110 | /// will be directly attached. If there is a gap, the gap will be filled in by linear interpolation. |
---|
| 2111 | /// \param func a c2_function with a defined domain to be appended to the collection |
---|
[1230] | 2112 | void append_function(const c2_function<float_type> &func) throw (c2_exception); |
---|
[807] | 2113 | protected: |
---|
[1230] | 2114 | std::vector<c2_const_ptr<float_type> > functions; |
---|
[807] | 2115 | mutable int lastKLow; |
---|
| 2116 | }; |
---|
| 2117 | |
---|
| 2118 | #include "c2_function.icc" |
---|
| 2119 | |
---|
| 2120 | #endif |
---|