#======================================================================

#                    . I N T E R P . P O L Y 

#                    doc: Thu Nov 23 21:30:25 2000

#                    dlm: Tue Aug  5 14:20:19 2008

#                    (c) 2000 A.M. Thurnherr

#                    uE-Info: 58 0 NIL 0 0 72 10 2 4 NIL ofnI

#======================================================================

# polynomial interpolation

# HISTORY:

# 	Nov 23, 2000: - created

#	Nov 24, 2000: - added -c)enter x-val in window

#	Dec 13, 2002: - replaced -n by -o

#	Oct 22, 2003: - changed -e to -r (conflicts with -e option of [resample])

#	Jul 28, 2006: - removed debugging code

#				  - added xf to ISInit() args

#				  - BUG: sufficient-input data test was buggy

#   Aug 22, 2006: - adapted to work with [match]

#	Sep 11, 2006: - BUG: sufficient-input ermesg was wrong

#   Aug  5, 2008: - added idr param to IS_init()

# see [.interp.linear] for documentation of interface

# NOTES:

#	- $opt_o <poly order> is required

#	- $opt_r is used to return the error estimates instead of the

#	  interpolated values

#	- -o 0 is almost like using -s subsample, the only difference

#	  being that subsample returns the real x values while -o 0

#	  returns the subsampled x values

#	- on -c the window for each interpolation estimate is centered

#	  around the interpolated x value; this implies window-sliding

#	  between tabulated x values resulting in functional disconti-

#	  nuities

require	"$ANTS/polint.pl";

$IS_opts = "co:r";

$IS_optsUsage = "[e-r)ror estimates] [-c)enter x-val] -o <poly order>";

sub IS_usage()

{

	die("$0 (.interp.poly): ERROR! -o required\n")

		unless (defined($opt_o) && $opt_o >= 0);

	&antsInfo("WARNING: -c generates discontinuous function")

		if ($opt_c);

	&antsInfo("WARNING: high-order polynomial interpolation deprecated")

		if ($opt_o > 5);

	&antsInfo("order-1 polynomial; -s linear is more efficient")

		if ($opt_o == 1);

	&antsInfo("giving error estimates") if ($opt_r);

}

sub IS_init($$$$)

{

	my($bR,$idR,$f,$xf) = @_;

	die(sprintf("$0 (.interp.poly): ERROR! " .

				"not enough data points (%d) for order-$opt_o poly\n",

				 $#{$bR}+1))

		if ($opt_o > $#{$bR});

}

sub findWindow($$$$$)							# find window of size n=opt_o+1

{

	my($bR,$f,$v,$i,$n) = @_;

	my($cn);									# current win size

	return ($bR->[$i+1][$f] - $v < $v - $bR->[$i][$f]) ? $i+1 : $i

		if ($n == 1);							# nearest neighbor

	if ($opt_c) {								# center window around x value

		for ($cn=2; $cn<$n; $cn++) {			# grow window

			$i-- if (($i>0 && $bR->[$i+$cn-1][$f] - $v >= $v - $bR->[$i][$f]) ||

					 ($i+$cn-1 == $#{$bR}));

		}

	} else {

		$i -= int(($n-1)/2);						# ideally

		$i  = 0 if ($i < 0);						# off leftmost column

		$i  = $#{$bR}+1 - $n if ($i+$n > $#{$bR}+1);# off rightmost column

	}

	return $i;

}

sub IS_interpolate($$$$$$)

{

	my($bR,$idR,$xf,$xv,$xi,$f) = @_;

	return $xv if ($xf == $f);

	$xi = &findWindow($bR,$xf,$xv,$xi,$opt_o+1);

	my($y,$dy) = &polint($bR,$xf,$xv,$xi,$opt_o+1,$f);

	return $opt_r ? $dy : $y;

}