function P1TestScript()
% Checks out MySpline on a small random-data example

% Set up a random example
a = 1;
b = 5;
n = 10;
x = linspace(a,b,n)';
y = randn(n,1);
xVals = linspace(a,b,500)';
clc

% First derivative end conditions
sL = 2;
sR = -3;
alphaCVL = MySplineCVL(a,b,n,y,1,[sL sR]);
yValsCVL = BsplineEvalCVL3(a,b,n,alphaCVL,xVals);
alphaYou = MySpline(a,b,n,y,1,[sL sR]);
yValsYou = BsplineEvalCVL3(a,b,n,alphaYou,xVals);
errFirstDerAlpha = norm(alphaCVL - alphaYou);
errFirstDerVal = norm(yValsCVL - yValsYou);
fprintf('First derivative end conditions...\n\n')
fprintf(' || alphaCVL - alphaYou || = %10.3e\n',errFirstDerAlpha)
fprintf(' || yValsCVL - yValsYou || = %10.3e\n',errFirstDerVal)

% Second derivative end conditions
s2L = 2;
s2R = -3;
alphaCVL = MySplineCVL(a,b,n,y,2,[sL sR]);
yValsCVL = BsplineEvalCVL3(a,b,n,alphaCVL,xVals);
alphaYou = MySpline(a,b,n,y,2,[sL sR]);
yValsYou = BsplineEvalCVL3(a,b,n,alphaYou,xVals);
errSecondDerAlpha = norm(alphaCVL - alphaYou);
errSecondDerVal = norm(yValsCVL - yValsYou);
fprintf('\nSecond derivative end conditions...\n\n')
fprintf(' || alphaCVL - alphaYou || = %10.3e\n',errSecondDerAlpha)
fprintf(' || yValsCVL - yValsYou || = %10.3e\n',errSecondDerVal)

% Not-a-Knot spline
alphaCVL = MySplineCVL(a,b,n,y,3);
yValsCVL = BsplineEvalCVL3(a,b,n,alphaCVL,xVals);
alphaYou = MySpline(a,b,n,y,3);
yValsYou = BsplineEvalCVL3(a,b,n,alphaYou,xVals);
errNotAKnotAlpha = norm(alphaCVL - alphaYou);
errNotAKnotVal = norm(yValsCVL - yValsYou);
fprintf('\nNot-a-Knot end conditions...\n\n')
fprintf(' || alphaCVL - alphaYou || = %10.3e\n',errNotAKnotAlpha)
fprintf(' || yValsCVL - yValsYou || = %10.3e\n',errNotAKnotVal)

% Periodic end conditions
alphaCVL = MySplineCVL(a,b,n,y,4);
yValsCVL = BsplineEvalCVL3(a,b,n,alphaCVL,xVals);
alphaYou = MySpline(a,b,n,y,4);
yValsYou = BsplineEvalCVL3(a,b,n,alphaYou,xVals);
errPeriodicAlpha = norm(alphaCVL - alphaYou);
errPeriodicVal = norm(yValsCVL - yValsYou);
fprintf('\nPeriodic end conditions...\n\n')
fprintf(' || alphaCVL - alphaYou || = %10.3e\n',errPeriodicAlpha)
fprintf(' || yValsCVL - yValsYou || = %10.3e\n',errPeriodicVal)

function alpha = MySplineCVL(a,b,n,y,type,e)
% a and b are scalars with a<b.
% n is a positive integer >=4.
% y is a column n-vector
% alpha is a column (n+2)-vector with the property that if
%
% s(z) = alpha(1)B_{0}(z) + ... + alpha_{n+2}B_{n+1}(z)
%
% s the B-spline representation, then s(x(k)) = y(k) k=1:n.
%
% If type==1 then e must be a 2-vector and s'(x(1)) = e(1) and s'(x(n)) = e(2).
% If type==2 then e must be a 2-vector and s''(x(1)) = e(1) and s''(x(n)) = e(2).
% If type==3 then s'''(z) is continuous at z = x(2) and z = x(n-1).
% If type==4 then s'(x(1)) = s'(x(n)) and s''(x(1)) = s''(x(n)).

x = linspace(a,b,n)';
h = (b-a)/(n-1);
A = zeros(n+2,n+2);
b = zeros(n+2,1);
% Set up the n equations based on s(x(k)) = y(k)), k=1:n ...
for k=2:n+1
b(k) = y(k-1);
A(k,k) = 2/3;
A(k,k-1) = 1/6;
A(k,k+1) = 1/6;
end
b(2:n+1) = y;
% Now deal with the end conditions...
% In the following, B1_{k} first derivative of B_{k}
% In the following, B2_{k} second derivative of B_{k}
% In the following, B3_{k} third derivative of B_{k}
if type==1
% First derivative end conditions...
% At x(1):
% e(1) = alpha(1)B1_{0}(x(1)) + alpha(2)B1_{1}(x(1)) + alpha(3)B1_{2}(x(1))
% = (alpha(1)*(-3) + alpha(2)*0 + alpha(3)*3)/(6h)
A(1,1:3) = [-3 0 3]/(6*h);
b(1) = e(1);
% At x(n):
% e(2) = alpha(n)B1_{n-1}(x(n)) + alpha(n+1)B1_{n}(x(n)) + alpha(n+2)B1_{n+1}(x(n))
% = (alpha(n)*(-3) + alpha(n+1)*0 + alpha(n+2)*3)/(6h)
A(n+2,n:n+2) = [-3 0 3]/(6*h);
b(n+2) = e(2);
elseif type==2
% Second derivative end conditions...
% At x(1):
% e(1) = alpha(1)B2_{0}(x(1)) + alpha(2)B2_{1}(x(1)) + alpha(3)B2_{2}(x(1))
% = (alpha(1) + alpha(2)*(-2) + alpha(3))/h^2
A(1,1:3) = [1 -2 1]/h^2;
b(1) = e(1);
% At x(n):
% e(2) = alpha(n)B2_{n-1}(x(n)) + alpha(n+1)B2_{n}(x(n)) + alpha(n+2)B2_{n+1}(x(n))
% = (alpha(n) + alpha(n+1)*(-2) + alpha(n+2)*3)/h^2
A(n+2,n:n+2) = [1 -2 1]/h^2;
b(n+2) = e(2);
elseif type==3
% Not-A-Knot
% On [x(1),x(2)]
% s'''(z) = alpha(1)B3_{0}(x(2)) + alpha(2)B3_{1}(x(2)) + alpha(3)B3_{2}(x(2)) + alpha(4)B3_{3}(x(2))
% = (alpha(1)*(-1) + alpha(2)*3 + alpha(3)(-3) + alpha(4)*1)/h^3
% On [x(2),x(3)]
% s'''(z) = alpha(2)B3_{1}(x(2)) + alpha(3)B3_{2}(x(2)) + alpha(4)B3_{3}(x(2)) + alpha(5)B3_{4}(x(2))
% = (alpha(2)*(-1) + alpha(3)*3 + alpha(4)(-3) + alpha(5)*1)/h^3
% Equating we get
% -alpha(1) +4*alpha(2) -6*alpha(3) + 4*alpha(4) - alpha(5) = 0
A(1,1:5) = [-1 4 -6 4 -1];
% On [x(n-2),x(n-1)]
% s'''(z) = alpha(n-2)B3_{n-3}(x(n-1)) + alpha(n-1)B3_{n-2}(x(n-1)) + alpha(n)B3_{n-1}(x(n-1)) + alpha(n+1)B3_{n}(x(n-1))
% = (alpha(n-2)*(-1) + alpha(n-1)*3 + alpha(n)(-3) + alpha(n+1)*1)/h^3
% On [x(n-1),x(n)]
% s'''(z) = alpha(n-1)B3_{n-2}(x(n-1)) + alpha(n)B3_{n-1}(x(n-1)) + alpha(n+1)B3_{n}(x(n-1)) + alpha(n+2)B3_{n+1}(x(n-1))
% = (alpha(n-1)*(-1) + alpha(n)*3 + alpha(n+1)(-3) + alpha(n+2)*1)/h^3
% Equating we get
% -alpha(n-2) +4*alpha(n-1) -6*alpha(n) + 4*alpha(n+1) - alpha(n+2) = 0
A(n+2,n-2:n+2) = [-1 4 -6 4 -1];
else
% Periodic Spline
% First derivative continuity:
% s'(x(1)) = alpha(1)B1_{0}(x(1)) + alpha(2)B1_{1}(x(1)) + alpha(3)B1_{2}(x(1))
% = (alpha(1)(-3) + alpha(2)*0 + alpha(3)*3)/h
% s'(x(n)) = alpha(n)B1_{n-1}(x(n)) + alpha(n+1)B1_{n}(x(n)) + alpha(n+2)B1_{n+1}(x(n))
% = (alpha(n)*(-3) + alpha(n+1)*0 + alpha(n+2)*3)/(6h)
% and so
% -3alpha(1) + 3alpha(3) +3alpha(n) - 3alpha(n+2) = 0
A(1,[1 3 n n+2]) = [-3 3 3 -3];
% Second derivative continuity:
% s''(x(1)) = alpha(1)B2_{0}(x(1)) + alpha(2)B2_{1}(x(1)) + alpha(3)B2_{2}(x(1))
% = (alpha(1) + alpha(2)*(-2) + alpha(3))/h^2
% s''(x(n)) = alpha(n)B2_{n-1}(x(n)) + alpha(n+1)B2_{n}(x(n)) + alpha(n+2)B2_{n+1}(x(n))
% = (alpha(n) + alpha(n+1)*(-2) + alpha(n+2))/h^2
% and so
% alpha(1) -2alpha(2) + alpha(3) -alpha(n) + 2alpha(n+1) -alpha(n+2) = 0
A(n+2,[1 2 3 n n+1 n+2]) = [1 -2 1 -1 2 -1];
end
% Solve the linear system
alpha = A\b;

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