Parent Directory
| 
Revision Log
Copied PG modules over from old system. -sam
1 #!/usr/bin/perl -w 2 3 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 15 #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 20 $obj -> initialize; # calculates the approximation 21 22 #get methods 23 $rf_function = $obj -> rf_f; 24 $rf_function_derivative = $obj -> rf_fp; 25 $rf_function_2nd_derivative = $obj -> rf_fpp; 26 27 $rh_critical_points =$obj -> rh_critical_points 28 $rh_inflection_points =$obj -> rh_inflection_points 29 30 31 32 33 =head1 DESCRIPTION 34 35 This module defines an object containing a Hermite spline approximation to a function. 36 The approximation 37 consists of a piecewise cubic polynomial which agrees with the original 38 function and its first derivative at 39 the node points. 40 41 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 the number of extraneous critical points 45 introduced. Too many points will introduce many small variations and 46 a large number of extraneous critical points. 47 48 There are even more extraneous inflections points. 49 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 89 $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 97 $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 112 # define the minimum points 113 114 } 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 123 } 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 $self ->{rf_f}; 144 } 145 sub rf_fp { 146 my $self = shift; 147 my $rf_fp =shift; 148 $self ->{rf_fp} = $rf_fp if defined( $rf_fp ); 149 $self ->{rf_fp}; 150 } 151 sub rf_fpp { 152 my $self = shift; 153 my $rf_fpp =shift; 154 $self ->{rf_fpp} = $rf_fpp if defined( $rf_fpp ); 155 $self ->{rf_fpp}; 156 } 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 204 # coefficients for the approximating polynomial 205 206 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 220 $rh_roots -> {$z1} = &$rf_function($z1); 221 } 222 } 223 } else { 224 my $discriminent = (4*$a[2]**2 - 12*$a[1]*$a[3]); 225 226 227 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 } 237 } 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 243 # coefficients for the approximating polynomial 244 245 my @a = ( $l, 246 $lp, 247 -$lp/2 + (3*(-$lp - 2*($l - $r) - $rp))/2 + $rp/2, 248 $lp + 2*($l - $r) + $rp 249 ); 250 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 257 } 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 $a[0] = $y0; 269 $a[1] = $yp0; 270 $a[2] = -3*$y0 - 2*$yp0 +3*$y1 -$yp1; 271 $a[3] = 2*$y0 + $yp0 - 2*$y1 +$yp1; 272 273 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 279 }; 280 $f; 281 } 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 $a[0] = $y0; 291 $a[1] = $yp0; 292 $a[2] = -3*$y0 - 2*$yp0 +3*$y1 -$yp1; 293 $a[3] = 2*$y0 + $yp0 - 2*$y1 +$yp1; 294 295 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 303 304 $fp; 305 306 } 307 308 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 $a[0] = $y0; 316 $a[1] = $yp0; 317 $a[2] = -3*$y0 - 2*$yp0 +3*$y1 -$yp1; 318 $a[3] = 2*$y0 + $yp0 - 2*$y1 +$yp1; 319 320 321 322 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 329 $fpp; 330 } 331 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 345 while (@xvals) { 346 $x1 = shift @xvals; 347 $y1 = shift @yvals; 348 $yp1 = shift @ypvals; 349 350 push @polys, cubic_hermite($x0, $y0, $yp0 , $x1, $y1 , $yp1); 351 $x0 = $x1; 352 $y0 = $y1; 353 $yp0 = $yp1; 354 } 355 356 357 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 365 while (@xvals && $x > $xvals[0]) { # find the function for this range of x 366 shift(@xvals); 367 $fun = shift(@fns); 368 } 369 370 # 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 400 401 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 409 while (@xvals && $x > $xvals[0]) { # find the function for this range of x 410 shift(@xvals); 411 $fun = shift(@fns); 412 } 413 414 # 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 443 444 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 452 while (@xvals && $x > $xvals[0]) { # find the function for this range of x 453 shift(@xvals); 454 $fun = shift(@fns); 455 } 456 457 # 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 |