"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Libtmp/Slatec/slatec.pd" between
PDL-2.077.tar.gz and PDL-2.078.tar.gz

About: PDL (Perl Data Language) aims to turn perl into an efficient numerical language for scientific computing (similar to IDL and MatLab).

slatec.pd  (PDL-2.077):slatec.pd  (PDL-2.078)
skipping to change at line 1497 skipping to change at line 1497
=cut =cut
', ',
'Piecewise Cubic Hermite function to B-Spline converter.' 'Piecewise Cubic Hermite function to B-Spline converter.'
); );
################################################################## ##################################################################
################################################################## ##################################################################
#
# This version of polfit accepts bad values and also allows broadcasting # This version of polfit accepts bad values and also allows broadcasting
#
#
# indices: # indices:
# n runs across input points; # n runs across input points;
# foo runs across wacky Slatec buffer size; # foo runs across wacky Slatec buffer size;
# bar runs across polynomial coefficients. # bar runs across polynomial coefficients.
#
pp_def('polfit', pp_def('polfit',
Pars => 'x(n); y(n); w(n); longlong maxdeg(); longlong [o]ndeg(); [o]eps(); [o ]r(n); longlong [o]ierr(); [o]a(foo); [o]coeffs(bar);[t]xtmp(n);[t]ytmp(n);[t]wt mp(n);[t]rtmp(n)', Pars => 'x(n); y(n); w(n); longlong maxdeg(); longlong [o]ndeg(); [o]eps(); [o ]r(n); longlong [o]ierr(); [o]a(foo); [o]coeffs(bar);[t]xtmp(n);[t]ytmp(n);[t]wt mp(n);[t]rtmp(n)',
OtherPars => '',
Code => '
PDL_LongLong maxord, ord, k, ns = $SIZE(n);
$GENERIC() zero = 0;
PDL_FORTRAN($TFD(polfit,dpolft))(&ns,$P(x),$P(y),$P(w),$P(maxdeg),$P(
ndeg),$P(eps),$P(r),$P(ierr),$P(a));
maxord = ($P(a))[0]+0.5;
ord = ($P(a))[maxord * 3 + 2];
if(ord >= $maxdeg()) {
ord = $maxdeg();
}
PDL_FORTRAN($TFD(p,dp)coef)( &(ord), &(zero), $P(coeffs), $P(a));
for(k=ord+1; k<=$maxdeg(); k++)
($P(coeffs))[k] = 0;
',
GenericTypes => ['F','D'], GenericTypes => ['F','D'],
HandleBad => 1, HandleBad => 1,
NoBadifNaN => 1, NoBadifNaN => 1,
BadCode => 'PDL_LongLong ns = $SIZE(n), i, j = 0; Code => 'PDL_LongLong howmany = PDL_IF_BAD(0,$SIZE(n)), i;
if(ns<$maxdeg()) if($SIZE(n)<$maxdeg())
$CROAK("Must have at least <n> points to fit <n> coefficients"); $CROAK("Must have at least <n> points to fit <n> coefficients");
PDL_IF_BAD(
for (i=0;i<ns;i++) { /* get rid of bad values. Call polfit with loop(n) %{ /* get rid of bad values. Call polfit with [xyw]tmp
[xyw]tmp instead of [xyz]. */ instead of [xyz]. */
if ($ISGOOD(y(n=>i)) && $ISGOOD(x(n=>i)) && $ISGOOD(w(n=>i))) { if ($ISGOOD(y()) && $ISGOOD(x()) && $ISGOOD(w())) {
$xtmp(n=>j) = $x(n=>i); $xtmp(n=>howmany) = $x();
$ytmp(n=>j) = $y(n=>i); $ytmp(n=>howmany) = $y();
$wtmp(n=>j) = $w(n=>i); $wtmp(n=>howmany) = $w();
j++; howmany++;
} }
} %}
if (j <= $maxdeg()) { if (howmany <= $maxdeg()) {
/* Not enough good datapoints -- set this whole row to BAD. */ /* Not enough good datapoints -- set this whole row to BAD. */
for (i=0;i<ns;i++) { loop(n) %{ $SETBAD(r()); %}
$SETBAD(r(n=>i));
}
$ierr() = 2; $ierr() = 2;
} else { } else,) {
/* Found enough good datapoints for a fit -- do the fit */ /* Found enough good datapoints for a fit -- do the fit */
PDL_LongLong k;
PDL_LongLong ord;
PDL_LongLong maxord;
$GENERIC() zero = 0; $GENERIC() zero = 0;
/* Do the fit */ /* Do the fit */
PDL_FORTRAN($TFD(polfit,dpolft)) PDL_FORTRAN($TFD(polfit,dpolft))(
(&j,$P(xtmp),$P(ytmp),$P(wtmp),$P(maxdeg),$P(ndeg),$P(eps),$ &howmany,
P(rtmp),$P(ierr),$P(a)); PDL_IF_BAD($P(xtmp),$P(x)),
PDL_IF_BAD($P(ytmp),$P(y)),
maxord = ($P(a))[0]+0.5; PDL_IF_BAD($P(wtmp),$P(w)),
ord = ($P(a))[maxord * 3 + 2]; $P(maxdeg),$P(ndeg),$P(eps),
if(ord >= $maxdeg()) { PDL_IF_BAD($P(rtmp),$P(r)),
ord = $maxdeg(); $P(ierr),$P(a)
} );
PDL_LongLong maxord = $a(foo=>0)+0.5;
PDL_LongLong ord = PDLMIN($maxdeg(),$a(foo=>(maxord*3)+2));
/* Extract the polynomial coefficients into coeffs -- used for ba d values */ /* Extract the polynomial coefficients into coeffs -- used for ba d values */
PDL_FORTRAN($TFD(,d)pcoef)(&(ord), &(zero), $P(coeffs), $P(a)); PDL_FORTRAN($TFD(,d)pcoef)(&(ord), &(zero), $P(coeffs), $P(a));
for(k=ord+1; k<=$maxdeg(); k++) for(i=ord+1; i<=$maxdeg(); i++)
($P(coeffs))[k] = 0; $coeffs(bar=>i) = 0;
j=0; PDL_IF_BAD(
for (i=0;i<ns;i++) { /* replace bad values */ howmany=0;
if ($ISGOOD(y(n=>i))) { loop(n) %{ /* replace bad values */
$r(n=>i) = $rtmp(n=>j); if ($ISGOOD(y())) {
j++; $r() = $rtmp(n=>howmany);
} else if($ISGOOD(x(n=>i))) { howmany++;
} else if($ISGOOD(x())) {
/* Bad values are omitted from the call to polfit, so we mus t reconstitute them on return */ /* Bad values are omitted from the call to polfit, so we mus t reconstitute them on return */
/* (just because a value is bad in y, does not mean the fit itself is bad there) */ /* (just because a value is bad in y, does not mean the fit itself is bad there) */
/* */ /* */
/* The number in ord is not the number of coefficients in t he polynomial, it is the highest */ /* The number in ord is not the number of coefficients in t he polynomial, it is the highest */
/* order coefficient -- so 3 for a cubic, which has 4 coeff icients. */ /* order coefficient -- so 3 for a cubic, which has 4 coeff icients. */
/* --CED */
PDL_LongLong ii; PDL_LongLong ii;
$GENERIC() acc = 0; $GENERIC() acc = 0;
for( ii=ord; ii>0; ii-- ) { for( ii=ord; ii>0; ii-- ) {
acc += $coeffs(bar=>ii); acc += $coeffs(bar=>ii);
acc *= $x(n=>i); acc *= $x();
} }
acc += $coeffs(bar=>0); acc += $coeffs(bar=>0);
$r(n=>i) = acc; $r() = acc;
} else { } else {
/* x and y are bad here... */ /* x and y are bad here... */
$SETBAD(r(n=>i)); $SETBAD(r());
} }
} %},)
}', }',
Doc => 'Fit discrete data in a least squares sense by polynomials Doc => 'Fit discrete data in a least squares sense by polynomials
in one variable. C<x()>, C<y()> and C<w()> must be of the same type. in one variable. C<x()>, C<y()> and C<w()> must be of the same type.
This version handles bad values appropriately', This version handles bad values appropriately',
); );
################################################################## ##################################################################
################################################################## ##################################################################
pp_addpm(<<'EOD'); pp_addpm(<<'EOD');
=head1 AUTHOR =head1 AUTHOR
Copyright (C) 1997 Tuomas J. Lukka. Copyright (C) 1997 Tuomas J. Lukka.
Copyright (C) 2000 Tim Jenness, Doug Burke. Copyright (C) 2000 Tim Jenness, Doug Burke.
All rights reserved. There is no warranty. You are allowed All rights reserved. There is no warranty. You are allowed
to redistribute this software / documentation under certain to redistribute this software / documentation under certain
conditions. For details, see the file COPYING in the PDL conditions. For details, see the file COPYING in the PDL
distribution. If this file is separated from the PDL distribution, distribution. If this file is separated from the PDL distribution,
the copyright notice should be included in the file. the copyright notice should be included in the file.
=cut =cut
EOD EOD
pp_done(); pp_done();
 End of changes. 23 change blocks. 
66 lines changed or deleted 38 lines changed or added

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)