function P4Grade()
% Checks out the quality of Kappa1Est as a 1-norm condition estimator.
% Reports on the ratio of Kappa1Est(T)*norm(T,1)/condest(T)

mode = 3;
nVals = [100 200 400 800]; % Results essential the same for larger n
kappaVals = [1000 1000000 1000000000] % Results not sensitive to kappa
clc
fprintf('Set mode = %1d in randsvd\n',mode)
fprintf('"Ratio" means Kappa1Est(T)*norm(T,1)/condest(T)\n\n')
disp(' n kappa2(T) minRatio maxRatio aveRatio')
disp('------------------------------------------------------------')
for n = nVals
for kappa = kappaVals
% Initializations for the experiments for this particular n and
% kappa...
nRepeat = 3200/n;
sumRatio = 0;
maxRatio = 0;
minRatio = inf;
for k=1:nRepeat
T = gallery('randsvd', n, kappa, mode, 0, n-1);
% Estimate mu = cond(T,1) two ways...
mu_condest = condest(T);
[y,d] = Kappa1EstCVL(T);
mu_Kappa1Est = norm(y,1)*norm(T,1);
ratio = mu_Kappa1Est/mu_condest;
minRatio = min(ratio,minRatio);
maxRatio = max(ratio,maxRatio);
sumRatio = sumRatio+ratio;
end
fprintf('%4d %10.2e %11.2f %11.2f %11.2f\n',n,kappa,minRatio,maxRatio,sumRatio/nRepeat)
end
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

% Not a precise formula for grading, but here are some guidelines
% -1 for no repetitionof experiments for a particular n and kappa.
% Kappaa1est is a heuristic method and you can't draw conclusions from one
% problem.

% Not fussy about your range of n values and kappa values, but there had to
% be a range.

% -1 for terrible displays of results. A table has columns and headings.
% Not hard to do. Lines of output like
% n = 50 kappa = 1000 my_ratio = .1234 matlabratio = .6789
% are not presentable. Why not a table with headings n, kappa, etc?

% Given how condition estimation is used and how important it is not to
% underestimate the value, you should do more than just print averages.
%
%