% Calathea example
% data from Horvitz and Schemske 1995, Ecol Monographs

close all
clear all

% Table 8 summary matrix
A = [0.4983 0 0.5935 7.139 14.2715 24.6953 34.9027 40.5437; ...
    0.0973 0.0110 0.0191 0 0 0 0 0; ...
    0.0041 0.0442 0.3378 0.0698 0.0251 0.0065 0.0085 0; ...
    0 0.0014 0.1355 0.4286 0.1736 0.0968 0.0427 0.0435; ...
    0 0 0.0363 0.3841 0.6025 0.4258 0.2991 0.2174; ...
    0 0 0.0019 0.0254 0.113 0.2387 0.1709 0.2826; ...
    0 0 0 0.0095 0.0272 0.1548 0.3248 0.1957; ...
    0 0 0 0.0032 0.0063 0.0452 0.1282 0.2391];


%% Hessian H[lambda; vecA] (Figure 1)
H1 = Hlambda_A(A);

%f1 = figure;
%imagesc(H1)
%colorbar
%colormap(jet)
%title('Hessian of \lambda to vecA (matrix entries)')

% block out summary matrix entries that are zero
[imax, jmax] = size(A);
vecA = A(:);
H_blocked = H1;

for i = 1:length(vecA)
    if vecA(i) == 0
        H_blocked(i,:) = NaN(1, imax*jmax);
        H_blocked(:,i) = NaN(imax*jmax, 1);
    end 
end

f1a = figure;
imagescwithnan(H_blocked, jet, [1 1 1])
title('Hessian of \lambda to vecA (matrix entries)')

f1b = figure;
stem(sort(H_blocked(:)))
title('Hessian of \lambda to vecA (matrix entries)')
    

%% Hessian H[lambda; sigma] (Figure 2)
[U, F] = UFdecomp(A);
H2 = Hlambda_sigma(U, F);

f2a = figure;
imagesc(H2)
colorbar
colormap(jet)
title('Hessian of \lambda to \sigma (stage specific survival)')

f2b = figure;
stem(sort(H2(:)))
title('Hessian of \sigma (stage specific survival)')


%% Hessian H[lambda; sigma] bar graphs (Figure 3)
f3 = figure;
subplot(1,3,1)
bar(diag(H2))
title('Figure 3a')

subplot(1,3,2)
bar(H2(1,:))
title('Figure 3b')

subplot(1,3,3)
bar(H2(2,:))
title('Figure 3c')