[807] | 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 |
---|