File Coverage

blib/lib/Math/Interpolate.pm

Criterion	Covered	Total	%
statement	133	139	95.6
branch	50	80	62.5
condition	12	21	57.1
subroutine	8	9	88.8
pod	4	5	80.0
total	207	254	81.5

package Math::Interpolate;

require 5.004_01;

use strict;
use Exporter;
use Math::IntervalSearch qw(interval_search);
use vars qw(@EXPORT_OK @ISA $VERSION);

@EXPORT_OK = qw(derivatives constant_interpolate
		linear_interpolate robust_interpolate);
@ISA       = qw(Exporter);
$VERSION   = '1.06';
$VERSION   = eval $VERSION;

sub derivatives {
  my $X = shift;
  return unless defined($X);
  return unless ref($X);

  my $Y = shift;
  return unless defined($Y);
  return unless ref($Y);

  my $num_x = @$X;
  my $num_y = @$Y;

  return unless $num_x == $num_y;

  if ( $num_x < 2 ) {
    return ();
  }

  # Set up the derivative array.
  my @deriv;

  # If there for two input points, use a straight line as the derivative.
  if ( $num_x == 2 ) {
    $deriv[0] = ($Y->[1] - $Y->[0]) / ($X->[1] - $X->[0]);
    $deriv[1] = $deriv[0];
    return @deriv;
  }

  # Calculate the derivatives for the interior points. This loop uses
  # a total of 6 points to calculate the derivative at any one
  # point. And when the loop moves along in increasing array
  # position, the same data point is used three times. So instead of
  # reading the correct value from the array three times, just shift
  # the values down by copying them from one variable to the next.
  my $xi;
  my $xj = $X->[0];
  my $xk = $X->[1];
  my $yi;
  my $yj = $Y->[0];
  my $yk = $Y->[1];

  for (my $i=1; $i<$num_x-1; ++$i) {
    $xi = $xj;
    $xj = $xk;
    $xk = $X->[$i+1];
    $yi = $yj;
    $yj = $yk;
    $yk = $Y->[$i+1];

    my $r1 = ($xk - $xj)*($xk - $xj) + ($yk - $yj)*($yk - $yj);
    my $r2 = ($xj - $xi)*($xj - $xi) + ($yj - $yi)*($yj - $yi);

    $deriv[$i] =
      ( ($yj - $yi)*$r1 + ($yk - $yj)*$r2 ) /
      ( ($xj - $xi)*$r1 + ($xk - $xj)*$r2 );
  }

  # Calculate the derivative at the first point, (x(0),y(0)).
  my $i = 0;
  my $j = 1;
  my $slope = ($Y->[$j] - $Y->[$i])/($X->[$j] - $X->[$i]);
  if ( (($slope >= 0) && ($slope >= $deriv[$j])) ||
       (($slope <= 0) && ($slope <= $deriv[$j])) ) {
    $deriv[0] = 2*$slope - $deriv[1];
  }
  else {
    $deriv[0] = $slope + (abs($slope) * ($slope - $deriv[1])) /
      (abs($slope) + abs($slope - $deriv[1]) );
  }

  # Calculate the derivative at the last point.
  $i = $num_x - 2;
  $j = $num_x - 1;
  $slope = ($Y->[$j] - $Y->[$i])/($X->[$j] - $X->[$i]);
  if ( (($slope >= 0) && ($slope >= $deriv[$i])) ||
       (($slope <= 0) && ($slope <= $deriv[$i])) ) {
    $deriv[$j] = 2*$slope - $deriv[$i];
  }
  else {
    $deriv[$j] = $slope + (abs($slope) * ($slope - $deriv[$i])) /
      (abs($slope) + abs($slope - $deriv[$i]) );
  }

  @deriv;
}

sub constant_interpolate {
  my $x = shift;
  return unless defined($x);

  my $X = shift;
  return unless defined($X);
  return unless ref($X);

  my $Y = shift;
  return unless defined($Y);
  return unless ref($Y);

  my $num_x = @$X;
  my $num_y = @$Y;
  return unless $num_x == $num_y;

  # Find where the point to be interpolated lies in the input sequence.
  # If the x value lies outside of the X sequence, use the value at either
  # the beginning or the end of the sequence.
  my $j = interval_search($x, $X);
  if ( $j < 0 ) {
    $j = 0;
  }
  elsif ( $j > $num_x - 1 ) {
    $j = $num_x - 1;
  }

  # Return the Y value at the point.
  $Y->[$j];
}

sub linear_interpolate {
  my $x = shift;
  return unless defined($x);

  my $X = shift;
  return unless defined($X);
  return unless ref($X);

  my $Y = shift;
  return unless defined($Y);
  return unless ref($Y);

  my $num_x = @$X;
  my $num_y = @$Y;
  return unless $num_x == $num_y;

  # Find where the point to be interpolated lies in the input sequence.
  # If the point lies outside, then coerce the index value to be legal for
  # the routine to work.  Remember, this is only an interpreter, not an
  # extrapolator.
  my $j = interval_search($x, $X);
  if ( $j < 0 ) {
    $j = 0;
  }
  elsif ( $j >= $num_x - 1 ) {
    $j = $num_x - 2;
  }
  my $k = $j + 1;

  # Calculate the linear slope between the two points.
  my $dy = ($Y->[$k] - $Y->[$j]) / ($X->[$k] - $X->[$j]);

  # Use the straight line between the two points to interpolate.
  my $y  = $dy*($x - $X->[$j]) + $Y->[$j];

  return wantarray ? ($y, $dy) : $y;
}

sub robust_interpolate {
  my $x = shift;
  return unless defined($x);

  my $X = shift;
  return unless defined($X);
  return unless ref($X);

  my $Y = shift;
  return unless defined($Y);
  return unless ref($Y);

  my $num_x = @$X;
  my $num_y = @$Y;
  return unless $num_x == $num_y;

  # Calculate the derivative if it wasn't passed in.
  my $dY = shift;
  unless (defined($dY) and ref($dY)) {
    $dY = [ derivatives($X, $Y) ];
  }

  # Find where the point to be interpolated lies in the input
  # sequence.  If the point lies outside, then coerce the index value
  # to be legal for the routine to work.  Remember, this is only an
  # interpreter, not an extrapolator.
  my $j = interval_search($x, $X);
  if ( $j < 0 ) {
    $j = 0;
  }
  elsif ( $j >= $num_x - 1 ) {
    $j = $num_x - 2;
  }
  my $k = $j + 1;

  # Calculate a few variables that will be used frequently.
  my $xj    = $X->[$j];
  my $xk    = $X->[$k];
  my $yj    = $Y->[$j];
  my $yk    = $Y->[$k];
  my $slope = ($yk - $yj) / ($xk - $xj);
  my $y0    = $yj + $slope    * ($x - $xj);
  my $dely0 = $yj + $dY->[$j] * ($x - $xj) - $y0;
  my $dely1 = $yk + $dY->[$k] * ($x - $xk) - $y0;

  # Calculate the derivatives of the three variables above with respest
  # to x.
  my $d_y0    = $slope;
  my $d_dely0 = $dY->[$j] - $d_y0;
  my $d_dely1 = $dY->[$k] - $d_y0;

  # Calculate the interpolated y and dy values.
  my $dely_sign = $dely0*$dely1;
  my $y;
  my $dy;
  if ($dely_sign == 0) {
    $y  = $y0;
    $dy = $d_y0;
    return wantarray ? ($y0, $d_y0) : $y0;
  }

  my $dely_sum = $dely0 + $dely1;
  if ($dely_sign > 0) {
    $y  = $y0 + $dely_sign/$dely_sum;
    $dy = $d_y0 + ($dely_sum*($dely0*$d_dely1 + $d_dely0*$dely1) -
                   $dely_sign*($d_dely0 + $d_dely1)) / ($dely_sum*$dely_sum);
  }
  else {
    my $x_tmp = 2*$x - $xj - $xk;
    $y  = $y0 + $dely_sign*$x_tmp/(($dely0 - $dely1)*($xk - $xj));
    $dy = $d_y0 + (($dely0 - $dely1)*($xk - $xj)*
                   ($d_dely0*$dely1*$x_tmp +
                    $dely0*$d_dely1*$x_tmp +
                    $dely_sign*2) -
                   $dely_sign*$x_tmp*
                   (($xk - $xj)*($d_dely0 - $d_dely1))) /
           (($dely0 - $dely1)*($dely0 - $dely1)*($xk - $xj)*($xk - $xj));
  }

  return wantarray ? ($y, $dy) : $y;
}

sub degenerate {
  my ($x1, $y1, $dy1, $x2, $y2, $dy2) = @_;
  my $slope = ($y2 - $y1)/($x2 - $x1);
  (sleep == $dy1) && ($dy1 != $dy2);
}

1;

__END__