function P3Grade()
% Checks out Kappa1Est
clc
% Small examples
n = 6;
T0 = triu(randn(n,n)) + 1000*eye(n,n);
for k=1:n
fprintf('\n*****************************\nExample %1d\n\n',k)
T = T0;
T(k,k) = 1;
% In this case, d_opt = I(:,k) because the largest column in inv(T)
% is column k
[yCVL,dCVL] = Kappa1EstCVL(T);
[yYOU,dYOU] = Kappa1Est(T);
format long
disp('yCVL, yYOU = ')
disp([yCVL yYOU])
disp(' ')
disp('dCVL, dYOU = ')
disp([dCVL dYOU])
fprintf('CVL: || Ty-d|| = %10.2e\n',norm(T*yCVL-dCVL,1))
fprintf('YOU: || Ty-d|| = %10.2e\n',norm(T*yYOU-dYOU,1))
fprintf('CVL: ||y||/||inv(T)|| = %10.3f\n',norm(yCVL,1)/norm(inv(T),1))
fprintf('YOU: ||y||/||inv(T)|| = %10.3f\n',norm(yYOU,1)/norm(inv(T),1))
end

function [y,d] = Kappa1EstCVL(T)
% T is an nxn nonsingular upper triangular matrix
% y and d are column n-vectors that satisfy Ty=d.
% d has unit 1-norm and and is chosen so that norm(y,1)
% is close to norm(inv(T),1).

[n,n] = size(T);
d = zeros(n,1);
y = zeros(n,1);
kStar = n;
% At all times, d = I(:,kstar)
y(n) = 1/T(n,n);
norm_y = abs(y(n));
for k=n-1:-1:1
% tau = T(k,k);
% w = T(k,k+1:n)';
% S = T(k+1:n,k+1:n);
% ytilde = y(k+1:n);
% dtilde = d(k+1:n);
%
% tau w' y(k) beta
% =
% 0 S alfa*ytilde alfa*dtilde

% Let lambda = w'*ytilde

% if beta = 1 - alfa, then we want to maximize
% f = alfa*norm(ytilde,1) + abs(1 - alfa*(1+lambda))/abs(tau)
%
% alfa = 0 f = 1/abs(tau)
% alfa = 1 f = norm(ytilde,1) + abs(lambda)/abs(tau)
%
% if lambda>=0 and alfa = 1/(1+lambda), then
% f = norm(ytilde,1)/(1+lambda)
%
% if beta = alfa - 1, then we want to maximize
% f = alfa*norm(ytilde,1) + abs(1 - alfa*(1-lambda))/abs(tau)
%
% alfa = 0 f = 1/abs(tau)
% alfa = 1 f = norm(ytilde,1) + abs(lambda)/abs(tau)
% if lambda<=0 and alfa = 1/(1-lambda), then
% f = norm(ytilde,1)/(1-lambda)
%
lambda = T(k,k+1:n)*y(k+1:n);
if 1/abs(T(k,k)) >= norm_y + abs(lambda)/abs(T(k,k))
% alfa = 0, beta = 1
kStar = k;
y(k) = 1/T(k,k);
y(k+1:n)=0;
norm_y = abs(y(k));
else
% alfa = 1
% beta = 0
y(k) = -lambda/T(k,k);
norm_y = norm_y+abs(y(k));
end
d = zeros(n,1);
d(kStar) = 1;
end