function P3TestScript()
% Checks out BsplineLS
% We will generate some random data based on the function
% f(x) = exp(-x/5)sin(2x) over the interval [a,b]
a = 0;
b = 15;
m = 300;
xVals = [a;a+(b-a)*sort(rand(m-2,1));b];
yVals = exp(-.2*xVals).*sin(2*xVals) + .1*randn(m,1);
plot(xVals,yVals)
% Next we will approximate with b-splines B_{0},...,B_{n+1} where
n = 20;
alphaCVL = BsplineLSCVL(xVals,yVals,n);
BsplineValsCVL = BsplineEvalCVL2(a,b,n,alphaCVL,xVals);
alphaYou = BsplineLS(xVals,yVals,n);
BsplineValsYou = BsplineEvalCVL2(a,b,n,alphaYou,xVals);
hold on
plot(xVals,BsplineValsYou,'r','linewidth',3)
hold off
title(sprintf('|| alphaCVL - alphaYou || = %10.3e',norm(alphaCVL-alphaYou)),'fontsize',14)
shg

function alpha = BsplineLSCVL(zVals,yVals,n)
% zVals and yVals are column m-vectors with a = zVals(1) < ... < zVals(m) = b
% n is a positive integer that satisfies 4<=n<<m. Let h = (b-a)/(n-1).
% alpha is a column (n+2)-vector with the property that if x(k) = a + (k-1)h
% and
% s(z) = alpha(1)B_{0}(z) + ... + alpha_{n+2}B_{n+1}(z)
%
% with B_{k}(z) = Bstar((z-x(k))/h), then
%
% (s(zVals(1)) - yVals(1))^2 + ... + (s(zVals(m)) - yVals(m))^2
% is minimized.
m = length(zVals);
a = zVals(1);
b = zVals(m);
h = (b-a)/(n-1);
% Set up a length-(n+2) struct array where B(k+1).iFirst <= B(k+1).iLast
% are integers with the property that zVals(iFirst:iLast)) is the vector
% of z values in the subinterval [a+(k-3)h,a+(k+1)h] and B(k+1).vals is
% the vector of values obtained when basis spline B_{k} is evaluated at
% zVals(iFirst:iLast)).
for k=0:n+1
xk = a+(k-1)*h;
idx = find(xk-2*h <= zVals & zVals <=xk+2*h);
BkVals = Bstar((zVals(idx)-xk)/h);
B(k+1) = struct('iFirst',idx(1),'iLast',idx(length(idx)),'vals',BkVals);
end
% Set up the matrix M...
M = zeros(n+2,n+2);
for i=0:n+1
for j=max(0,i-3):i
% Compute M(i+1,j+1), the innerproduct associated with B_{i} and B_{j}
% Note that we can abbreviate the j-loop since B_{i}(x)B_{j}(x) is always
% zero whenever |i-j|>3.
iFirst = max(B(i+1).iFirst,B(j+1).iFirst);
iLast = min(B(i+1).iLast,B(j+1).iLast);
range = iFirst:iLast;
% B_{i} and B_{j} are nonzero on zVals(iFirst:iLast);
M(i+1,j+1) = B(i+1).vals(range-B(i+1).iFirst+1)'*B(j+1).vals(range-B(j+1).iFirst+1);
M(j+1,i+1) = M(i+1,j+1);
end
end
% Set up the rhs
rhs = zeros(n+2,1);
for i=0:n+1
rhs(i+1) = B(i+1).vals'*yVals(B(i+1).iFirst:B(i+1).iLast);
end
alpha = M\rhs;

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;