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 | // |
---|
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 | * |
---|
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 |
---|
37 | */ |
---|
38 | |
---|
39 | #ifndef __has_c2_function_hh |
---|
40 | #define __has_c2_function_hh 1 |
---|
41 | |
---|
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 | |
---|
53 | #include <cmath> |
---|
54 | #include <vector> |
---|
55 | #include <utility> |
---|
56 | #include <string> |
---|
57 | #include <stdexcept> |
---|
58 | #include <typeinfo> |
---|
59 | #include <sstream> |
---|
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. |
---|
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 | */ |
---|
94 | |
---|
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 | |
---|
115 | /** |
---|
116 | \brief the parent class for all c2_functions. |
---|
117 | \ingroup abstract_classes |
---|
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, |
---|
120 | composed via c2_composed_function_, |
---|
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 |
---|
130 | inverse_integrated_density_function(). |
---|
131 | |
---|
132 | For a discussion of memory management, see \ref memory_management |
---|
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 |
---|
140 | "c2_function.hh,v 1.238 2008/05/22 12:45:19 marcus Exp"; |
---|
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 |
---|
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 | } |
---|
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. |
---|
202 | /// \see \ref rootfinder_subsec "Root finding sample" |
---|
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. |
---|
224 | /// \return sum of partial integrals, which is the definite integral from the first value in \a xgrid to the last. |
---|
225 | float_type partial_integrals(std::vector<float_type> xgrid, std::vector<float_type> *partials = 0, |
---|
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); |
---|
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. |
---|
238 | /// It will automatically be sized appropriately, if provided, to contain \a n - 1 elements where \a n is the length of \a xgrid |
---|
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. |
---|
247 | /// \return sum of partial integrals, which is the definite integral from the first value in \a xgrid to the last. |
---|
248 | float_type integral(float_type xmin, float_type xmax, std::vector<float_type> *partials = 0, |
---|
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); |
---|
251 | |
---|
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 | |
---|
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 | /// |
---|
302 | /// It is not maintained automatically in general! The root finder, integrator, and sampler do increment it. |
---|
303 | /// \return number of evaluations logged since last reset. |
---|
304 | volatile size_t get_evaluations() const { return evaluations; } |
---|
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 |
---|
315 | bool check_monotonicity(const std::vector<float_type> &data, const char message[]) const throw(c2_exception); |
---|
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 | |
---|
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 | |
---|
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 |
---|
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 ; |
---|
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 |
---|
344 | /// \param [in,out] grid the grid to refine in place |
---|
345 | /// \param refinement the number of new steps for each old step |
---|
346 | void refine_sampling_grid(std::vector<float_type> &grid, size_t refinement) const; |
---|
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. |
---|
353 | c2_function<float_type> &normalized_function(float_type xmin, float_type xmax, float_type norm=1.0) const throw(c2_exception); |
---|
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. |
---|
359 | c2_function<float_type> &square_normalized_function(float_type xmin, float_type xmax, float_type norm=1.0) const throw(c2_exception); |
---|
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. |
---|
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); |
---|
369 | |
---|
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; } |
---|
401 | |
---|
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; } |
---|
423 | |
---|
424 | protected: |
---|
425 | c2_function(const c2_function<float_type> &src) : sampling_grid(0), |
---|
426 | no_overwrite_grid(false), |
---|
427 | fXMin(src.fXMin), fXMax(src.fXMax), root_info(0), owner_count(0) |
---|
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()), |
---|
432 | fXMax(std::numeric_limits<float_type>::max()), root_info(0), owner_count(0) |
---|
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 | |
---|
443 | std::vector<float_type> * sampling_grid; |
---|
444 | bool no_overwrite_grid; |
---|
445 | |
---|
446 | float_type fXMin, fXMax; |
---|
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 | }; |
---|
469 | |
---|
470 | |
---|
471 | /// \brief structure used to pass information recursively in integrator. |
---|
472 | /// |
---|
473 | /// the \a abs_tol is scaled by a factor of two at each division. |
---|
474 | /// Everything else is just passed down. |
---|
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; |
---|
481 | }; |
---|
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 | }; |
---|
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. |
---|
507 | float_type integrate_step(struct c2_integrate_recur &rb) const throw(c2_exception); |
---|
508 | |
---|
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, |
---|
517 | /// to avoid the necessity of evaluating the function on the brackets every time |
---|
518 | /// if the brackets have not been changed. |
---|
519 | /// it is declared as a pointer, since many c2_functions may never need one allocated |
---|
520 | mutable struct c2_root_info *root_info; |
---|
521 | |
---|
522 | mutable size_t owner_count; |
---|
523 | }; |
---|
524 | |
---|
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 | */ |
---|
765 | /// \brief a container into which any other c2_function can be dropped, to allow expressions |
---|
766 | /// with replacable components. |
---|
767 | /// \ingroup containers |
---|
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 |
---|
773 | /// be used, for example, to patch badly-behaved functions with c2_piecewise_function_p by |
---|
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. |
---|
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> { |
---|
783 | public: |
---|
784 | /// \brief construct the container with no function |
---|
785 | c2_plugin_function_p() : c2_function<float_type>(), func() {} |
---|
786 | /// \brief construct the container with a pre-defined function |
---|
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()); |
---|
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 | { |
---|
800 | if(!func.valid()) throw c2_exception("c2_plugin_function_p called uninitialized"); |
---|
801 | return func->value_with_derivatives(x, yprime, yprime2); |
---|
802 | } |
---|
803 | /// \brief destructor |
---|
804 | ~c2_plugin_function_p() { } |
---|
805 | |
---|
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 | } |
---|
815 | protected: |
---|
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() { } |
---|
835 | |
---|
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(); } |
---|
839 | }; |
---|
840 | |
---|
841 | /// \brief Provides support for c2_function objects which are constructed from two other c2_function |
---|
842 | /// objects. |
---|
843 | /// \ingroup abstract_classes |
---|
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 | { |
---|
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); |
---|
853 | } |
---|
854 | |
---|
855 | /// \brief destructor releases ownership of member functions |
---|
856 | /// |
---|
857 | virtual ~c2_binary_function() { } |
---|
858 | |
---|
859 | protected: |
---|
860 | /// \brief construct the binary function |
---|
861 | /// \param combiner pointer to the function which actualy knows how to execute the binary |
---|
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 |
---|
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) |
---|
869 | { |
---|
870 | set_domain( |
---|
871 | (left.xmin() > right.xmin()) ? left.xmin() : right.xmin(), |
---|
872 | (left.xmax() < right.xmax()) ? left.xmax() : right.xmax() |
---|
873 | ); |
---|
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 |
---|
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) { } |
---|
882 | |
---|
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 | |
---|
892 | }; |
---|
893 | |
---|
894 | /// \brief Create a very lightweight method to return a scalar multiple of another function. |
---|
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> { |
---|
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 |
---|
904 | c2_scaled_function_p(const c2_function<float_type> &outer, float_type scale) : |
---|
905 | c2_function<float_type>(), func(outer), yscale(scale) { } |
---|
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: |
---|
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; |
---|
926 | float_type yscale; |
---|
927 | }; |
---|
928 | |
---|
929 | /// \brief A container into which any other c2_function can be dropped. |
---|
930 | /// \ingroup containers |
---|
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 | /// |
---|
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> { |
---|
937 | public: |
---|
938 | /// \brief construct the container |
---|
939 | /// |
---|
940 | /// \param f the function to be cached |
---|
941 | c2_cached_function_p(const c2_function<float_type> &f) : c2_function<float_type>(), |
---|
942 | func(f), init(false) {} |
---|
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: |
---|
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; |
---|
962 | mutable bool init; |
---|
963 | mutable float_type x0, y, yp, ypp; |
---|
964 | |
---|
965 | }; |
---|
966 | |
---|
967 | /// \brief Provides function composition (nesting) |
---|
968 | /// \ingroup arithmetic_functions |
---|
969 | /// This allows evaluation of \a f(g(x)) where \a f and \a g are c2_function objects. |
---|
970 | /// |
---|
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> { |
---|
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 |
---|
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()); } |
---|
981 | /// \brief Create a stub just for the combiner to avoid statics. |
---|
982 | c2_composed_function_p() : c2_binary_function<float_type>(combine) {} |
---|
983 | |
---|
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) |
---|
987 | { |
---|
988 | float_type y0, y1; |
---|
989 | if(yprime || yprime2) { |
---|
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; |
---|
995 | } else { |
---|
996 | y0=right(x); |
---|
997 | y1=left(y0); |
---|
998 | } |
---|
999 | return y1; |
---|
1000 | } |
---|
1001 | }; |
---|
1002 | |
---|
1003 | /// \brief create a c2_function which is the sum of two other c2_function objects. |
---|
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> { |
---|
1007 | public: |
---|
1008 | /// \brief construct \a left + \a right |
---|
1009 | /// \param left the left function |
---|
1010 | /// \param right the right function |
---|
1011 | c2_sum_p(const c2_function<float_type> &left, const c2_function<float_type> &right) : c2_binary_function<float_type>(combine, left, right) {} |
---|
1012 | /// \brief Create a stub just for the combiner to avoid statics. |
---|
1013 | c2_sum_p() : c2_binary_function<float_type>(combine) {} ; // create a stub just for the combiner to avoid statics |
---|
1014 | |
---|
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) |
---|
1018 | { |
---|
1019 | float_type y0, y1; |
---|
1020 | if(yprime || yprime2) { |
---|
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; |
---|
1026 | } else { |
---|
1027 | y0=left(x); |
---|
1028 | y1=right(x); |
---|
1029 | } |
---|
1030 | return y0+y1; |
---|
1031 | } |
---|
1032 | }; |
---|
1033 | |
---|
1034 | |
---|
1035 | /// \brief create a c2_function which is the difference of two other c2_functions. |
---|
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> { |
---|
1039 | public: |
---|
1040 | /// \brief construct \a left - \a right |
---|
1041 | /// \param left the left function |
---|
1042 | /// \param right the right function |
---|
1043 | c2_diff_p(const c2_function<float_type> &left, const c2_function<float_type> &right) : c2_binary_function<float_type>(combine, left, right) {} |
---|
1044 | /// \brief Create a stub just for the combiner to avoid statics. |
---|
1045 | c2_diff_p() : c2_binary_function<float_type>(combine) {} ; // create a stub just for the combiner to avoid statics |
---|
1046 | |
---|
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) |
---|
1050 | { |
---|
1051 | float_type y0, y1; |
---|
1052 | if(yprime || yprime2) { |
---|
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; |
---|
1058 | } else { |
---|
1059 | y0=left(x); |
---|
1060 | y1=right(x); |
---|
1061 | } |
---|
1062 | return y0-y1; |
---|
1063 | } |
---|
1064 | }; |
---|
1065 | |
---|
1066 | |
---|
1067 | /// \brief create a c2_function which is the product of two other c2_functions. |
---|
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> { |
---|
1071 | public: |
---|
1072 | /// \brief construct \a left * \a right |
---|
1073 | /// \param left the left function |
---|
1074 | /// \param right the right function |
---|
1075 | c2_product_p(const c2_function<float_type> &left, const c2_function<float_type> &right) : c2_binary_function<float_type>(combine, left, right) {} |
---|
1076 | /// \brief Create a stub just for the combiner to avoid statics. |
---|
1077 | c2_product_p() : c2_binary_function<float_type>(combine) {} ; // create a stub just for the combiner to avoid statics |
---|
1078 | |
---|
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) |
---|
1082 | { |
---|
1083 | float_type y0, y1; |
---|
1084 | if(yprime || yprime2) { |
---|
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; |
---|
1090 | } else { |
---|
1091 | y0=left(x); |
---|
1092 | y1=right(x); |
---|
1093 | } |
---|
1094 | return y0*y1; |
---|
1095 | } |
---|
1096 | }; |
---|
1097 | |
---|
1098 | |
---|
1099 | /// \brief create a c2_function which is the ratio of two other c2_functions. |
---|
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> { |
---|
1103 | public: |
---|
1104 | /// \brief construct \a left / \a right |
---|
1105 | /// \param left the left function |
---|
1106 | /// \param right the right function |
---|
1107 | c2_ratio_p(const c2_function<float_type> &left, const c2_function<float_type> &right) : c2_binary_function<float_type>(combine, left, right) {} |
---|
1108 | /// \brief Create a stub just for the combiner to avoid statics. |
---|
1109 | c2_ratio_p() : c2_binary_function<float_type>(combine) {} ; // create a stub just for the combiner to avoid statics |
---|
1110 | |
---|
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) |
---|
1114 | { |
---|
1115 | float_type y0, y1; |
---|
1116 | if(yprime || yprime2) { |
---|
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); |
---|
1122 | } else { |
---|
1123 | y0=left(x); |
---|
1124 | y1=right(x); |
---|
1125 | } |
---|
1126 | return y0/y1; |
---|
1127 | } |
---|
1128 | |
---|
1129 | }; |
---|
1130 | |
---|
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> { |
---|
1136 | public: |
---|
1137 | c2_constant_p(float_type x) : c2_function<float_type>(), value(x) {} |
---|
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 | |
---|
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 | |
---|
1345 | /** |
---|
1346 | \brief create a cubic spline interpolation of a set of (x,y) pairs |
---|
1347 | \ingroup interpolators |
---|
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 |
---|
1361 | c2_ptr<double> c2p; |
---|
1362 | c2_factory<double> c2; |
---|
1363 | std::vector<double> xvals(10), yvals(10); |
---|
1364 | // < fill in xvals and yvals > |
---|
1365 | c2p myfunc=c2.interpolating_function().load(xvals, yvals,true,0,true,0); |
---|
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 |
---|
1372 | |
---|
1373 | The factory function c2_factory::interpolating_function() creates *new interpolating_function_p() |
---|
1374 | */ |
---|
1375 | |
---|
1376 | template <typename float_type=double> class interpolating_function_p : public c2_function<float_type> { |
---|
1377 | public: |
---|
1378 | /// \brief an empty linear-linear cubic-spline interpolating_function_p |
---|
1379 | /// |
---|
1380 | /// lots to say here, but see Numerical Recipes for a discussion of cubic splines. |
---|
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. |
---|
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 |
---|
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); |
---|
1441 | |
---|
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); |
---|
1478 | |
---|
1479 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception); |
---|
1480 | |
---|
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 | |
---|
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. |
---|
1491 | /// The vectors will have their sizes set correctly on return. |
---|
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 | |
---|
1512 | /// \brief create a new interpolating_function_p which is the \a source |
---|
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 |
---|
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; |
---|
1519 | |
---|
1520 | /// \brief create a new interpolating_function_p which is the parent interpolating_function_p |
---|
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. |
---|
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 |
---|
1529 | ) const; |
---|
1530 | /// \brief produce a newly resampled interpolating_function_p which is the specified sum. |
---|
1531 | /// \param rhs the function to add, pointwise |
---|
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. |
---|
1536 | /// \param rhs the function to subtract, pointwise |
---|
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. |
---|
1541 | /// \param rhs the function to multiply, pointwise |
---|
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. |
---|
1546 | /// \param rhs the function to divide, pointwise |
---|
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; |
---|
1559 | |
---|
1560 | protected: |
---|
1561 | /// \brief create the spline coefficients |
---|
1562 | void spline( |
---|
1563 | bool lowerSlopeNatural, float_type lowerSlope, |
---|
1564 | bool upperSlopeNatural, float_type upperSlope |
---|
1565 | ) throw(c2_exception); |
---|
1566 | |
---|
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 | |
---|
1570 | std::vector<float_type> Xraw, X, F, y2; |
---|
1571 | c2_const_ptr<float_type> sampler_function; |
---|
1572 | bool xInverted; |
---|
1573 | mutable size_t lastKLow; |
---|
1574 | }; |
---|
1575 | |
---|
1576 | /// \brief A spline with X transformed into log space. |
---|
1577 | /// \ingroup interpolators |
---|
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. |
---|
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> { |
---|
1583 | public: |
---|
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>(); } |
---|
1591 | }; |
---|
1592 | |
---|
1593 | |
---|
1594 | /// \brief A spline with Y transformed into log space. |
---|
1595 | /// \ingroup interpolators |
---|
1596 | /// Most useful for functions looking like y=exp(x) |
---|
1597 | /// |
---|
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> { |
---|
1600 | public: |
---|
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>(); } |
---|
1608 | }; |
---|
1609 | |
---|
1610 | |
---|
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. |
---|
1614 | /// |
---|
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>(); } |
---|
1625 | }; |
---|
1626 | |
---|
1627 | |
---|
1628 | /// \brief A spline with X in reciprocal space and Y transformed in log space. |
---|
1629 | /// \ingroup interpolators |
---|
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. |
---|
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> { |
---|
1635 | public: |
---|
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>(); } |
---|
1643 | }; |
---|
1644 | |
---|
1645 | /// \brief compute sin(x) with its derivatives. |
---|
1646 | /// \ingroup math_functions |
---|
1647 | /// |
---|
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> { |
---|
1650 | public: |
---|
1651 | /// \brief constructor. |
---|
1652 | c2_sin_p() : c2_function<float_type>() {} |
---|
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 |
---|
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; |
---|
1662 | }; |
---|
1663 | |
---|
1664 | /// \brief compute cos(x) with its derivatives. |
---|
1665 | /// \ingroup math_functions |
---|
1666 | /// |
---|
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> { |
---|
1669 | public: |
---|
1670 | /// \brief constructor. |
---|
1671 | c2_cos_p() : c2_sin_p<float_type>() {} |
---|
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 | }; |
---|
1676 | |
---|
1677 | /// \brief compute tan(x) with its derivatives. |
---|
1678 | /// \ingroup math_functions |
---|
1679 | /// |
---|
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> { |
---|
1682 | public: |
---|
1683 | /// \brief constructor. |
---|
1684 | c2_tan_p() : c2_function<float_type>() {} |
---|
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 | }; |
---|
1695 | |
---|
1696 | /// \brief compute log(x) with its derivatives. |
---|
1697 | /// \ingroup math_functions |
---|
1698 | /// |
---|
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> { |
---|
1701 | public: |
---|
1702 | /// \brief constructor. |
---|
1703 | c2_log_p() : c2_function<float_type>() {} |
---|
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 | }; |
---|
1708 | |
---|
1709 | /// \brief compute exp(x) with its derivatives. |
---|
1710 | /// \ingroup math_functions |
---|
1711 | /// |
---|
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> { |
---|
1714 | public: |
---|
1715 | /// \brief constructor. |
---|
1716 | c2_exp_p() : c2_function<float_type>() {} |
---|
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 | }; |
---|
1721 | |
---|
1722 | /// \brief compute sqrt(x) with its derivatives. |
---|
1723 | /// \ingroup math_functions |
---|
1724 | /// |
---|
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> { |
---|
1727 | public: |
---|
1728 | /// \brief constructor. |
---|
1729 | c2_sqrt_p() : c2_function<float_type>() {} |
---|
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 | }; |
---|
1734 | |
---|
1735 | /// \brief compute scale/x with its derivatives. |
---|
1736 | /// \ingroup parametric_functions |
---|
1737 | /// |
---|
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> { |
---|
1740 | public: |
---|
1741 | /// \brief constructor. |
---|
1742 | c2_recip_p(float_type scale) : c2_function<float_type>(), rscale(scale) {} |
---|
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 | }; |
---|
1758 | |
---|
1759 | /// \brief compute x with its derivatives. |
---|
1760 | /// \ingroup math_functions |
---|
1761 | /// |
---|
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> { |
---|
1764 | public: |
---|
1765 | /// \brief constructor. |
---|
1766 | c2_identity_p() : c2_function<float_type>() {} |
---|
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 |
---|
1774 | \ingroup parametric_functions |
---|
1775 | for example, given a c2_function \a f |
---|
1776 | \code |
---|
1777 | c2_function<double> &F=c2_linear<double>(1.2, 2.0, 3.0)(f); |
---|
1778 | \endcode |
---|
1779 | produces a new c2_function F=2.0+3.0*(\a f - 1.2) |
---|
1780 | |
---|
1781 | The factory function c2_factory::linear() creates *new c2_linear_p |
---|
1782 | */ |
---|
1783 | template <typename float_type=double> class c2_linear_p : public c2_function<float_type> { |
---|
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 |
---|
1789 | c2_linear_p(float_type x0, float_type y0, float_type slope) : |
---|
1790 | c2_function<float_type>(), xint(x0), intercept(y0), m(slope) {} |
---|
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: |
---|
1802 | c2_linear_p() {} // do not allow naked construction... it is usually an accident. |
---|
1803 | }; |
---|
1804 | |
---|
1805 | /** |
---|
1806 | \brief create a quadratic mapping of another function |
---|
1807 | \ingroup parametric_functions |
---|
1808 | for example, given a c2_function \a f |
---|
1809 | \code |
---|
1810 | c2_function<double> &F=c2_quadratic<double>(1.2, 2.0, 3.0, 4.0)(f); |
---|
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 | |
---|
1816 | The factory function c2_factory::quadratic() creates *new c2_quadratic_p |
---|
1817 | */ |
---|
1818 | template <typename float_type=double> class c2_quadratic_p : public c2_function<float_type> { |
---|
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 |
---|
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 |
---|
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: |
---|
1839 | c2_quadratic_p() {} // do not allow naked construction... it is usually an accident. |
---|
1840 | }; |
---|
1841 | |
---|
1842 | /** |
---|
1843 | \brief create a power law mapping of another function |
---|
1844 | \ingroup parametric_functions |
---|
1845 | for example, given a c2_function \a f |
---|
1846 | \code |
---|
1847 | c2_power_law_p<double> PLaw(1.2, 2.5); |
---|
1848 | c2_composed_function_p<double> &F=PLaw(f); |
---|
1849 | \endcode |
---|
1850 | produces a new c2_function F=1.2 * f^2.5 |
---|
1851 | |
---|
1852 | The factory function c2_factory::power_law() creates *new c2_power_law_p |
---|
1853 | */ |
---|
1854 | template <typename float_type=double> class c2_power_law_p : public c2_function<float_type> { |
---|
1855 | public: |
---|
1856 | /// \brief Construct the operator |
---|
1857 | /// \param scale the multipler |
---|
1858 | /// \param power the exponent |
---|
1859 | c2_power_law_p(float_type scale, float_type power) : |
---|
1860 | c2_function<float_type>(), a(scale), b(power) {} |
---|
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: |
---|
1871 | c2_power_law_p() {} // do not allow naked construction... it is usually an accident. |
---|
1872 | }; |
---|
1873 | |
---|
1874 | /** |
---|
1875 | \brief create the formal inverse function of another function |
---|
1876 | \ingroup containers |
---|
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. |
---|
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> { |
---|
1893 | public: |
---|
1894 | /// \brief Construct the operator |
---|
1895 | /// \param source the function to be inverted |
---|
1896 | c2_inverse_function_p(const c2_function<float_type> &source); |
---|
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 |
---|
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 | |
---|
1936 | protected: |
---|
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; |
---|
1941 | }; |
---|
1942 | |
---|
1943 | /** |
---|
1944 | \brief |
---|
1945 | An interpolating_function_p which is the cumulative integral of a histogram. |
---|
1946 | \ingroup interpolators |
---|
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 | |
---|
1953 | template <typename float_type=double> class accumulated_histogram : public interpolating_function_p <float_type> { |
---|
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 | /** |
---|
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 |
---|
1983 | */ |
---|
1984 | |
---|
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); |
---|
1989 | |
---|
1990 | /** |
---|
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); |
---|
2010 | |
---|
2011 | /// \brief create a c2_function which smoothly connects two other c2_functions. |
---|
2012 | /// \ingroup parametric_functions |
---|
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. |
---|
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> { |
---|
2027 | public: |
---|
2028 | /// \brief construct the container from two functions |
---|
2029 | /// \param x0 the point at which to match \a f1 and its derivatives |
---|
2030 | /// \param f0 the function on the left side to be connected |
---|
2031 | /// \param x2 the point at which to match \a f2 and its derivatives |
---|
2032 | /// \param f2 the function on the right side to be connected |
---|
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 |
---|
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, |
---|
2037 | bool auto_center, float_type y1); |
---|
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 | |
---|
2066 | /// \brief destructor |
---|
2067 | virtual ~c2_connector_function_p(); |
---|
2068 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw (c2_exception); |
---|
2069 | protected: |
---|
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; |
---|
2077 | }; |
---|
2078 | |
---|
2079 | /// \brief create a c2_function which is a piecewise assembly of other c2_functions. |
---|
2080 | /// \ingroup containers |
---|
2081 | /// The functions must have increasing, non-overlapping domains. Any empty space |
---|
2082 | /// between functions will be filled with a linear interpolation. |
---|
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 | /// |
---|
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 | /// |
---|
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> { |
---|
2099 | public: |
---|
2100 | /// \brief construct the container |
---|
2101 | c2_piecewise_function_p(); |
---|
2102 | /// \brief destructor |
---|
2103 | virtual ~c2_piecewise_function_p(); |
---|
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 |
---|
2112 | void append_function(const c2_function<float_type> &func) throw (c2_exception); |
---|
2113 | protected: |
---|
2114 | std::vector<c2_const_ptr<float_type> > functions; |
---|
2115 | mutable int lastKLow; |
---|
2116 | }; |
---|
2117 | |
---|
2118 | #include "c2_function.icc" |
---|
2119 | |
---|
2120 | #endif |
---|