Skip to content

Commit

Permalink
make PCHIP routines skip parameter be [o] at end of Pars
Browse files Browse the repository at this point in the history
  • Loading branch information
mohawk2 committed Dec 26, 2024
1 parent 637be67 commit 14ed9a7
Show file tree
Hide file tree
Showing 3 changed files with 20 additions and 26 deletions.
14 changes: 4 additions & 10 deletions lib/PDL/Func.pm
Original file line number Diff line number Diff line change
Expand Up @@ -761,13 +761,7 @@ sub _interp_hermite {
my ( $self, $xi, $x, $y ) = @_;
# get gradient
my $g = $self->_get_value( 'g' );
my @highest_dims;
for my $param ($x,$y,$g,$xi) {
my @this_dims = $param->dims;
@highest_dims = @this_dims if @this_dims > @highest_dims;
}
shift @highest_dims; # don't care about bottom one
my ( $yi, $ierr ) = PDL::Primitive::pchip_chfe( $x, $y, $g, PDL->ones(PDL::long(),@highest_dims), $xi );
my ( $yi, $ierr ) = PDL::Primitive::pchip_chfe( $x, $y, $g, $xi );
$self->{flags}{routine} = "pchip_chfe";
$self->_set_value( err => $ierr );
if ( all $ierr == 0 ) {
Expand Down Expand Up @@ -814,7 +808,7 @@ sub gradient {
# get values in one go
my ( $x, $y, $g ) = $self->_get_value( qw( x y g ) );

my ( $yi, $gi, $ierr ) = PDL::Primitive::pchip_chfd( $x, $y, $g, 0, $xi );
my ( $yi, $gi, $ierr ) = PDL::Primitive::pchip_chfd( $x, $y, $g, $xi );
$self->{flags}{routine} = "pchip_chfd";
$self->_set_value( err => $ierr );

Expand Down Expand Up @@ -903,7 +897,7 @@ sub integrate {
my ( $ans, $ierr );

if ( $type eq "x" ) {
( $ans, $ierr ) = PDL::Primitive::pchip_chia( $x, $y, $g, 0, $lo, $hi );
( $ans, $ierr ) = PDL::Primitive::pchip_chia( $x, $y, $g, $lo, $hi );
$self->{flags}{routine} = "pchip_chia";

if ( all $ierr == 0 ) {
Expand All @@ -918,7 +912,7 @@ sub integrate {
}

} else {
( $ans, $ierr ) = PDL::Primitive::pchip_chid( $x, $y, $g, 0, $lo, $hi );
( $ans, $ierr ) = PDL::Primitive::pchip_chid( $x, $y, $g, $lo, $hi );
$self->{flags}->{routine} = "pchip_chid";

if ( all $ierr == 0 ) {
Expand Down
16 changes: 8 additions & 8 deletions lib/PDL/Primitive.pd
Original file line number Diff line number Diff line change
Expand Up @@ -5463,8 +5463,8 @@ EOF
);

pp_def('pchip_chfd',
Pars => 'x(n); f(n); d(n); int [io] skip(); xe(ne);
[o] fe(ne); [o] de(ne); indx [o] ierr()',
Pars => 'x(n); f(n); d(n); xe(ne);
[o] fe(ne); [o] de(ne); indx [o] ierr(); int [o] skip()',
RedoDimsCode => <<'EOF',
if ($SIZE(n) < 2) $CROAK("NUMBER OF DATA POINTS LESS THAN TWO");
if ($SIZE(ne) < 1) $CROAK("NUMBER OF EVALUATION POINTS LESS THAN ONE");
Expand Down Expand Up @@ -5542,8 +5542,8 @@ EOF
);

pp_def('pchip_chfe',
Pars => 'x(n); f(n); d(n); int [io] skip(); xe(ne);
[o] fe(ne); indx [o] ierr()',
Pars => 'x(n); f(n); d(n); xe(ne);
[o] fe(ne); indx [o] ierr(); int [o] skip()',
RedoDimsCode => <<'EOF',
if ($SIZE(n) < 2) $CROAK("NUMBER OF DATA POINTS LESS THAN TWO");
if ($SIZE(ne) < 1) $CROAK("NUMBER OF EVALUATION POINTS LESS THAN ONE");
Expand Down Expand Up @@ -5631,8 +5631,8 @@ EOF
);

pp_def('pchip_chia',
Pars => 'x(n); f(n); d(n); int [io]skip(); la(); lb();
[o]ans(); indx [o]ierr()',
Pars => 'x(n); f(n); d(n); la(); lb();
[o]ans(); indx [o]ierr(); int [o]skip()',
GenericTypes => $F,
RedoDimsCode => <<'EOF',
if ($SIZE(n) < 2) $CROAK("NUMBER OF DATA POINTS LESS THAN TWO");
Expand Down Expand Up @@ -5810,8 +5810,8 @@ EOF

pp_def('pchip_chid',
Pars => 'x(n); f(n); d(n);
int [io]skip(); indx ia(); indx ib();
[o]ans(); indx [o]ierr()',
indx ia(); indx ib();
[o]ans(); indx [o]ierr(); int [o]skip()',
RedoDimsCode => <<'EOF',
if ($SIZE(n) < 2) $CROAK("NUMBER OF DATA POINTS LESS THAN TWO");
EOF
Expand Down
16 changes: 8 additions & 8 deletions t/primitive-interpolate.t
Original file line number Diff line number Diff line change
Expand Up @@ -60,13 +60,13 @@ subtest PCHIP => sub {
is_pdl $g, pdl('0 1.5 3.75 5.8333333 7.875 9.9 11.916667 13.928571 15.9375 18; 0 1.75 10.230769 25.107143 46.061224 73.039474 106.02752 145.02027 190.01554 241'), 'pchip_chic';
($g) = pchip_chim($x, $y);
is_pdl $g, pdl('0 1.5 3.75 5.8333333 7.875 9.9 11.916667 13.928571 15.9375 18; 0 1.75 10.230769 25.107143 46.061224 73.039474 106.02752 145.02027 190.01554 241'), 'pchip_chim';
my ($yi) = pchip_chfe($x, $y, $g, [1,1], $xi);
my ($yi) = pchip_chfe($x, $y, $g, $xi);
is_pdl $yi, my $yi_exp = pdl('48.56375 54.173375 61.777925 71.38055 82.98225; -10.973827 12.780893 56.345513 125.71307 226.88177'), 'pchip_chfe';
($yi) = pchip_chfd($x, $y, $g, [0,0], $xi);
($yi) = pchip_chfd($x, $y, $g, $xi);
is_pdl $yi, $yi_exp, 'pchip_chfd';
my ($integral) = pchip_chia($x, $y, $g, [0,0], 3, 4);
my ($integral) = pchip_chia($x, $y, $g, 3, 4);
is_pdl $integral, pdl(55.6298611111111,20.7538265306122), 'pchip_chia';
($integral) = pchip_chid($x, $y, $g, [0,0], 3, 4);
($integral) = pchip_chid($x, $y, $g, 3, 4);
is_pdl $integral, pdl(55.6298611111111,20.7538265306122), 'pchip_chid';
my ($t, $bcoef) = pchip_chbs($x, $y, $g, 0);
is_pdl $t, pdl('0 0 0 0 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 9 9; 0 0 0 0 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 9 9'), 'pchip_chbs t';
Expand All @@ -93,13 +93,13 @@ subtest PCHIP => sub {
is_pdl $d3, $d->double, {atol=>2};
my $xe = float( pdl( 4 .. 8 ) + 0.5 );
my ( $fe, $de );
( $fe, $de, $err ) = pchip_chfd( $x, $f, $d, 1, $xe );
( $fe, $de, $err ) = pchip_chfd( $x, $f, $d, $xe );
is_pdl $err, indx(0);
$answer = $xe*$xe*$xe + 425.42352;
is_pdl $fe, $answer, {rtol=>1e-5};
$answer = 3.0*$xe*$xe;
is_pdl $de, $answer, {rtol=>2e-2};
( $fe, $err ) = pchip_chfe( $x, $f, $d, 1, $xe );
( $fe, $err ) = pchip_chfe( $x, $f, $d, $xe );
is_pdl $fe, $xe*$xe*$xe + 425.42352, {rtol=>1e-3};
is_pdl $err, indx 0;
$x = float( 1, 2, 3, 5, 6, 7 );
Expand All @@ -110,13 +110,13 @@ subtest PCHIP => sub {
$f = $x * $x;
( $d, $err ) = pchip_chim($x, $f);
my $ans = pdl( 9.0**3, (8.0**3-1.0**3) ) / 3.0;
( my $int, $err ) = pchip_chia($x, $f, $d, my $skip=zeroes(2), pdl(0.0,1.0), pdl(9.0,8.0));
( my $int, $err ) = pchip_chia($x, $f, $d, pdl(0.0,1.0), pdl(9.0,8.0));
is_pdl $err, indx 0,0;
is_pdl $int, $ans, {atol=>4e-2};
my $hi = pdl( $x->at(9), $x->at(7) );
my $lo = pdl( $x->at(0), $x->at(1) );
$ans = ($hi**3 - $lo**3) / 3;
( $int, $err ) = pchip_chid( $x, $f, $d, $skip=zeroes(2), pdl(0,1), pdl(9,7) );
( $int, $err ) = pchip_chid( $x, $f, $d, pdl(0,1), pdl(9,7) );
is_pdl $err, indx 0,0;
is_pdl $int, $ans, {atol=>6e-2};
};
Expand Down

0 comments on commit 14ed9a7

Please sign in to comment.