chbs more PDL-idiomatic interface

PDLPorters · Dec 17, 2024 · 49f4121 · 49f4121
1 parent b26fd1f
commit 49f4121
Show file tree

Hide file tree

Showing 2 changed files with 19 additions and 34 deletions.
diff --git a/lib/PDL/Primitive.pd b/lib/PDL/Primitive.pd
@@ -5902,33 +5902,27 @@ EOF
 );
 
 pp_def('pchip_chbs',
-  Pars => 'x(n); f(n); d(n); longlong knotyp(); longlong [io]nknots(); [io]t(tsize=CALC(2*$SIZE(n)+4));
-    [o]bcoef(bsize=CALC(2*$SIZE(n))); longlong [o]ndim(); longlong [o]kord(); indx [o]ierr()',
+  Pars => 'x(n); f(n); d(n); sbyte knotyp();
+    [o]t(nknots=CALC(2*$SIZE(n)+4));
+    [o]bcoef(ndim=CALC(2*$SIZE(n))); indx [o]ierr()',
   GenericTypes => $F,
   Code => pp_line_numbers(__LINE__, <<'EOF'),
 /* Local variables */
-PDL_Indx n = $SIZE(n), ndim = $ndim() = n << 1;
-$kord() = 4;
+PDL_Indx n = $SIZE(n), ndim = $SIZE(ndim);
 $ierr() = 0;
 /*  Check argument validity.  Set up knot sequence if OK. */
 if ($knotyp() > 2) {
   $ierr() = -1;
   $CROAK("KNOTYP GREATER THAN 2");
 }
-if ($knotyp() < 0) {
-  if ($nknots() != $ndim() + 4) {
-    $ierr() = -2;
-    $CROAK("KNOTYP.LT.0 AND NKNOTS.NE.(2*N+4)");
-  }
-} else {
+if ($knotyp() >= 0) {
 /*      Set up knot sequence. */
-  $nknots() = $ndim() + 4;
   do { /* inline dpchkt */
 /*  Set interior knots. */
     PDL_Indx j = 0;
     loop (n) %{
       j += 2;
-      $t(tsize=>j+1) = $t(tsize=>j) = $x();
+      $t(nknots=>j+1) = $t(nknots=>j) = $x();
     %}
 /*   Assertion:  At this point T(3),...,T(NDIM+2) have been set and */
 /*         J=NDIM+1. */
@@ -5937,31 +5931,31 @@ if ($knotyp() < 0) {
     $GENERIC() hend = $x(n=>n-1) - $x(n=>n-2);
     if ($knotyp() == 1) {
 /*      Extrapolate. */
-      $t(tsize=>1) = $x(n=>0) - hbeg;
-      $t(tsize=>ndim+2) = $x(n=>n-1) + hend;
+      $t(nknots=>1) = $x(n=>0) - hbeg;
+      $t(nknots=>ndim+2) = $x(n=>n-1) + hend;
     } else if ($knotyp() == 2) {
 /*      Periodic. */
-      $t(tsize=>1) = $x(n=>0) - hend;
-      $t(tsize=>ndim+2) = $x(n=>n-1) + hbeg;
+      $t(nknots=>1) = $x(n=>0) - hend;
+      $t(nknots=>ndim+2) = $x(n=>n-1) + hbeg;
     } else {
 /*      Quadruple end knots. */
-      $t(tsize=>1) = $x(n=>0);
-      $t(tsize=>ndim+2) = $x(n=>n-1);
+      $t(nknots=>1) = $x(n=>0);
+      $t(nknots=>ndim+2) = $x(n=>n-1);
     }
-    $t(tsize=>0) = $t(tsize=>1);
-    $t(tsize=>ndim+3) = $t(tsize=>ndim+2);
+    $t(nknots=>0) = $t(nknots=>1);
+    $t(nknots=>ndim+3) = $t(nknots=>ndim+2);
   } while (0); /* end inline dpchkt */
 }
 /*  Compute B-spline coefficients. */
-$GENERIC() hnew = $t(tsize=>2) - $t(tsize=>0);
+$GENERIC() hnew = $t(nknots=>2) - $t(nknots=>0);
 loop (n) %{
   $GENERIC() hold = hnew;
 /*      The following requires mixed mode arithmetic. */
   $GENERIC() dov3 = $d() / 3;
-  $bcoef(bsize=>2*n) = $f() - hold * dov3;
+  $bcoef(ndim=>2*n) = $f() - hold * dov3;
 /*      The following assumes T(2*K+1) = X(K). */
-  hnew = $t(tsize=>2*n+4) - $t(tsize=>2*n+2);
-  $bcoef(bsize=>2*n+1) = $f() + hnew * dov3;
+  hnew = $t(nknots=>2*n+4) - $t(nknots=>2*n+2);
+  $bcoef(ndim=>2*n+1) = $f() + hnew * dov3;
 %}
 EOF
   Doc => <<'EOF',
@@ -5989,13 +5983,9 @@ Restrictions/assumptions:
    1. N.GE.2 .  (not checked)
    2. X(i).LT.X(i+1), i=1,...,N .  (not checked)
    4. KNOTYP.LE.2 .  (error return if not)
-  *5. NKNOTS = NDIM+4 = 2*N+4 .  (error return if not)
   *6. T(2*k+1) = T(2*k) = X(k), k=1,...,N .  (not checked)
      * Indicates this applies only if KNOTYP.LT.0 .
 
-Array sizes: C<tsize = 2*n + 4>, C<bsize = 2*n>,
-and C<ndim = 2*n>.
-
 C<f> is the array of dependent variable values.
 C<f(I)> is the value corresponding to C<x(I)>.
 

diff --git a/t/primitive-interpolate.t b/t/primitive-interpolate.t
@@ -58,14 +58,9 @@ subtest PCHIP => sub {
   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';
+  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';
   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';