diff --git a/lib/PDL/Primitive.pd b/lib/PDL/Primitive.pd index 8042e2dd0..838948cfd 100644 --- 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, C, -and C. - C is the array of dependent variable values. C is the value corresponding to C. diff --git a/t/primitive-interpolate.t b/t/primitive-interpolate.t index 1596c595d..58a4094c1 100644 --- 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';