[system] / trunk / pg / lib / Hermite.pm Repository:
ViewVC logotype

Annotation of /trunk/pg/lib/Hermite.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1079 - (view) (download) (as text)

1 : sh002i 1050
2 : apizer 1079
3 : sh002i 1050 use strict;
4 :     use Carp;
5 :    
6 :     =head1 NAME
7 :    
8 :     Hermite.pm
9 :    
10 :     =head1 SYNPOSIS
11 :    
12 :     Usage:
13 :     $obj = new Hermit(\@x_values, \y_valuses \@yp_values);
14 : apizer 1079
15 : sh002i 1050 #get and set methods
16 :     $ra_x_values = $obj -> ra_x(\@x_values);
17 :     $ra_y_values = $obj -> ra_y;
18 :     $ra_yp_values = $obj -> ra_yp;
19 : apizer 1079
20 : sh002i 1050 $obj -> initialize; # calculates the approximation
21 : apizer 1079
22 : sh002i 1050 #get methods
23 :     $rf_function = $obj -> rf_f;
24 :     $rf_function_derivative = $obj -> rf_fp;
25 :     $rf_function_2nd_derivative = $obj -> rf_fpp;
26 : apizer 1079
27 : sh002i 1050 $rh_critical_points =$obj -> rh_critical_points
28 :     $rh_inflection_points =$obj -> rh_inflection_points
29 :    
30 :    
31 :    
32 : apizer 1079
33 : sh002i 1050 =head1 DESCRIPTION
34 :    
35 :     This module defines an object containing a Hermite spline approximation to a function.
36 :     The approximation
37 : apizer 1079 consists of a piecewise cubic polynomial which agrees with the original
38 : sh002i 1050 function and its first derivative at
39 :     the node points.
40 :    
41 : apizer 1079 This is useful for creating on the fly graphics. Care must be taken to use a small
42 :     number of points spaced reasonably far apart, preferably
43 :     points with alternating slopes, since this will minimize
44 : sh002i 1050 the number of extraneous critical points
45 : apizer 1079 introduced. Too many points will introduce many small variations and
46 : sh002i 1050 a large number of extraneous critical points.
47 :    
48 : apizer 1079 There are even more extraneous inflections points.
49 : sh002i 1050 This parameter is probably not very useful. A different
50 :     approximation scheme needs to be used.
51 :    
52 :     =cut
53 :    
54 :    
55 :     # an object that contains a cubic hermite spline
56 :    
57 :     package Hermite;
58 :    
59 :    
60 :     sub new {
61 :     my $class = shift;
62 :     my ($ra_x,$ra_y,$ra_yp) = @_;
63 :     my $self = {
64 :     ra_x => [],
65 :     ra_y => [],
66 :     ra_yp => [],
67 :     rf_f => 0, # refers to a subroutine for the function
68 :     rf_fp => 0, # refers to a subroutine for the derivative of the function
69 :     rf_fpp => 0, # refers to a subroutine for the second derivative of the function
70 :     rh_critical_points => {},
71 :     rh_inflection_points => {},
72 :     rh_maximum_points => {},
73 :     rh_minimum_points => {},
74 :     };
75 :     bless $self, $class;
76 :     $self->define($ra_x,$ra_y,$ra_yp) if ref($ra_x) eq 'ARRAY' and ref($ra_y) eq 'ARRAY' and ref($ra_yp) eq 'ARRAY';
77 :     $self;
78 :     }
79 :    
80 :     sub define{
81 :     my $self = shift;
82 :     my $ra_x = shift;
83 :     my $ra_y = shift;
84 :     my $ra_yp = shift;
85 :     $self->ra_x($ra_x);
86 :     $self->ra_y($ra_y);
87 :     $self->ra_yp($ra_yp);
88 : apizer 1079
89 : sh002i 1050 $self -> initialize();
90 :     ($self->{x_ref},$self->{y_ref}, $self->{yp_ref} )
91 :     }
92 :    
93 :     sub initialize {
94 :     my $self = shift;
95 :     # define the functions rf_f
96 : apizer 1079
97 : sh002i 1050 $self -> {rf_f} = hermite_spline($self -> {ra_x}, $self -> {ra_y}, $self -> {ra_yp} );
98 :     # # define the function rf_fp
99 :     $self -> {rf_fp} = hermite_spline_p($self -> {ra_x}, $self -> {ra_y}, $self -> {ra_yp} );
100 :     # # define the function rf_fp
101 :     $self -> {rf_fpp} = hermite_spline_pp($self -> {ra_x}, $self -> {ra_y}, $self -> {ra_yp} );
102 :     # # define the critical points
103 :     %{$self -> {rh_critical_points}} = critical_points($self -> {ra_x}, $self -> {ra_y}, $self -> {ra_yp},
104 :     $self -> { rf_f },
105 :     );
106 :     # # define the inflection_points
107 :     %{$self -> {rh_inflection_points}} = inflection_points($self -> {ra_x}, $self -> {ra_y}, $self -> {ra_yp} ,
108 :     $self -> { rf_f },
109 :     );
110 :     # # define the maximum points
111 : apizer 1079
112 : sh002i 1050 # define the minimum points
113 : apizer 1079
114 : sh002i 1050 }
115 :    
116 :     # fix up these accesses to do more checking
117 :     sub ra_x {
118 :     my $self = shift;
119 :     my $ra_x = shift;
120 :     @{$self->{ra_x}} = @$ra_x if ref($ra_x) eq 'ARRAY';
121 :     $self->{ra_x};
122 : apizer 1079
123 : sh002i 1050 }
124 :    
125 :     sub ra_y {
126 :     my $self = shift;
127 :     my $ra_y = shift;
128 :     @{$self->{ra_y}} = @$ra_y if ref($ra_y) eq 'ARRAY';
129 :     $self->{ra_y};
130 :     }
131 :    
132 :     sub ra_yp {
133 :     my $self = shift;
134 :     my $ra_yp = shift;
135 :     @{$self->{ra_yp}} = @$ra_yp if ref($ra_yp) eq 'ARRAY';
136 :     $self->{ra_yp};
137 :     }
138 :    
139 :     sub rf_f {
140 :     my $self = shift;
141 :     my $rf_f =shift;
142 :     $self ->{rf_f} = $rf_f if defined( $rf_f );
143 : apizer 1079 $self ->{rf_f};
144 : sh002i 1050 }
145 :     sub rf_fp {
146 :     my $self = shift;
147 :     my $rf_fp =shift;
148 :     $self ->{rf_fp} = $rf_fp if defined( $rf_fp );
149 : apizer 1079 $self ->{rf_fp};
150 : sh002i 1050 }
151 :     sub rf_fpp {
152 :     my $self = shift;
153 :     my $rf_fpp =shift;
154 :     $self ->{rf_fpp} = $rf_fpp if defined( $rf_fpp );
155 : apizer 1079 $self ->{rf_fpp};
156 : sh002i 1050 }
157 :    
158 :     sub rh_critical_points {
159 :     my $self = shift;
160 :     $self -> { rh_critical_points};
161 :     }
162 :    
163 :     sub rh_inflection_points {
164 :     my $self = shift;
165 :     $self -> { rh_inflection_points};
166 :     }
167 :    
168 :     ##Internal subroutines
169 :    
170 :    
171 :     sub critical_points{
172 :     my $ra_x = shift;
173 :     my $ra_y =shift;
174 :     my $ra_yp = shift;
175 :     my $rf_hermite_fun =shift;
176 :     my %critical_points = ();
177 :     my $last_index = @$ra_x -2;
178 :     foreach my $i (0 .. $last_index ) {
179 :     internal_critical_points($ra_x->[$i], $ra_y->[$i], $ra_yp->[$i],
180 :     $ra_x->[$i+1], $ra_y->[$i+1], $ra_yp->[$i+1],
181 :     \%critical_points,$rf_hermite_fun);
182 :     }
183 :     %critical_points;
184 :     }
185 :    
186 :     sub inflection_points{
187 :     my $ra_x = shift;
188 :     my $ra_y =shift;
189 :     my $ra_yp = shift;
190 :     my $rf_hermite_fun =shift;
191 :     my %inflection_points = ();
192 :     my $last_index = @$ra_x -2;
193 :     foreach my $i (0 .. $last_index ) {
194 :     internal_inflection_points($ra_x->[$i], $ra_y->[$i], $ra_yp->[$i],
195 :     $ra_x->[$i+1], $ra_y->[$i+1], $ra_yp->[$i+1],
196 :     \%inflection_points,$rf_hermite_fun);
197 :     }
198 :     %inflection_points;
199 :     }
200 :     sub internal_critical_points{
201 :     my ($x0,$l,$lp, $x1,$r,$rp, $rh_roots ,$rf_function) = @_;
202 :     #data for one segment of the hermite spline
203 : apizer 1079
204 : sh002i 1050 # coefficients for the approximating polynomial
205 : apizer 1079
206 : sh002i 1050 my @a = ( $l,
207 :     $lp,
208 :     -$lp/2 + (3*(-$lp - 2*($l - $r) - $rp))/2 + $rp/2,
209 :     $lp + 2*($l - $r) + $rp
210 :     );
211 :    
212 :     my ($root1, $root2, $z1,$z2);
213 :     if ( $a[3] == 0 ) {
214 :     if ( $a[2] == 0 ) {
215 :     } else {
216 :     $root1 = -$a[1]/( 2*$a[2] );
217 :     if ( 0 <= $root1 and $root1 < 1) {
218 :     $z1 = $root1*($x1-$x0) + $x0;
219 : apizer 1079
220 :     $rh_roots -> {$z1} = &$rf_function($z1);
221 : sh002i 1050 }
222 :     }
223 :     } else {
224 :     my $discriminent = (4*$a[2]**2 - 12*$a[1]*$a[3]);
225 :    
226 : apizer 1079
227 : sh002i 1050 if ( $discriminent >= 0 ) {
228 :     $discriminent = $discriminent**0.5;
229 :     $root1 = (-2*$a[2] - $discriminent )/( 6*$a[3] );
230 :     $root2 = (-2*$a[2] + $discriminent )/( 6*$a[3] );
231 :     $z1 = $root1*($x1-$x0) + $x0;
232 :     $z2 = $root2*($x1-$x0) + $x0;
233 :     $rh_roots -> {$z1} = &$rf_function($z1) if 0 <= $root1 and $root1 < 1;
234 :     $rh_roots -> {$z2} = &$rf_function($z1) if 0 <= $root2 and $root2 < 1;
235 :     }
236 : apizer 1079 }
237 : sh002i 1050 }
238 :    
239 :     sub internal_inflection_points {
240 :     my ($x0,$l,$lp,,$x1,$r,$rp,$rh_roots,$rf_function) = @_;
241 :     #data for one segment of the hermite spline
242 : apizer 1079
243 : sh002i 1050 # coefficients for the approximating polynomial
244 : apizer 1079
245 : sh002i 1050 my @a = ( $l,
246 :     $lp,
247 :     -$lp/2 + (3*(-$lp - 2*($l - $r) - $rp))/2 + $rp/2,
248 :     $lp + 2*($l - $r) + $rp
249 : apizer 1079 );
250 : sh002i 1050 if ($a[3] == 0 ) {
251 :     } else {
252 :     my $root1 = -$a[2]/( 3*$a[3] );
253 :     my $z1 = $root1*($x1-$x0) + $x0;
254 :     $rh_roots -> {$z1} = &$rf_function($z1) if 0 <= $root1 and $root1 < 1;
255 :     }
256 : apizer 1079
257 : sh002i 1050 }
258 :    
259 :    
260 :    
261 :     sub cubic_hermite {
262 :     my ( $x0, $y0,$yp0, $x1, $y1, $yp1 ) = @_;
263 :     my @a;
264 :     my $width = $x1 - $x0;
265 :     $yp0 = $yp0*$width; # normalize to unit width
266 :     $yp1 = $yp1*$width;
267 :    
268 : apizer 1079 $a[0] = $y0;
269 : sh002i 1050 $a[1] = $yp0;
270 :     $a[2] = -3*$y0 - 2*$yp0 +3*$y1 -$yp1;
271 :     $a[3] = 2*$y0 + $yp0 - 2*$y1 +$yp1;
272 : apizer 1079
273 : sh002i 1050 my $f = sub {
274 :     my $x = shift;
275 :     #normalize to unit width
276 :     $x = ( $x - $x0 )/$width;
277 :     ( ($a[3]*$x + $a[2]) * $x + $a[1] )*$x + $a[0];
278 : apizer 1079
279 : sh002i 1050 };
280 : apizer 1079 $f;
281 : sh002i 1050 }
282 :    
283 :     sub cubic_hermite_p {
284 :     my ( $x0, $y0,$yp0, $x1, $y1, $yp1 )=@_;
285 :     my @a;
286 :     my $width = $x1 - $x0;
287 :     $yp0 = $yp0*$width; # normalize to unit width
288 :     $yp1 = $yp1*$width;
289 :    
290 : apizer 1079 $a[0] = $y0;
291 : sh002i 1050 $a[1] = $yp0;
292 :     $a[2] = -3*$y0 - 2*$yp0 +3*$y1 -$yp1;
293 :     $a[3] = 2*$y0 + $yp0 - 2*$y1 +$yp1;
294 : apizer 1079
295 : sh002i 1050 my $fp = sub {
296 :     my $x = shift;
297 :     #normalize to unit width
298 :     $x = ( $x - $x0 )/$width;
299 :     ( (3*$a[3]*$x + 2*$a[2]) * $x + $a[1] )/$width ;
300 :     };
301 :    
302 : apizer 1079
303 :    
304 : sh002i 1050 $fp;
305 :    
306 : apizer 1079 }
307 :    
308 : sh002i 1050 sub cubic_hermite_pp {
309 :     my ( $x0, $y0,$yp0, $x1, $y1, $yp1 ) = @_;
310 :     my @a;
311 :     my $width = $x1 - $x0;
312 :     $yp0 = $yp0*$width; # normalize to unit width
313 :     $yp1 = $yp1*$width;
314 :    
315 : apizer 1079 $a[0] = $y0;
316 : sh002i 1050 $a[1] = $yp0;
317 :     $a[2] = -3*$y0 - 2*$yp0 +3*$y1 -$yp1;
318 :     $a[3] = 2*$y0 + $yp0 - 2*$y1 +$yp1;
319 :    
320 : apizer 1079
321 :    
322 : sh002i 1050 my $fpp = sub {
323 :     my $x = shift;
324 :     #normalize to unit width
325 :     $x = ( $x - $x0 )/$width;
326 :     (6*$a[3]*$x + 2*$a[2])/$width**2 ;
327 :     };
328 : apizer 1079
329 : sh002i 1050 $fpp;
330 : apizer 1079 }
331 : sh002i 1050
332 :    
333 :    
334 :     sub hermite_spline {
335 :     my ($xref, $yref, $ypref) = @_;
336 :     my @xvals = @$xref;
337 :     my @yvals = @$yref;
338 :     my @ypvals = @$ypref;
339 :     my $x0 = shift @xvals;
340 :     my $y0 = shift @yvals;
341 :     my $yp0 = shift @ypvals;
342 :     my ($x1,$y1,$yp1);
343 :     my @polys; #calculate a hermite polynomial evaluator for each region
344 : apizer 1079
345 : sh002i 1050 while (@xvals) {
346 :     $x1 = shift @xvals;
347 :     $y1 = shift @yvals;
348 :     $yp1 = shift @ypvals;
349 : apizer 1079
350 : sh002i 1050 push @polys, cubic_hermite($x0, $y0, $yp0 , $x1, $y1 , $yp1);
351 :     $x0 = $x1;
352 :     $y0 = $y1;
353 :     $yp0 = $yp1;
354 :     }
355 : apizer 1079
356 :    
357 : sh002i 1050 my $hermite_spline_function = sub {
358 :     my $x = shift;
359 :     my $y;
360 :     my $fun;
361 :     my @xvals = @$xref;
362 :     my @fns = @polys;
363 :     return &{$fns[0]} ($x) if $x == $xvals[0]; #handle left most endpoint
364 : apizer 1079
365 : sh002i 1050 while (@xvals && $x > $xvals[0]) { # find the function for this range of x
366 :     shift(@xvals);
367 :     $fun = shift(@fns);
368 :     }
369 : apizer 1079
370 : sh002i 1050 # now that we have the left hand of the input
371 :     #check first that x isn't out of range to the left or right
372 :     if (@xvals && defined($fun) ) {
373 :     $y =&$fun($x);
374 :     }
375 :     $y;
376 :     };
377 :     $hermite_spline_function;
378 :     }
379 :    
380 :     sub hermite_spline_p {
381 :     my ($xref, $yref, $ypref) = @_;
382 :     my @xvals = @$xref;
383 :     my @yvals = @$yref;
384 :     my @ypvals = @$ypref;
385 :     my $x0 = shift @xvals;
386 :     my $y0 = shift @yvals;
387 :     my $yp0 = shift @ypvals;
388 :     my ($x1,$y1,$yp1);
389 :     my @polys; #calculate a hermite polynomial evaluator for each region
390 :     while (@xvals) {
391 :     $x1 = shift @xvals;
392 :     $y1 = shift @yvals;
393 :     $yp1 = shift @ypvals;
394 :     push @polys, cubic_hermite_p($x0, $y0, $yp0 , $x1, $y1 , $yp1);
395 :     $x0 = $x1;
396 :     $y0 = $y1;
397 :     $yp0 = $yp1;
398 :     }
399 : apizer 1079
400 :    
401 : sh002i 1050 my $hermite_spline_function_p = sub {
402 :     my $x = shift;
403 :     my $y;
404 :     my $fun;
405 :     my @xvals = @$xref;
406 :     my @fns = @polys;
407 :     return $y=&{$fns[0]} ($x) if $x == $xvals[0]; #handle left most endpoint
408 : apizer 1079
409 : sh002i 1050 while (@xvals && $x > $xvals[0]) { # find the function for this range of x
410 :     shift(@xvals);
411 :     $fun = shift(@fns);
412 :     }
413 : apizer 1079
414 : sh002i 1050 # now that we have the left hand of the input
415 :     #check first that x isn't out of range to the left or right
416 :     if (@xvals && defined($fun) ) {
417 :     $y =&$fun($x);
418 :     }
419 :     $y;
420 :     };
421 :     $hermite_spline_function_p;
422 :     }
423 :     sub hermite_spline_pp {
424 :     my ($xref, $yref, $ypref) = @_;
425 :     my @xvals = @$xref;
426 :     my @yvals = @$yref;
427 :     my @ypvals = @$ypref;
428 :     my $x0 = shift @xvals;
429 :     my $y0 = shift @yvals;
430 :     my $yp0 = shift @ypvals;
431 :     my ($x1,$y1,$yp1);
432 :     my @polys; #calculate a hermite polynomial evaluator for each region
433 :     while (@xvals) {
434 :     $x1 = shift @xvals;
435 :     $y1 = shift @yvals;
436 :     $yp1 = shift @ypvals;
437 :     push @polys, cubic_hermite_pp($x0, $y0, $yp0 , $x1, $y1 , $yp1);
438 :     $x0 = $x1;
439 :     $y0 = $y1;
440 :     $yp0 = $yp1;
441 :     }
442 : apizer 1079
443 :    
444 : sh002i 1050 my $hermite_spline_function_pp = sub {
445 :     my $x = shift;
446 :     my $y;
447 :     my $fun;
448 :     my @xvals = @$xref;
449 :     my @fns = @polys;
450 :     return $y=&{$fns[0]} ($x) if $x == $xvals[0]; #handle left most endpoint
451 : apizer 1079
452 : sh002i 1050 while (@xvals && $x > $xvals[0]) { # find the function for this range of x
453 :     shift(@xvals);
454 :     $fun = shift(@fns);
455 :     }
456 : apizer 1079
457 : sh002i 1050 # now that we have the left hand of the input
458 :     #check first that x isn't out of range to the left or right
459 :     if (@xvals && defined($fun) ) {
460 :     $y =&$fun($x);
461 :     }
462 :     $y;
463 :     };
464 :     $hermite_spline_function_pp;
465 :     }
466 :    
467 :    
468 :     1;

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9