function P2_TestScript()
% Generate a large random s(z) = alpha(1)B_{0} + ... + alpha(n+2)B_{n+1}..
n = 400;
alpha = randn(n+2,1);
% Scatter m z-values across the interval [a,b] where
m = 4000; a = 1; b = 5;
% Include the breakpoints..
zVals = a+ (b-a)*rand(m-n,1);
zVals = sort([zVals;linspace(a,b,n)']);
% Compare our three methods with yours...
tic
yVals1 = BsplineEvalCVL1(a,b,n,alpha,zVals);
t1 = toc;
tic
yVals2 = BsplineEvalCVL2(a,b,n,alpha,zVals);
t2 = toc;
tic
yVals3 = BsplineEvalCVL3(a,b,n,alpha,zVals);
t3 = toc;
tic
yValsYou = BsplineEval(a,b,n,alpha,zVals);
tYou = toc;
clc
fprintf('n = %1d\nm = %1d\n\n',n,m)
fprintf('Time CVL1 = %6.3f\n',t1)
fprintf('Time CVL2 = %6.3f\n',t2)
fprintf('Time CVL3 = %6.3f\n',t3)
fprintf('Time You = %6.3f\n',tYou)
fprintf('\n max| yValsCVL(k) - yValsYou(k) | = %10.3e\n',max(abs(yVals1-yValsYou)))

function yVals = BsplineEvalCVL1(a,b,n,alpha,zVals)
% a and b are scalars with a<b.
% n is a positive integer >=4.
% alpha is a column (n+2)-vector,
% zVals is a column m-vector with a <= zVals(1) <= ... <= zVals(m) <= b
% yVals is a column m-vector with the property that if
%
% s(z) = alpha(1)B_{0}(z) + ... + alpha_{n+2}B_{n+1}(z)
%
% then yVals(k) = s(zVals(k)), k=1:m

m = length(zVals);
yVals = zeros(m,1);
% The breakpoints..
h = (b-a)/(n-1);
for k=1:m
for j=0:n+1
xj = a + (j-1)*h;
yVals(k) = yVals(k) + alpha(j+1)*Bstar((zVals(k)-xj)/h);
end
end

function yVals = BsplineEvalCVL2(a,b,n,alpha,zVals)
% a and b are scalars with a<b.
% n is a positive integer >=4.
% alpha is a column (n+2)-vector,
% zVals is a column m-vector with a <= zVals(1) <= ... <= zVals(m) <= b
% yVals is a column m-vector with the property that if
%
% s(z) = alpha(1)B_{0}(z) + ... + alpha_{n+2}B_{n+1}(z)
%
% then yVals(k) = s(zVals(k)), k=1:m

m = length(zVals);
yVals = zeros(m,1);
% The breakpoints..
h = (b-a)/(n-1);
for k=0:n+1
% add in the contribution of B_{k}
xk = a+(k-1)*h;
idx = find(xk-2*h <= zVals & zVals<=xk+2*h);
yVals(idx) = yVals(idx) + alpha(k+1)*MyBstar((zVals(idx)-xk)/h);
end

function y = MyBstar(z)
% z is a column m-vector and y is a column m-vector with the property that
% y(k) = B_{*}(x(k)), k=1:m
m = length(z);
y = zeros(m,1);

i1 = -2<=z & z<=-1;

% z-values between -2 and -1 ...
idx = find(-2<=z & z <=-1);
y(idx) = ((2+z(idx)).^3)/6;

% z-values between -1 and 0 ...
idx = find(-1<=z & z <=0);
t = (1+z(idx));
y(idx) = ((-.5*t+.5).*t+.5).*t+(1/6);

% z-values between 0 and 1 ...
idx = find(0<=z & z <=1);
t = (1-z(idx));
y(idx) = ((-.5*t+.5).*t+.5).*t+(1/6);

% z-values between 1 and 2 ...
idx = find(1<=z & z <=2);
y(idx) = ((2-z(idx)).^3)/6;

function yVals = BsplineEvalCVL3(a,b,n,alpha,zVals)
% a and b are scalars with a<b.
% n is a positive integer >=4.
% alpha is a column (n+2)-vector,
% zVals is a column m-vector with a <= zVals(1) <= ... <= zVals(m) <= b
% yVals is a column m-vector with the property that if
%
% s(z) = alpha(1)B_{0}(z) + ... + alpha_{n+2}B_{n+1}(z)
%
% then yVals(k) = s(zVals(k)), k=1:m

m = length(zVals);
yVals = zeros(m,1);
% The breakpoint definitions.
% Let x(k) = a + (k-1)h and refer to [x(k),x(k+1)] as the kth interval
% where
h = (b-a)/(n-1);
% If
idx = min(floor((zVals-a)/h)+1,n-1);
% then zVals(k) is in interval idx(k)
% We say [x(k),x(k+1)] is a "serious interval if it contains at least
% one z-value. Here is the vector of serious interval indices..
SeriousIntervals = [1;1+find(diff(idx)>0)];
for i=1:length(SeriousIntervals)
% The ith serious interval is interval k where
k = idx(SeriousIntervals(i));
% zVals(iL:iR) houses the z-values that are in interval k where
iL = SeriousIntervals(i);
if i<length(SeriousIntervals)
iR = SeriousIntervals(i+1)-1;
else
iR = m;
end
zk = zVals(iL:iR);
% The basis function that are nonzero on the kth interval
% are B_{k-1}, B_{k}, B_{k+1}, and B_{k+2}
wk = (zk-a)/h;
tau = k-wk;
% Add in contribution from B_{k-1}...
yVals(iL:iR) = yVals(iL:iR) + (alpha(k)/6)*tau.^3;
% Add in contribution from B_{k}...
yVals(iL:iR) = yVals(iL:iR) + (alpha(k+1)/2)*(( (-tau + 1).*tau + 1).*tau + 1/3);
% Add in contribution from B_{k+1}...
tau = 1-tau;
yVals(iL:iR) = yVals(iL:iR) + (alpha(k+2)/2)*(( (-tau + 1).*tau + 1).*tau + 1/3);
% Add in contribution from B_{k+2}...
yVals(iL:iR) = yVals(iL:iR) + (alpha(k+3)/6)*tau.^3;
end

function y = Bstar(z)
% z is a column m-vector and y is a column m-vector with the property that
% y(k) = B_{*}(x(k)), k=1:m
m = length(z);
y = zeros(m,1);
for k=1:m
t = z(k);
if -2<=t && t<=-1
y(k) = (2+t)^3/6;
elseif -1<=t && t<=0
y(k) = (-3*(1+t)^3 + 3*(1+t)^2 + 3*(1+t) + 1)/6;
elseif 0<=t && t<=1
y(k) = (-3*(1-t)^3 + 3*(1-t)^2 + 3*(1-t) + 1)/6;
elseif 1<=t && t<=2
y(k) = (2-t)^3/6;
end
end
% This can be made much more efficient by using the Matlab function
% find and vectorizing.