Hello,
Just in case anyone is using the cumuintover function I sent out a few
months ago, I found a small bug. The fourth element of the results was
always zero, due to a failure to use decimals to indicate floating
point constants in one line of the PP definition. The bug did not
affect any other element. Here's the fixed PP def.
no PDL::NiceSlice;
__DATA__
__Pdlpp__
# from equation 4.1.14, Press et al. 2nd Ed chapter 4.1
pp_def('cumuintover',
Pars => 'f(n); int+ [o]g(n);',
Code => '
int i;
int ns = $SIZE(n);
$GENERIC(g) tmp = 0;
threadloop %{
tmp = (3.0/8.0) * $f(n=>0) +
(7.0/6.0) * $f(n=>1) +
(23.0/24.0) * $f(n=>2);
for (i=5;i<ns;i++) { $g(n=>i) = tmp +
(23./24.) * $f(n=>i-2) +
(7./6.) * $f(n=>i-1) +
(3./8.) * $f(n=>i);
tmp += $f(n=>i-2);
}
$g(n=>0) = 0;
$g(n=>1) = ($f(n=>0) + $f(n=>1))/2.0;
$g(n=>2) = ($f(n=>0) + 4. * $f(n=>1) + $f(n=>2) ) / 3.0;
$g(n=>3) = (3./8.) * ($f(n=>0) + $f(n=>3) + 3. * ( $f(n=>1) + $f(n=>2)
) );
$g(n=>4) = (1./45.) * ( 14. * ( $f(n=>0) + $f(n=>4) )
+ 64. * ( $f(n=>1) + $f(n=>3) )
+ 24. * $f(n=>2)
);
%}
'
);
On 15:53, Thu 15 Mar 07, Roban Hultman Kramer wrote:
> Hello all,
>
> Just in case it's of interest to anyone else, I've implemented a
> "cumuintover" function. Analogously to "cumusumover" and
> "cumuprodover" it returns a piddle with the same dimensions as the
> input containing the integral over the first dimension.
>
> Note that it will crash if you feed it a piddle with a first dimension
> of less than 6. If it's of use to anyone, I could modify it to behave
> like intover and work for dims < 6.
>
> This is my first use of PP, so I welcome any feedback folks might
> have.
>
> Attached is a demo of the code using inline, along with the pp_def.
>
> Here's what I get:
>
> perl demo_cumuintover.pl
> x | num | ana | rel err
> 0 0 0 nan
> 1 2.13 2.13 8.43e-08
> 2 14.1 14.1 3.68e-08
> 3 56.5 56.5 1.91e-08
> 4 162 162 1.23e-08
> 5 377 377 9.35e-09
> 6 758 758 7.95e-09
> 7 1.38e+03 1.38e+03 7.37e-09
> 8 2.33e+03 2.33e+03 7.3e-09
> 9 3.7e+03 3.7e+03 7.62e-09
> 10 5.63e+03 5.63e+03 8.3e-09
>
> -roban
>
>
>
>
> #!/usr/bin/perl -w
> use strict;
> use PDL;
> use PDL::NiceSlice;
> use Inline qw(Pdlpp);
>
> my $dx = 0.1;
> my $x = $dx * pdl [0..100];
> my $y = $x**2 + 2*$x**3 + exp($x/2);
> my $g_num1 = $dx * PDL::cumuintover($y);
> my $g_ana = $x**3/3 + 2*$x**4/4 + 2 * exp($x/2) - 2;
> print " x | num | ana | rel err\n";
> wcols "%8.3g ", $x(0:-1:10), $g_num1(0:-1:10), $g_ana(0:-1:10),
> ($g_num1(0:-1:10) - $g_ana(0:-1:10))/$g_ana(0:-1:10);
>
>
> no PDL::NiceSlice;
>
> __DATA__
>
> __Pdlpp__
>
>
> # from equation 4.1.14, Press et al. 2nd Ed chapter 4.1
> pp_def('cumuintover',
> Pars => 'f(n); int+ [o]g(n);',
> Code => '
> int i;
> int ns = $SIZE(n);
> $GENERIC(g) tmp = 0;
> threadloop %{
> tmp = (3.0/8.0) * $f(n=>0) +
> (7.0/6.0) * $f(n=>1) +
> (23.0/24.0) * $f(n=>2);
> for (i=5;i<ns;i++) { $g(n=>i) = tmp +
> (23./24.) * $f(n=>i-2) +
> (7./6.) * $f(n=>i-1) +
> (3./8.) * $f(n=>i);
> tmp += $f(n=>i-2);
> }
> $g(n=>0) = 0;
> $g(n=>1) = ($f(n=>0) + $f(n=>1))/2.0;
> $g(n=>2) = ($f(n=>0) + 4. * $f(n=>1) + $f(n=>2) ) / 3.0;
> $g(n=>3) = (3/8) * ($f(n=>0) + $f(n=>3) + 3. * ( $f(n=>1) + $f(n=>2) ) );
> $g(n=>4) = (1./45.) * ( 14. * ( $f(n=>0) + $f(n=>4) )
> + 64. * ( $f(n=>1) + $f(n=>3) )
> + 24. * $f(n=>2)
> );
> %}
> '
> );
>
> _______________________________________________
> Perldl mailing list
> [email protected]
> http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
_______________________________________________
Perldl mailing list
[email protected]
http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
