source: trunk/examples/extended/electromagnetic/TestEm7/include/c2_function.icc@ 1036

Last change on this file since 1036 was 807, checked in by garnier, 17 years ago

update
File size: 41.3 KB
Line 
1/**
2 * \file
3 * \brief Provides code 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.cc,v 1.43 2007/11/12 20:22:54 marcus Exp
10 */
11
12#include <iostream>
13#include <vector>
14#include <algorithm>
15#include <cstdlib>
16#include <numeric>
17#include <functional>
18#include <iterator>
19#include <cmath>
20#include <limits>
21#include <sstream>
22
23template <typename float_type> const std::string c2_function<float_type>::cvs_file_vers() const
24{ return "c2_function.cc,v 1.43 2007/11/12 20:22:54 marcus Exp"; }
25
26// find a pre-bracketed root of a c2_function, which is a MUCH easier job than general root finding
27// since the derivatives are known exactly, and smoothness is guaranteed.
28// this searches for f(x)=value, to make life a little easier than always searching for f(x)=0
29template <typename float_type> float_type c2_function<float_type>::find_root(float_type lower_bracket, float_type upper_bracket,
30 float_type start, float_type value, int *error,
31 float_type *final_yprime, float_type *final_yprime2) const throw(c2_exception)
32{
33 // find f(x)=value within the brackets, using the guarantees of smoothness associated with a c2_function
34 // can use local f(x)=a*x**2 + b*x + c and solve quadratic to find root, then iterate
35 reset_evaluations();
36
37 float_type yp, yp2; // we will make unused pointers point here, to save null checks later
38 if (!final_yprime) final_yprime=&yp;
39 if (!final_yprime2) final_yprime2=&yp2;
40
41 float_type ftol=5*(std::numeric_limits<float_type>::epsilon()*std::abs(value)+std::numeric_limits<float_type>::min());
42 float_type xtol=5*(std::numeric_limits<float_type>::epsilon()*(std::abs(upper_bracket)+std::abs(lower_bracket))+std::numeric_limits<float_type>::min());
43
44 float_type root=start; // start looking in the middle
45 if(error) *error=0; // start out with error flag set to OK, if it is expected
46
47 float_type c, b;
48
49 // this new logic is to keep track of where we were before, and lower the number of
50 // function evaluations if we are searching inside the same bracket as before.
51 // Since this root finder has, very often, the bracket of the entire domain of the function,
52 // this makes a big difference, especially to c2_inverse_function
53 if(!rootInitialized || upper_bracket != lastRootUpperX || lower_bracket != lastRootLowerX) {
54 lastRootUpperY=value_with_derivatives(upper_bracket, final_yprime, final_yprime2);
55 increment_evaluations();
56 lastRootUpperX=upper_bracket;
57
58 lastRootLowerY=value_with_derivatives(lower_bracket, final_yprime, final_yprime2);
59 increment_evaluations();
60 lastRootLowerX=lower_bracket;
61 rootInitialized=true;
62 }
63
64 float_type clower=lastRootLowerY-value;
65 float_type cupper=lastRootUpperY-value;
66 if(clower*cupper >0) {
67 // argh, no sign change in here!
68 if(error) { *error=1; return 0.0; }
69 else {
70 std::ostringstream outstr;
71 outstr << "unbracketed root in find_root at xlower= " << lower_bracket << ", xupper= " << upper_bracket;
72 outstr << ", value= " << value << ": bailing";
73 throw c2_exception(outstr.str().c_str());
74 }
75 }
76
77 float_type delta=upper_bracket-lower_bracket; // first error step
78 c=value_with_derivatives(root, final_yprime, final_yprime2)-value; // compute initial values
79 b=*final_yprime; // make a local copy for readability
80 increment_evaluations();
81
82 while(
83 std::abs(delta) > xtol && // absolute x step check
84 std::abs(c) > ftol && // absolute y tolerance
85 std::abs(c) > xtol*std::abs(b) // comparison to smallest possible Y step from derivative
86 )
87 {
88 float_type a=(*final_yprime2)/2; // second derivative is 2*a
89 float_type disc=b*b-4*a*c;
90 // std::cout << std::endl << "find_root_debug a,b,c,d " << a << " " << b << " " << c << " " << disc << std::endl;
91
92 if(disc >= 0) {
93 float_type q=-0.5*((b>=0)?(b+std::sqrt(disc)):(b-std::sqrt(disc)));
94 if(q*q > std::abs(a*c)) delta=c/q; // since x1=q/a, x2=c/q, x1/x2=q^2/ac, this picks smaller step
95 else delta=q/a;
96 root+=delta;
97 }
98
99 if(disc < 0 || root<lower_bracket || root>upper_bracket ||
100 std::abs(delta) >= 0.5*(upper_bracket-lower_bracket)) {
101 // if we jump out of the bracket, or aren't converging well, bisect
102 root=0.5*(lower_bracket+upper_bracket);
103 delta=upper_bracket-lower_bracket;
104 }
105 c=value_with_derivatives(root, final_yprime, final_yprime2)-value; // compute initial values
106 b=*final_yprime; // make a local copy for readability
107 increment_evaluations();
108
109 // now, close in bracket on whichever side this still brackets
110 if(c*clower < 0.0) {
111 cupper=c;
112 upper_bracket=root;
113 } else {
114 clower=c;
115 lower_bracket=root;
116 }
117 // std::cout << "find_root_debug x, y, dx " << root << " " << c << " " << delta << std::endl;
118 }
119 return root;
120}
121
122/* def partial_integrals(self, xgrid):
123 Return the integrals of a function between the sampling points xgrid. The sum is the definite integral.
124 This method uses an exact integration of the polynomial which matches the values and derivatives at the
125 endpoints of a segment. Its error scales as h**6, if the input functions really are smooth.
126 This could very well be used as a stepper for adaptive Romberg integration.
127 For InterpolatingFunctions, it is likely that the Simpson's rule integrator is sufficient.
128 #the weights come from an exact mathematica solution to the 5th order polynomial with the given values & derivatives
129 #yint = (y0+y1)*dx/2 + dx^2*(yp0-yp1)/10 + dx^3 * (ypp0+ypp1)/120 )
130*/
131
132// the recursive part of the integrator is agressively designed to minimize copying of data... lots of pointers
133template <typename float_type> float_type c2_function<float_type>::integrate_step(c2_integrate_recur &rb) const
134{
135 struct c2_integrate_fblock *fbl[3]={rb.f0, rb.f1, rb.f2};
136 struct c2_integrate_fblock f1; // will hold new middle values
137 float_type retvals[2]={0.0,0.0};
138 float_type lr[2];
139
140 // std::cout << "entering with " << rb.f0->x << " " << rb.f1->x << " " << rb.f2->x << std::endl;
141
142 int depth=rb.depth; // save this from the recursion block
143 float_type abs_tol=rb.abs_tol; // this is the value we will pass down
144 float_type *rblr=rb.lr; // save pointer to our parent's lr[2] array since it will get trampled in recursion
145
146 if(!depth) {
147 switch(rb.derivs) {
148 case 0:
149 rb.eps_scale=0.1; rb.extrap_coef=16; break;
150 case 1:
151 rb.eps_scale=0.1; rb.extrap_coef=64; break;
152 case 2:
153 rb.eps_scale=0.02; rb.extrap_coef=1024; break;
154 default:
155 throw c2_exception("derivs must be 0, 1 or 2 in partial_integrals");
156 }
157
158 rb.extrap2=1.0/(rb.extrap_coef-1.0);
159 }
160
161 for (int i=0; i<(depth==0?1:2); i++) { // handle left and right intervals, but only left one for depth=0
162 struct c2_integrate_fblock *f0=fbl[i], *f2=fbl[i+1];
163 f1.x=0.5*(f0->x + f2->x); // center of interval
164 float_type dx=f2->x - f0->x;
165 float_type dx2 = 0.5*dx;
166 float_type total;
167
168 f1.y=value_with_derivatives(f1.x, &(f1.yp), &(f1.ypp));
169 increment_evaluations();
170
171 // check for underflow on step size, which prevents us from achieving specified accuracy.
172 if(std::abs(dx) < std::abs(f1.x)*rb.rel_tol) {
173 std::ostringstream outstr;
174 outstr << "Step size underflow in adaptive_partial_integrals at depth=" << depth << ", x= " << f1.x;
175 throw c2_exception(outstr.str().c_str());
176 }
177
178 if(!depth) { // top level, total has not been initialized yet
179 switch(rb.derivs) { // create estimate of next lower order for first try
180 case 0:
181 total=0.5*(f0->y+f2->y)*dx; break;
182 case 1:
183 total=(f0->y+4.0*f1.y+f2->y)*dx/6.0; break;
184 case 2:
185 total=( (14*f0->y + 32*f1.y + 14*f2->y) + dx * (f0->yp - f2->yp) ) * dx /60.; break;
186 default:
187 total=0.0; // just to suppress missing default warnings
188 }
189 } else total=rblr[i]; // otherwise, get it from previous level
190
191 float_type left, right;
192
193 switch(rb.derivs) {
194 case 2:
195 // use ninth-order estimates for each side, from full set of all values (!) (Thanks, Mathematica!)
196 left= ( ( (169*f0->ypp + 1024*f1.ypp - 41*f2->ypp)*dx2 +
197 (2727*f0->yp - 5040*f1.yp + 423*f2->yp) )*dx2 +
198 (17007*f0->y + 24576*f1.y - 1263*f2->y) )* (dx2/40320.0);
199 right= ( ( (169*f2->ypp + 1024*f1.ypp - 41*f0->ypp)*dx2 -
200 (2727*f2->yp - 5040*f1.yp + 423*f0->yp) )*dx2 +
201 (17007*f2->y + 24576*f1.y - 1263*f0->y) )* (dx2/40320.0);
202 // std::cout << f0->x << " " << f1.x << " " << f2->x << std::endl ;
203 // std::cout << f0->y << " " << f1.y << " " << f2->y << " " << left << " " << right << " " << total << std::endl ;
204 break;
205 case 1:
206 left= ( (202*f0->y + 256*f1.y + 22*f2->y) + dx*(13*f0->yp - 40*f1.yp - 3*f2->yp) ) * dx /960.;
207 right= ( (202*f2->y + 256*f1.y + 22*f0->y) - dx*(13*f2->yp - 40*f1.yp - 3*f0->yp) ) * dx /960.;
208 break;
209 case 0:
210 left= (5*f0->y + 8*f1.y - f2->y)*dx/24.;
211 right= (5*f2->y + 8*f1.y - f0->y)*dx/24.;
212 break;
213 default:
214 left=right=0.0; // suppress warnings about missing default
215 break;
216 }
217
218 lr[0]= left; // left interval
219 lr[1]= right; // right interval
220 float_type lrsum=left+right;
221
222 float_type eps=std::abs(total-lrsum)*rb.eps_scale;
223 if(rb.extrapolate) eps*=rb.eps_scale;
224
225 if(!rb.adapt || eps < abs_tol || eps < std::abs(total)*rb.rel_tol) {
226 if(depth==0 || !rb.extrapolate) retvals[i]=lrsum;
227 else {
228 retvals[i]=(rb.extrap_coef*lrsum - total)*rb.extrap2;
229 // std::cout << "extrapolating " << lrsum << " " << total << " " << retvals[i] << std::endl;
230
231 }
232 } else {
233 rb.depth=depth+1; // increment depth counter
234 rb.lr=lr; // point to our left-right values array for recursion
235 rb.abs_tol=abs_tol*0.5; // each half has half the error budget
236 rb.f0=f0; rb.f1=&f1; rb.f2=f2; // insert pointers to data into our recursion block
237 // std::cout << "recurring with " << f0->x << " " << f1.x << " " << f2->x << std::endl ;
238 retvals[i]=integrate_step(rb); // and recur
239 }
240 }
241 return retvals[0]+retvals[1];
242}
243
244template <typename float_type> bool c2_function<float_type>::check_monotonicity(
245 const std::vector<float_type> &data, const char message[]) throw(c2_exception)
246{
247 size_t np=data.size();
248 if(np < 2) return false; // one point has no direction!
249
250 bool rev=(data[1] < data[0]); // which way do data point?
251 size_t i;
252
253 if(!rev) for(i = 2; i < np && (data[i-1] < data[i]) ; i++);
254 else for(i = 2; i < np &&(data[i-1] > data[i]) ; i++);
255
256 if(i != np) throw c2_exception(message);
257
258 return rev;
259}
260
261template <typename float_type> void c2_function<float_type>::set_sampling_grid(const std::vector<float_type> &grid) throw(c2_exception)
262{
263 bool rev=check_monotonicity(grid, "set_sampling_grid: sampling grid not monotonic");
264
265 if(!sampling_grid || no_overwrite_grid) sampling_grid=new std::vector<float_type>;
266 (*sampling_grid)=grid; no_overwrite_grid=0;
267
268 if(rev) std::reverse(sampling_grid->begin(), sampling_grid->end()); // make it increasing
269}
270
271template <typename float_type> std::vector<float_type> &c2_function<float_type>::get_sampling_grid(float_type xmin, float_type xmax) const
272{
273 std::vector<float_type> *result=new std::vector<float_type>;
274
275 if( !(sampling_grid) || !(sampling_grid->size()) || (xmax <= sampling_grid->front()) || (xmin >= sampling_grid->back()) ) {
276 // nothing is known about the function in this region, return xmin and xmax
277 result->push_back(xmin);
278 result->push_back(xmax);
279 } else {
280 std::vector<float_type> &sg=*sampling_grid; // just a shortcut
281 int np=sg.size();
282 int klo=0, khi=np-1, firstindex=0, lastindex=np-1;
283
284 result->push_back(xmin);
285
286 if(xmin > sg.front() ) {
287 // hunt through table for position bracketing starting point
288 while(khi-klo > 1) {
289 int km=(khi+klo)/2;
290 if(sg[km] > xmin) khi=km;
291 else klo=km;
292 }
293 khi=klo+1;
294 // khi now points to first point definitively beyond our first point, or last point of array
295 firstindex=khi;
296 khi=np-1; // restart upper end of search
297 }
298
299 if(xmax < sg.back()) {
300 // hunt through table for position bracketing starting point
301 while(khi-klo > 1) {
302 int km=(khi+klo)/2;
303 if(sg[km] > xmax) khi=km;
304 else klo=km;
305 }
306 khi=klo+1;
307 // khi now points to first point definitively beyond our last point, or last point of array
308 lastindex=klo;
309 }
310
311 int initsize=result->size();
312 result->resize(initsize+(lastindex-firstindex+2));
313 std::copy(sg.begin()+firstindex, sg.begin()+lastindex+1, result->begin()+initsize);
314 result->back()=xmax;
315
316 // this is the unrefined sampling grid... now check for very close points on front & back and fix if needed.
317 preen_sampling_grid(result);
318 }
319 return *result;
320}
321
322template <typename float_type> void c2_function<float_type>::preen_sampling_grid(std::vector<float_type> *result) const
323{
324 // this is the unrefined sampling grid... now check for very close points on front & back and fix if needed.
325 if(result->size() > 2) { // may be able to prune dangerously close points near the ends if there are at least 3 points
326 bool deleteit=false;
327 float_type x0=(*result)[0], x1=(*result)[1];
328 float_type dx1=x1-x0;
329
330 float_type ftol=10.0*(std::numeric_limits<float_type>::epsilon()*(std::abs(x0)+std::abs(x1))+std::numeric_limits<float_type>::min());
331 if(dx1 < ftol) deleteit=true;
332 float_type dx2=(*result)[2]-x0;
333 if(dx1/dx2 < 0.1) deleteit=true; // endpoint is very close to internal interesting point
334
335 if(deleteit) result->erase(result->begin()+1); // delete redundant interesting point
336 }
337
338 if(result->size() > 2) { // may be able to prune dangerously close points near the ends if there are at least 3 points
339 bool deleteit=false;
340 int pos=result->size()-3;
341 float_type x0=(*result)[pos+1], x1=(*result)[pos+2];
342 float_type dx1=x1-x0;
343
344 float_type ftol=10.0*(std::numeric_limits<float_type>::epsilon()*(std::abs(x0)+std::abs(x1))+std::numeric_limits<float_type>::min());
345 if(dx1 < ftol) deleteit=true;
346 float_type dx2=x1-(*result)[pos];
347 if(dx1/dx2 < 0.1) deleteit=true; // endpoint is very close to internal interesting point
348
349 if(deleteit) result->erase(result->end()-2); // delete redundant interesting point
350 }
351}
352
353template <typename float_type> std::vector<float_type> &c2_function<float_type>::
354 refine_sampling_grid(const std::vector<float_type> &grid, size_t refinement) const
355{
356 size_t np=grid.size();
357 size_t count=(np-1)*refinement + 1;
358 float_type dxscale=1.0/refinement;
359
360 std::vector<float_type> *result=new std::vector<float_type>(count);
361
362 for(size_t i=0; i<(np-1); i++) {
363 float_type x=grid[i];
364 float_type dx=(grid[i+1]-x)*dxscale;
365 for(size_t j=0; j<refinement; j++, x+=dx) (*result)[i*refinement+j]=x;
366 }
367 (*result)[count-1]=grid.back();
368 return *result;
369}
370
371template <typename float_type> float_type c2_function<float_type>::integral(float_type xmin, float_type xmax, std::vector<float_type> *partials,
372 float_type abs_tol, float_type rel_tol, int derivs, bool adapt, bool extrapolate) const
373{
374 std::vector<float_type> &grid=get_sampling_grid(xmin, xmax);
375 float_type intg=partial_integrals(grid, partials, abs_tol, rel_tol, adapt, extrapolate);
376 delete &grid;
377 return intg;
378}
379
380template <typename float_type> c2_function<float_type> &c2_function<float_type>::normalized_function(float_type xmin, float_type xmax, float_type norm)
381{
382 float_type intg=integral(xmin, xmax);
383 return *new c2_scaled_function<float_type>(*this, norm/intg);
384}
385
386template <typename float_type> c2_function<float_type> &c2_function<float_type>::square_normalized_function(float_type xmin, float_type xmax, float_type norm)
387{
388 c2_quadratic<float_type> q(0., 0., 0., 1.);
389 c2_composed_function<float_type> mesquared(q,*this);
390
391 std::vector<float_type> grid(get_sampling_grid(xmin, xmax));
392 float_type intg=mesquared.partial_integrals(grid);
393
394 return *new c2_scaled_function<float_type>(*this, std::sqrt(norm/intg));
395}
396
397template <typename float_type> c2_function<float_type> &c2_function<float_type>::square_normalized_function(
398 float_type xmin, float_type xmax, const c2_function<float_type> &weight, float_type norm)
399{
400 c2_quadratic<float_type> q(0., 0., 0., 1.);
401 c2_composed_function<float_type> mesquared(q,*this);
402 c2_product<float_type> weighted(mesquared, weight);
403
404 std::vector<float_type> grid(get_sampling_grid(xmin, xmax));
405 float_type intg=weighted.partial_integrals(grid);
406
407 return *new c2_scaled_function<float_type>(*this, std::sqrt(norm/intg));
408}
409
410template <typename float_type> float_type c2_function<float_type>::partial_integrals(
411 std::vector<float_type> xgrid, std::vector<float_type> *partials,
412 float_type abs_tol, float_type rel_tol, int derivs, bool adapt, bool extrapolate) const
413{
414 int np=xgrid.size();
415
416 struct c2_integrate_fblock f0, f2;
417 struct c2_integrate_recur rb;
418 rb.rel_tol=rel_tol;
419 rb.extrapolate=extrapolate;
420 rb.adapt=adapt;
421 rb.derivs=derivs;
422
423 reset_evaluations(); // counter returns with total evaluations needed for this integral
424
425 if(partials) partials->resize(np-1);
426
427 float_type sum=0.0;
428
429 f2.x=xgrid[0];
430 f2.y=value_with_derivatives(f2.x, &f2.yp, &f2.ypp);
431 increment_evaluations();
432
433 for(int i=0; i<np-1; i++) {
434 f0=f2; // copy upper bound to lower before computing new upper bound
435
436 f2.x=xgrid[i+1];
437 f2.y=value_with_derivatives(f2.x, &f2.yp, &f2.ypp);
438 increment_evaluations();
439
440 rb.depth=0;
441 rb.abs_tol=abs_tol;
442 rb.f0=&f0; rb.f1=&f2; rb.f2=&f2; // we are really only using the left half for the top level
443 rb.lr=0; // pointer is meaningless; will be filled in in recursion
444 float_type ps=integrate_step(rb);
445 sum+=ps;
446 if(partials) (*partials)[i]=ps;
447 }
448 return sum;
449}
450
451// declare singleton functions for most common c2_function instances
452#define c2_singleton(X) template <typename float_type> const c2_##X<float_type> c2_##X<float_type>::X=c2_##X();
453c2_singleton(sin)
454c2_singleton(cos)
455c2_singleton(tan)
456c2_singleton(log)
457c2_singleton(exp)
458c2_singleton(sqrt)
459c2_singleton(identity)
460
461// reciprocal is actually parametric (a/x), but make singleton 1/x
462template <typename float_type> const c2_recip<float_type> c2_recip<float_type>::recip=c2_recip(1.0);
463
464#undef c2_singleton
465
466// generate a sampling grid at points separated by dx=5, which is intentionally
467// incommensurate with pi and 2*pi so grid errors are somewhat randomized
468template <typename float_type> std::vector<float_type> &c2_sin<float_type>::get_sampling_grid(float_type xmin, float_type xmax)
469{
470 std::vector<float_type> *result=new std::vector<float_type>;
471
472 for(; xmin < xmax; xmin+=5.0) result->push_back(xmin);
473 result->push_back(xmax);
474 this->preen_sampling_grid(result);
475 return *result;
476}
477
478template <typename float_type> float_type Identity(float_type x) { return x; } // a useful function
479template <typename float_type> float_type f_one(float_type) { return 1.0; } // the first derivative of identity
480template <typename float_type> float_type f_zero(float_type) { return 0.0; } // the second derivative of identity
481
482// The constructor
483template <typename float_type> void interpolating_function<float_type>::init(
484 const std::vector<float_type> &x, const std::vector<float_type> &f,
485 bool lowerSlopeNatural, float_type lowerSlope,
486 bool upperSlopeNatural, float_type upperSlope,
487 float_type (*inputXConversion)(float_type),
488 float_type (*inputXConversionPrime)(float_type),
489 float_type (*inputXConversionDPrime)(float_type),
490 float_type (*inputYConversion)(float_type),
491 float_type (*inputYConversionPrime)(float_type),
492 float_type (*inputYConversionDPrime)(float_type),
493 float_type (*outputYConversion)(float_type)
494 ) throw(c2_exception)
495{
496 X= x;
497 F= f;
498
499 // Xraw is useful in some of the arithmetic operations between interpolating functions
500 Xraw=x;
501
502 set_domain(std::min(Xraw.front(), Xraw.back()),std::max(Xraw.front(), Xraw.back()));
503
504 fXin=inputXConversion;
505 fXinPrime=inputXConversionPrime;
506 fXinDPrime=inputXConversionDPrime;
507 fYin=inputYConversion;
508 fYinPrime=inputYConversionPrime;
509 fYinDPrime=inputYConversionDPrime;
510 fYout=outputYConversion;
511
512 if(x.size() != f.size()) {
513 throw c2_exception("interpolating_function::init() -- x & y inputs are of different size");
514 }
515
516 size_t np=X.size(); // they are the same now, so lets take a short cut
517
518 if(np < 2) {
519 throw c2_exception("interpolating_function::init() -- input < 2 elements ");
520 }
521
522 bool xraw_rev=check_monotonicity(Xraw,
523 "interpolating_function::init() non-monotonic raw x input"); // which way does raw X point? sampling grid MUST be increasing
524
525 if(!xraw_rev) { // we can use pointer to raw X values if they are in the right order
526 set_sampling_grid_pointer(Xraw); // our intial grid of x values is certainly a good guess for 'interesting' points
527 } else {
528 set_sampling_grid(Xraw); // make a copy of it, and assure it is in right order
529 }
530
531 if(fXin) { // check if X scale is nonlinear, and if so, do transform
532 if(!lowerSlopeNatural) lowerSlope /= fXinPrime(X[0]);
533 if(!upperSlopeNatural) upperSlope /= fXinPrime(X[np-1]);
534 for(size_t i=0; i<np; i++) X[i]=fXin(X[i]);
535 } else {
536 fXin=Identity<float_type>;
537 fXinPrime=f_one<float_type>;
538 fXinDPrime=f_zero<float_type>;
539 }
540
541 if(inputYConversion) { // check if Y scale is nonlinear, and if so, do transform
542 if(!lowerSlopeNatural) lowerSlope *= fYinPrime(F[0]);
543 if(!upperSlopeNatural) upperSlope *= fYinPrime(F[np-1]);
544 for(size_t i=0; i<np; i++) F[i]=inputYConversion(F[i]);
545 } else {
546 fYin=Identity<float_type>;
547 fYinPrime=f_one<float_type>;
548 fYinDPrime=f_zero<float_type>;
549 fYout=Identity<float_type>;
550 }
551
552 xInverted=check_monotonicity(X,
553 "interpolating_function::init() non-monotonic transformed x input");
554
555 // construct spline tables here.
556 // this code is a re-translation of the pythonlabtools spline algorithm from pythonlabtools.sourceforge.net
557
558 std::vector<float_type> u(np), dy(np-1), dx(np-1), dxi(np-1), dx2i(np-2), siga(np-2), dydx(np-1);
559
560 std::transform(X.begin()+1, X.end(), X.begin(), dx.begin(), std::minus<float_type>() ); // dx=X[1:] - X [:-1]
561 for(size_t i=0; i<dxi.size(); i++) dxi[i]=1.0/dx[i]; // dxi = 1/dx
562 for(size_t i=0; i<dx2i.size(); i++) dx2i[i]=1.0/(X[i+2]-X[i]);
563
564 std::transform(F.begin()+1, F.end(), F.begin(), dy.begin(), std::minus<float_type>() ); // dy = F[i+1]-F[i]
565 std::transform(dx2i.begin(), dx2i.end(), dx.begin(), siga.begin(), std::multiplies<float_type>()); // siga = dx[:-1]*dx2i
566 std::transform(dxi.begin(), dxi.end(), dy.begin(), dydx.begin(), std::multiplies<float_type>()); // dydx=dy/dx
567
568 // u[i]=(y[i+1]-y[i])/float(x[i+1]-x[i]) - (y[i]-y[i-1])/float(x[i]-x[i-1])
569 std::transform(dydx.begin()+1, dydx.end(), dydx.begin(), u.begin()+1, std::minus<float_type>() ); // incomplete rendering of u = dydx[1:]-dydx[:-1]
570
571 y2.resize(np,0.0);
572
573 if(lowerSlopeNatural) {
574 y2[0]=u[0]=0.0;
575 } else {
576 y2[0]= -0.5;
577 u[0]=(3.0*dxi[0])*(dy[0]*dxi[0] -lowerSlope);
578 }
579
580 for(size_t i=1; i < np -1; i++) { // the inner loop
581 float_type sig=siga[i-1];
582 float_type p=sig*y2[i-1]+2.0;
583 y2[i]=(sig-1.0)/p;
584 u[i]=(6.0*u[i]*dx2i[i-1] - sig*u[i-1])/p;
585 }
586
587 float_type qn, un;
588
589 if(upperSlopeNatural) {
590 qn=un=0.0;
591 } else {
592 qn= 0.5;
593 un=(3.0*dxi[dxi.size()-1])*(upperSlope- dy[dy.size()-1]*dxi[dxi.size()-1] );
594 }
595
596 y2[np-1]=(un-qn*u[np-2])/(qn*y2[np-2]+1.0);
597 for (size_t k=np-1; k != 0; k--) y2[k-1]=y2[k-1]*y2[k]+u[k-1];
598
599 lastKLow=-1; // flag new X search required for next evaluation
600}
601
602// This function is the reason for this class to exist
603// it computes the interpolated function, and (if requested) its proper first and second derivatives including all coordinate transforms
604template <typename float_type> float_type interpolating_function<float_type>::value_with_derivatives(
605 float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
606{
607 if(x < this->xmin() || x > this->xmax()) {
608 std::ostringstream outstr;
609 outstr << "Interpolating function argument " << x << " out of range " << this->xmin() << " -- " << this ->xmax() << ": bailing";
610 throw c2_exception(outstr.str().c_str());
611 }
612
613 float_type xraw=x;
614
615 // template here is impossible! if(fXin && fXin != (Identity<float_type>) )
616 x=fXin(x); // save time by explicitly testing for identity function here
617
618 int klo=0, khi=X.size()-1;
619
620 if(!xInverted) { // select search depending on whether transformed X is increasing or decreasing
621 if(lastKLow >=0 && (X[lastKLow] <= x) && (X[lastKLow+1] >= x) ) { // already bracketed
622 klo=lastKLow;
623 } else if(lastKLow >=0 && (X[lastKLow+1] <= x) && (X[lastKLow+2] > x)) { // in next bracket to the right
624 klo=lastKLow+1;
625 } else if(lastKLow > 0 && (X[lastKLow-1] <= x) && (X[lastKLow] > x)) { // in next bracket to the left
626 klo=lastKLow-1;
627 } else { // not bracketed, not close, start over
628 // search for new KLow
629 while(khi-klo > 1) {
630 int km=(khi+klo)/2;
631 if(X[km] > x) khi=km;
632 else klo=km;
633 }
634 }
635 } else {
636 if(lastKLow >=0 && (X[lastKLow] >= x) && (X[lastKLow+1] <= x) ) { // already bracketed
637 klo=lastKLow;
638 } else if(lastKLow >=0 && (X[lastKLow+1] >= x) && (X[lastKLow+2] < x)) { // in next bracket to the right
639 klo=lastKLow+1;
640 } else if(lastKLow > 0 && (X[lastKLow-1] >= x) && (X[lastKLow] < x)) { // in next bracket to the left
641 klo=lastKLow-1;
642 } else { // not bracketed, not close, start over
643 // search for new KLow
644 while(khi-klo > 1) {
645 int km=(khi+klo)/2;
646 if(X[km] < x) khi=km;
647 else klo=km;
648 }
649 }
650 }
651
652 khi=klo+1;
653 lastKLow=klo;
654
655 float_type h=X[khi]-X[klo];
656
657 float_type a=(X[khi]-x)/h;
658 float_type b=1.0-a;
659 float_type ylo=F[klo], yhi=F[khi], y2lo=y2[klo], y2hi=y2[khi];
660 float_type y=a*ylo+b*yhi+((a*a*a-a)*y2lo+(b*b*b-b)*y2hi)*(h*h)/6.0;
661
662 // template here is impossible! if(fYin && fYin != Identity)
663 y=fYout(y); // save time by explicitly testing for identity function here
664
665 if(yprime || yprime2) {
666 float_type fpi=1.0/fYinPrime(y);
667 float_type gp=fXinPrime(xraw);
668 float_type yp0=(yhi-ylo)/h+((3*b*b-1)*y2hi-(3*a*a-1)*y2lo)*h/6.0; // the derivative in interpolating table coordinates
669
670 // from Mathematica Dt[InverseFunction[f][y[g[x]]], x]
671 if(yprime) *yprime=gp*yp0*fpi; // the real derivative of the inverse transformed output
672 if(yprime2) {
673 float_type ypp0=b*y2hi+a*y2lo;
674 float_type fpp=fYinDPrime(y);
675 float_type gpp=fXinDPrime(xraw);
676 // also from Mathematica Dt[InverseFunction[f][y[g[x]]], {x,2}]
677 if(yprime2) *yprime2=(gp*gp*ypp0 + yp0*gpp - gp*gp*yp0*yp0*fpp*fpi*fpi)*fpi;
678 }
679 }
680
681 return y;
682}
683
684template <typename float_type> void interpolating_function<float_type>::set_lower_extrapolation(float_type bound)
685{
686 int kl = 0 ;
687 int kh=kl+1;
688 float_type xx=fXin(bound);
689 float_type h0=X[kh]-X[kl];
690 float_type h1=xx-X[kl];
691 float_type yextrap=F[kl]+((F[kh]-F[kl])/h0 - h0*(y2[kl]+2.0*y2[kh])/6.0)*h1+y2[kl]*h1*h1/2.0;
692
693 X.insert(X.begin(), xx);
694 F.insert(F.begin(), yextrap);
695 y2.insert(y2.begin(), y2.front()); // duplicate first or last element
696 Xraw.insert(Xraw.begin(), bound);
697 if (bound < this->fXMin) this->fXMin=bound; // check for reversed data
698 else this->fXMax=bound;
699
700 //printf("%10.4f %10.4f %10.4f %10.4f %10.4f\n", bound, xx, h0, h1, yextrap);
701 //for(int i=0; i<X.size(); i++) printf("%4d %10.4f %10.4f %10.4f %10.4f \n", i, Xraw[i], X[i], F[i], y2[i]);
702}
703
704template <typename float_type> void interpolating_function<float_type>::set_upper_extrapolation(float_type bound)
705{
706 int kl = X.size()-2 ;
707 int kh=kl+1;
708 float_type xx=fXin(bound);
709 float_type h0=X[kh]-X[kl];
710 float_type h1=xx-X[kl];
711 float_type yextrap=F[kl]+((F[kh]-F[kl])/h0 - h0*(y2[kl]+2.0*y2[kh])/6.0)*h1+y2[kl]*h1*h1/2.0;
712
713 X.insert(X.end(), xx);
714 F.insert(F.end(), yextrap);
715 y2.insert(y2.end(), y2.back()); // duplicate first or last element
716 Xraw.insert(Xraw.end(), bound);
717 if (bound < this->fXMin) this->fXMin=bound; // check for reversed data
718 else this->fXMax=bound;
719 //printf("%10.4f %10.4f %10.4f %10.4f %10.4f\n", bound, xx, h0, h1, yextrap);
720 //for(int i=0; i<X.size(); i++) printf("%4d %10.4f %10.4f %10.4f %10.4f \n", i, Xraw[i], X[i], F[i], y2[i]);
721}
722
723// move derivatives into our internal coordinates (use splint to go the other way!)
724template <typename float_type> void interpolating_function<float_type>::localize_derivatives(
725 float_type xraw, float_type y, float_type yp, float_type ypp, float_type *y0, float_type *yprime, float_type *yprime2) const
726{
727 float_type fp=fYinPrime(y);
728 float_type gp=fXinPrime(xraw);
729 float_type fpp=fYinDPrime(y);
730 float_type gpp=fXinDPrime(xraw);
731
732 if(y0) *y0=fYin(y);
733 if(yprime) *yprime=yp*fp/gp; // Mathematica Dt[f[y[InverseFunction[g][x]]], x]
734 if(yprime2) *yprime2=( yp*yp*fpp - fp*yp*gpp/gp + fp*ypp )/(gp*gp) ; // Mathematica Dt[f[y[InverseFunction[g][x]]], {x,2}]
735}
736
737// return a new interpolating_function which is the unary function of an existing interpolating_function
738// can also be used to generate a resampling of another c2_function on a different grid
739// by creating a=interpolating_function(x,x)
740// and doing b=a.unary_operator(c) where c is a c2_function (probably another interpolating_function)
741
742template <typename float_type> interpolating_function<float_type>&
743 interpolating_function<float_type>::unary_operator(const c2_function<float_type> &source) const
744{
745 size_t np=X.size();
746 std::vector<float_type>yv(np);
747 c2_composed_function<float_type> comp(source, *this);
748 float_type yp0, yp1, ypp;
749
750 for(size_t i=0; i<np; i++) {
751 yv[i]=source(fYout(F[i])); // copy pointwise the function of our data values
752 }
753
754 comp(Xraw.front(), &yp0, &ypp); // get derivative at front
755 comp(Xraw.back(), &yp1, &ypp); // get derivative at back
756
757 return *new interpolating_function(Xraw, yv, false, yp0, false, yp1,
758 fXin, fXinPrime, fXinDPrime,
759 fYin, fYinPrime, fYinDPrime, fYout);
760}
761
762template <typename float_type> void
763interpolating_function<float_type>::get_data(std::vector<float_type> &xvals, std::vector<float_type> &yvals) const throw()
764{
765
766 xvals=Xraw;
767 yvals.resize(F.size());
768
769 for(size_t i=0; i<F.size(); i++) yvals[i]=fYout(F[i]);
770}
771
772template <typename float_type> interpolating_function<float_type> &
773 interpolating_function<float_type>::binary_operator(const c2_function<float_type> &rhs,
774 c2_binary_function<float_type> *combining_stub) const
775{
776 size_t np=X.size();
777 std::vector<float_type> yv(np);
778 c2_constant<float_type> fval;
779 c2_constant<float_type> yval;
780 float_type yp0, yp1, ypp;
781
782 for(size_t i=0; i<np; i++) {
783 fval.reset(fYout(F[i])); // update the constant function pointwise
784 yval.reset(rhs(Xraw[i]));
785 yv[i]=(*combining_stub).combine(fval, yval, Xraw[i], (float_type *)0, (float_type *)0); // compute rhs & combine without derivatives
786 }
787
788 (*combining_stub).combine(*this, rhs, Xraw.front(), &yp0, &ypp); // get derivative at front
789 (*combining_stub).combine(*this, rhs, Xraw.back(), &yp1, &ypp); // get derivative at back
790
791 delete combining_stub;
792
793 return *new interpolating_function(Xraw, yv, false, yp0, false, yp1,
794 fXin, fXinPrime, fXinDPrime,
795 fYin, fYinPrime, fYinDPrime, fYout);
796}
797
798template <typename float_type> float_type c2_f_logprime(float_type x) { return 1.0/x; } // the derivative of log(x)
799template <typename float_type> float_type c2_f_logprime2(float_type x) { return -1.0/(x*x); } // the second derivative of log(x)
800
801template <typename float_type> log_lin_interpolating_function<float_type>::log_lin_interpolating_function(
802 const std::vector<float_type> &x, const std::vector<float_type> &f,
803 bool lowerSlopeNatural, float_type lowerSlope,
804 bool upperSlopeNatural, float_type upperSlope)
805 : interpolating_function<float_type>()
806{
807 init(x, f, lowerSlopeNatural, lowerSlope, upperSlopeNatural, upperSlope,
808 (float_type (*)(float_type)) (std::log), c2_f_logprime, c2_f_logprime2, 0, 0, 0, 0);
809}
810
811template <typename float_type> lin_log_interpolating_function<float_type>::lin_log_interpolating_function(
812 const std::vector<float_type> &x, const std::vector<float_type> &f,
813 bool lowerSlopeNatural, float_type lowerSlope,
814 bool upperSlopeNatural, float_type upperSlope)
815 : interpolating_function<float_type>()
816{
817 init(x, f, lowerSlopeNatural, lowerSlope, upperSlopeNatural, upperSlope,
818 0, 0, 0,
819 (float_type (*)(float_type)) (std::log), c2_f_logprime, c2_f_logprime2,
820 (float_type (*)(float_type)) (std::exp) );
821}
822
823template <typename float_type> log_log_interpolating_function<float_type>::log_log_interpolating_function(
824 const std::vector<float_type> &x, const std::vector<float_type> &f,
825 bool lowerSlopeNatural, float_type lowerSlope,
826 bool upperSlopeNatural, float_type upperSlope)
827 : interpolating_function<float_type>()
828{
829 init(x, f, lowerSlopeNatural, lowerSlope, upperSlopeNatural, upperSlope,
830 (float_type (*)(float_type)) (std::log), c2_f_logprime, c2_f_logprime2,
831 (float_type (*)(float_type)) (std::log), c2_f_logprime, c2_f_logprime2,
832 (float_type (*)(float_type)) (std::exp) );
833}
834
835template <typename float_type> float_type c2_f_recip(float_type x) { return 1.0/x; }
836template <typename float_type> float_type c2_f_recipprime(float_type x) { return -1.0/(x*x); } // the derivative of 1/x
837template <typename float_type> float_type c2_f_recipprime2(float_type x) { return 2.0/(x*x*x); } // the second derivative of 1/x
838
839template <typename float_type> arrhenius_interpolating_function<float_type>::arrhenius_interpolating_function(
840 const std::vector<float_type> &x, const std::vector<float_type> &f,
841 bool lowerSlopeNatural, float_type lowerSlope,
842 bool upperSlopeNatural, float_type upperSlope)
843 : interpolating_function<float_type>()
844{
845 init(x, f, lowerSlopeNatural, lowerSlope, upperSlopeNatural, upperSlope,
846 c2_f_recip, c2_f_recipprime, c2_f_recipprime2,
847 (float_type (*)(float_type)) (std::log), c2_f_logprime, c2_f_logprime2,
848 (float_type (*)(float_type)) (std::exp) );
849}
850
851template <typename float_type> c2_inverse_function<float_type>::c2_inverse_function(const c2_function<float_type> &source)
852 : c2_plugin_function<float_type>(source)
853{
854 float_type l=source.xmin();
855 float_type r=source.xmax();
856 start_hint=(l+r)*0.5; // guess that we start in the middle
857 // compute our domain assuming the function is monotonic so its values on its domain boundaries are our domain
858 float_type ly=source(l);
859 float_type ry=source(r);
860 if (ly > ry) {
861 float_type t=ly; ly=ry; ry=t;
862 }
863 set_domain(ly, ry);
864}
865
866template <typename float_type> float_type c2_inverse_function<float_type>::value_with_derivatives(
867 float_type x, float_type *yprime, float_type *yprime2
868 ) const throw(c2_exception)
869{
870 float_type l=this->func->xmin();
871 float_type r=this->func->xmax();
872 float_type yp, ypp;
873 float_type y=this->func->find_root(l, r, get_start_hint(x), x, 0, &yp, &ypp);
874 start_hint=y;
875 if(yprime) *yprime=1.0/yp;
876 if(yprime2) *yprime2=-ypp/(yp*yp*yp);
877 return y;
878}
879
880//accumulated_histogram starts with binned data, generates the integral, and generates a piecewise linear interpolating_function
881//If drop_zeros is true, it merges empty bins together before integration
882//Note that the resulting interpolating_function is guaranteed to be increasing (if drop_zeros is false)
883// or stricly increasing (if drop_zeros is true)
884//If inverse_function is true, it drop zeros, integrates, and returns the inverse function which is useful
885// for random number generation based on the input distribution.
886//If normalize is true, the big end of the integral is scaled to 1.
887//If the data are passed in reverse order (large X first), the integral is carried out from the big end,
888// and then the data are reversed to the result in in increasing X order.
889template <typename float_type> accumulated_histogram<float_type>::accumulated_histogram(
890 const std::vector<float_type>binedges, const std::vector<float_type> binheights,
891 bool normalize, bool inverse_function, bool drop_zeros)
892{
893
894 int np=binheights.size();
895
896 std::vector<float_type> be, bh;
897 if(drop_zeros || inverse_function) { //inverse functions cannot have any zero bins or they have vertical sections
898 if(binheights[0] || !inverse_function) { // conserve lower x bound if not an inverse function
899 be.push_back(binedges[0]);
900 bh.push_back(binheights[0]);
901 }
902 for(int i=1; i<np-1; i++) {
903 if(binheights[i]) {
904 be.push_back(binedges[i]);
905 bh.push_back(binheights[i]);
906 }
907 }
908 if(binheights[np-1] || !inverse_function) {
909 bh.push_back(binheights[np-1]);
910 be.push_back(binedges[np-1]);
911 be.push_back(binedges[np]); // push both sides of the last bin if needed
912 }
913 np=bh.size(); // set np to compressed size of bin array
914 } else {
915 be=binedges;
916 bh=binheights;
917 }
918 std::vector<float_type> cum(np+1, 0.0);
919 for(int i=1; i<=np; i++) cum[i]=bh[i]*(be[i]-be[i-1])+cum[i-1]; // accumulate bins, leaving bin 0 as 0
920 if(be[1] < be[0]) { // if bins passed in backwards, reverse them
921 std::reverse(be.begin(), be.end());
922 std::reverse(cum.begin(), cum.end());
923 for(unsigned int i=0; i<cum.size(); i++) cum[i]*=-1; // flip sign on reversed data
924 }
925 if(normalize) {
926 float_type m=1.0/std::max(cum[0], cum[np]);
927 for(int i=0; i<=np; i++) cum[i]*=m;
928 }
929 if(inverse_function) interpolating_function<float_type>(cum, be); // use cum as x axis in inverse function
930 else interpolating_function<float_type>(be, cum); // else use lower bin edge as x axis
931 std::fill(this->y2.begin(), this->y2.end(), 0.0); // clear second derivatives, to we are piecewise linear
932}
933
934template <typename float_type> c2_piecewise_function<float_type>::c2_piecewise_function()
935: c2_function<float_type>(), lastKLow(-1)
936{
937 this->sampling_grid=new std::vector<float_type>; // this always has a smapling grid
938}
939
940template <typename float_type> c2_piecewise_function<float_type>::~c2_piecewise_function()
941{
942 size_t np=functions.size();
943 for(size_t i=0; i<np; i++) if(owns[i]) delete functions[i];
944}
945
946template <typename float_type> float_type c2_piecewise_function<float_type>::value_with_derivatives(
947 float_type x, float_type *yprime, float_type *yprime2
948 ) const throw(c2_exception)
949{
950
951 size_t np=functions.size();
952 if(!np) throw c2_exception("attempting to evaluate an empty piecewise function");
953
954 if(x < this->xmin() || x > this->xmax()) {
955 std::ostringstream outstr;
956 outstr << "piecewise function argument " << x << " out of range " << this->xmin() << " -- " << this->xmax();
957 throw c2_exception(outstr.str().c_str());
958 }
959
960 int klo=0;
961
962 if(lastKLow >= 0 && functions[lastKLow]->xmin() <= x && functions[lastKLow]->xmax() > x) {
963 klo=lastKLow;
964 } else {
965 int khi=np;
966 while(khi-klo > 1) {
967 int km=(khi+klo)/2;
968 if(functions[km]->xmin() > x) khi=km;
969 else klo=km;
970 }
971 }
972 lastKLow=klo;
973 return functions[klo]->value_with_derivatives(x, yprime, yprime2);
974}
975
976template <typename float_type> void c2_piecewise_function<float_type>::append_function(
977 c2_function<float_type> &func, bool pass_ownership) throw(c2_exception)
978{
979 if(functions.size()) { // check whether there are any gaps to fill, etc.
980 c2_function<float_type> &tail=*(functions.back());
981 float_type x0=tail.xmax();
982 float_type x1=func.xmin();
983 if(x0 < x1) {
984 // must insert a connector if x0 < x1
985 float_type y0=tail(x0);
986 float_type y1=func(x1);
987 c2_function<float_type> *connector=new c2_linear<float_type>(x0, y0, (y1-y0)/(x1-x0));
988 connector->set_domain(x0,x1);
989 functions.push_back(connector);
990 owns.push_back(true);
991 this->sampling_grid->push_back(x1);
992 } else if(x0>x1) throw c2_exception("function domains not increasing in c2_piecewise_function");
993 }
994 functions.push_back(&func);
995 owns.push_back(pass_ownership);
996 // extend our domain to include all known functions
997 this->set_domain(functions.front()->xmin(), functions.back()->xmax());
998 // extend our sampling grid with the new function's grid, with the first point dropped to avoid duplicates
999 std::vector<float_type> &newgrid=func.get_sampling_grid(func.xmin(), func.xmax());
1000 this->sampling_grid->insert(this->sampling_grid->end(), newgrid.begin()+1, newgrid.end());
1001 delete &newgrid;
1002}
1003
1004template <typename float_type> c2_connector_function<float_type>::c2_connector_function(
1005 const c2_function<float_type> &f1, const c2_function<float_type> &f2, float_type x0, float_type x2,
1006 bool auto_center, float_type y1)
1007
1008: c2_function<float_type>()
1009{
1010 float_type y0, yp0, ypp0, y2, yp2, ypp2;
1011 fdx=(x2-x0)/2.0;
1012 fhinv=1.0/fdx;
1013 fx1=(x0+x2)/2.0;
1014
1015 y0=f1.value_with_derivatives(x0, &yp0, &ypp0); // get left wall values from conventional computation
1016 y2=f2.value_with_derivatives(x2, &yp2, &ypp2); // get right wall values from conventional computation
1017
1018 // scale derivs to put function on [-1,1] since mma solution is done this way
1019 yp0*=fdx;
1020 yp2*=fdx;
1021 ypp0*=fdx*fdx;
1022 ypp2*=fdx*fdx;
1023
1024 float_type ff0=(8*(y0 + y2) + 5*(yp0 - yp2) + ypp0 + ypp2)/16.0;
1025 if(auto_center) y1=ff0; // forces ff to be 0 if we are auto-centering
1026
1027 // y[x_] = y1 + x (a + b x) + (x-1) x (x+1) (c + d x + e x^2 + f x^3)
1028 fy1=y1;
1029 fa=-(y0 - y2)/2.;
1030 fb=(y0 - 2*y1 + y2)/2.;
1031 fc=(7*(y0 - y2 + yp0 + yp2) + ypp0 - ypp2)/16.;
1032 fd=(32*y1 - 16*(y2 + y0) + 9*(yp2 - yp0) - ypp0 - ypp2)/16.;
1033 fe=(3*(y2 - y0 - yp0 - yp2) - ypp0 + ypp2)/16.;
1034 ff=(ff0 - y1);
1035 // y'[x] = a + 2 b x + (3x^2 - 1) (c + d x + e x^2 + f x^3) + (x-1) x (x+1) (d + 2 e x + 3 f x^2 )
1036 // y''[x] = 2b + (x-1) x (x+1) (2 e + 6 f x) + 2 (3 x^2 -1) (d + 2 e x + 3 f x^2 ) + 6 x (c + d x + e x^2 + f x^3)
1037 this->set_domain(x0,x2); // this is where the function is valid
1038}
1039
1040template <typename float_type> c2_connector_function<float_type>::~c2_connector_function()
1041{
1042}
1043
1044template <typename float_type> float_type c2_connector_function<float_type>::value_with_derivatives(
1045 float_type x, float_type *yprime, float_type *yprime2
1046 ) const throw(c2_exception)
1047{
1048
1049 float_type dx=(x-fx1)*fhinv;
1050 float_type q1=fc + dx*(fd + dx*(fe + dx*ff));
1051 float_type xp1=(dx-1)*(dx+1)*dx;
1052
1053 float_type y= fy1 + dx*(fa+fb*dx) + xp1*q1;
1054 if(yprime || yprime2) {
1055 float_type q2 =fd + dx*(2*fe + dx*3*ff);
1056 float_type q3=2*fe+6*ff*dx;
1057 float_type xp2=(3*dx*dx-1);
1058 if(yprime) *yprime=(fa + 2*fb*dx + xp2*q1 + xp1*q2)*fhinv;
1059 if(yprime2) *yprime2=(2*fb+xp1*q3+2*xp2*q2+6*dx*q1)*fhinv*fhinv;
1060 }
1061 return y;
1062}
Note: See TracBrowser for help on using the repository browser.