close all; clear; clc;

%% Make random A
N = 50; % Size of A
A = rand(N,N); A = A*A'; % Random positive semi-definite matrix

%% Eigenvalues with eig
tic
[V,D] = eig(A);
toc
D_eig = diag(D)';

%% Eigenvalues by optimisation
plt = false; % Allow disabling plot for timing
options = optimset('Display','none');%,'MaxIter',1e4,'MaxFunEvals',1e4);
cm = plasma(N);
tic
for j = 1:size(A,1)
if j > 1
    lossfun = @(L) lossfun(L) + log10( (L-d(j-1)).^-2 );
else
    lossfun = @(L) log10( det(A - (L)*eye(size(A,1))).^2 );
end
d(j) = fminsearch(lossfun,0,options);
%
if plt
    Ls = linspace(0,20,1e3+1);
    for v1i = 1:numel(Ls)
        L_surf(v1i) = lossfun(Ls(v1i));
    end
    LL(j,:) = L_surf;
    if j == 1
        figure(1); clf;
        grid on; box on; hold on; axis square;
        xlabel('L'); ylabel('loss function');
    end
    plot(Ls,LL(j,:),'color',cm(j,:)); drawnow;
    end
end
toc
fprintf('Optimisation is *not* faster than ''eig'', but it does work\n');
D_optim = sort(d);

%%
delta_D = abs(D_eig-D_optim);
fprintf('median delta = %.1f eps\n',median(delta_D./eps));
fprintf('max delta = %.1f eps\n',max(delta_D./eps));