PDL::Interpolate::Slatec - simple interface to SLATEC interpolation routines

PDL documentation Contained in the PDL distribution.

Index

Code Index:

NAME

Top

PDL::Interpolate::Slatec - simple interface to SLATEC interpolation routines

SYNOPSIS

Top

 use PDL::Interpolate::Slatec;
 use PDL::Math;

 # somewhat pointless way to estimate cos and sin,
 # but is shows that you can thread if you want to
 #
 my $x   = pdl( 0 .. 45 ) * 4 * 3.14159 / 180;
 my $y   = cat( sin($x), cos($x) );
 #
 my $obj = new PDL::Interpolate::Slatec( x => $x, y = $y );
 #
 my $xi  = pdl( 0.5, 1.5, 2.5 );
 my $yi  = $obj->interpolate( $xi );
 #
 print "cos( $xi ) equals ", $yi->slice(':,(0)'), "\n";
 cos( [0.5 1.5 2.5] ) equals  [0.87759844 0.070737667 -0.80115622]
 #
 print "sin( $xi ) equals ", $yi->slice(':,(1)'), "\n";
 sin( [0.5 1.5 2.5] ) equals  [ 0.4794191 0.99768655 0.59846449]
 #
 print cos($xi), "\n", sin($xi), "\n";
 [0.87758256 0.070737202 -0.80114362]
 [0.47942554 0.99749499 0.59847214]

DESCRIPTION

Top

Use the interface defined by PDL::Interpolate to provide a simple way to use the SLATEC interpolation routines (e.g. see PDL::Slatec (PDL::Slatec)). Hence the name for this library - as returned by the library method - is "Slatec".

Currently, only the piecewise cubic Hermite polynomial routines (Piecewise_cubic_Hermite_interpol in PDL::Slatec) are available (type == "pch").

Attributes

The following changes are made to the attributes of PDL::Interpolate:

 Attribute  Flag  Description
 bc         sgr   boundary conditions
 g          g     estimated gradient at x positions

 Attribute  Default    Allowed values
 bc         "simple"   see Boundary conditions section
 type       "pch"

Given the initial set of points (x,y), the "pch" library estimates the gradient at these points using the given boundary conditions (as specified by the bc attribute). The estimated gradient can be obtained using

 $gradient = $obj->get( 'g' );

As described in the interpolate method, the "pch" routines can also estimate the gradient, as well as the function value, for a set of $xi.

Boundary conditions for the pch routines

If your data is monotonic, and you are not too bothered about edge effects, then the default value of bc of "simple" is for you. Otherwise, take a look at the description of PDL::Slatec::chic (chic in PDL::Slatec) and use a hash reference for the bc attribute, with the following keys:

monotonic

0 if the interpolant is to be monotonic in each interval (so the gradient will be 0 at each switch point), otherwise the gradient is calculated using a 3-point difference formula at switch points. If > 0 then the interpolant is forced to lie close to the data, if < 0 no such control is imposed. Default = 0.

start

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 -5 .. 5, as given for the ic parameter of chic (chic in PDL::Slatec). The second element, only used if options 2, 1, -1, or 2 are chosen, contains the value of the vc parameter. Default = [ 0 ].

end

As for start, but for the end of the data.

An example would be

 $obj->set( bc => { start => [ 1, 0 ], end => [ 1, -1 ] }

which sets the first derivative at the first point to 0, and at the last point to -1.

Errors

The status method provides a simple mechanism to check if the previous method was successful. The err attribute contains the $ierr piddle returned by the Slatec routine if a more precise diagnostic is required. To find out which routine was called, use the routine method.

interpolate

 my $yi          = $obj->interpolate( $xi );
 my ( $yi, $gi ) = $obj->interpolate( $xi );

Returns the interpolated function and derivative at a given set of points.

If evaluated in scalar mode, it returns only the interpolated function values.

integrate

 my $ans = $obj->integrate( index => pdl( 2, 5 ) );
 my $ans = $obj->integrate( x => pdl( 2.3, 4.5 ) );

Integrate the function stored in the PDL::Interpolate::Slatec object.

The integration can either be between points of the original x array (index), or arbitrary x values (x). For both cases, a two element piddle should be given, to specify the start and end points of the integration.

index

The values given refer to the indices of the points in the x array.

x

The array contains the actual values to integrate between.

If the status method returns a value of -1, then one or both of the integration limits did not lie inside the x array. Caveat emptor with the result in such a case.

The reason for using piddles, rather than arrays, is that it allows for threading.

AUTHOR

Top

Copyright (C) 2000 Doug Burke (burke@ifa.hawaii.edu). 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.

SEE ALSO

Top

PDL::Interpolate, PDL, perltoot(1).

PDL documentation Contained in the PDL distribution. 
package PDL::Interpolate::Slatec;

use PDL::Interpolate;
use PDL::Slatec;
use Carp;

use strict;
use vars qw( @ISA );

@ISA     = qw ( PDL::Interpolate );

####################################################################
#

####################################################################
# 
## Public routines:

sub new {
    my $this  = shift;
    my $class = ref($this) || $this;
    my $self  = $class->SUPER::new();  # note: do not send in values

    # change from PDL::Interpolate to PDL::Interpolate::Slatec
    bless ($self, $class);

    # change class attributes
    $self->_change_attr( 
			 bc => { required => 1, settable => 1 }, # already gettable
			 );
    $self->_set_value( bc => "simple", type => "pch" );
    $self->_add_attr( 
		      g => { gettable => 1 },
		     );

    $self->{flags}->{library} = "Slatec";
    $self->{flags}->{routine} = "none";

    # set variables
    $self->set( @_ );

    return $self;         
                                               
} # sub: new

####################################################################

# set up the interpolation
#
sub _initialise {
    my $self = shift;

    # set up error flags
    $self->{flags}->{status} = 0;
    $self->{flags}->{routine} = "none";

    # get values in one go
    my ( $x, $y, $g, $bc ) = $self->_get_value( qw( x y g bc ) );

    # check 1st dimention of x and y are the same
    #  ie allow the possibility of threading
    my $xdim = $x->getdim( 0 );
    my $ydim = $y->getdim( 0 );
    croak "ERROR: x and y piddles must have the same first dimension.\n"
	unless $xdim == $ydim;

    # if a gradient has been specified, then we don't need to do anything
    # - other than check the dimensions
    if ( defined $g ) {
	croak "ERROR: gradient piddle must have the same first dimension as x and y.\n"
	    unless $g->getdim( 0 ) == $xdim;
	$self->{flags}->{status} = 1;
	return;
    }

    my $ierr;
    if ( ref($bc) eq "HASH" ) {
	my $monotonic = $bc->{monotonic} || 0;
	my $start     = $bc->{start}     || [ 0 ];
	my $end       = $bc->{end}       || [ 0 ];

	my $ic = $x->short( $start->[0], $end->[0] );
	my $vc = $x->float( 0, 0 ); # it will get promoted if required

	if ( $#$start == 1 ) { $vc->set( 0, $start->[1] ); }
	if ( $#$end   == 1 ) { $vc->set( 1, $end->[1] ); }

	my $wk = $x->zeroes( $x->float, 2*$xdim );
	( $g, $ierr ) = chic( $ic, $vc, $monotonic, $x, $y, $wk );

	$self->{flags}->{routine} = "chic";

    } elsif ( $bc eq "simple" ) {
	# chim
	( $g, $ierr ) = chim( $x, $y );

	$self->{flags}->{routine} = "chim";

    } else {
	# Unknown boundary condition
	croak "ERROR: unknown boundary condition <$bc>.\n";
	# return;
    }

    $self->_set_value( g => $g, err => $ierr );

    if ( all $ierr == 0 ) {
	# everything okay
	$self->{flags}->{status} = 1;
    } elsif ( any $ierr < 0 ) {
	# a problem
	$self->{flags}->{status} = 0;
    } else {
	# there were switches in monotonicity
	$self->{flags}->{status} = -1;
    }

} # sub: _initialise

####################################################################


sub interpolate {
    my $self = shift;
    my $xi   = shift;

    croak 'Usage: $obj->interpolate( $xi )' . "\n"
	unless defined $xi;

    # check everything is fine
    $self->_check_attr();

    # get values in one go
    my ( $x, $y, $g ) = $self->_get_value( qw( x y g ) );

    my ( $yi, $gi, $ierr );

    if ( wantarray ) {
	( $yi, $gi, $ierr ) = chfd( $x, $y, $g, 0, $xi );
	$self->{flags}->{routine} = "chfd";
    } else {
	( $yi, $ierr ) = chfe( $x, $y, $g, 0, $xi );
	$self->{flags}->{routine} = "chfe";
    }

    # set err/status info
    $self->_set_value( err => $ierr );

    if ( all $ierr == 0 ) {
	# everything okay
	$self->{flags}->{status} = 1;
    } elsif ( all $ierr > 0 ) {
	# extrapolation was required
	$self->{flags}->{status} = -1;
    } else {
	# a problem
	$self->{flags}->{status} = 0;
    }
	
    return wantarray ? ( $yi, $gi ) : $yi;

} # sub: interpolate


sub integrate {
    my $self = shift;

    croak 'Usage: $obj->integrate( $type => $limits )' . "\n"
	unless $#_ == 1;

    # check everything is fine
    $self->_check_attr();

    $self->{flags}->{status} = 0;
    $self->{flags}->{routine} = "none";

    my ( $type, $indices ) = ( @_ );

    croak "Unknown type ($type) sent to integrate method.\n"
	unless $type eq "x" or $type eq "index";

    my $fdim = $indices->getdim(0);
    croak "Indices must have a first dimension of 2, not $fdim.\n"
	unless $fdim == 2;

    my $lo = $indices->slice('(0)');
    my $hi = $indices->slice('(1)');

    my ( $x, $y, $g ) = $self->_get_value( qw( x y g ) );
    my ( $ans, $ierr );

    if ( $type eq "x" ) {
	( $ans, $ierr ) = chia( $x, $y, $g, 0, $lo, $hi );
	$self->{flags}->{routine} = "chia";

	if ( all $ierr == 0 ) {
	    # everything okay
	    $self->{flags}->{status} = 1;
	} elsif ( any $ierr < 0 ) {
	    # a problem
	    $self->{flags}->{status} = 0;
	} else {
	    # out of range
	    $self->{flags}->{status} = -1;
	}

    } else {
	( $ans, $ierr ) = chid( $x, $y, $g, 0, $lo, $hi );
	$self->{flags}->{routine} = "chid";

	if ( all $ierr == 0 ) {
	    # everything okay
	    $self->{flags}->{status} = 1;
	} elsif ( all $ierr != -4 ) {
	    # a problem
	    $self->{flags}->{status} = 0;
	} else {
	    # out of range (ierr == -4)
	    $self->{flags}->{status} = -1;
	}

    }

    $self->_set_value( err => $ierr );
    return $ans;
}


####################################################################
# End with a true
1;