1 | /** |
---|
2 | * \file |
---|
3 | * \brief Provides the headers for the general c2_function algebra which supports |
---|
4 | * fast, flexible operations on piecewise-twice-differentiable functions |
---|
5 | * |
---|
6 | * \author Created by R. A. Weller and Marcus H. Mendenhall on 7/9/05. |
---|
7 | * \author Copyright 2005 __Vanderbilt University__. All rights reserved. |
---|
8 | * |
---|
9 | * \version c2_function.hh,v 1.53 2007/11/12 13:58:57 marcus Exp |
---|
10 | */ |
---|
11 | |
---|
12 | #ifndef __has_C2Functions_c2_h |
---|
13 | #define __has_C2Functions_c2_h 1 |
---|
14 | |
---|
15 | #include <cmath> |
---|
16 | #include <vector> |
---|
17 | #include <string> |
---|
18 | |
---|
19 | /// \brief the exception class for c2_function operations. |
---|
20 | class c2_exception : public std::exception { |
---|
21 | public: |
---|
22 | /// \brief construct the exception with an error message |
---|
23 | /// \param msgcode the message |
---|
24 | c2_exception(const char msgcode[]) : info(msgcode) { } |
---|
25 | virtual ~c2_exception() throw() { } |
---|
26 | /** Returns a C-style character string describing the general cause |
---|
27 | * of the current error. */ |
---|
28 | virtual const char* what() const throw() { return info.c_str(); } |
---|
29 | private: |
---|
30 | std::string info; |
---|
31 | }; |
---|
32 | |
---|
33 | // put these forward references here, and with a bogus typename to make swig happy. |
---|
34 | template <typename float_type> class c2_composed_function; |
---|
35 | template <typename float_type> class c2_sum; |
---|
36 | template <typename float_type> class c2_diff; |
---|
37 | template <typename float_type> class c2_product; |
---|
38 | template <typename float_type> class c2_ratio; |
---|
39 | |
---|
40 | /** |
---|
41 | \brief the parent class for all c2_functions. |
---|
42 | |
---|
43 | c2_functions know their value, first, and second derivative at almost every point. |
---|
44 | They can be efficiently combined with binary operators, via c2_binary_function, |
---|
45 | composed via c2_composed_function, |
---|
46 | have their roots found via find_root(), |
---|
47 | and be adaptively integrated via partial_integrals() or integral(). |
---|
48 | They also can carry information with them about how to find 'interesting' points on the function. |
---|
49 | This information is set with set_sampling_grid() and extracted with get_sampling_grid(). |
---|
50 | |
---|
51 | Particularly important subclasses are the interpolating functions classes, |
---|
52 | interpolating_function , lin_log_interpolating_function, log_lin_interpolating_function, |
---|
53 | log_log_interpolating_function, and arrhenius_interpolating_function, |
---|
54 | as well as the template functions |
---|
55 | inverse_integrated_density(). |
---|
56 | |
---|
57 | \warning |
---|
58 | The composite flavors of c2_functions (c2_sum, c2_composed_function, c2_binary_function, e.g.) make no effort to manage |
---|
59 | deletion of their component functions. |
---|
60 | These are just container classes, and the user (along with normal automatic variable semantics) |
---|
61 | is responsible for the lifetime of components. |
---|
62 | Inappropriate attention to this can cause massive memory leaks. |
---|
63 | However, in most cases these do exactly what is intended. |
---|
64 | The classes will be left this way since the only other option is to use copy constructors on everything, |
---|
65 | which would make this all very slow. |
---|
66 | |
---|
67 | */ |
---|
68 | template <typename float_type=double> class c2_function { |
---|
69 | public: |
---|
70 | /// \brief get versioning information for the header file |
---|
71 | /// \return the CVS Id string |
---|
72 | const std::string cvs_header_vers() const { return |
---|
73 | "c2_function.hh,v 1.53 2007/11/12 13:58:57 marcus Exp"; |
---|
74 | } |
---|
75 | |
---|
76 | /// \brief get versioning information for the source file |
---|
77 | /// \return the CVS Id string |
---|
78 | const std::string cvs_file_vers() const ; |
---|
79 | |
---|
80 | public: |
---|
81 | /// \brief destructor |
---|
82 | virtual ~c2_function() { if(sampling_grid && !no_overwrite_grid) delete sampling_grid; } |
---|
83 | |
---|
84 | /// \brief get the value and derivatives. |
---|
85 | /// |
---|
86 | /// There is required checking for null pointers on the derivatives, |
---|
87 | /// and most implementations should operate faster if derivatives are not needed. |
---|
88 | /// \param[in] x the point at which to evaluate the function |
---|
89 | /// \param[out] yprime the first derivative (if pointer is non-null) |
---|
90 | /// \param[out] yprime2 the second derivative (if pointer is non-null) |
---|
91 | /// \return the value of the function |
---|
92 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) =0 ; // { return 0; }; |
---|
93 | |
---|
94 | /// \brief evaluate the function in the classic way, ignoring derivatives. |
---|
95 | /// \param x the point at which to evaluate |
---|
96 | /// \return the value of the function |
---|
97 | inline float_type operator () (float_type x) const throw(c2_exception) |
---|
98 | { return value_with_derivatives(x, (float_type *)0, (float_type *)0); } |
---|
99 | |
---|
100 | /// \brief compose this function outside another. |
---|
101 | /// \param inner the inner function |
---|
102 | /// \return the composed function |
---|
103 | c2_composed_function<float_type> & operator ()(const c2_function<float_type> &inner) const |
---|
104 | { return *new c2_composed_function<float_type>((*this), inner); } |
---|
105 | |
---|
106 | /// \brief get the value and derivatives. |
---|
107 | /// |
---|
108 | /// \param[in] x the point at which to evaluate the function |
---|
109 | /// \param[out] yprime the first derivative (if pointer is non-null) |
---|
110 | /// \param[out] yprime2 the second derivative (if pointer is non-null) |
---|
111 | /// \return the value of the function |
---|
112 | inline float_type operator () (float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
113 | { return value_with_derivatives(x, yprime, yprime2); } |
---|
114 | |
---|
115 | /// \brief solve f(x)==value very efficiently, with explicit knowledge of derivatives of the function |
---|
116 | /// |
---|
117 | /// find_root solves by iterated inverse quadratic extrapolation for a solution to f(x)=y. It |
---|
118 | /// includes checks against bad convergence, so it should never be able to fail. Unlike typical |
---|
119 | /// secant method or fancier Brent's method finders, this does not depend in any strong wasy on the |
---|
120 | /// brackets, unless the finder has to resort to successive approximations to close in on a root. |
---|
121 | /// Often, it is possible to make the brackets equal to the domain of the function, if there is |
---|
122 | /// any clue as to where the root lies, as given by the parameter \a start. |
---|
123 | /// \param lower_bracket the lower bound for the search |
---|
124 | /// \param upper_bracket the upper bound for the search. Function sign must be |
---|
125 | /// opposite to that at \a lower_bracket |
---|
126 | /// \param start starting value for the search |
---|
127 | /// \param value the value of the function being sought (solves f(x) = \a value) |
---|
128 | /// \param[out] error If pointer is zero, errors raise exception. Otherwise, returns error here. |
---|
129 | /// \param[out] final_yprime If pointer is not zero, return derivative of function at root |
---|
130 | /// \param[out] final_yprime2 If pointer is not zero, return second derivative of function at root |
---|
131 | /// \return the position of the root. |
---|
132 | float_type find_root(float_type lower_bracket, float_type upper_bracket, float_type start, |
---|
133 | float_type value, int *error=0, |
---|
134 | float_type *final_yprime=0, float_type *final_yprime2=0 ) const throw(c2_exception) ; // solve f(x)=value |
---|
135 | |
---|
136 | /// \brief for points in xgrid, adaptively return Integral[f(x),{x,xgrid[i],xgrid[i+1]}] and return in vector, along with sum |
---|
137 | /// |
---|
138 | /// partial_integrals uses a method with an error O(dx**10) with full information from the derivatives, |
---|
139 | /// and falls back to lower order methods if informed of incomplete derivatives. |
---|
140 | /// It uses exact midpoint splitting of the intervals for recursion, resulting in no recomputation of the function |
---|
141 | /// during recursive descent at previously computed points. |
---|
142 | /// \param xgrid points between which to evaluate definite integrals. |
---|
143 | /// \param partials if non-NULL, a vector in which to receive the partial integrals. |
---|
144 | /// It will automatically be sized apprpropriately, if provided, to contain \a n - 1 elements where \a n is the length of \a xgrid |
---|
145 | /// \param abs_tol the absolute error bound for each segment |
---|
146 | /// \param rel_tol the fractional error bound for each segment. |
---|
147 | /// If the error is smaller than either the relative or absolute tolerance, the integration step is finished. |
---|
148 | /// \param derivs number of derivatives to trust, which sets the order of the integrator. The order |
---|
149 | /// is 3*\a derivs + 4. \a derivs can be 0, 1, or 2. |
---|
150 | /// \param adapt if true, use recursive adaptation, otherwise do simple evaluation on the grid provided |
---|
151 | /// with no error checking. |
---|
152 | /// \param extrapolate if true, use simple Richardson extrapolation on the final 2 steps to reduce the error. |
---|
153 | /// \return sum of partial integrals, whcih is the definite integral from the first value in \a xgrid to the last. |
---|
154 | float_type partial_integrals(std::vector<float_type> xgrid, std::vector<float_type> *partials = 0, |
---|
155 | float_type abs_tol=1e-12, float_type rel_tol=1e-12, int derivs=2, bool adapt=true, bool extrapolate=true) const; |
---|
156 | |
---|
157 | /// \brief a fully-automated integrator which uses the information provided by the get_sampling_grid() function |
---|
158 | /// to figure out what to do. |
---|
159 | /// |
---|
160 | /// It returns the integral of the function over the domain requested |
---|
161 | /// with error tolerances as specified. It is just a front-end to partial_integrals() |
---|
162 | /// |
---|
163 | /// \param xmin lower bound of the domain for integration |
---|
164 | /// \param xmax upper bound of the domain for integration |
---|
165 | /// \param partials if non-NULL, a vector in which to receive the partial integrals. |
---|
166 | /// It will automatically be sized apprpropriately, if provided, to contain \a n - 1 elements where \a n is the length of \a xgrid |
---|
167 | /// \param abs_tol the absolute error bound for each segment |
---|
168 | /// \param rel_tol the fractional error bound for each segment. |
---|
169 | /// If the error is smaller than either the relative or absolute tolerance, the integration step is finished. |
---|
170 | /// \param derivs number of derivatives to trust, which sets the order of the integrator. The order |
---|
171 | /// is 3*\a derivs + 4. \a derivs can be 0, 1, or 2. |
---|
172 | /// \param adapt if true, use recursive adaptation, otherwise do simple evaluation on the grid provided |
---|
173 | /// with no error checking. |
---|
174 | /// \param extrapolate if true, use simple Richardson extrapolation on the final 2 steps to reduce the error. |
---|
175 | /// \return sum of partial integrals, whcih is the definite integral from the first value in \a xgrid to the last. |
---|
176 | float_type integral(float_type xmin, float_type xmax, std::vector<float_type> *partials = 0, |
---|
177 | float_type abs_tol=1e-12, float_type rel_tol=1e-12, int derivs=2, bool adapt=true, bool extrapolate=true) const; |
---|
178 | |
---|
179 | /// \brief return the lower bound of the domain for this function as set by set_domain() |
---|
180 | inline float_type xmin() const { return fXMin; } |
---|
181 | /// \brief return the upper bound of the domain for this function as set by set_domain() |
---|
182 | inline float_type xmax() const { return fXMax; } |
---|
183 | /// \brief set the domain for this function. |
---|
184 | void set_domain(float_type xmin, float_type xmax) { fXMin=xmin; fXMax=xmax; } |
---|
185 | |
---|
186 | /// \brief this is a counter owned by the function but which can be used to monitor efficiency of algorithms. |
---|
187 | /// |
---|
188 | /// It is not maintained automatically in general! |
---|
189 | /// The adaptive integrator and root finder do clear it at the start and update it for performance checking. |
---|
190 | /// \return number of evaluations logged since last reset. |
---|
191 | volatile int get_evaluations() const { return evaluations; } |
---|
192 | /// \brief reset the counter |
---|
193 | void reset_evaluations() const { evaluations=0; } // evaluations are 'invisible' to constant |
---|
194 | /// \brief count evaluations |
---|
195 | inline void increment_evaluations() const { evaluations++; } |
---|
196 | |
---|
197 | /// \brief check that a vector is monotonic, throw an exception if not, and return a flag if it is reversed |
---|
198 | /// |
---|
199 | /// \param data a vector of data points which are expected to be monotonic. |
---|
200 | /// \param message an informative string to include in an exception if this throws c2_exception |
---|
201 | /// \return true if in decreasing order, false if increasing |
---|
202 | bool check_monotonicity(const std::vector<float_type> &data, const char message[]) throw(c2_exception); |
---|
203 | |
---|
204 | /// \brief establish a grid of 'interesting' points on the function. |
---|
205 | /// |
---|
206 | /// The sampling grid describes a reasonable initial set of points to look at the function. |
---|
207 | /// this should generally be set at a scale which is quite coarse, and sufficient for initializing |
---|
208 | /// adaptive integration or possibly root bracketing. For sampling a function to build a new interpolating |
---|
209 | /// function, one may want to refine this for accuracy. However, interpolating_functions themselves |
---|
210 | /// return their original X grid by default, so refining the grid in this case might be a bad idea. |
---|
211 | /// \param grid a vector of abscissas. The contents is copied into an internal vector, so the \a grid can be discarded after passingin. |
---|
212 | virtual void set_sampling_grid(const std::vector<float_type> &grid) throw(c2_exception); |
---|
213 | |
---|
214 | /// \brief return the grid of 'interesting' points along this function which lie in the region requested |
---|
215 | /// |
---|
216 | /// if a sampling grid is defined, work from there, otherwise return vector of (xmin, xmax) |
---|
217 | /// \param xmin the lower bound for which the function is to be sampled |
---|
218 | /// \param xmax the upper bound for which the function is to be sampled |
---|
219 | /// \return a new vector containing the samplng grid. |
---|
220 | virtual std::vector<float_type> &get_sampling_grid(float_type xmin, float_type xmax) const ; |
---|
221 | |
---|
222 | /// \brief clean up endpoints on a grid of points |
---|
223 | /// |
---|
224 | /// \param[in,out] result the sampling grid with excessively closely space endpoints removed. |
---|
225 | /// The grid is modified in place. |
---|
226 | void preen_sampling_grid(std::vector<float_type> *result) const; |
---|
227 | /// \brief refine a grid by splitting each interval into more intervals |
---|
228 | /// \param grid the grid to refine |
---|
229 | /// \param refinement the number of new steps for each old step |
---|
230 | /// \return a new sampling grid with more points. |
---|
231 | std::vector<float_type> & refine_sampling_grid(const std::vector<float_type> &grid, size_t refinement) const; |
---|
232 | |
---|
233 | /// \brief create a new c2_function from this one which is normalized on the interval |
---|
234 | /// \param xmin lower bound of the domain for integration |
---|
235 | /// \param xmax upper bound of the domain for integration |
---|
236 | /// \param norm the desired integral for the function over the region |
---|
237 | /// \return a new c2_function with the desired \a norm. |
---|
238 | c2_function<float_type> &normalized_function(float_type xmin, float_type xmax, float_type norm=1.0); |
---|
239 | /// \brief create a new c2_function from this one which is square-normalized on the interval |
---|
240 | /// \param xmin lower bound of the domain for integration |
---|
241 | /// \param xmax upper bound of the domain for integration |
---|
242 | /// \param norm the desired integral for the function over the region |
---|
243 | /// \return a new c2_function with the desired \a norm. |
---|
244 | c2_function<float_type> &square_normalized_function(float_type xmin, float_type xmax, float_type norm=1.0); |
---|
245 | /// \brief create a new c2_function from this one which is square-normalized with the provided \a weight on the interval |
---|
246 | /// \param xmin lower bound of the domain for integration |
---|
247 | /// \param xmax upper bound of the domain for integration |
---|
248 | /// \param weight a c2_function providing the weight |
---|
249 | /// \param norm the desired integral for the function over the region |
---|
250 | /// \return a new c2_function with the desired \a norm. |
---|
251 | c2_function<float_type> &square_normalized_function(float_type xmin, float_type xmax, const c2_function<float_type> &weight, float_type norm=1.0); |
---|
252 | |
---|
253 | /// \brief factory function to create a c2_sum from an regular algebraic expression. |
---|
254 | /// \note |
---|
255 | /// be very wary of ownership issues if this is used in a complex expression. |
---|
256 | c2_sum<float_type> &operator + (const c2_function<float_type> &rhs) { return *new c2_sum<float_type>(*this, rhs); } |
---|
257 | /// \brief factory function to create a c2_diff from an regular algebraic expression. |
---|
258 | /// \note |
---|
259 | /// be very wary of ownership issues if this is used in a complex expression. |
---|
260 | c2_diff<float_type> &operator - (const c2_function<float_type> &rhs) { return *new c2_diff<float_type>(*this, rhs); } |
---|
261 | /// \brief factory function to create a c2_product from an regular algebraic expression. |
---|
262 | /// \note |
---|
263 | /// be very wary of ownership issues if this is used in a complex expression. |
---|
264 | c2_product<float_type> &operator * (const c2_function<float_type> &rhs) { return *new c2_product<float_type>(*this, rhs); } |
---|
265 | /// \brief factory function to create a c2_ratio from an regular algebraic expression. |
---|
266 | /// \note |
---|
267 | /// be very wary of ownership issues if this is used in a complex expression. |
---|
268 | c2_ratio<float_type> &operator / (const c2_function<float_type> &rhs) { return *new c2_ratio<float_type>(*this, rhs); } |
---|
269 | |
---|
270 | |
---|
271 | |
---|
272 | std::vector<float_type> *sampling_grid; |
---|
273 | bool no_overwrite_grid; |
---|
274 | |
---|
275 | protected: |
---|
276 | c2_function(const c2_function<float_type> &src) : sampling_grid(0), |
---|
277 | no_overwrite_grid(false), |
---|
278 | fXMin(src.fXMin), fXMax(src.fXMax), rootInitialized(false) |
---|
279 | {} // copy constructor only copies domain, and is only for internal use |
---|
280 | c2_function() : |
---|
281 | sampling_grid(0), no_overwrite_grid(0), |
---|
282 | fXMin(-std::numeric_limits<float_type>::max()), |
---|
283 | fXMax(std::numeric_limits<float_type>::max()), |
---|
284 | rootInitialized(false) |
---|
285 | {} // prevent accidental naked construction (impossible any since this has pure virtual methods) |
---|
286 | |
---|
287 | // this should only be called very early on, by a constructor, before anyone else |
---|
288 | // sets a sampling grid, or it will leak memory |
---|
289 | virtual void set_sampling_grid_pointer(std::vector<float_type> &grid) |
---|
290 | { |
---|
291 | if (sampling_grid && !no_overwrite_grid) delete sampling_grid; // grid was ours, lose it. |
---|
292 | sampling_grid=&grid; no_overwrite_grid=1; |
---|
293 | } |
---|
294 | |
---|
295 | float_type fXMin, fXMax; |
---|
296 | mutable int evaluations; |
---|
297 | |
---|
298 | private: |
---|
299 | /// \brief structure used for recursion in adaptive integrator. |
---|
300 | /// |
---|
301 | /// Contains all the information for the function at one point. |
---|
302 | struct c2_integrate_fblock { float_type x, y, yp, ypp; }; |
---|
303 | /// \brief structure used to pass information recursively. |
---|
304 | /// |
---|
305 | /// the \a abs_tol is scaled by a factor of two at each division. |
---|
306 | /// Everything else is just passed down. |
---|
307 | struct c2_integrate_recur { struct c2_integrate_fblock *f0, *f1, *f2; |
---|
308 | float_type abs_tol, rel_tol, *lr, eps_scale, extrap_coef, extrap2; |
---|
309 | int depth, derivs; |
---|
310 | bool adapt, extrapolate; |
---|
311 | }; |
---|
312 | |
---|
313 | /// \brief Carry out the recursive subdivision and integration. |
---|
314 | /// |
---|
315 | /// This passes information recursively through the \a recur block pointer |
---|
316 | /// to allow very efficient recursion. |
---|
317 | /// \param rb a pointer to the recur struct. |
---|
318 | float_type integrate_step(struct c2_integrate_recur &rb) const; |
---|
319 | |
---|
320 | /// these carry a memory of the last root bracketing, |
---|
321 | /// to avoid the necessity of evaluating the function on the brackets every time |
---|
322 | /// if the brackets have not been changed. |
---|
323 | mutable float_type lastRootLowerX, lastRootUpperX, lastRootLowerY, lastRootUpperY; |
---|
324 | mutable int rootInitialized; |
---|
325 | |
---|
326 | }; |
---|
327 | |
---|
328 | /// \brief a container into which any other c2_function can be dropped, to allow expressions |
---|
329 | /// with replacable components. |
---|
330 | /// |
---|
331 | ///It is useful for plugging different InterpolatingFunctions into a c2_function expression. |
---|
332 | ///It saves a lot of effort in other places with casting away const declarations. |
---|
333 | /// |
---|
334 | /// It is also useful as a wrapper for a function if it is necessary to have a copy of a function |
---|
335 | /// which has a different domain or sampling grid than the parent function. This can be |
---|
336 | /// be used, for example, to patch badly-behaved functions with c2_piecewise_function by |
---|
337 | /// taking the parent function, creating two plugins of it with domains on each side of the |
---|
338 | /// nasty bit, and then inserting a nice function in the hole. |
---|
339 | template <typename float_type=double> class c2_plugin_function : public c2_function<float_type> { |
---|
340 | public: |
---|
341 | /// \brief construct the container with no function |
---|
342 | c2_plugin_function() : c2_function<float_type>(), func(0), owns(false) {} |
---|
343 | /// \brief construct the container with a pre-defined function |
---|
344 | c2_plugin_function(const c2_function<float_type> &f) : c2_function<float_type>(), func(0), owns(false) { set_function(f); } |
---|
345 | /// \brief fill the container with a new function |
---|
346 | void set_function(const c2_function<float_type> &f) |
---|
347 | { |
---|
348 | if(owns && func) delete func; |
---|
349 | func=&f; set_domain(f.xmin(), f.xmax()); |
---|
350 | set_sampling_grid_pointer(*f.sampling_grid); |
---|
351 | owns=false; |
---|
352 | } |
---|
353 | /// \copydoc c2_function::value_with_derivatives |
---|
354 | /// Uses the internal function pointer set by set_function(). |
---|
355 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
356 | { |
---|
357 | if(!func) throw c2_exception("c2_plugin_function<float_type> called uninitialized"); |
---|
358 | return this->func->value_with_derivatives(x, yprime, yprime2); |
---|
359 | } |
---|
360 | /// \brief clear the function |
---|
361 | /// |
---|
362 | /// Any attempt to use this c2_plugin_function throws an exception if the saved function is cleared. |
---|
363 | void unset_function(void) { if(owns && func) delete func; func=0; owns=false; } |
---|
364 | /// \brief destructor |
---|
365 | ~c2_plugin_function() { if(owns && func) delete func; } |
---|
366 | /// \brief tell us we should delete the function at destruction. NOT sticky when function is reset |
---|
367 | void set_ownership() { this->owns=true; } |
---|
368 | |
---|
369 | protected: |
---|
370 | const c2_function<float_type> *func; |
---|
371 | bool owns; |
---|
372 | |
---|
373 | }; |
---|
374 | |
---|
375 | /// \brief Provides support for c2_function objects which are constructed from two other c2_function |
---|
376 | /// objects. |
---|
377 | /// |
---|
378 | /// It provides a very primitive ownership concept, so that the creator can tag a function |
---|
379 | /// as being owned by this function, so that when this function is deleted, the owned function will be deleted, too. |
---|
380 | /// This allows a piece of code to create various c2_function objects, combine them with binary operators, |
---|
381 | /// appropriately mark wich ones have no other possible owners, and return the final function with |
---|
382 | /// reasonable faith that everything will get cleaned up when the final function is deleted. Note that |
---|
383 | /// none of this marking is automatic, to keep this class very lightweight. |
---|
384 | template <typename float_type=double> class c2_binary_function : public c2_function<float_type> { |
---|
385 | public: |
---|
386 | |
---|
387 | |
---|
388 | /// \brief function to manage the binary operation, used by c2_binary_function::value_with_derivatives() |
---|
389 | /// |
---|
390 | /// Normally not used alone, but can be used to combine functions in other contexts. |
---|
391 | /// See interpolating_function::binary_operator() for an example. |
---|
392 | /// \param left the function on the left of the binary operator or outside the composition |
---|
393 | /// \param right the function to the right of the operator or inside the composition |
---|
394 | /// \param[in] x the point at which to evaluate the function |
---|
395 | /// \param[out] yprime the first derivative (if pointer is non-null) |
---|
396 | /// \param[out] yprime2 the second derivative (if pointer is non-null) |
---|
397 | /// \return the value of the function |
---|
398 | virtual float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, |
---|
399 | float_type x, float_type *yprime, float_type *yprime2) const =0; |
---|
400 | |
---|
401 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw (c2_exception) |
---|
402 | { |
---|
403 | return this->combine(this->Left, this->Right, x, yprime, yprime2); |
---|
404 | } |
---|
405 | |
---|
406 | /// \brief allow c2_binary_function to remember ownership of contained functions for automatic cleanup |
---|
407 | /// |
---|
408 | /// upon destruction, this will cause disposal of the left member function |
---|
409 | void set_left_ownership(void) { leftown=true; } |
---|
410 | /// \brief allow c2_binary_function to remember ownership of contained functions for automatic cleanup |
---|
411 | /// |
---|
412 | /// upon destruction, this will cause disposal of the right member function |
---|
413 | void set_right_ownership(void) { rightown=true; } |
---|
414 | /// \brief allow c2_binary_function to remember ownership of contained functions for automatic cleanup |
---|
415 | /// |
---|
416 | /// upon destruction, this will cause disposal of both member functions |
---|
417 | void set_ownership(void) { leftown=rightown=true; } |
---|
418 | /// \brief destructor executes disposal of member functions if flagged |
---|
419 | /// |
---|
420 | /// depends on judicious use of set_ownership(), set_right_ownership(), or set_left_ownership() |
---|
421 | virtual ~c2_binary_function() { |
---|
422 | if(leftown) delete &Left; |
---|
423 | if(rightown) delete &Right; |
---|
424 | } |
---|
425 | |
---|
426 | protected: |
---|
427 | /// \brief construct the binary function |
---|
428 | /// \param left the c2_function to be used in the left side of the binary relation |
---|
429 | /// \param right the c2_function to be used in the right side of the binary relation |
---|
430 | c2_binary_function( const c2_function<float_type> &left, const c2_function<float_type> &right) : |
---|
431 | c2_function<float_type>(), |
---|
432 | Left(left), Right(right), leftown(false), rightown(false) |
---|
433 | { |
---|
434 | set_domain( |
---|
435 | (left.xmin() > right.xmin()) ? left.xmin() : right.xmin(), |
---|
436 | (left.xmax() < right.xmax()) ? left.xmax() : right.xmax() |
---|
437 | ); |
---|
438 | } |
---|
439 | |
---|
440 | /// \brief construct a 'stub' c2_binary_function, which provides access to the combine() function |
---|
441 | /// \note Do not evaluate a 'stub' ever. It is only used so that combine() can be called |
---|
442 | c2_binary_function() : c2_function<float_type>(), |
---|
443 | Left(*((c2_function<float_type> *)0)), Right(*((c2_function<float_type> *)0)) {} |
---|
444 | |
---|
445 | const c2_function<float_type> &Left, &Right; |
---|
446 | bool leftown, rightown; |
---|
447 | |
---|
448 | }; |
---|
449 | |
---|
450 | /// \brief Create a very lightweight method to return a scalar multiple of another function. |
---|
451 | /// |
---|
452 | template <typename float_type=double> class c2_scaled_function : public c2_plugin_function<float_type> { |
---|
453 | public: |
---|
454 | /// \brief construct the function with its scale factor. |
---|
455 | /// |
---|
456 | /// \param outer the function to be scaled |
---|
457 | /// \param scale the multiplicative scale factor |
---|
458 | c2_scaled_function(const c2_function<float_type> &outer, float_type scale) : |
---|
459 | c2_plugin_function<float_type>(outer), yscale(scale) { } |
---|
460 | |
---|
461 | /// \brief set a new scale factor |
---|
462 | /// \param scale the new factor |
---|
463 | void reset(float_type scale) { yscale=scale; } |
---|
464 | |
---|
465 | /// \copydoc c2_function::value_with_derivatives |
---|
466 | /// |
---|
467 | /// provide our own value_with_derivatives which bypasses the combiner for quicker operation |
---|
468 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw (c2_exception) |
---|
469 | { |
---|
470 | float_type y=this->func->value_with_derivatives(x, yprime, yprime2); |
---|
471 | if(yprime) (*yprime)*=yscale; |
---|
472 | if(yprime2) (*yprime2)*=yscale; |
---|
473 | return y*yscale; |
---|
474 | } |
---|
475 | |
---|
476 | protected: |
---|
477 | c2_scaled_function<float_type>() {} // hide default constructor, since its use is almost always an error. |
---|
478 | float_type yscale; |
---|
479 | }; |
---|
480 | |
---|
481 | /// \brief A container into which any other c2_function can be dropped. |
---|
482 | /// |
---|
483 | /// It allows a function to be pre-evaluated at a point, and used at multiple places in an expression |
---|
484 | /// efficiently. If it is re-evaluated at the previous point, it returns the remembered values; |
---|
485 | /// otherwise, it re-evauates the function at the new point. |
---|
486 | /// |
---|
487 | template <typename float_type=double> class c2_cached_function : public c2_plugin_function<float_type> { |
---|
488 | public: |
---|
489 | /// \brief construct the container |
---|
490 | /// |
---|
491 | /// \param f the function to be cached |
---|
492 | c2_cached_function(const c2_function<float_type> &f) : c2_plugin_function<float_type>(f), init(false) {} |
---|
493 | /// \copydoc c2_function::value_with_derivatives |
---|
494 | /// |
---|
495 | /// Checks to see if the function is being re-evaluated at the previous point, and |
---|
496 | /// returns remembered values if so. |
---|
497 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
498 | { |
---|
499 | if(!init || x != x0) { |
---|
500 | y=this->func->value_with_derivatives(x, &yp, &ypp); |
---|
501 | x0=x; |
---|
502 | init=true; |
---|
503 | } |
---|
504 | if(yprime) *yprime=yp; |
---|
505 | if(yprime2) *yprime2=ypp; |
---|
506 | return y; |
---|
507 | } |
---|
508 | |
---|
509 | protected: |
---|
510 | c2_cached_function() {} // hide default constructor, since its use is almost always an error. |
---|
511 | mutable bool init; |
---|
512 | mutable float_type x0, y, yp, ypp; |
---|
513 | |
---|
514 | }; |
---|
515 | |
---|
516 | /// \brief Provides function composition (nesting) |
---|
517 | /// |
---|
518 | /// This allows evaluation of \a f(g(x)) where \a f and \a g are c2_function objects. |
---|
519 | /// |
---|
520 | /// \note See c2_binary_function for discussion of ownership. |
---|
521 | template <typename float_type=double> class c2_composed_function : public c2_binary_function<float_type> { |
---|
522 | public: |
---|
523 | |
---|
524 | /// \brief construct \a outer( \a inner (x)) |
---|
525 | /// \note See c2_binary_function for discussion of ownership. |
---|
526 | /// \param outer the outer function |
---|
527 | /// \param inner the inner function |
---|
528 | c2_composed_function(const c2_function<float_type> &outer, const c2_function<float_type> &inner) : c2_binary_function<float_type>(outer, inner) {} |
---|
529 | /// \brief Create a stub just for the combiner to avoid statics. |
---|
530 | c2_composed_function() : c2_binary_function<float_type>() {} ; |
---|
531 | |
---|
532 | virtual float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, |
---|
533 | float_type x, float_type *yprime, float_type *yprime2) const |
---|
534 | { |
---|
535 | float_type y0, yp0, ypp0, y1, yp1, ypp1; |
---|
536 | float_type *yp0p, *ypp0p, *yp1p, *ypp1p; |
---|
537 | if(yprime || yprime2) { |
---|
538 | yp0p=&yp0; ypp0p=&ypp0; yp1p=&yp1; ypp1p=&ypp1; |
---|
539 | } else { |
---|
540 | yp0p=ypp0p=yp1p=ypp1p=0; |
---|
541 | } |
---|
542 | |
---|
543 | y0=right.value_with_derivatives(x, yp0p, ypp0p); |
---|
544 | y1=left.value_with_derivatives(y0, yp1p, ypp1p); |
---|
545 | if(yprime) *yprime=yp1*yp0; |
---|
546 | if(yprime2) *yprime2=ypp0*yp1+yp0*yp0*ypp1; |
---|
547 | return y1; |
---|
548 | } |
---|
549 | }; |
---|
550 | |
---|
551 | /// \brief create a c2_function which is the sum of two other c2_function objects. |
---|
552 | /// |
---|
553 | /// \note See c2_binary_function for discussion of ownership. |
---|
554 | template <typename float_type=double> class c2_sum : public c2_binary_function<float_type> { |
---|
555 | public: |
---|
556 | /// \brief construct \a left + \a right |
---|
557 | /// \note See c2_binary_function for discussion of ownership. |
---|
558 | /// \param left the left function |
---|
559 | /// \param right the right function |
---|
560 | c2_sum(const c2_function<float_type> &left, const c2_function<float_type> &right) : c2_binary_function<float_type>(left, right) {} |
---|
561 | /// \brief Create a stub just for the combiner to avoid statics. |
---|
562 | c2_sum() : c2_binary_function<float_type>() {} ; // create a stub just for the combiner to avoid statics |
---|
563 | |
---|
564 | // function to do derivative arithmetic for sums |
---|
565 | virtual float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, |
---|
566 | float_type x, float_type *yprime, float_type *yprime2) const |
---|
567 | { |
---|
568 | float_type y0, yp0, ypp0, y1, yp1, ypp1; |
---|
569 | float_type *yp0p, *ypp0p, *yp1p, *ypp1p; |
---|
570 | if(yprime || yprime2) { |
---|
571 | yp0p=&yp0; ypp0p=& ypp0; yp1p=&yp1; ypp1p=&ypp1; |
---|
572 | } else { |
---|
573 | yp0p=ypp0p=yp1p=ypp1p=0; |
---|
574 | } |
---|
575 | y0=left.value_with_derivatives(x, yp0p, ypp0p); |
---|
576 | y1=right.value_with_derivatives(x, yp1p, ypp1p); |
---|
577 | if(yprime) *yprime=yp0+yp1; |
---|
578 | if(yprime2) *yprime2=ypp0+ypp1; |
---|
579 | return y0+y1; |
---|
580 | } |
---|
581 | }; |
---|
582 | |
---|
583 | |
---|
584 | /// \brief create a c2_function which is the difference of two other c2_functions. |
---|
585 | /// |
---|
586 | /// \note See c2_binary_function for discussion of ownership. |
---|
587 | template <typename float_type=double> class c2_diff : public c2_binary_function<float_type> { |
---|
588 | public: |
---|
589 | /// \brief construct \a left - \a right |
---|
590 | /// \note See c2_binary_function for discussion of ownership. |
---|
591 | /// \param left the left function |
---|
592 | /// \param right the right function |
---|
593 | c2_diff(const c2_function<float_type> &left, const c2_function<float_type> &right) : c2_binary_function<float_type>(left, right) {} |
---|
594 | /// \brief Create a stub just for the combiner to avoid statics. |
---|
595 | c2_diff() : c2_binary_function<float_type>() {} ; // create a stub just for the combiner to avoid statics |
---|
596 | |
---|
597 | // function to do derivative arithmetic for diffs |
---|
598 | virtual float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, |
---|
599 | float_type x, float_type *yprime, float_type *yprime2) const |
---|
600 | { |
---|
601 | float_type y0, yp0, ypp0, y1, yp1, ypp1; |
---|
602 | float_type *yp0p, *ypp0p, *yp1p, *ypp1p; |
---|
603 | if(yprime || yprime2) { |
---|
604 | yp0p=&yp0; ypp0p=&ypp0; yp1p=&yp1; ypp1p=&ypp1; |
---|
605 | } else { |
---|
606 | yp0p=ypp0p=yp1p=ypp1p=0; |
---|
607 | } |
---|
608 | y0=left.value_with_derivatives(x, yp0p, ypp0p); |
---|
609 | y1=right.value_with_derivatives(x, yp1p, ypp1p); |
---|
610 | if(yprime) *yprime=yp0-yp1; |
---|
611 | if(yprime2) *yprime2=ypp0-ypp1; |
---|
612 | return y0-y1; |
---|
613 | } |
---|
614 | }; |
---|
615 | |
---|
616 | |
---|
617 | /// \brief create a c2_function which is the product of two other c2_functions. |
---|
618 | /// |
---|
619 | /// \note See c2_binary_function for discussion of ownership. |
---|
620 | template <typename float_type=double> class c2_product : public c2_binary_function<float_type> { |
---|
621 | public: |
---|
622 | /// \brief construct \a left * \a right |
---|
623 | /// \note See c2_binary_function for discussion of ownership. |
---|
624 | /// \param left the left function |
---|
625 | /// \param right the right function |
---|
626 | c2_product(const c2_function<float_type> &left, const c2_function<float_type> &right) : c2_binary_function<float_type>(left, right) {} |
---|
627 | /// \brief Create a stub just for the combiner to avoid statics. |
---|
628 | c2_product() : c2_binary_function<float_type>() {} ; // create a stub just for the combiner to avoid statics |
---|
629 | |
---|
630 | virtual float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, |
---|
631 | float_type x, float_type *yprime, float_type *yprime2) const |
---|
632 | { |
---|
633 | float_type y0, yp0, ypp0, y1, yp1, ypp1; |
---|
634 | float_type *yp0p, *ypp0p, *yp1p, *ypp1p; |
---|
635 | if(yprime || yprime2) { |
---|
636 | yp0p=&yp0; ypp0p=&ypp0; yp1p=&yp1; ypp1p=&ypp1; |
---|
637 | } else { |
---|
638 | yp0p=ypp0p=yp1p=ypp1p=0; |
---|
639 | } |
---|
640 | y0=left.value_with_derivatives(x, yp0p, ypp0p); |
---|
641 | y1=right.value_with_derivatives(x, yp1p, ypp1p); |
---|
642 | if(yprime) *yprime=y1*yp0+y0*yp1; |
---|
643 | if(yprime2) *yprime2=ypp0*y1+2.0*yp0*yp1+ypp1*y0; |
---|
644 | return y0*y1; |
---|
645 | } |
---|
646 | }; |
---|
647 | |
---|
648 | |
---|
649 | /// \brief create a c2_function which is the ratio of two other c2_functions. |
---|
650 | /// |
---|
651 | /// \note See c2_binary_function for discussion of ownership. |
---|
652 | template <typename float_type=double> class c2_ratio : public c2_binary_function<float_type> { |
---|
653 | public: |
---|
654 | /// \brief construct \a left / \a right |
---|
655 | /// \note See c2_binary_function for discussion of ownership. |
---|
656 | /// \param left the left function |
---|
657 | /// \param right the right function |
---|
658 | c2_ratio(const c2_function<float_type> &left, const c2_function<float_type> &right) : c2_binary_function<float_type>(left, right) {} |
---|
659 | /// \brief Create a stub just for the combiner to avoid statics. |
---|
660 | c2_ratio() : c2_binary_function<float_type>() {} ; // create a stub just for the combiner to avoid statics |
---|
661 | |
---|
662 | virtual float_type combine(const c2_function<float_type> &left, const c2_function<float_type> &right, |
---|
663 | float_type x, float_type *yprime, float_type *yprime2) const |
---|
664 | { |
---|
665 | float_type y0, yp0, ypp0, y1, yp1, ypp1; |
---|
666 | float_type *yp0p, *ypp0p, *yp1p, *ypp1p; |
---|
667 | if(yprime || yprime2) { |
---|
668 | yp0p=&yp0; ypp0p=&ypp0; yp1p=&yp1; ypp1p=&ypp1; |
---|
669 | } else { |
---|
670 | yp0p=ypp0p=yp1p=ypp1p=0; |
---|
671 | } |
---|
672 | y0=left.value_with_derivatives(x, yp0p, ypp0p); |
---|
673 | y1=right.value_with_derivatives(x, yp1p, ypp1p); |
---|
674 | if(yprime) *yprime=(yp0*y1-y0*yp1)/(y1*y1); // first deriv of ratio |
---|
675 | if(yprime2) *yprime2=(y1*y1*ypp0+y0*(2*yp1*yp1-y1*ypp1)-2*y1*yp0*yp1)/(y1*y1*y1); // second deriv of ratio |
---|
676 | return y0/y1; |
---|
677 | } |
---|
678 | |
---|
679 | }; |
---|
680 | |
---|
681 | /// \brief a c2_function which is constant : can do interpolating_function f2=f1 + c2_constant(11.3) |
---|
682 | template <typename float_type> class c2_constant : public c2_function<float_type> { |
---|
683 | public: |
---|
684 | c2_constant(float_type x=0.0) : c2_function<float_type>(), value(x) {} |
---|
685 | void reset(float_type val) { value=val; } |
---|
686 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
687 | { if(yprime) *yprime=0; if(yprime2) *yprime2=0; return value; } |
---|
688 | |
---|
689 | private: |
---|
690 | float_type value; |
---|
691 | }; |
---|
692 | |
---|
693 | /** |
---|
694 | \brief create a cubic spline interpolation of a set of (x,y) pairs |
---|
695 | |
---|
696 | This is one of the main reasons for c2_function objects to exist. |
---|
697 | |
---|
698 | It provides support for cubic spline interpolation of data provides from tables of \a x, \a y pairs. |
---|
699 | It supports automatic, transparent linearization of the data before storing in its tables (through |
---|
700 | subclasses such as |
---|
701 | log_lin_interpolating_function, lin_log_interpolating_function, and |
---|
702 | log_log_interpolating_function) to permit very high accuracy representations of data which have a suitable |
---|
703 | structure. It provides utility functions LinearInterpolatingGrid() and LogLogInterpolatingGrid() |
---|
704 | to create grids for mapping other functions onto a arithmetic or geometric grid. |
---|
705 | |
---|
706 | In its simplest form, an untransformed cubic spline of a data set, using natural boundary conditions |
---|
707 | (vanishing second derivative), is created as: \n |
---|
708 | \code |
---|
709 | std::vector<double> xvals(10), yvals(10); |
---|
710 | // < fill in xvals and yvals > |
---|
711 | interpolating_function<double> myfunc(xvals, yvals); |
---|
712 | // and it can be evaluated at a point for its value only by: |
---|
713 | double y=myfunc(x); |
---|
714 | // or it can be evaluated with its derivatives by |
---|
715 | double yprime, yprime2; |
---|
716 | double y=myfunc(x,&yprime, &yprime2); |
---|
717 | \endcode |
---|
718 | */ |
---|
719 | |
---|
720 | template <typename float_type=double> class interpolating_function : public c2_function<float_type> { |
---|
721 | public: |
---|
722 | /// \brief create the most general interpolating_function which defaults to linear-linear space |
---|
723 | /// |
---|
724 | /// lots to say here, but see Numerical Recipes for a discussion of cubic splines. |
---|
725 | /// \param x the list of abscissas. Must be either strictly increasing or strictly decreasing. |
---|
726 | /// Strictly increasing is preferred, as less memory is used since a copy is not required for the sampling grid. |
---|
727 | /// \param f the list of function values. |
---|
728 | /// \param lowerSlopeNatural if true, set y''(first point)=0, otherwise compute it from \a lowerSope |
---|
729 | /// \param lowerSlope derivative of the function at the lower bound, used only if \a lowerSlopeNatural is false |
---|
730 | /// \param upperSlopeNatural if true, set y''(last point)=0, otherwise compute it from \a upperSope |
---|
731 | /// \param upperSlope derivative of the function at the upper bound, used only if \a upperSlopeNatural is false |
---|
732 | /// \param inputXConversion a function (not a c2_function) which converts \a x into the internal space. \n |
---|
733 | /// If this is NULL, use linear space and ignore inputXConversionPrime, inputXConversionDPrime |
---|
734 | /// \param inputYConversion a function (not a c2_function) which converts \a y into the internal space. \n |
---|
735 | /// If this is NULL, use linear space and ignore inputYConversionPrime, inputYConversionDPrime, outputYConversion |
---|
736 | /// \param outputYConversion a function (not a c2_function) which converts \a y out of the internal space |
---|
737 | /// \param inputXConversionPrime the derivative of \a inputXConversion |
---|
738 | /// \param inputYConversionPrime the derivative of \a inputYConversion |
---|
739 | /// \param inputXConversionDPrime the second derivative of \a inputXConversion |
---|
740 | /// \param inputYConversionDPrime the second derivative of \a inputYConversion |
---|
741 | interpolating_function(const std::vector<float_type> &x, const std::vector<float_type> &f, |
---|
742 | bool lowerSlopeNatural=true, float_type lowerSlope=0.0, |
---|
743 | bool upperSlopeNatural=true, float_type upperSlope=0.0, |
---|
744 | float_type (*inputXConversion)(float_type)=0, |
---|
745 | float_type (*inputYConversion)(float_type)=0, |
---|
746 | float_type (*outputYConversion)(float_type)=0, |
---|
747 | float_type (*inputXConversionPrime)(float_type)=0, |
---|
748 | float_type (*inputYConversionPrime)(float_type)=0, |
---|
749 | float_type (*inputXConversionDPrime)(float_type)=0, |
---|
750 | float_type (*inputYConversionDPrime)(float_type)=0 |
---|
751 | ) throw(c2_exception) : c2_function<float_type>() |
---|
752 | { init(x, f, lowerSlopeNatural, lowerSlope, upperSlopeNatural, upperSlope, |
---|
753 | inputXConversion, inputYConversion, outputYConversion, |
---|
754 | inputXConversionPrime, inputYConversionPrime, |
---|
755 | inputXConversionDPrime, inputYConversionDPrime |
---|
756 | ); |
---|
757 | } |
---|
758 | |
---|
759 | /// \brief copy constructor |
---|
760 | /// \param rhs interpolating_function to copy from |
---|
761 | interpolating_function(const interpolating_function <float_type> &rhs) : c2_function<float_type>(rhs), |
---|
762 | Xraw(rhs.Xraw), X(rhs.X), F(rhs.F), y2(rhs.y2), |
---|
763 | fXin(rhs.fXin), fYin(rhs.fYin), fYout(rhs.fYout), |
---|
764 | fXinPrime(rhs.fXinPrime), fYinPrime(rhs.fYinPrime), |
---|
765 | fXinDPrime(rhs.fXinDPrime), fYinDPrime(rhs.fYinDPrime) , |
---|
766 | xInverted(rhs.xInverted), lastKLow(-1) |
---|
767 | { set_sampling_grid_pointer(Xraw); } |
---|
768 | |
---|
769 | virtual ~interpolating_function() { } // just to suppress warnings about no virtual destructor |
---|
770 | |
---|
771 | /// \brief retrieve copies of the x & y tables from which this was built |
---|
772 | /// |
---|
773 | /// This is often useful in the creation of new interpolating functions with transformed data. |
---|
774 | /// The vectorswill have their sizes set correctly on return. |
---|
775 | /// \param [in, out] xvals the abscissas |
---|
776 | /// \param [in, out] yvals the ordinates |
---|
777 | void get_data(std::vector<float_type> &xvals, std::vector<float_type> &yvals) const throw() ; |
---|
778 | |
---|
779 | /// \brief enable extrapolation of the function below the tabulated data. |
---|
780 | /// |
---|
781 | /// This allows the interpolator to be extrapolated outside the bounds of the data, |
---|
782 | /// using whatever derivatives it already had at the lower bound. |
---|
783 | /// \param bound the abscissa to which the function should be extended. |
---|
784 | void set_lower_extrapolation(float_type bound); |
---|
785 | /// \brief enable extrapolation of the function above the tabulated data. |
---|
786 | /// |
---|
787 | /// This allows the interpolator to be extrapolated outside the bounds of the data, |
---|
788 | /// using whatever derivatives it already had at the upper bound. |
---|
789 | /// \param bound the abscissa to which the function should be extended. |
---|
790 | void set_upper_extrapolation(float_type bound); |
---|
791 | |
---|
792 | // these functions correctly combine the interpolating function with another interpolating function |
---|
793 | // preserving the X bounds and mapping functions of the host (left hand) function. |
---|
794 | |
---|
795 | /// \brief create a new interpolating_function which is the \a source |
---|
796 | /// function applied to every point in the interpolating tables |
---|
797 | /// |
---|
798 | /// This carefully manages the derivative of the composed function at the two ends. |
---|
799 | /// \param source the function to apply |
---|
800 | /// \return a new interpolating_function with the same mappings for x and y |
---|
801 | interpolating_function <float_type> & unary_operator(const c2_function<float_type> &source) const; |
---|
802 | |
---|
803 | /// \brief create a new interpolating_function which is the parent interpolating_function |
---|
804 | /// combined with \a rhs using \a combiner at every point in the interpolating tables |
---|
805 | /// |
---|
806 | /// This carefully manages the derivative of the composed function at the two ends. |
---|
807 | /// \param rhs the function to apply |
---|
808 | /// \param combining_stub a function which defines which binary operation to use. |
---|
809 | /// \return a new interpolating_function with the same mappings for x and y |
---|
810 | interpolating_function <float_type> & binary_operator(const c2_function<float_type> &rhs, |
---|
811 | c2_binary_function<float_type> *combining_stub |
---|
812 | ) const; |
---|
813 | |
---|
814 | // InterpolatingFunctions override the c2_function operators, since they explicitly re-generate the interpolation table |
---|
815 | // when they are applied. If this is not desired, these operators are not virtual, so the interpolating_function |
---|
816 | // can be upcast back to a c2_function to produce unprocessed binaries. |
---|
817 | |
---|
818 | /// \brief produce a newly resampled interpolating_function which is the specified sum. |
---|
819 | /// \param rhs the function to add, pointwise |
---|
820 | /// \return a new interpolating_function |
---|
821 | /// \note |
---|
822 | /// InterpolatingFunctions override the c2_function operators, since they explicitly re-generate the interpolation table |
---|
823 | /// when they are applied. If this is not desired, these operators are not virtual, so the interpolating_function |
---|
824 | /// can be upcast back to a c2_function to produce unprocessed binaries. |
---|
825 | interpolating_function <float_type> & operator + (const c2_function<float_type> &rhs) const { |
---|
826 | return binary_operator(rhs, new c2_sum<float_type>()); } |
---|
827 | /// \brief produce a newly resampled interpolating_function which is the specified difference. |
---|
828 | /// \param rhs the function to subtract, pointwise |
---|
829 | /// \return a new interpolating_function |
---|
830 | /// \note |
---|
831 | /// InterpolatingFunctions override the c2_function operators, since they explicitly re-generate the interpolation table |
---|
832 | /// when they are applied. If this is not desired, these operators are not virtual, so the interpolating_function |
---|
833 | /// can be upcast back to a c2_function to produce unprocessed binaries. |
---|
834 | interpolating_function <float_type> & operator - (const c2_function<float_type> &rhs) const { |
---|
835 | return binary_operator(rhs, new c2_diff<float_type>()); } |
---|
836 | /// \brief produce a newly resampled interpolating_function which is the specified product. |
---|
837 | /// \param rhs the function to multiply, pointwise |
---|
838 | /// \return a new interpolating_function |
---|
839 | /// \note |
---|
840 | /// InterpolatingFunctions override the c2_function operators, since they explicitly re-generate the interpolation table |
---|
841 | /// when they are applied. If this is not desired, these operators are not virtual, so the interpolating_function |
---|
842 | /// can be upcast back to a c2_function to produce unprocessed binaries. |
---|
843 | interpolating_function <float_type> & operator * (const c2_function<float_type> &rhs) const { |
---|
844 | return binary_operator(rhs, new c2_product<float_type>()); } |
---|
845 | /// \brief produce a newly resampled interpolating_function which is the specified ratio. |
---|
846 | /// \param rhs the function to divide, pointwise |
---|
847 | /// \return a new interpolating_function |
---|
848 | /// \note |
---|
849 | /// InterpolatingFunctions override the c2_function operators, since they explicitly re-generate the interpolation table |
---|
850 | /// when they are applied. If this is not desired, these operators are not virtual, so the interpolating_function |
---|
851 | /// can be upcast back to a c2_function to produce unprocessed binaries. |
---|
852 | interpolating_function <float_type> & operator / (const c2_function<float_type> &rhs) const { |
---|
853 | return binary_operator(rhs, new c2_ratio<float_type>()); } |
---|
854 | /// \brief produce a newly resampled interpolating_function which is the specified sum. |
---|
855 | /// \param rhs a constant to add, pointwise |
---|
856 | /// \return a new interpolating_function |
---|
857 | /// \note |
---|
858 | /// InterpolatingFunctions override the c2_function operators, since they explicitly re-generate the interpolation table |
---|
859 | /// when they are applied. If this is not desired, these operators are not virtual, so the interpolating_function |
---|
860 | /// can be upcast back to a c2_function to produce unprocessed binaries. |
---|
861 | interpolating_function <float_type> & operator + (float_type rhs) const { return (*this)+c2_constant<float_type>(rhs); } |
---|
862 | /// \brief produce a newly resampled interpolating_function which is the specified difference. |
---|
863 | /// \param rhs a constant to subtract, pointwise |
---|
864 | /// \return a new interpolating_function |
---|
865 | /// \note |
---|
866 | /// InterpolatingFunctions override the c2_function operators, since they explicitly re-generate the interpolation table |
---|
867 | /// when they are applied. If this is not desired, these operators are not virtual, so the interpolating_function |
---|
868 | /// can be upcast back to a c2_function to produce unprocessed binaries. |
---|
869 | interpolating_function <float_type> & operator - (float_type rhs) const { return (*this)-c2_constant<float_type>(rhs); } |
---|
870 | /// \brief produce a newly resampled interpolating_function which is the specified product. |
---|
871 | /// \param rhs a constant to multiply, pointwise |
---|
872 | /// \return a new interpolating_function |
---|
873 | /// \note |
---|
874 | /// InterpolatingFunctions override the c2_function operators, since they explicitly re-generate the interpolation table |
---|
875 | /// when they are applied. If this is not desired, these operators are not virtual, so the interpolating_function |
---|
876 | /// can be upcast back to a c2_function to produce unprocessed binaries. |
---|
877 | interpolating_function <float_type> & operator * (float_type rhs) const { return (*this)*c2_constant<float_type>(rhs); } |
---|
878 | /// \brief produce a newly resampled interpolating_function which is the specified ratio. |
---|
879 | /// \param rhs a constant to divide, pointwise |
---|
880 | /// \return a new interpolating_function |
---|
881 | /// \note |
---|
882 | /// InterpolatingFunctions override the c2_function operators, since they explicitly re-generate the interpolation table |
---|
883 | /// when they are applied. If this is not desired, these operators are not virtual, so the interpolating_function |
---|
884 | /// can be upcast back to a c2_function to produce unprocessed binaries. |
---|
885 | interpolating_function <float_type> & operator / (float_type rhs) const { return (*this)/c2_constant<float_type>(rhs); } |
---|
886 | |
---|
887 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception); |
---|
888 | |
---|
889 | /// \brief move value & derivatives into our internal coordinates (use splint to go the other way!) |
---|
890 | /// \note why? |
---|
891 | void localize_derivatives(float_type xraw, float_type y, float_type yprime, float_type yprime2, float_type *y0, float_type *yp0, float_type *ypp0) const; |
---|
892 | |
---|
893 | protected: |
---|
894 | |
---|
895 | interpolating_function() : c2_function<float_type>() { } // default constructor is never used, prevent accidents by protecting it. |
---|
896 | |
---|
897 | /// \brief do the dirty work of constructing the spline. See interpolating_function constructor for details. |
---|
898 | void init(const std::vector<float_type> &, const std::vector<float_type> &, |
---|
899 | bool lowerSlopeNatural, float_type lowerSlope, |
---|
900 | bool upperSlopeNatural, float_type upperSlope, |
---|
901 | float_type (*inputXConversion)(float_type)=0, |
---|
902 | float_type (*inputXConversionPrime)(float_type)=0, |
---|
903 | float_type (*inputXConversionDPrime)(float_type)=0, |
---|
904 | float_type (*inputYConversion)(float_type)=0, |
---|
905 | float_type (*inputYConversionPrime)(float_type)=0, |
---|
906 | float_type (*inputYConversionDPrime)(float_type)=0, |
---|
907 | float_type (*outputYConversion)(float_type)=0 |
---|
908 | ) throw(c2_exception) ; |
---|
909 | |
---|
910 | std::vector<float_type> Xraw, X, F, y2; |
---|
911 | |
---|
912 | float_type (*fXin)(float_type), (*fYin)(float_type), (*fYout)(float_type); |
---|
913 | float_type (*fXinPrime)(float_type), (*fYinPrime)(float_type); |
---|
914 | float_type (*fXinDPrime)(float_type), (*fYinDPrime)(float_type); |
---|
915 | |
---|
916 | int xInverted; |
---|
917 | mutable int lastKLow; |
---|
918 | }; |
---|
919 | |
---|
920 | /// \brief An interpolatingFunction with X transformed into log space. |
---|
921 | /// |
---|
922 | /// Most useful for functions looking like y=log(x) or any other function with a huge X dynamic range, |
---|
923 | /// and a slowly varying Y. |
---|
924 | template <typename float_type=double> class log_lin_interpolating_function : public interpolating_function <float_type> { |
---|
925 | public: |
---|
926 | /// \brief Construct the function. |
---|
927 | /// \param x the list of abscissas. Must be either strictly increasing or strictly decreasing. |
---|
928 | /// Strictly increasing is preferred, as less memory is used since a copy is not required for the sampling grid. |
---|
929 | /// \param f the list of function values. |
---|
930 | /// \param lowerSlopeNatural if true, set y''(first point)=0 in LogLin space, otherwise compute it from \a lowerSope |
---|
931 | /// \param lowerSlope derivative of the function at the lower bound, used only if \a lowerSlopeNatural is false |
---|
932 | /// \param upperSlopeNatural if true, set y''(last point)=0 in LogLin space, otherwise compute it from \a upperSope |
---|
933 | /// \param upperSlope derivative of the function at the upper bound, used only if \a upperSlopeNatural is false |
---|
934 | log_lin_interpolating_function(const std::vector<float_type> &x, const std::vector<float_type> &f, |
---|
935 | bool lowerSlopeNatural=true, float_type lowerSlope=0.0, |
---|
936 | bool upperSlopeNatural=true, float_type upperSlope=0.0); |
---|
937 | protected: |
---|
938 | log_lin_interpolating_function() {} // do not allow naked construction... it is usually an accident. |
---|
939 | }; |
---|
940 | |
---|
941 | |
---|
942 | /// \brief An interpolatingFunction with Y transformed into log space. |
---|
943 | /// |
---|
944 | /// Most useful for functions looking like y=exp(x) |
---|
945 | template <typename float_type=double> class lin_log_interpolating_function : public interpolating_function <float_type> { |
---|
946 | public: |
---|
947 | /// \brief Construct the function. |
---|
948 | /// \param x the list of abscissas. Must be either strictly increasing or strictly decreasing. |
---|
949 | /// Strictly increasing is preferred, as less memory is used since a copy is not required for the sampling grid. |
---|
950 | /// \param f the list of function values. |
---|
951 | /// \param lowerSlopeNatural if true, set y''(first point)=0 in LinLog space, otherwise compute it from \a lowerSope |
---|
952 | /// \param lowerSlope derivative of the function at the lower bound, used only if \a lowerSlopeNatural is false |
---|
953 | /// \param upperSlopeNatural if true, set y''(last point)=0 in LinLog space, otherwise compute it from \a upperSope |
---|
954 | /// \param upperSlope derivative of the function at the upper bound, used only if \a upperSlopeNatural is false |
---|
955 | lin_log_interpolating_function(const std::vector<float_type> &x, const std::vector<float_type> &f, |
---|
956 | bool lowerSlopeNatural=true, float_type lowerSlope=0.0, |
---|
957 | bool upperSlopeNatural=true, float_type upperSlope=0.0); |
---|
958 | protected: |
---|
959 | lin_log_interpolating_function() {} // do not allow naked construction... it is usually an accident. |
---|
960 | }; |
---|
961 | |
---|
962 | |
---|
963 | /// \brief An interpolatingFunction with X and Y transformed into log space. |
---|
964 | /// |
---|
965 | /// Most useful for functions looking like y=x^n or any other function with a huge X and Y dynamic range. |
---|
966 | template <typename float_type=double> class log_log_interpolating_function : public interpolating_function <float_type> { |
---|
967 | public: |
---|
968 | /// \brief Construct the function. |
---|
969 | /// \param x the list of abscissas. Must be either strictly increasing or strictly decreasing. |
---|
970 | /// Strictly increasing is preferred, as less memory is used since a copy is not required for the sampling grid. |
---|
971 | /// \param f the list of function values. |
---|
972 | /// \param lowerSlopeNatural if true, set y''(first point)=0 in LogLog space, otherwise compute it from \a lowerSope |
---|
973 | /// \param lowerSlope derivative of the function at the lower bound, used only if \a lowerSlopeNatural is false |
---|
974 | /// \param upperSlopeNatural if true, set y''(last point)=0 in LogLog space, otherwise compute it from \a upperSope |
---|
975 | /// \param upperSlope derivative of the function at the upper bound, used only if \a upperSlopeNatural is false |
---|
976 | log_log_interpolating_function(const std::vector<float_type> &x, const std::vector<float_type> &f, |
---|
977 | bool lowerSlopeNatural=true, float_type lowerSlope=0.0, |
---|
978 | bool upperSlopeNatural=true, float_type upperSlope=0.0); |
---|
979 | protected: |
---|
980 | log_log_interpolating_function() {} // do not allow naked construction... it is usually an accident. |
---|
981 | }; |
---|
982 | |
---|
983 | |
---|
984 | /// \brief An interpolating_function with X in reciprocal space and Y transformed in log space. |
---|
985 | /// |
---|
986 | /// Most useful for thermodynamic types of data where Y is roughly A*exp(-B/x). |
---|
987 | /// Typical examples are reaction rate data, and thermistor calibration data. |
---|
988 | template <typename float_type=double> class arrhenius_interpolating_function : public interpolating_function <float_type> { |
---|
989 | public: |
---|
990 | /// \brief Construct the function. |
---|
991 | /// \param x the list of abscissas. Must be either strictly increasing or strictly decreasing. |
---|
992 | /// Strictly increasing is preferred, as less memory is used since a copy is not required for the sampling grid. |
---|
993 | /// \param f the list of function values. |
---|
994 | /// \param lowerSlopeNatural if true, set y''(first point)=0 in Arrhenius space, otherwise compute it from \a lowerSope |
---|
995 | /// \param lowerSlope derivative of the function at the lower bound, used only if \a lowerSlopeNatural is false |
---|
996 | /// \param upperSlopeNatural if true, set y''(last point)=0 in Arrhenius space, otherwise compute it from \a upperSope |
---|
997 | /// \param upperSlope derivative of the function at the upper bound, used only if \a upperSlopeNatural is false |
---|
998 | arrhenius_interpolating_function(const std::vector<float_type> &x, const std::vector<float_type> &f, |
---|
999 | bool lowerSlopeNatural=true, float_type lowerSlope=0.0, |
---|
1000 | bool upperSlopeNatural=true, float_type upperSlope=0.0); |
---|
1001 | protected: |
---|
1002 | arrhenius_interpolating_function() {} // do not allow naked construction... it is usually an accident. |
---|
1003 | }; |
---|
1004 | |
---|
1005 | /** |
---|
1006 | \brief create a linear-linear interpolating grid with both x & y set to |
---|
1007 | (xmin, xmin+dx, ... xmin + (count-1)*dx ) |
---|
1008 | |
---|
1009 | very useful for transformaiton with other functions e.g. |
---|
1010 | \code |
---|
1011 | f=c2_sin<double>::sin(LinearInterpolatingGrid(-0.1,0.1, 65)) |
---|
1012 | \endcode |
---|
1013 | creates a spline table of sin(x) slightly beyond the first period |
---|
1014 | \param xmin the starting point for the grid |
---|
1015 | \param dx the step size for the grid |
---|
1016 | \param count the number of points in the grid |
---|
1017 | \return an identity interpolating_function with the requested grid |
---|
1018 | */ |
---|
1019 | template <typename float_type> interpolating_function <float_type> &linear_interpolating_grid(float_type xmin, float_type dx, int count) { |
---|
1020 | std::vector<float_type> x(count); |
---|
1021 | for(int i=0; i<count; i++) x[i]=xmin + i * dx; |
---|
1022 | return *new interpolating_function <float_type>(x,x); |
---|
1023 | } |
---|
1024 | |
---|
1025 | /** |
---|
1026 | \brief create a log-log interpolating grid with both x & y set to |
---|
1027 | (xmin, xmin*dx, ... xmin * dx^(count-1) ) |
---|
1028 | |
---|
1029 | very useful for transformaiton with other functions e.g. |
---|
1030 | \code |
---|
1031 | f=c2_log<double>::log(LogLogInterpolatingGrid(2, 1.1, 65)) |
---|
1032 | \endcode |
---|
1033 | creates a spline table of log(x) |
---|
1034 | \param xmin the starting point for the grid |
---|
1035 | \param dx the ratio between points |
---|
1036 | \param count the number of points in the grid |
---|
1037 | \return an identity log_log_interpolating_function with the requested grid |
---|
1038 | */ |
---|
1039 | template <typename float_type> log_log_interpolating_function <float_type> &log_log_interpolating_grid(float_type xmin, float_type dx, int count) { |
---|
1040 | std::vector<float_type> x(count); |
---|
1041 | x[0]=xmin; |
---|
1042 | for(int i=1; i<count; i++) x[i]=dx*x[i-1]; |
---|
1043 | return *new log_log_interpolating_function<float_type>(x,x); |
---|
1044 | } |
---|
1045 | |
---|
1046 | /// \brief compute sin(x) with its derivatives. |
---|
1047 | /// |
---|
1048 | /// Creates a singleton instance c2_sin::sin of itself for convenient access. |
---|
1049 | template <typename float_type=double> class c2_sin : public c2_function<float_type> { |
---|
1050 | public: |
---|
1051 | /// \brief constructor. There is alread a singleton c2_sin::sin, which usually suffices. |
---|
1052 | c2_sin() {} |
---|
1053 | |
---|
1054 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
1055 | { float_type q=std::sin(x); if(yprime) *yprime=std::cos(x); if(yprime2) *yprime2=-q; return q; } |
---|
1056 | |
---|
1057 | /// \brief return a grid dynamically, suitable for use with trig functions with period 2*pi |
---|
1058 | /// \param xmin the lower bound for the grid |
---|
1059 | /// \param xmax upper bound for the grid |
---|
1060 | /// \return a new sampling grid. |
---|
1061 | virtual std::vector<float_type> &get_sampling_grid(float_type xmin, float_type xmax); |
---|
1062 | /// \brief the static singleton |
---|
1063 | static const c2_sin sin; |
---|
1064 | }; |
---|
1065 | /// \brief compute cos(x) with its derivatives. |
---|
1066 | /// |
---|
1067 | /// Creates a singleton instance c2_cos::cos of itself for convenient access. |
---|
1068 | template <typename float_type=double> class c2_cos : public c2_sin<float_type> { |
---|
1069 | public: |
---|
1070 | /// \brief constructor. There is already a singleton c2_cos::cos, which usually suffices. |
---|
1071 | c2_cos() {} |
---|
1072 | |
---|
1073 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
1074 | { float_type q=std::cos(x); if(yprime) *yprime=-std::sin(x); if(yprime2) *yprime2=-q; return q; } |
---|
1075 | /// \brief the static singleton |
---|
1076 | static const c2_cos cos; |
---|
1077 | }; |
---|
1078 | /// \brief compute tan(x) with its derivatives. |
---|
1079 | /// |
---|
1080 | /// Creates a singleton instance c2_tan::tan of itself for convenient access. |
---|
1081 | template <typename float_type=double> class c2_tan : public c2_function<float_type> { |
---|
1082 | public: |
---|
1083 | /// \brief constructor. There is already a singleton c2_tan::tan, which usually suffices. |
---|
1084 | c2_tan() {} |
---|
1085 | |
---|
1086 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
1087 | { |
---|
1088 | float_type c=std::cos(x), s=std::sin(x); |
---|
1089 | float_type t=s/c; |
---|
1090 | float_type yp=1/(c*c); |
---|
1091 | if(yprime) *yprime=yp; if(yprime2) *yprime2=2*t*yp; |
---|
1092 | return t; |
---|
1093 | } |
---|
1094 | /// \brief the static singleton |
---|
1095 | static const c2_tan tan; |
---|
1096 | }; |
---|
1097 | /// \brief compute log(x) with its derivatives. |
---|
1098 | /// |
---|
1099 | /// Creates a singleton instance c2_log::log of itself for convenient access. |
---|
1100 | template <typename float_type=double> class c2_log : public c2_function<float_type> { |
---|
1101 | public: |
---|
1102 | /// \brief constructor. There is already a singleton c2_log::log, which usually suffices. |
---|
1103 | c2_log() {} |
---|
1104 | |
---|
1105 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
1106 | { if(yprime) *yprime=1.0/x; if(yprime2) *yprime2=-1.0/(x*x); return std::log(x); } |
---|
1107 | /// \brief the static singleton |
---|
1108 | static const c2_log log; |
---|
1109 | }; |
---|
1110 | /// \brief compute exp(x) with its derivatives. |
---|
1111 | /// |
---|
1112 | /// Creates a singleton instance c2_exp::exp of itself for convenient access. |
---|
1113 | template <typename float_type=double> class c2_exp : public c2_function<float_type> { |
---|
1114 | public: |
---|
1115 | /// \brief constructor. There is already a singleton c2_exp::exp, which usually suffices. |
---|
1116 | c2_exp() {} |
---|
1117 | |
---|
1118 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
1119 | { float_type q=std::exp(x); if(yprime) *yprime=q; if(yprime2) *yprime2=q; return q; } |
---|
1120 | /// \ brief the static singleton |
---|
1121 | static const c2_exp exp; |
---|
1122 | }; |
---|
1123 | /// \brief compute sqrt(x) with its derivatives. |
---|
1124 | /// |
---|
1125 | /// Creates a singleton instance c2_sqrt::sqrt of itself for convenient access. |
---|
1126 | template <typename float_type=double> class c2_sqrt : public c2_function<float_type> { |
---|
1127 | public: |
---|
1128 | /// \brief constructor. There is already a singleton c2_sqrt::sqrt, which usually suffices. |
---|
1129 | c2_sqrt() {} |
---|
1130 | |
---|
1131 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
1132 | { float_type q=std::sqrt(x); if(yprime) *yprime=0.5/q; if(yprime2) *yprime2=-0.25/(x*q); return q; } |
---|
1133 | /// \brief the static singleton |
---|
1134 | static const c2_sqrt sqrt; |
---|
1135 | }; |
---|
1136 | /// \brief compute scale/x with its derivatives. |
---|
1137 | /// |
---|
1138 | /// Creates a singleton instance c2_recip:recip of itself for convenient access. |
---|
1139 | template <typename float_type=double> class c2_recip : public c2_function<float_type> { |
---|
1140 | public: |
---|
1141 | /// \brief constructor. There is already a singleton c2_recip::recip, which usually suffices. |
---|
1142 | c2_recip(float_type scale) : rscale(scale) {} |
---|
1143 | |
---|
1144 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
1145 | { |
---|
1146 | float_type q=1.0/x; |
---|
1147 | float_type y=rscale*q; |
---|
1148 | if(yprime) *yprime=-y*q; |
---|
1149 | if(yprime2) *yprime2=2*y*q*q; |
---|
1150 | return y; |
---|
1151 | } |
---|
1152 | /// \brief reset the scale factor |
---|
1153 | /// \param scale the new numerator |
---|
1154 | void reset(float_type scale) { rscale=scale; } |
---|
1155 | /// \brief the static singleton |
---|
1156 | static const c2_recip recip; |
---|
1157 | private: |
---|
1158 | float_type rscale; |
---|
1159 | }; |
---|
1160 | /// \brief compute x with its derivatives. |
---|
1161 | /// |
---|
1162 | /// Creates a singleton instance c2_identity::identity of itself for convenient access. |
---|
1163 | template <typename float_type=double> class c2_identity : public c2_function<float_type> { |
---|
1164 | public: |
---|
1165 | /// \brief constructor. There is already a singleton c2_identity::identity, which usually suffices. |
---|
1166 | c2_identity() {} |
---|
1167 | |
---|
1168 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
1169 | { if(yprime) *yprime=1.0; if(yprime2) *yprime2=0; return x; } |
---|
1170 | /// \brief the static singleton |
---|
1171 | static const c2_identity identity; |
---|
1172 | }; |
---|
1173 | |
---|
1174 | /** |
---|
1175 | \brief create a linear mapping of another function |
---|
1176 | |
---|
1177 | for example, given a c2_function \a f |
---|
1178 | \code |
---|
1179 | c2_linear<double> L(1.2, 2.0, 3.0); |
---|
1180 | c2_composed_function<double> &F=L(f); |
---|
1181 | \endcode |
---|
1182 | produces a new c2_function F=2.0+3.0*(\a f - 1.2) |
---|
1183 | */ |
---|
1184 | template <typename float_type=double> class c2_linear : public c2_function<float_type> { |
---|
1185 | public: |
---|
1186 | /// \brief Construct the operator f=y0 + slope * (x-x0) |
---|
1187 | /// \param x0 the x offset |
---|
1188 | /// \param y0 the y-intercept i.e. f(x0) |
---|
1189 | /// \param slope the slope of the mapping |
---|
1190 | c2_linear(float_type x0, float_type y0, float_type slope) : xint(x0), intercept(y0), m(slope) {} |
---|
1191 | /// \brief Change the slope and intercepts after construction. |
---|
1192 | /// \param x0 the x offset |
---|
1193 | /// \param y0 the y-intercept |
---|
1194 | /// \param slope the slope of the mapping |
---|
1195 | void reset(float_type x0, float_type y0, float_type slope) { xint=x0; intercept=y0; m=slope; } |
---|
1196 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
1197 | { if(yprime) *yprime=m; if(yprime2) *yprime2=0; return m*(x-xint)+intercept; } |
---|
1198 | |
---|
1199 | private: |
---|
1200 | float_type xint, intercept, m; |
---|
1201 | protected: |
---|
1202 | c2_linear() {} // do not allow naked construction... it is usually an accident. |
---|
1203 | }; |
---|
1204 | |
---|
1205 | /** |
---|
1206 | \brief create a quadratic mapping of another function |
---|
1207 | |
---|
1208 | for example, given a c2_function \a f |
---|
1209 | \code |
---|
1210 | c2_quadratic<double> Q(1.2, 2.0, 3.0, 4.0); |
---|
1211 | c2_composed_function<double> &F=Q(f); |
---|
1212 | \endcode |
---|
1213 | produces a new c2_function F=2.0 + 3.0*(f-1.2) + 4.0*(f-1.2)^2 |
---|
1214 | |
---|
1215 | note that the parameters are overdetermined, but allows the flexibility of two different representations |
---|
1216 | |
---|
1217 | */ |
---|
1218 | template <typename float_type=double> class c2_quadratic : public c2_function<float_type> { |
---|
1219 | public: |
---|
1220 | /// \brief Construct the operator |
---|
1221 | /// \param x0 the center around which the powers are computed |
---|
1222 | /// \param y0 the value of the function at \a x = \a x0 |
---|
1223 | /// \param xcoef the scale on the (\a x - \a x0) term |
---|
1224 | /// \param x2coef the scale on the (\a x - \a x0)^2 term |
---|
1225 | c2_quadratic(float_type x0, float_type y0, float_type xcoef, float_type x2coef) : intercept(y0), center(x0), a(x2coef), b(xcoef) {} |
---|
1226 | /// Modify the mapping after construction |
---|
1227 | /// \param x0 the new center around which the powers are computed |
---|
1228 | /// \param y0 the new value of the function at \a x = \a x0 |
---|
1229 | /// \param xcoef the new scale on the (\a x - \a x0) term |
---|
1230 | /// \param x2coef the new scale on the (\a x - \a x0)^2 term |
---|
1231 | void reset(float_type x0, float_type y0, float_type xcoef, float_type x2coef) { intercept=y0; center=x0; a=x2coef; b=xcoef; } |
---|
1232 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
1233 | { float_type dx=x-center; if(yprime) *yprime=2*a*dx+b; if(yprime2) *yprime2=2*a; return a*dx*dx+b*dx+intercept; } |
---|
1234 | |
---|
1235 | private: |
---|
1236 | float_type intercept, center, a, b; |
---|
1237 | protected: |
---|
1238 | c2_quadratic() {} // do not allow naked construction... it is usually an accident. |
---|
1239 | }; |
---|
1240 | |
---|
1241 | /** |
---|
1242 | \brief create a power law mapping of another function |
---|
1243 | |
---|
1244 | for example, given a c2_function \a f |
---|
1245 | \code |
---|
1246 | c2_power_law<double> PLaw(1.2, 2.5); |
---|
1247 | c2_composed_function<double> &F=PLaw(f); |
---|
1248 | \endcode |
---|
1249 | produces a new c2_function F=1.2 * f^2.5 |
---|
1250 | |
---|
1251 | */ |
---|
1252 | template <typename float_type=double> class c2_power_law : public c2_function<float_type> { |
---|
1253 | public: |
---|
1254 | /// \brief Construct the operator |
---|
1255 | /// \param scale the multipler |
---|
1256 | /// \param power the exponent |
---|
1257 | c2_power_law(float_type scale, float_type power) : a(scale), b(power) {} |
---|
1258 | /// \brief Modify the mapping after construction |
---|
1259 | /// \param scale the new multipler |
---|
1260 | /// \param power the new exponent |
---|
1261 | void reset(float_type scale, float_type power) { a=scale; b=power; } |
---|
1262 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) |
---|
1263 | { 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; } |
---|
1264 | |
---|
1265 | private: |
---|
1266 | float_type a, b; |
---|
1267 | protected: |
---|
1268 | c2_power_law() {} // do not allow naked construction... it is usually an accident. |
---|
1269 | }; |
---|
1270 | |
---|
1271 | /** |
---|
1272 | \brief create the formal inverse function of another function |
---|
1273 | |
---|
1274 | for example, given a c2_function \a f |
---|
1275 | \code |
---|
1276 | c2_inverse_function<double> inv(f); |
---|
1277 | a=f(x); |
---|
1278 | x1=inv(a); |
---|
1279 | \endcode |
---|
1280 | will return x1=x to machine precision. The important part of this |
---|
1281 | is that the resulting function is a first-class c2_function, so it knows its |
---|
1282 | derivatives, too, unlike the case of a simple root-finding inverse. This means |
---|
1283 | it can be integrated (for example) quite efficiently. |
---|
1284 | |
---|
1285 | Note that it is a subclass of c2_scaled_function only to manage ownership of another c2_function. |
---|
1286 | */ |
---|
1287 | template <typename float_type=double> class c2_inverse_function : public c2_plugin_function<float_type> { |
---|
1288 | public: |
---|
1289 | /// \brief Construct the operator |
---|
1290 | /// \param source the function to be inverted |
---|
1291 | c2_inverse_function(const c2_function<float_type> &source); |
---|
1292 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception); |
---|
1293 | |
---|
1294 | /// \brief give the function a hint as to where to look for its inverse |
---|
1295 | /// \param hint the likely value of the inverse, which defaults to whatever the evaluation returned. |
---|
1296 | void set_start_hint(float_type hint) const { start_hint=hint; } |
---|
1297 | |
---|
1298 | /// \brief get the starting hint. |
---|
1299 | /// |
---|
1300 | /// This is virtual so if there is a better way, this can be easily overridden. |
---|
1301 | /// It is used in value_with_derivatives() to guess where to start the root finder. |
---|
1302 | /// \param x the abscissa for which an estimate is needed |
---|
1303 | virtual float_type get_start_hint(float_type x) const { return start_hint; } |
---|
1304 | |
---|
1305 | protected: |
---|
1306 | c2_inverse_function() {} // do not allow naked construction... it is usually an accident. |
---|
1307 | mutable float_type start_hint; |
---|
1308 | }; |
---|
1309 | |
---|
1310 | /** |
---|
1311 | \brief |
---|
1312 | An interpolating_function which is the cumulative integral of a histogram. |
---|
1313 | |
---|
1314 | Note than binedges should be one element longer than binheights, since the lower & upper edges are specified. |
---|
1315 | Note that this is a malformed spline, since the second derivatives are all zero, so it has less continuity. |
---|
1316 | Also, note that the bin edges can be given in backwards order to generate the |
---|
1317 | reversed accumulation (starting at the high end) |
---|
1318 | */ |
---|
1319 | |
---|
1320 | template <typename float_type=double> class accumulated_histogram : public interpolating_function <float_type> { |
---|
1321 | public: |
---|
1322 | /// \brief Construct the integrated histogram |
---|
1323 | /// \param binedges the edges of the bins in \a binheights. It should have one more element than \a binheights |
---|
1324 | /// \param binheights the number of counts in each bin. |
---|
1325 | /// \param normalize if true, normalize integral to 1 |
---|
1326 | /// \param inverse_function if true, drop zero channels, and return inverse function for random generation |
---|
1327 | /// \param drop_zeros eliminate null bins before integrating, so integral is strictly monotonic. |
---|
1328 | accumulated_histogram(const std::vector<float_type>binedges, const std::vector<float_type> binheights, |
---|
1329 | bool normalize=false, bool inverse_function=false, bool drop_zeros=true); |
---|
1330 | |
---|
1331 | }; |
---|
1332 | |
---|
1333 | /** |
---|
1334 | \brief Construct a function useful for generation of random numbers from the given distribution |
---|
1335 | |
---|
1336 | inverse_integrated_density<InterpolatingFunctionFlavor>() starts with a probability density c2_function, generates the integral, |
---|
1337 | and generates an interpolating_function which, when evaluated using a uniform random on [0,1] returns values |
---|
1338 | with a density distribution equal to the input distribution |
---|
1339 | If the data are passed in reverse order (large X first), the integral is carried out from the big end. |
---|
1340 | |
---|
1341 | \sa template <typename Intermediate, typename Final> Final inverse_integrated_density(const std::vector, c2_function &) |
---|
1342 | |
---|
1343 | \param bincenters points at which to sample the c2_function \a binheights |
---|
1344 | \param binheights a c2_function which describes the random number distribution to be produced. |
---|
1345 | \return an interpolating_function of the type requested in the template which, |
---|
1346 | if evaluated randomly with a uniform variate on [0,1) produces numbers |
---|
1347 | distributed according to \a binheights |
---|
1348 | */ |
---|
1349 | |
---|
1350 | template <typename float_type, typename Final > |
---|
1351 | Final &inverse_integrated_density(const std::vector<float_type> &bincenters, c2_function<float_type> &binheights) |
---|
1352 | { |
---|
1353 | std::vector<float_type> integral; |
---|
1354 | |
---|
1355 | // integrate from first to last bin in original order, leaving results in integral |
---|
1356 | // ask for relative error of 1e-6 on each bin, with absolute error set to 0 (since we don't know the data scale). |
---|
1357 | float_type sum=binheights.partial_integrals(bincenters, &integral, 0.0, 1e-6); |
---|
1358 | // the integral vector now has partial integrals... it must be accumulated by summing |
---|
1359 | integral.insert(integral.begin(), 0.0); // integral from start to start is 0 |
---|
1360 | float_type scale=1.0/sum; |
---|
1361 | for(size_t i=1; i<integral.size(); i++) integral[i]=integral[i]*scale + integral[i-1]; |
---|
1362 | integral.back()=1.0; // force exact value on boundary |
---|
1363 | |
---|
1364 | return *new Final(integral, bincenters, |
---|
1365 | false, 1.0/(scale*binheights(bincenters.front() )), |
---|
1366 | false, 1.0/(scale*binheights(bincenters.back() )) |
---|
1367 | ); // use integral as x axis in inverse function |
---|
1368 | } |
---|
1369 | |
---|
1370 | /** |
---|
1371 | \brief Construct a function useful for generation of random numbers from the given distribution |
---|
1372 | |
---|
1373 | \code |
---|
1374 | template <typename Intermediate, typename Final> |
---|
1375 | Final & inverse_integrated_density(const std::vector &bincenters, const std::vector &binheights) |
---|
1376 | \endcode |
---|
1377 | is a variant of \code |
---|
1378 | template <typename Final> |
---|
1379 | Final & inverse_integrated_density(const std::vector &bincenters, c2_function &binheights) |
---|
1380 | \endcode |
---|
1381 | which takes two std::vectors and generates the intermediate interpolating_function required for |
---|
1382 | inverse_integrated_density(), and then calls it. |
---|
1383 | |
---|
1384 | \param bincenters points at which \a binheights are defined |
---|
1385 | \param binheights an std::vector which describes the random number distribution to be produced. |
---|
1386 | \return an interpolating_function of the type requested in the template which, |
---|
1387 | if evaluated randomly with a uniform variate on [0,1) produces numbers |
---|
1388 | distributed according to \a binheights |
---|
1389 | */ |
---|
1390 | |
---|
1391 | template <typename float_type, typename Intermediate, typename Final> Final |
---|
1392 | &inverse_integrated_density(const std::vector<float_type> &bincenters, const std::vector<float_type> &binheights) |
---|
1393 | { |
---|
1394 | std::vector<float_type> be(bincenters), bh(binheights); |
---|
1395 | |
---|
1396 | if(be[1] < be[0]) { // reverse data for interpolator if x axis passed in backwards |
---|
1397 | std::reverse(be.begin(), be.end()); |
---|
1398 | std::reverse(bh.begin(), bh.end()); |
---|
1399 | } |
---|
1400 | |
---|
1401 | Intermediate temp(be, bh); // create a temporary interpolating_function to integrate |
---|
1402 | Final &result=inverse_integrated_density<Final>(bincenters, temp); |
---|
1403 | |
---|
1404 | return result; |
---|
1405 | } |
---|
1406 | |
---|
1407 | /// \brief create a c2_function which smoothly connects two other c2_functions. |
---|
1408 | /// |
---|
1409 | /// This takes two points and generates a polynomial which matches two c2_function arguments |
---|
1410 | /// at those two points, with two derivatives at each point, and an arbitrary value at the center of the |
---|
1411 | /// region. It is useful for splicing together functions over rough spots (0/0, for example). |
---|
1412 | /// |
---|
1413 | /// If \a auto_center is true, the value at the midpoint is computed so that the resulting polynomial is |
---|
1414 | /// of order 5. If \a auto_center is false, the value \a y1 is used at the midpoint, resulting in a |
---|
1415 | /// polynomial of order 6. |
---|
1416 | template <typename float_type=double> class c2_connector_function : public c2_function<float_type> { |
---|
1417 | public: |
---|
1418 | /// \brief construct the container |
---|
1419 | /// \param f1 the function on the left side to be connected |
---|
1420 | /// \param f2 the function on the right side to be connected |
---|
1421 | /// \param x0 the point at which to match \a f1 and its derivatives |
---|
1422 | /// \param x2 the point at which to match \a f2 and its derivatives |
---|
1423 | /// \param auto_center if true, no midpoint value is specified. If false, match the value \a y1 at the midpoint |
---|
1424 | /// \param y1 the value to match at the midpoint, if \a auto_center is false |
---|
1425 | /// \return a c2_function with domain (\a x0,\a x2) which smoothly connects \a f1 and \a f2 |
---|
1426 | c2_connector_function(const c2_function<float_type> &f1, const c2_function<float_type> &f2, float_type x0, float_type x2, |
---|
1427 | bool auto_center, float_type y1); |
---|
1428 | /// \brief destructor |
---|
1429 | virtual ~c2_connector_function(); |
---|
1430 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw (c2_exception); |
---|
1431 | protected: |
---|
1432 | float_type fx1, fhinv, fdx, fy1, fa, fb, fc, fd, fe, ff; |
---|
1433 | }; |
---|
1434 | |
---|
1435 | |
---|
1436 | |
---|
1437 | /// \brief create a c2_function which is a piecewise assembly of other c2_functions. |
---|
1438 | /// |
---|
1439 | /// The functions must have increasing, non-overlapping domains. Any empty space |
---|
1440 | /// between functions will be filled with a linear interpolation. |
---|
1441 | /// \note The creation of the container results in the creation of an explicit sampling grid. |
---|
1442 | /// If this is used with functions with a large domain, or which generate very dense sampling grids, |
---|
1443 | /// it could eat a lot of memory. Do not abuse this by using functions which can generate gigantic grids. |
---|
1444 | /// |
---|
1445 | /// See c2_plugin_function for a discussion of how this might be used. |
---|
1446 | template <typename float_type=double> class c2_piecewise_function : public c2_function<float_type> { |
---|
1447 | public: |
---|
1448 | /// \brief construct the container |
---|
1449 | c2_piecewise_function(); |
---|
1450 | /// \brief destructor |
---|
1451 | virtual ~c2_piecewise_function(); |
---|
1452 | virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw (c2_exception); |
---|
1453 | /// \brief append a new function to the sequence |
---|
1454 | /// |
---|
1455 | /// This takes a c2_function, and appends it onto the end of the piecewise collection. |
---|
1456 | /// The domain of the function (which MUST be set) specifies the place it will be used in |
---|
1457 | /// the final function. If the domain exactly abuts the domain of the previous function, it |
---|
1458 | /// will be directly attached. If there is a gap, the gap will be filled in by linear interpolation. |
---|
1459 | /// If the function being appended is to be deleted automatically when this container is deleted, set the pass_ownership flag. |
---|
1460 | /// \param func a c2_function with a defined domain to be appended to the collection |
---|
1461 | /// \param pass_ownership if set, \a func will be deleted when the container is destroyed |
---|
1462 | void append_function(c2_function<float_type> &func, bool pass_ownership) throw (c2_exception); |
---|
1463 | protected: |
---|
1464 | std::vector<c2_function<float_type> *> functions; |
---|
1465 | std::vector<bool> owns; |
---|
1466 | mutable int lastKLow; |
---|
1467 | }; |
---|
1468 | |
---|
1469 | #include "c2_function.icc" |
---|
1470 | |
---|
1471 | #endif |
---|