pp_addpm({At=>Top},<<'EOD');
=head1 NAME
PDL::Slatec - PDL interface to the slatec numerical programming library
=head1 SYNOPSIS
use PDL::Slatec;
($ndeg, $r, $ierr, $a) = polyfit($x, $y, $w, $maxdeg, $eps);
=head1 DESCRIPTION
This module serves the dual purpose of providing an interface to
parts of the slatec library and showing how to interface PDL
to an external library.
Using this library requires a fortran compiler; the source for the routines
is provided for convenience.
Currently available are routines to:
manipulate matrices; calculate FFT's;
fit data using polynomials;
and interpolate/integrate data using piecewise cubic Hermite interpolation.
=head2 Piecewise cubic Hermite interpolation (PCHIP)
PCHIP is the slatec package of routines to perform piecewise cubic
Hermite interpolation of data.
It features software to produce a monotone and "visually pleasing"
interpolant to monotone data.
According to Fritsch & Carlson ("Monotone piecewise
cubic interpolation", SIAM Journal on Numerical Analysis
17, 2 (April 1980), pp. 238-246),
such an interpolant may be more reasonable than a cubic spline if
the data contains both "steep" and "flat" sections.
Interpolation of cumulative probability distribution functions is
another application.
These routines are cryptically named (blame FORTRAN),
beginning with 'ch', and accept either float or double piddles.
Most of the routines require an integer parameter called C<check>;
if set to 0, then no checks on the validity of the input data are
made, otherwise these checks are made.
The value of C<check> can be set to 0 if a routine
such as L<chim|/chim> has already been successfully called.
=over 4
=item *
If not known, estimate derivative values for the points
using the L<chim|/chim>, L<chic|/chic>, or L<chsp|/chsp> routines
(the following routines require both the function (C<f>)
and derivative (C<d>) values at a set of points (C<x>)).
=item *
Evaluate, integrate, or differentiate the resulting PCH
function using the routines:
L<chfd|/chfd>; L<chfe|/chfe>; L<chia|/chia>; L<chid|/chid>.
=item *
If desired, you can check the monotonicity of your
data using L<chcm|/chcm>.
=back
=cut
EOD
# ' un-confuse emacs
# if define chbs, then add something like the following to point 3:
#
# or use L<chbs|/chbs> to convert a PCH function into B-representation
# for use with the B-spline routines of slatec
# (although no interface to them currently exist).
#
# add function definitions after finishing the first pp_addpm(), since this
# adds a '=head1 FUNCTIONS' line at the end of the text
pp_addpm(<<'END');
=head2 eigsys
=for ref
Eigenvalues and eigenvectors of a real positive definite symmetric matrix.
=for usage
($eigvals,$eigvecs) = eigsys($mat)
Note: this function should be extended to calculate only eigenvalues if called
in scalar context!
=head2 matinv
=for ref
Inverse of a square matrix
=for usage
($inv) = matinv($mat)
=head2 polyfit
Convenience wrapper routine about the C<polfit> C<slatec> function.
Separates supplied arguments and return values.
=for ref
Fit discrete data in a least squares sense by polynomials
in one variable. Handles threading correctly--one can pass in a 2D PDL (as C<$y>)
and it will pass back a 2D PDL, the rows of which are the polynomial regression
results (in C<$r> corresponding to the rows of $y.
=for usage
($ndeg, $r, $ierr, $a, $coeffs, $rms) = polyfit($x, $y, $w, $maxdeg, [$eps]);
$coeffs = polyfit($x,$y,$w,$maxdeg,[$eps]);
where on input:
C<$x> and C<$y> are the values to fit to a polynomial.
C<$w> are weighting factors
C<$maxdeg> is the maximum degree of polynomial to use and
C<$eps> is the required degree of fit.
and the output switches on list/scalar context.
In list context:
C<$ndeg> is the degree of polynomial actually used
C<$r> is the values of the fitted polynomial
C<$ierr> is a return status code, and
C<$a> is some working array or other (preserved for historical purposes)
C<$coeffs> is the polynomial coefficients of the best fit polynomial.
C<$rms> is the rms error of the fit.
In scalar context, only $coeffs is returned.
Historically, C<$eps> was modified in-place to be a return value of the
rms error. This usage is deprecated, and C<$eps> is an optional parameter now.
It is still modified if present.
C<$a> is a working array accessible to Slatec - you can feed it to several
other Slatec routines to get nice things out. It does not thread
correctly and should probably be fixed by someone. If you are
reading this, that someone might be you.
=for bad
This version of polyfit handles bad values correctly. Bad values in
$y are ignored for the fit and give computed values on the fitted
curve in the return. Bad values in $x or $w are ignored for the fit and
result in bad elements in the output.
=head2 polycoef
Convenience wrapper routine around the C<pcoef> C<slatec> function.
Separates supplied arguments and return values.
=for ref
Convert the C<polyfit>/C<polfit> coefficients to Taylor series form.
=for usage
$tc = polycoef($l, $c, $a);
=head2 polyvalue
Convenience wrapper routine around the C<pvalue> C<slatec> function.
Separates supplied arguments and return values.
For multiple input x positions, a corresponding y position is calculated.
The derivatives PDL is one dimensional (of size C<nder>) if a single x
position is supplied, two dimensional if more than one x position is
supplied.
=for ref
Use the coefficients generated by C<polyfit> (or C<polfit>) to evaluate
the polynomial fit of degree C<l>, along with the first C<nder> of its
derivatives, at a specified point.
=for usage
($yfit, $yp) = polyvalue($l, $nder, $x, $a);
=head2 detslatec
=for ref
compute the determinant of an invertible matrix
=for example
$mat = zeroes(5,5); $mat->diagonal(0,1) .= 1; # unity matrix
$det = detslatec $mat;
Usage:
=for usage
$determinant = detslatec $matrix;
=for sig
Signature: detslatec(mat(n,m); [o] det())
C<detslatec> computes the determinant of an invertible matrix and barfs if
the matrix argument provided is non-invertible. The matrix threads as usual.
This routine was previously known as C<det> which clashes now with
L<det|PDL::MatrixOps/det> which is provided by L<PDL::MatrixOps>.
=head2 fft
=for ref
Fast Fourier Transform
=for example
$v_in = pdl(1,0,1,0);
($azero,$a,$b) = PDL::Slatec::fft($v_in);
C<PDL::Slatec::fft> is a convenience wrapper for L<ezfftf|ezfftf>, and
performs a Fast Fourier Transform on an input vector C<$v_in>. The
return values are the same as for L<ezfftf|ezfftf>.
=head2 rfft
=for ref
reverse Fast Fourier Transform
=for example
$v_out = PDL::Slatec::rfft($azero,$a,$b);
print $v_in, $vout
[1 0 1 0] [1 0 1 0]
C<PDL::Slatec::rfft> is a convenience wrapper for L<ezfftb|ezfftb>,
and performs a reverse Fast Fourier Transform. The input is the same
as the output of L<PDL::Slatec::fft|/PDL::Slatec::fft>, and the output
of C<rfft> is a data vector, similar to what is input into
L<PDL::Slatec::fft|/PDL::Slatec::fft>.
=cut
END
use strict;
# for MDim, ld[str] is interpreted as "leading dimension of ..."
# Making the BAD BAD BAD assumption that PDL_Long == int
# in fortran. BAD BAD BAD XXX (I'm going to regret this)
my %ftypes = (S => 'F', D => 'D');
sub firstpar {
$_[0] =~ /^\(([^),]+)[),]/ or die "Can't find first par from $_[0]";
$1
}
# whether or not to append undercores
my $uscore = (-e "f77_underscore" ? "_" : "");
# used in defslatec()
#my %ignore_ppar = ( Incfd => 1, CheckFlag => 1 );
my %ignore_ppar = ( Incfd => 1 );
my %prototype = ( F => "float", D => "double" );
# an alternative is to declare the function in the Code section
# of pp_def(), using something like:
#
# my $codeproto = "\$T".(join '',map {$_->[0]} @talts)."(".
# (join ',',map {$_->[1].$uscore} @talts).") ();";
# if ( defined $fpar ) {
# $codeproto = "\$T".(join '',map {$_->[0]} @talts)."(float,double) $codeproto";
# }
# $codeproto = "extern $codeproto";
#
# and then add `$codeproto . "\n" .' to the beginning of the Code
# section.
#
# this then gets rid of the need of the prototype file.
#
open( PROTOS, "> SlatecProtos.h" );
# defslatec( $pdlname, $funcnames, $argstr, $docstring, $funcref )
#
# $pdlname is the name of the PDL function to be created
# $funcnames is a reference to a hash array, whose keys define
# the single (S), and double precision (D) names of the
# SLATEC routines to be linked to.
#
# $argstr is a list of arguments expected by the SLATEC routine
# - some of the allowed type names are:
# FuncRet
# - specifies that this is a function, not a subroutine, and
# that the output of the function should be stored in this
# variable
# Incfd
# - used in the PCHIP functions to specify the INCFD argument
# that we force to be 1, so the user never has to specify it
# (this allows the PCHIP routines to use 2D data, but as it's
# done in FORTRAN array order, and PDL has a much richer way
# of accessing parts of an array we force the data to be 1D).
# CheckFlag
# - the PCHIP routined may change the value from 0 to 1 if an
# error occurs but the checks were successful. As this complicates
# things we copy the user value to a temporary variable,
# so that the sent in value is not changed.
# FortranIndex
# - pchid()/dpchid() require FORTRAN array indices, so
# this type flags that we should add 1 onto the input values
# before sending to the slatec function
#
# $docstring gives the text to be used as the function dicumentation
#
# $funcref gets placed within a '=for ref' pod statement at the
# start of the documentation - ie it is placed before the
# text within $docstring. This string gets printed out
# in the perldl or pdl2 shell after a '?? string' command
#
sub defslatec {
my $debug = 0; # print out calls to pp_def
my($pname,$fnames,$argstr,$docstring,$funcref) = @_;
my @args = map {/^\s*$/ ? () : $_} split ';', $argstr;
my @args2 = map {
/^\s*([a-zA-Z]+)\s+ # "Type name"
((?:\[[^]]*\])?)\s* # Options
([a-zA-Z]+)\s* # Par name
((?:\([^)]*\))?)\s*$ # Dims
/x or die("Invalid slatec par $_");
[$1,$2,$3,$4]} @args;
# is this for a function (Type name eq "FuncRet")
# or a subroutine?
my $fpar;
foreach ( @args2 ) {
next unless $_->[0] eq "FuncRet";
die "Only one FuncRet allowed in pars list.\n" if defined $fpar;
$fpar = "\$$_->[2]()";
}
my @ppars = map {
if($_->[0] =~ /^M?Dim$/ or defined $ignore_ppar{$_->[0]} ) {
()
} else {
(($_->[0] eq "Mat" or $_->[0] eq "FuncRet")
and join '',@{$_}[1,2,3]) or
(($_->[0] eq "IntFlag" or $_->[0] eq "FortranIndex" or $_->[0] eq "CheckFlag")
and "int ".join '',@{$_}[1,2,3]) or
die "Invalid ppars ",(join ',',@$_),"\n";
}
} @args2;
# uncomment the following line to see what perl thinks the input pars are
##print "Pars: ",(join ';',@ppars),"\n";
my @talts = map {
defined $ftypes{$_} or die "FTYPE $_ NOT THERE\n";
[$ftypes{$_},$fnames->{$_}]
} sort keys %$fnames;
my $func = "\$T".(join '',map {$_->[0]} @talts) . "(" .
(join ',',map {$_->[1].$uscore} @talts).")";
if ( defined $fpar ) { $func = "$fpar = $func"; }
my %lds = map {
($_->[0] eq "Mat" and $_->[3] ne "()") ?
("ld".$_->[2] => "&\$PRIV(__".firstpar($_->[3])."_size)")
: ()
} @args2;
my @funcpars;
foreach ( @args2 ) {
next if $_->[0] eq "FuncRet";
if ( $_->[0] eq "Mat" or $_->[0] eq "IntFlag" ) {
push @funcpars, "\$P($_->[2])";
} elsif ( $_->[0] eq "Dim" ) {
push @funcpars, "&\$PRIV(__$_->[2]_size)";
} elsif ( $_->[0] eq "MDim" ) {
push @funcpars, $lds{$_->[2]};
} elsif ( $_->[0] eq "Incfd" or $_->[0] eq "CheckFlag" ) {
push @funcpars, "&_" . lc($_->[0]);
} elsif ( $_->[0] eq "FortranIndex" ) {
push @funcpars, "&_$_->[2]";
} else {
die "Invalid args2";
}
}
# _incfd = 1 makes sure PCHIP code treats piddle as 1D
# _checkflag - copy input data to a temporary variable, in case
# the PCHIP routine decides to change it
#
my @ifincode;
foreach ( @args2 ) {
if ( $_->[0] eq "Incfd" ) {
push @ifincode, "int _" . lc($_->[0]) . " = 1;";
} elsif ( $_->[0] eq "CheckFlag" ) {
push @ifincode, "int _" . lc($_->[0]) . " = \$$_->[2]();";
} elsif ( $_->[0] eq "FortranIndex" ) {
# convert from C to F77 index
push @ifincode, "int _$_->[2] = \$$_->[2]() + 1;"
}
}
foreach ( @talts ) {
my $codeproto = "extern ";
if ( defined $fpar ) { $codeproto .= "$prototype{$_->[0]} "; }
else { $codeproto .= "int "; }
$codeproto .= "$_->[1]$uscore ();";
print PROTOS $codeproto . "\n";
}
# add on the function reference, if supplied, to the start of
# the doc string
if ( defined $docstring ) {
$docstring = "\n=for ref\n\n$funcref\n\n$docstring" if defined $funcref;
} else {
$docstring = '';
}
# If debug flag set, then print out pp_def call for each call to defslatec
if ($debug) {
my $pars = (join ';',@ppars);
my $code = (join '',@ifincode) . "\n " . $func . " (". (join ',',@funcpars) . ");\n";
my $generictypes = "[" . join (", ", map {$_->[0]} @talts) . "],\n";
print <<"ENDDBG";
pp_def($pname,
Pars => $pars,
OtherPars => '',
Code => $code,
GenericTypes => $generictypes,
Doc => $docstring
);
ENDDBG
}
pp_def($pname,
Pars => (join ';',@ppars),
OtherPars => '',
Code => (join '',@ifincode) . "\n " .
$func . " (". (join ',',@funcpars) . ");\n",
# . (join '',@ifoutcode),
GenericTypes => [map {$_->[0]} @talts],
Doc => $docstring
# %$opts,
);
} # sub: defslatec()
pp_addhdr(qq|
#include "SlatecProtos.h"
void MAIN__ () {
/* Cheat to define MAIN__ symbol */
croak("This should never happen");
}
void slatecbarf$uscore() {
croak("slatec called halt");
}
|);
pp_add_exported('',"eigsys matinv polyfit polycoef polyvalue");
pp_addpm(<<'END');
use PDL::Core;
use PDL::Basic;
use PDL::Primitive;
use PDL::Ufunc;
use strict;
# Note: handles only real symmetric positive-definite.
*eigsys = \&PDL::eigsys;
sub PDL::eigsys {
my($h) = @_;
$h = float($h);
rs($h,
(my $eigval=PDL->null),
(long (pdl (1))),(my $eigmat=PDL->null),
(my $fvone = PDL->null),(my $fvtwo = PDL->null),
(my $errflag=PDL->null)
);
# print $covar,$eigval,$eigmat,$fvone,$fvtwo,$errflag;
if(sum($errflag) > 0) {
barf("Non-positive-definite matrix given to eigsys: $h\n");
}
return ($eigval,$eigmat);
}
*matinv = \&PDL::matinv;
sub PDL::matinv {
my($m) = @_;
my(@dims) = $m->dims;
# Keep from dumping core (FORTRAN does no error checking)
barf("matinv requires a 2-D square matrix")
unless( @dims >= 2 && $dims[0] == $dims[1] );
$m = $m->copy(); # Make sure we don't overwrite :(
gefa($m,(my $ipvt=null),(my $info=null));
if(sum($info) > 0) {
barf("Uninvertible matrix given to inv: $m\n");
}
gedi($m,$ipvt,(pdl 0,0),(null),(long( pdl (1))));
$m;
}
*detslatec = \&PDL::detslatec;
sub PDL::detslatec {
my($m) = @_;
$m = $m->copy(); # Make sure we don't overwrite :(
gefa($m,(my $ipvt=null),(my $info=null));
if(sum($info) > 0) {
barf("Uninvertible matrix given to inv: $m\n");
}
gedi($m,$ipvt,(my $det=null),(null),(long( pdl (10))));
return $det->slice('(0)')*10**$det->slice('(1)');
}
sub prepfft {
my($n) = @_;
my $tmp = PDL->zeroes(float(),$n*3+15);
$n = pdl $n;
ezffti($n,$tmp);
return $tmp;
}
sub fft (;@) {
my($v) = @_;
my $ws = prepfft($v->getdim(0));
ezfftf($v,(my $az = PDL->null), (my $a = PDL->null),
(my $b = PDL->null), $ws);
return ($az,$a,$b);
}
sub rfft {
my($az,$a,$b) = @_;
my $ws = prepfft($a->getdim(0));
my $v = $a->copy();
ezfftb($v,$az,$a,$b,$ws);
return $v;
}
# polynomial fitting routines
# simple wrappers around the SLATEC implementations
*polyfit = \&PDL::polyfit;
sub PDL::polyfit {
barf 'Usage: polyfit($x, $y, $w, $maxdeg, [$eps]);'
unless (@_ == 5 || @_==4 );
my ($x_in, $y_in, $w_in, $maxdeg_in, $eps_in) = @_;
# if $w_in does not match the data vectors ($x_in, $y_in), then we can resize
# it to match the size of $y_in :
$w_in = $w_in + $y_in->zeros;
# Create the output arrays
my $r = PDL->null;
# A array needs some work space
my $sz = ((3 * $x_in->getdim(0)) + (3*$maxdeg_in) + 3); # Buffer size called for by Slatec
my @otherdims = $_[0]->dims; shift @otherdims; # Thread dims
my $a = PDL::new_from_specification('PDL',$x_in->type,$sz,@otherdims);
my $coeffs = PDL::new_from_specification('PDL',$x_in->type, $maxdeg_in + 1, @otherdims);
my $ierr = PDL->null;
my $ndeg = PDL->null;
# Now call polfit
my $rms = pdl($eps_in);
polfit($x_in, $y_in, $w_in, $maxdeg_in, $ndeg, $rms, $r, $ierr, $a, $coeffs);
# Preserve historic compatibility by flowing rms error back into the argument
if( UNIVERSAL::isa($eps_in,'PDL') ){
$eps_in .= $rms;
}
# Return the arrays
if(wantarray) {
return ($ndeg, $r, $ierr, $a, $coeffs, $rms );
} else {
return $coeffs;
}
}
*polycoef = \&PDL::polycoef;
sub PDL::polycoef {
barf 'Usage: polycoef($l, $c, $a);'
unless @_ == 3;
# Allocate memory for return PDL
# Simply l + 1 but cant see how to get PP to do this - TJ
# Not sure the return type since I do not know
# where PP will get the information from
my $tc = PDL->zeroes( abs($_[0]) + 1 );
# Run the slatec routine
pcoef($_[0], $_[1], $tc, $_[2]);
# Return results
return $tc;
}
*polyvalue = \&PDL::polyvalue;
sub PDL::polyvalue {
barf 'Usage: polyvalue($l, $nder, $x, $a);'
unless @_ == 4;
# Two output arrays
my $yfit = PDL->null;
# This one must be preallocated and must take into account
# the size of $x if greater than 1
my $yp;
if ($_[2]->getdim(0) == 1) {
$yp = $_[2]->zeroes($_[1]);
} else {
$yp = $_[2]->zeroes($_[1], $_[2]->getdim(0));
}
# Run the slatec function
pvalue($_[0], $_[2], $yfit, $yp, $_[3]);
# Returns
return ($yfit, $yp);
}
END
defslatec(
'svdc',{S => 'ssvdc'},
'Mat x (n,p);
MDim ldx;
Dim n;
Dim p;
Mat [o] s (p);
Mat [o] e (p);
Mat [o] u (n,p);
MDim ldu;
Mat [o] v (p,p);
MDim ldv;
Mat [o] work (n);
IntFlag job ();
IntFlag [o] info ();
',
'singular value decomposition of a matrix'
);
defslatec(
'poco',{S => 'spoco', D => 'dpoco'},
'Mat a (n,n);
MDim lda;
Dim n;
Mat rcond ();
Mat [o] z (n);
IntFlag [o] info ();
',
'Factor a real symmetric positive definite matrix
and estimate the condition number of the matrix.'
);
defslatec(
'geco',{S => 'sgeco', D => 'dgeco'},
'Mat a (n,n);
MDim lda;
Dim n;
IntFlag [o] ipvt (n);
Mat [o] rcond ();
Mat [o] z (n);
',
'Factor a matrix using Gaussian elimination and estimate
the condition number of the matrix.'
);
defslatec(
'gefa',{S => 'sgefa', D => 'dgefa'},
'Mat a (n,n);
MDim lda;
Dim n;
IntFlag [o] ipvt (n);
IntFlag [o] info ();
',
'Factor a matrix using Gaussian elimination.'
);
# XXX Ensure two == 2!!
#
# pofa and sqrdc aren't (yet?) implemented
#
defslatec(
'podi',{S => 'spodi', D => 'dpodi'},
'Mat a (n,n);
MDim lda;
Dim n;
Mat [o] det (two=2);
IntFlag job ();
',
'Compute the determinant and inverse of a certain real
symmetric positive definite matrix using the factors
computed by L<poco|/poco>.'
);
defslatec(
'gedi',{S => 'sgedi', D => 'dgedi'},
'Mat a (n,n);
MDim lda;
Dim n;
IntFlag [o] ipvt (n);
Mat [o] det (two=2);
Mat [o] work (n);
IntFlag job ();
',
'Compute the determinant and inverse of a matrix using the
factors computed by L<geco|/geco> or L<gefa|/gefa>.'
);
defslatec(
'gesl',{S => 'sgesl', D => 'dgesl'},
'Mat a (lda,n);
MDim lda;
Dim n;
IntFlag ipvt (n);
Mat b (n);
IntFlag job ();
',
'Solve the real system C<A*X=B> or C<TRANS(A)*X=B> using the
factors computed by L<geco|/geco> or L<gefa|/gefa>.'
);
defslatec(
'rs', {S => 'rsfoo'},
'MDim lda;
Dim n;
Mat a (n,n);
Mat [o] w (n);
IntFlag matz ();
Mat [o] z (n,n);
Mat [t] fvone (n);
Mat [t] fvtwo (n);
IntFlag [o] ierr ();
',
'This subroutine calls the recommended sequence of
subroutines from the eigensystem subroutine package (EISPACK)
to find the eigenvalues and eigenvectors (if desired)
of a REAL SYMMETRIC matrix.'
);
# XXX wsave : at least 3n+15
defslatec(
'ezffti', {S => 'ezffti'},
'IntFlag n ();
Mat [o] wsave(foo);
',
'Subroutine ezffti initializes the work array C<wsave()>
which is used in both L<ezfftf|/ezfftf> and
L<ezfftb|/ezfftb>.
The prime factorization
of C<n> together with a tabulation of the trigonometric functions
are computed and stored in C<wsave()>.'
);
# XXX Correct for azero, a and b
defslatec(
'ezfftf', {S => 'ezfftf'},
'Dim n;
Mat r(n);
Mat [o] azero();
Mat [o] a(n);
Mat [o] b(n);
Mat wsave(foo);
'
);
defslatec(
'ezfftb', {S => 'ezfftb'},
'Dim n;
Mat [o] r(n);
Mat azero();
Mat a(n);
Mat b(n);
Mat wsave(foo);
'
);
##################################################################
##################################################################
defslatec(
'pcoef', {S => 'pcoef', D => 'dpcoef'},
'
IntFlag l ();
Mat c ();
Mat [o] tc (bar);
Mat a (foo);
',
'Convert the C<polfit> coefficients to Taylor series form.
C<c> and C<a()> must be of the same type.'
);
defslatec(
'pvalue', {S => 'pvalue', D => 'dp1vlu'},
'
IntFlag l ();
Dim nder;
Mat x ();
Mat [o] yfit ();
Mat [o] yp (nder);
Mat a (foo);
',
'Use the coefficients generated by C<polfit> to evaluate the
polynomial fit of degree C<l>, along with the first C<nder> of
its derivatives, at a specified point. C<x> and C<a> must be of the
same type.'
);
##################################################################
##################################################################
#
# PCHIP library
#
defslatec(
'chim', {S => 'pchim', D => 'dpchim'},
'Dim n;
Mat x (n);
Mat f (n);
Mat [o] d (n);
Incfd dummy;
IntFlag [o] ierr ();
',
'Calculate the derivatives at the given set of points (C<$x,$f>,
where C<$x> is strictly increasing).
The resulting set of points - C<$x,$f,$d>, referred to
as the cubic Hermite representation - can then be used in
other functions, such as L<chfe|/chfe>, L<chfd|/chfd>,
and L<chia|/chia>.
The boundary conditions are compatible with monotonicity,
and if the data are only piecewise monotonic, the interpolant
will have an extremum at the switch points; for more control
over these issues use L<chic|/chic>.
Error status returned by C<$ierr>:
=over 4
=item *
0 if successful.
=item *
E<gt> 0 if there were C<ierr> switches in the direction of
monotonicity (data still valid).
=item *
-1 if C<nelem($x) E<lt> 2>.
=item *
-3 if C<$x> is not strictly increasing.
=back
=cut
',
'Calculate the derivatives of (x,f(x)) using cubic Hermite interpolation.'
);
# switch has become mflag, since `switch' is a reserved word in
# C.
#
# can not say (nwk=2*n) --- the rhs has to equal a number
# -> could Basic/Gen/PP/Dims.pm be hacked to allow this?
#
# I didn't have much success with preceding wk by [t]
#
defslatec(
'chic', {S => 'pchic', D => 'dpchic'},
'IntFlag ic (two=2);
Mat vc (two=2);
Mat mflag ();
Dim n;
Mat x (n);
Mat f (n);
Mat [o] d (n);
Incfd dummy;
Mat wk (nwk);
Dim nwk;
IntFlag [o] ierr ();
',
'Calculate the derivatives at the given points (C<$x,$f>,
where C<$x> is strictly increasing).
Control over the boundary conditions is given by the
C<$ic> and C<$vc> piddles, and the value of C<$mflag> determines
the treatment of points where monotoncity switches
direction. A simpler, more restricted, interface is available
using L<chim|/chim>.
The first and second elements of C<$ic> determine the boundary
conditions at the start and end of the data respectively.
If the value is 0, then the default condition, as used by
L<chim|/chim>, is adopted.
If greater than zero, no adjustment for monotonicity is made,
otherwise if less than zero the derivative will be adjusted.
The allowed magnitudes for C<ic(0)> are:
=over 4
=item *
1 if first derivative at C<x(0)> is given in C<vc(0)>.
=item *
2 if second derivative at C<x(0)> is given in C<vc(0)>.
=item *
3 to use the 3-point difference formula for C<d(0)>.
(Reverts to the default b.c. if C<n E<lt> 3>)
=item *
4 to use the 4-point difference formula for C<d(0)>.
(Reverts to the default b.c. if C<n E<lt> 4>)
=item *
5 to set C<d(0)> so that the second derivative is
continuous at C<x(1)>.
(Reverts to the default b.c. if C<n E<lt> 4>)
=back
The values for C<ic(1)> are the same as above, except that
the first-derivative value is stored in C<vc(1)> for cases 1 and 2.
The values of C<$vc> need only be set if options 1 or 2 are chosen
for C<$ic>.
Set C<$mflag = 0> if interpolant is required to be monotonic in
each interval, regardless of the data. This causes C<$d> to be
set to 0 at all switch points. Set C<$mflag> to be non-zero to
use a formula based on the 3-point difference formula at switch
points. If C<$mflag E<gt> 0>, then the interpolant at swich points
is forced to not deviate from the data by more than C<$mflag*dfloc>,
where C<dfloc> is the maximum of the change of C<$f> on this interval
and its two immediate neighbours.
If C<$mflag E<lt> 0>, no such control is to be imposed.
The piddle C<$wk> is only needed for work space. However, I could
not get it to work as a temporary variable, so you must supply
it; it is a 1D piddle with C<2*n> elements.
Error status returned by C<$ierr>:
=over 4
=item *
0 if successful.
=item *
1 if C<ic(0) E<lt> 0> and C<d(0)> had to be adjusted for
monotonicity.
=item *
2 if C<ic(1) E<lt> 0> and C<d(n-1)> had to be adjusted
for monotonicity.
=item *
3 if both 1 and 2 are true.
=item *
-1 if C<n E<lt> 2>.
=item *
-3 if C<$x> is not strictly increasing.
=item *
-4 if C<abs(ic(0)) E<gt> 5>.
=item *
-5 if C<abs(ic(1)) E<gt> 5>.
=item *
-6 if both -4 and -5 are true.
=item *
-7 if C<nwk E<lt> 2*(n-1)>.
=back
=cut
',
'Calculate the derivatives of (x,f(x)) using cubic Hermite interpolation.'
);
# as above, have made wk an actual piddle, rather than a [t]
defslatec(
'chsp', {S => 'pchsp', D => 'dpchsp'},
'IntFlag ic (two=2);
Mat vc (two=2);
Dim n;
Mat x (n);
Mat f (n);
Mat [o] d (n);
Incfd dummy;
Mat wk (nwk);
Dim nwk;
IntFlag [o] ierr ();
',
'Calculate the derivatives, using cubic spline interpolation,
at the given points (C<$x,$f>), with the specified
boundary conditions.
Control over the boundary conditions is given by the
C<$ic> and C<$vc> piddles.
The resulting values - C<$x,$f,$d> - can
be used in all the functions which expect a cubic
Hermite function.
The first and second elements of C<$ic> determine the boundary
conditions at the start and end of the data respectively.
The allowed values for C<ic(0)> are:
=over 4
=item *
0 to set C<d(0)> so that the third derivative is
continuous at C<x(1)>.
=item *
1 if first derivative at C<x(0)> is given in C<vc(0>).
=item *
2 if second derivative at C<x(0>) is given in C<vc(0)>.
=item *
3 to use the 3-point difference formula for C<d(0)>.
(Reverts to the default b.c. if C<n E<lt> 3>.)
=item *
4 to use the 4-point difference formula for C<d(0)>.
(Reverts to the default b.c. if C<n E<lt> 4>.)
=back
The values for C<ic(1)> are the same as above, except that
the first-derivative value is stored in C<vc(1)> for cases 1 and 2.
The values of C<$vc> need only be set if options 1 or 2 are chosen
for C<$ic>.
The piddle C<$wk> is only needed for work space. However, I could
not get it to work as a temporary variable, so you must supply
it; it is a 1D piddle with C<2*n> elements.
Error status returned by C<$ierr>:
=over 4
=item *
0 if successful.
=item *
-1 if C<nelem($x) E<lt> 2>.
=item *
-3 if C<$x> is not strictly increasing.
=item *
-4 if C<ic(0) E<lt> 0> or C<ic(0) E<gt> 4>.
=item *
-5 if C<ic(1) E<lt> 0> or C<ic(1) E<gt> 4>.
=item *
-6 if both of the above are true.
=item *
-7 if C<nwk E<lt> 2*n>.
=item *
-8 in case of trouble solving the linear system
for the interior derivative values.
=back
=cut
',
'Calculate the derivatives of (x,f(x)) using cubic spline interpolation.'
);
defslatec(
'chfd', {S => 'pchfd', D => 'dpchfd'},
'Dim n;
Mat x (n);
Mat f (n);
Mat d (n);
Incfd dummy;
CheckFlag check ();
Dim ne;
Mat xe (ne);
Mat [o] fe (ne);
Mat [o] de (ne);
IntFlag [o] ierr ();
',
'Given a piecewise cubic Hermite function - such as from
L<chim|/chim> - evaluate the function (C<$fe>) and
derivative (C<$de>) at a set of points (C<$xe>).
If function values alone are required, use L<chfe|/chfe>.
Set C<check> to 0 to skip checks on the input data.
Error status returned by C<$ierr>:
=over 4
=item *
0 if successful.
=item *
E<gt>0 if extrapolation was performed at C<ierr> points
(data still valid).
=item *
-1 if C<nelem($x) E<lt> 2>
=item *
-3 if C<$x> is not strictly increasing.
=item *
-4 if C<nelem($xe) E<lt> 1>.
=item *
-5 if an error has occurred in a lower-level routine,
which should never happen.
=back
=cut
',
'Interpolate function and derivative values.'
);
defslatec(
'chfe', {S => 'pchfe', D => 'dpchfe'},
'Dim n;
Mat x (n);
Mat f (n);
Mat d (n);
Incfd dummy;
CheckFlag check ();
Dim ne;
Mat xe (ne);
Mat [o] fe (ne);
IntFlag [o] ierr ();
',
'Given a piecewise cubic Hermite function - such as from
L<chim|/chim> - evaluate the function (C<$fe>) at
a set of points (C<$xe>).
If derivative values are also required, use L<chfd|/chfd>.
Set C<check> to 0 to skip checks on the input data.
Error status returned by C<$ierr>:
=over 4
=item *
0 if successful.
=item *
E<gt>0 if extrapolation was performed at C<ierr> points
(data still valid).
=item *
-1 if C<nelem($x) E<lt> 2>
=item *
-3 if C<$x> is not strictly increasing.
=item *
-4 if C<nelem($xe) E<lt> 1>.
=back
=cut
',
'Interpolate function values.'
);
defslatec(
'chia', {S => 'pchia', D => 'dpchia'},
'Dim n;
Mat x (n);
Mat f (n);
Mat d (n);
Incfd dummy;
CheckFlag check ();
Mat a ();
Mat b ();
FuncRet [o] ans ();
IntFlag [o] ierr ();
',
'Evaluate the definite integral of a a piecewise
cubic Hermite function over an arbitrary interval,
given by C<[$a,$b]>.
See L<chid|/chid> if the integration limits are
data points.
Set C<check> to 0 to skip checks on the input data.
The values of C<$a> and C<$b> do not have
to lie within C<$x>, although the resulting integral
value will be highly suspect if they are not.
Error status returned by C<$ierr>:
=over 4
=item *
0 if successful.
=item *
1 if C<$a> lies outside C<$x>.
=item *
2 if C<$b> lies outside C<$x>.
=item *
3 if both 1 and 2 are true.
=item *
-1 if C<nelem($x) E<lt> 2>
=item *
-3 if C<$x> is not strictly increasing.
=item *
-4 if an error has occurred in a lower-level routine,
which should never happen.
=back
=cut
',
'Integrate (x,f(x)) over arbitrary limits.'
);
defslatec(
'chid', {S => 'pchid', D => 'dpchid'},
'Dim n;
Mat x (n);
Mat f (n);
Mat d (n);
Incfd dummy;
CheckFlag check ();
FortranIndex ia ();
FortranIndex ib ();
FuncRet [o] ans ();
IntFlag [o] ierr ();
',
'Evaluate the definite integral of a a piecewise
cubic Hermite function between C<x($ia)> and
C<x($ib)>.
See L<chia|/chia> for integration between arbitrary
limits.
Although using a fortran routine, the values of
C<$ia> and C<$ib> are zero offset.
Set C<check> to 0 to skip checks on the input data.
Error status returned by C<$ierr>:
=over 4
=item *
0 if successful.
=item *
-1 if C<nelem($x) E<lt> 2>.
=item *
-3 if C<$x> is not strictly increasing.
=item *
-4 if C<$ia> or C<$ib> is out of range.
=back
=cut
',
'Integrate (x,f(x)) between data points.'
);
defslatec(
'chcm', {S => 'pchcm', D => 'dpchcm'},
'Dim n;
Mat x (n);
Mat f (n);
Mat d (n);
Incfd dummy;
CheckFlag check ();
IntFlag [o] ismon (n);
IntFlag [o] ierr ();
',
'The outout piddle C<$ismon> indicates over
which intervals the function is monotonic.
Set C<check> to 0 to skip checks on the input data.
For the data interval C<[x(i),x(i+1)]>, the
values of C<ismon(i)> can be:
=over 4
=item *
-3 if function is probably decreasing
=item *
-1 if function is strictly decreasing
=item *
0 if function is constant
=item *
1 if function is strictly increasing
=item *
2 if function is non-monotonic
=item *
3 if function is probably increasing
=back
If C<abs(ismon(i)) == 3>, the derivative values are
near the boundary of the monotonicity region. A small
increase produces non-monotonicity, whereas a decrease
produces strict monotonicity.
The above applies to C<i = 0 .. nelem($x)-1>. The last element of
C<$ismon> indicates whether
the entire function is monotonic over $x.
Error status returned by C<$ierr>:
=over 4
=item *
0 if successful.
=item *
-1 if C<n E<lt> 2>.
=item *
-3 if C<$x> is not strictly increasing.
=back
=cut
',
'Check the given piecewise cubic Hermite function for monotonicity.'
);
=pod
ignore function for now although code is in slatec/ directory
=cut
# XXX tsize = 2*n+4
# bsize = 2*n
#
# ndim gets set to 2*n
#
# Changed by routine:
# nknots
# t
defslatec(
'chbs', {S => 'pchbs', D => 'dpchbs'},
'Dim n;
Mat x (n);
Mat f (n);
Mat d (n);
Incfd dummy;
IntFlag knotyp ();
IntFlag nknots ();
Mat t (tsize);
Mat [o] bcoef (bsize);
IntFlag [o] ndim ();
IntFlag [o] kord ();
IntFlag [o] ierr ();
',
'The resulting B-spline representation of the data
(i.e. C<nknots>, C<t>, C<bcoeff>, C<ndim>, and
C<kord>) can be evaluated by C<bvalu> (which is
currently not available).
Array sizes: C<tsize = 2*n + 4>, C<bsize = 2*n>,
and C<ndim = 2*n>.
C<knotyp> is a flag which controls the knot sequence.
The knot sequence C<t> is normally computed from C<$x>
by putting a double knot at each C<x> and setting the end knot pairs
according to the value of C<knotyp> (where C<m = ndim = 2*n>):
=over
=item *
0 - Quadruple knots at the first and last points.
=item *
1 - Replicate lengths of extreme subintervals:
C<t( 0 ) = t( 1 ) = x(0) - (x(1)-x(0))> and
C<t(m+3) = t(m+2) = x(n-1) + (x(n-1)-x(n-2))>
=item *
2 - Periodic placement of boundary knots:
C<t( 0 ) = t( 1 ) = x(0) - (x(n-1)-x(n-2))> and
C<t(m+3) = t(m+2) = x(n) + (x(1)-x(0))>
=item *
E<lt>0 - Assume the C<nknots> and C<t> were set in a previous call.
=back
C<nknots> is the number of knots and may be changed by the routine.
If C<knotyp E<gt>= 0>, C<nknots> will be set to C<ndim+4>,
otherwise it is an input variable, and an error will occur if its
value is not equal to C<ndim+4>.
C<t> is the array of C<2*n+4> knots for the B-representation
and may be changed by the routine.
If C<knotyp E<gt>= 0>, C<t> will be changed so that the
interior double knots are equal to the x-values and the
boundary knots set as indicated above,
otherwise it is assumed that C<t> was set by a
previous call (no check is made to verify that the data
forms a legitimate knot sequence).
Error status returned by C<$ierr>:
=over 4
=item *
0 if successful.
=item *
-4 if C<knotyp E<gt> 2>.
=item *
-5 if C<knotyp E<lt> 0> and C<nknots != 2*n + 4>.
=back
=cut
',
'Piecewise Cubic Hermite function to B-Spline converter.'
);
##################################################################
##################################################################
#
# This version of polfit accepts bad values and also allows threading
#
#
# indices:
# n runs across input points;
# foo runs across wacky Slatec buffer size;
# bar runs across polynomial coefficients.
#
pp_def('polfit',
Pars => 'x(n); y(n); w(n); int maxdeg(); int [o]ndeg(); [o]eps(); [o]r(n); int [o]ierr(); [o]a(foo); [o]coeffs(bar);[t]xtmp(n);[t]ytmp(n);[t]wtmp(n);[t]rtmp(n)',
OtherPars => '',
Code => '
int maxord;
int ord;
int k;
$GENERIC() zero = 0;
$TFD(polfit'.$uscore.',dpolft'.$uscore.') (&$PRIV(__n_size),$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();
}
$TFD(pcoef'.$uscore.',dpcoef'.$uscore.') ( &(ord), &(zero), $P(coeffs), $P(a));
for(k=ord+1; k<=$maxdeg(); k++)
($P(coeffs))[k] = 0;
',
GenericTypes => ['F','D'],
HandleBad => 1,
NoBadifNaN => 1,
BadCode => 'int ns = $SIZE(n);
int i;
int j = 0;
if($SIZE(n)<$maxdeg()) {
barf("polfit: Must have at least <n> points to fit <n> coefficients");
}
for (i=0;i<ns;i++) { /* get rid of bad values. Call polfit with [xyw]tmp instead of [xyz]. */
if ($ISGOOD(y(n=>i)) && $ISGOOD(x(n=>i)) && $ISGOOD(w(n=>i))) {
$xtmp(n=>j) = $x(n=>i);
$ytmp(n=>j) = $y(n=>i);
$wtmp(n=>j) = $w(n=>i);
j++;
}
}
if (j <= $maxdeg()) {
/* Not enough good datapoints -- set this whole row to BAD. */
for (i=0;i<ns;i++) {
$SETBAD(r(n=>i));
}
$ierr() = 2;
} else {
/* Found enough good datapoints for a fit -- do the fit */
int k;
int ord;
int maxord;
$GENERIC() zero = 0;
/* Do the fit */
$TFD(polfit'.$uscore.',dpolft'.$uscore.')
(&j,$P(xtmp),$P(ytmp),$P(wtmp),$P(maxdeg),$P(ndeg),$P(eps),$P(rtmp),$P(ierr),$P(a));
maxord = ($P(a))[0]+0.5;
ord = ($P(a))[maxord * 3 + 2];
if(ord >= $maxdeg()) {
ord = $maxdeg();
}
/* Extract the polynomial coefficients into coeffs -- used for bad values */
$TFD(pcoef'.$uscore.',dpcoef'.$uscore.') ( &(ord), &(zero), $P(coeffs), $P(a));
for(k=ord+1; k<=$maxdeg(); k++)
($P(coeffs))[k] = 0;
j=0;
for (i=0;i<ns;i++) { /* replace bad values */
if ($ISGOOD(y(n=>i))) {
$r(n=>i) = $rtmp(n=>j);
j++;
} else if($ISGOOD(x(n=>i))) {
/* Bad values are omitted from the call to polfit, so we must reconstitute them on return */
/* (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 the polynomial, it is the highest */
/* order coefficient -- so 3 for a cubic, which has 4 coefficients. */
/* --CED */
int ii;
$GENERIC() acc = 0;
for( ii=ord; ii>0; ii-- ) {
acc += $coeffs(bar=>ii);
acc *= $x(n=>i);
}
acc += $coeffs(bar=>0);
$r(n=>i) = acc;
} else {
/* x and y are bad here... */
$SETBAD(r(n=>i));
}
}
}',
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.
This version handles bad values appropriately',
);
#these two need to be done manually because we don't use defslatec for them
print PROTOS "extern int polfit_ ();\nextern int dpolft_ ();\n";
close( PROTOS );
##################################################################
##################################################################
pp_addpm(<<'EOD');
=head1 AUTHOR
Copyright (C) 1997 Tuomas J. Lukka.
Copyright (C) 2000 Tim Jenness, Doug Burke.
All rights reserved. There is no warranty. You are allowed
to redistribute this software / documentation under certain
conditions. For details, see the file COPYING in the PDL
distribution. If this file is separated from the PDL distribution,
the copyright notice should be included in the file.
=cut
EOD
pp_done();