% This tests my function interp1_quick for accuracy and speed.

% $Id:$
% Copyright J. V. Mansbridge, CSIRO, Thu May  7 11:48:53 EST 1998

n_rep = 125;
x = repmat([(0:.1:1)' 2*(-.1:.1:.9)' (20:30)' (-30:-20)'], 1, n_rep);
% x = flipud(repmat([(0:.1:1)' 2*(-.1:.1:.9)' (20:30)' (-30:-20)'], 1, n_rep));
y = repmat([(1:2:21)' (-3:2:17)' (20:30)' (-30:-20)'], 1, n_rep);

% xi as a scalar

xi = 0.22

tic
n_cols = size(x, 2);
yii = zeros(1, n_cols);
for ii = 1:n_cols
  yii(ii) = interp1(x(:, ii), y(:, ii), xi);
end
toc


tic
yi = interp1_quick(x, y, xi);
toc

range(yi - yii)
nani = isnan(yi);
nanii = isnan(yii);
range(nanii - nani)
num_nani = sum(nani)

% xi as a column vector

xi = rand(1, size(x, 2));

tic
n_cols = size(x, 2);
yii = zeros(1, n_cols);
for ii = 1:n_cols
  yii(ii) = interp1(x(:, ii), y(:, ii), xi(ii));
end
toc


tic
yi = interp1_quick(x, y, xi);
toc

range(yi - yii)
nani = isnan(yi);
nanii = isnan(yii);
range(nanii - nani)
num_nani = sum(nani)

% xi as a row vector

xi = rand(size(x, 2), 1);

tic
n_cols = size(x, 2);
yii = zeros(1, n_cols);
for ii = 1:n_cols
  yii(ii) = interp1(x(:, ii), y(:, ii), xi(ii));
end
toc


tic
yi = interp1_quick(x, y, xi);
toc

range(yi - yii)
nani = isnan(yi);
nanii = isnan(yii);
range(nanii - nani)
num_nani = sum(nani)
extra1 = y;
extra2 = y;

tic
[yi, extra1_i, extra2_i] = interp1_quick(x, y, xi, extra1, extra2);
time_2 = toc

range(yi - yii)
nani = isnan(yi);
nanii = isnan(yii);
range(nanii - nani)
num_nani = sum(nani)
range(yi - extra1_i)
range(yi - extra2_i)

% Check for non-monotonicity

x_bad = x;
x_bad(4,1) = 7;
tic
yi = interp1_quick(x_bad, y, xi);
toc
tic
[yi, extra1_i, extra2_i] = interp1_quick(x_bad, y, xi, extra1, extra2)
time_3 = toc

