port SLATEC PCHIP functions to C, add to Primitive

Basic/t/primitive-interpolate.t

@@ -3,6 +3,7 @@ use warnings;
 use Test::More;
 use Test::Exception;
 use PDL::LiteF;
+use Test::PDL;
 
 subtest interpol => sub {
 
@@ -30,7 +31,6 @@ subtest interpol => sub {
 };
 
 subtest interpND => sub {
-
     my $x     = xvals( 10, 10 ) + yvals( 10, 10 ) * 10;
     my $index = cat( 3 + xvals( 5, 5 ) * 0.25, 7 + yvals( 5, 5 ) * 0.25 )
       ->reorder( 2, 0, 1 );
@@ -40,4 +40,35 @@ subtest interpND => sub {
     ok !any( $y != $z ), "result";
 };
 
+subtest PCHIP => sub {
+  my $x = sequence(10);
+  my $y = pdl('43.3 44.3 47.3 52.3 59.3 68.3 79.3 92.3 107.3 124.3; -23 -22 -15 4 41 102 193 320 489 706');
+  my $xi = sequence(float,5) + 2.3;
+  my ($g) = pchip_chsp([0,0], [0,0], $x, $y);
+  is_pdl $g, pdl('0 2 4 6 8 10 12 14 16 18; 0 3 12 27 48 75 108 147 192 243'), 'pchip_chsp';
+  ($g) = pchip_chic([0,0], [0,0], 0, $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_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);
+  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);
+  is_pdl $yi, $yi_exp, 'pchip_chfd';
+  my ($integral) = pchip_chia($x, $y, $g, [0,0], 3, 4);
+  is_pdl $integral, pdl(55.6298611111111,20.7538265306122), 'pchip_chia';
+  ($integral) = pchip_chid($x, $y, $g, [0,0], 3, 4);
+  is_pdl $integral, pdl(55.6298611111111,20.7538265306122), 'pchip_chid';
+  my $nknots = zeroes(longlong, 2);
+  my $t = zeroes( 2*$x->dim(0)+4, 2 );
+  my ($bcoef, $ndim, $kord) = pchip_chbs($x, $y, $g, 0, $nknots, $t);
+  is_pdl $nknots, longlong(24,24), 'pchip_chbs nknots';
+  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';
+  is_pdl $bcoef, pdl('43.3 43.3 43.8 44.8 46.05 48.55 50.355556 54.244444 56.675 61.925 65 71.6 75.327778 83.272222 87.657143 96.942857 101.9875 112.6125 118.3 124.3; -23 -23 -22.583333 -21.416667 -18.410256 -11.589744 -4.3690476 12.369048 25.646259 56.353741 77.653509 126.34649 157.65749 228.34251 271.65991 368.34009 425.66149 552.33851 625.66667 706'), 'pchip_chbs bcoef';
+  is_pdl $ndim, longlong(20,20), 'pchip_chbs ndim';
+  is_pdl $kord, longlong(4,4), 'pchip_chbs kord';
+  my $x_slice = $x->slice('*1,0:-2'); # because calling with last value is out of range
+  my ($val) = pchip_bvalu($t, $bcoef, 0, $x_slice);
+  is_pdl $val->t, pdl('43.3 44.3 47.3 52.3 59.3 68.3 79.3 92.3 107.3; -23 -22 -15 4 41 102 193 320 489'), 'pchip_bvalu';
+};
+
 done_testing;

Changes

@@ -32,6 +32,7 @@
 - split PDL::Slatec out to separate distro
 - split PDL::Fit out to separate distro
 - split PDL::Opt::Simplex out to separate distro
+- add Primitive::pchip_{chsp,chic,chim,chfe,chfd,chia,chid,chbs,bvalu}
 
 2.095 2024-11-03
 - add PDL_GENTYPE_IS_{REAL,FLOATREAL,COMPLEX,SIGNED,UNSIGNED}_##ppsym (#502)

Libtmp/Func.pm

@@ -67,7 +67,6 @@ The valid attributes are:
 
 Use the piecewise cubic Hermite interpolation routines
 from the SLATEC library.
-Only available if L<PDL::Slatec> is installed.
 
 The valid attributes are:
 
@@ -94,7 +93,7 @@ C<interpolate>, C<gradient>, or C<integrate> is used.
 If your data is monotonic, and you are not too bothered about
 edge effects, then the default value of C<bc> of C<simple> is for you.
 Otherwise, take a look at the description of
-L<PDL::Slatec::chic|PDL::Slatec/chic> and use a hash reference
+L<PDL::Primitive::pchip_chic|PDL::Primitive/pchip_chic> and use a hash reference
 for the C<bc> attribute, with the following keys:
 
 =over 3
@@ -114,7 +113,7 @@ Default = B<0>.
 A perl list of one or two elements. The first element defines how the
 boundary condition for the start of the array is to be calculated;
 it has a range of C<-5 .. 5>, as given for the C<ic> parameter
-of L<chic|PDL::Slatec/chic>.
+of L<pchip_chic|PDL::Primitive/pchip_chic>.
 The second element, only used if options 2, 1, -1, or 2
 are chosen, contains the value of the C<vc> parameter.
 Default = B<[ 0 ]>.
@@ -136,7 +135,6 @@ and at the last point to -1.
 
 Use the cubic spline interpolation routines
 from the SLATEC library's PCHIP package.
-Only available if L<PDL::Slatec> is installed.
 
 The valid attributes are:
 
@@ -177,16 +175,6 @@ use parent qw(PDL::Exporter);
 our @EXPORT_OK = qw(pchip spline);
 our %EXPORT_TAGS = (Func=>[@EXPORT_OK]);
 
-####################################################################
-#
-# what modules are available ?
-#
-my %modules;
-BEGIN {
-    eval "use PDL::Slatec";
-    $modules{slatec} = ($@ ? 0 : 1);
-}
-
 ####################################################################
 ## Public routines:
 
@@ -284,11 +272,6 @@ sub _init_attr {
     croak "ERROR: Unknown interpolation scheme <$interpolate>.\n"
 	unless defined $attr{$interpolate};
 
-    # fall over if slatec library isn't present
-    # and asking for Hermite interpolation
-    croak "ERROR: Hermite interpolation is not available without PDL::Slatec.\n"
-        if $interpolate eq "Interpolate" and $modules{slatec} == 0;
-
     # clear out the old data (if it's not the first time through)
     $self->{attributes} = {};
     $self->{values}     = {};
@@ -407,26 +390,20 @@ sub _init_hermite {
 
     my ( $g, $ierr );
     if ( ref($bc) eq "HASH" ) {
-	croak "ERROR: Hermite interpolation is not available without PDL::Slatec.\n"
-	  if !$modules{slatec};
 	my $monotonic = $bc->{monotonic} || 0;
 	my $start     = $bc->{start}     || [ 0 ];
 	my $end       = $bc->{end}       || [ 0 ];
 
 	my $ic = [$start->[0], $end->[0]];
 	my $vc = [@$start == 2 ? $start->[1] : 0, @$end == 2 ? $end->[1] : 0];
 
-	( $g, $ierr ) = chic( $ic, $vc, $monotonic, $x, $y );
+	( $g, $ierr ) = PDL::Primitive::pchip_chic( $ic, $vc, $monotonic, $x, $y );
 
-	$self->{flags}{routine} = "chic";
+	$self->{flags}{routine} = "pchip_chic";
 
     } elsif ( $bc eq "simple" ) {
-	# chim
-	croak "ERROR: Hermite interpolation is not available without PDL::Slatec.\n"
-	  if $modules{slatec} == 0;
-	( $g, $ierr ) = chim( $x, $y );
-
-	$self->{flags}{routine} = "chim";
+	( $g, $ierr ) = PDL::Primitive::pchip_chim( $x, $y );
+	$self->{flags}{routine} = "pchip_chim";
 
     } else {
 	# Unknown boundary condition
@@ -449,8 +426,6 @@ sub _init_hermite {
 }
 
 sub _init_cspline {
-    croak "ERROR: CSpline interpolation is not available without PDL::Slatec\n"
-      if !$modules{slatec};
     my $self = shift;
     # set up error flags
     $self->{flags}{status} = 0;
@@ -468,8 +443,8 @@ sub _init_cspline {
 	my $end       = $bc->{end}       || [ 0 ];
 	my $ic = [$start->[0], $end->[0]];
 	my $vc = [@$start == 2 ? $start->[1] : 0, @$end == 2 ? $end->[1] : 0];
-	( $g, $ierr ) = chsp( $ic, $vc, $x, $y );
-	$self->{flags}{routine} = "chsp";
+	( $g, $ierr ) = PDL::Primitive::pchip_chsp( $ic, $vc, $x, $y );
+	$self->{flags}{routine} = "pchip_chsp";
     }
     $self->_set_value( g => $g, err => $ierr );
     if ( all $ierr == 0 ) {
@@ -792,8 +767,8 @@ sub _interp_hermite {
       @highest_dims = @this_dims if @this_dims > @highest_dims;
     }
     shift @highest_dims; # don't care about bottom one
-    my ( $yi, $ierr ) = chfe( $x, $y, $g, PDL->ones(PDL::long(),@highest_dims), $xi );
-    $self->{flags}{routine} = "chfe";
+    my ( $yi, $ierr ) = PDL::Primitive::pchip_chfe( $x, $y, $g, PDL->ones(PDL::long(),@highest_dims), $xi );
+    $self->{flags}{routine} = "pchip_chfe";
     $self->_set_value( err => $ierr );
     if ( all $ierr == 0 ) {
 	# everything okay
@@ -839,8 +814,8 @@ sub gradient {
     # get values in one go
     my ( $x, $y, $g ) = $self->_get_value( qw( x y g ) );
 
-    my ( $yi, $gi, $ierr ) = chfd( $x, $y, $g, 0, $xi );
-    $self->{flags}{routine} = "chfd";
+    my ( $yi, $gi, $ierr ) = PDL::Primitive::pchip_chfd( $x, $y, $g, 0, $xi );
+    $self->{flags}{routine} = "pchip_chfd";
     $self->_set_value( err => $ierr );
 
     if ( all $ierr == 0 ) {
@@ -928,8 +903,8 @@ sub integrate {
     my ( $ans, $ierr );
 
     if ( $type eq "x" ) {
-	( $ans, $ierr ) = chia( $x, $y, $g, 0, $lo, $hi );
-	$self->{flags}{routine} = "chia";
+	( $ans, $ierr ) = PDL::Primitive::pchip_chia( $x, $y, $g, 0, $lo, $hi );
+	$self->{flags}{routine} = "pchip_chia";
 
 	if ( all $ierr == 0 ) {
 	    # everything okay
@@ -943,8 +918,8 @@ sub integrate {
 	}
 
     } else {
-	( $ans, $ierr ) = chid( $x, $y, $g, 0, $lo, $hi );
-	$self->{flags}->{routine} = "chid";
+	( $ans, $ierr ) = PDL::Primitive::pchip_chid( $x, $y, $g, 0, $lo, $hi );
+	$self->{flags}->{routine} = "pchip_chid";
 
 	if ( all $ierr == 0 ) {
 	    # everything okay
@@ -1019,21 +994,15 @@ In the documentation, the methods are preceded by C<PDL::Func::>
 to avoid clashes with functions such as C<set> when using
 the C<help> or C<apropos> commands within I<perldl>.
 
-=head1 HISTORY
-
-Amalgamated C<PDL::Interpolate> and C<PDL::Interpolate::Slatec>
-to form C<PDL::Func>. Comments greatly appreciated on the
-current implementation, as it is not too sensible.
-
-Thanks to Robin Williams, Halldór Olafsson, and Vince McIntyre.
-
 =head1 AUTHOR
 
 Copyright (C) 2000,2001 Doug Burke ([email protected]).
 All rights reserved. There is no warranty.
 You are allowed to redistribute this software / documentation as
 described in the file COPYING in the PDL distribution.
 
+Thanks to Robin Williams, Halldór Olafsson, and Vince McIntyre.
+
 =cut
 
 ####################################################################

Libtmp/Func_demo.pm

@@ -2,7 +2,6 @@ package PDL::Demos::Func_demo;
 
 use PDL::Graphics::Simple;
 use PDL::Func;
-use PDL::Slatec; # no PCHIP without it
 
 sub info {('func', 'Interpolation etc (Req.: PDL::Graphics::Simple)')}
 

Libtmp/t/func.t

@@ -17,11 +17,6 @@ is_pdl $obj->get( 'err' ), long('1 0 0 0 0 0 0 0'), 'same error as direct';
 eval { $obj->gradient( $xi ); };
 like $@ , qr/can not call gradient/, 'calling unavailable method';
 
-if (!eval { require PDL::Slatec; PDL::Slatec->import; 1 }) {
-  done_testing;
-  exit;
-}
-
 $x = sequence(float,10);
 $y = $x*$x + 0.5;
 $obj->set( Interpolate => 'Hermite', x => $x, y => $y );
@@ -51,32 +46,4 @@ is_pdl $obj->interpolate( $xi ), $ans, {atol=>6, test_name=>'broadcasting non-si
 is_pdl pchip( $x, $y, $xi ), $ans, {atol=>6, test_name=>'pchip'};
 is_pdl spline( $x, $y, $xi ), $ans, {atol=>6, test_name=>'spline'};
 
-$x = sequence(10);
-$y = pdl('43.3 44.3 47.3 52.3 59.3 68.3 79.3 92.3 107.3 124.3; -23 -22 -15 4 41 102 193 320 489 706');
-my ($g) = chsp([0,0], [0,0], $x, $y);
-is_pdl $g, pdl('0 2 4 6 8 10 12 14 16 18; 0 3 12 27 48 75 108 147 192 243'), 'chsp';
-($g) = chic([0,0], [0,0], 0, $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'), 'chic';
-($g) = 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'), 'chim';
-my ($yi) = chfe($x, $y, $g, [1,1], $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'), 'chfe';
-($yi) = chfd($x, $y, $g, [0,0], $xi);
-is_pdl $yi, $yi_exp, 'chfd';
-my ($integral) = chia($x, $y, $g, [0,0], 3, 4);
-is_pdl $integral, pdl(55.6298611111111,20.7538265306122), 'chia';
-($integral) = chid($x, $y, $g, [0,0], 3, 4);
-is_pdl $integral, pdl(55.6298611111111,20.7538265306122), 'chid';
-my $nknots = zeroes(longlong, 2);
-my $t = zeroes( 2*$x->dim(0)+4, 2 );
-my ($bcoef, $ndim, $kord) = chbs($x, $y, $g, 0, $nknots, $t);
-is_pdl $nknots, longlong(24,24), 'chbs nknots';
-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'), 'chbs t';
-is_pdl $bcoef, pdl('43.3 43.3 43.8 44.8 46.05 48.55 50.355556 54.244444 56.675 61.925 65 71.6 75.327778 83.272222 87.657143 96.942857 101.9875 112.6125 118.3 124.3; -23 -23 -22.583333 -21.416667 -18.410256 -11.589744 -4.3690476 12.369048 25.646259 56.353741 77.653509 126.34649 157.65749 228.34251 271.65991 368.34009 425.66149 552.33851 625.66667 706'), 'chbs bcoef';
-is_pdl $ndim, longlong(20,20), 'chbs ndim';
-is_pdl $kord, longlong(4,4), 'chbs kord';
-my $x_slice = $x->slice('*1,0:-2'); # because calling with last value is out of range
-my ($val) = bvalu($t, $bcoef, 0, $x_slice);
-is_pdl $val->t, pdl('43.3 44.3 47.3 52.3 59.3 68.3 79.3 92.3 107.3; -23 -22 -15 4 41 102 193 320 489'), 'bvalu';
-
 done_testing;

MANIFEST

@@ -142,6 +142,7 @@ Basic/lib/PDL/PP/Dims.pm
 Basic/lib/PDL/PP/PDLCode.pm
 Basic/lib/PDL/PP/PdlParObj.pm
 Basic/lib/PDL/PP/Signature.pm
+Basic/lib/PDL/Primitive-pchip.c
 Basic/lib/PDL/Primitive-xoshiro256plus.c
 Basic/lib/PDL/Primitive.pd
 Basic/lib/PDL/QuickStart.pod