%computing limit of parametric path on a generalization of Levent's example
%when n=4 and sd=2

x11 = sym('x11','positive'); 
x33 = sym('x33','positive');
a = sym('a','positive');
datain = load('alogdetavec');
SD1 = load('SD1');
X1 = SD1.Xstore;
SD2 = load('SD2');
X2 = SD2.Xstore;
SD3 = load('SD3');
X3 = SD3.Xstore;
alogdetasd3 = datain.alogdetavec(:,1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%% sd=1 case
Xasd1 = [(1+3*a)/3 0 0 0
    0 (1+3*a)/3 0 0
    0 0 (1+3*a)/3 0
    0 0 0 a];
detXa = a*(1+3*a)^3/27;
alogdetasd1 = a*log(detXa);


%%%%%%%%%%%%%%%%%%%%%%%%%%%% sd=2 case
Xa = [x11 0 0 (1/2)*(x33-a)
    0 (1+3*a-x33-x11) 0 0
    0 0 x33 0
    (1/2)*(x33-a) 0 0 a];
detXa = det(Xa);
dx11 = diff(detXa,'x11');
x11val = simplify(solve(dx11,x11));
dx33 = simplify(diff(detXa,x33));
dx33new = simplify(subs(dx33,x11,x11val));
x33val = simplify(solve(dx33new,x33));
x33val = x33val(2);
x11final = simplify(subs(x11val,x33,x33val));

Xasd2 = subs(Xa,x33,x33val);
Xasd2 = subs(Xasd2,x11,x11final);

deta = simplify(subs(detXa,x11,x11final));
deta = simplify(subs(deta,x33,x33val));
alogdetasd2 = a*log(deta);


% %%%%%%%%%%%%%%%%%%%%%%%%%%%% sd=3 case
% x22 = 1+3*a-x33-x11;
% x14 = (x33-a)/2;
% x13 = (x22-a)/2;
% x34 = (1+2*a-x33-x11)*(x33-a)/(4*x11);
% Xa = [x11 0 x13 x14
%     0 x22 0 0 
%     x13 0 x33 x34
%     x14 0 x34 a];
% detXa = simplify(det(Xa));
% dx11 = simplify(diff(detXa,'x11'));
% x11val = simplify(solve(dx11,x11));
% dx33 = simplify(diff(detXa,x33));
% dx33new = simplify(subs(dx33,x11,x11val));
% x33val = simplify(solve(dx33new,x33));


%%%%%%%%%%%%%%%%%%%%%%%%%%%% plots
tts = logspace(-1,-14,100);
f1=zeros(1,length(tts));
f2=zeros(1,length(tts));
ind=0;
for tt=tts
	ind=ind+1;
	f1(ind)=subs(alogdetasd1,a,tt);
    f2(ind)=subs(alogdetasd2,a,tt);
end
figure(1)
plot(tts,f1,'k')
hold on
plot(tts,f2,'r')
%plot(tts,alogdetasd3,'b')
legend('sd1','sd2','Location','best')
title('Function value')
hold off


%%%%%%%%%%%%%%%%%%%%%%%%%%%% rate of convergence plots (function value)
ratesd1 = abs(f1(2:end))./abs(f1(1:end-1));
ratesd2 = abs(f2(2:end))./abs(f2(1:end-1));
ratesd3 = abs(alogdetasd3(2:end))./abs(alogdetasd3(1:end-1));
figure(2)
plot(ratesd1,'k')
hold on
plot(ratesd2,'r')
%plot(ratesd3(70:end),'b')
legend('sd1','sd2','Location','best')
title('Convergence rate of function value')
hold off



%%%%%%%%%%%%%%%%%%%%%%%%%%%% rate of convergence plots (iterates)
Xsd1 = eye(4)/3;
Xsd1(4,4) = 0;
Xsd2 = zeros(4);
Xsd2(1,1) = 0.6;
Xsd2(2,2) = 0.4;
Xsd3 = zeros(4);
Xsd3(1) = 1;

N = length(tts);
ratesd1 = zeros(N-1,1);
ratesd2 = zeros(N-1,1);
ratesd2adj = zeros(N-1,1);
for i=1:N-1
    ratesd1(i) = norm(subs(Xasd1,a,tts(i+1)) - Xsd1,'fro')/norm(subs(Xasd1,a,tts(i)) - Xsd1,'fro');
    ratesd2(i) = norm(subs(Xasd2,a,tts(i+1)) - Xsd2,'fro')/norm(subs(Xasd2,a,tts(i)) - Xsd2,'fro');
    ratesd2adj(i) = norm(subs(Xasd2,a,tts(i+1)) - Xsd2,'fro')/norm(subs(Xasd2,a,tts(i)) - Xsd2,'fro')^(2^(1/2));
end

figure(3)
clf
plot(ratesd1,'k')
hold on
plot(ratesd2,'r')
%plot(ratesd2adj,'b')
legend('sd1','sd2','Location','best')
title('Convergence Rate of Iterates')
hold off

%%%%%%%%%%%%%%%%%%%%%%%% log convergence for sd1
% lograte = zeros(N-2,1);
% for j=1:N-2
%     lograte(j) = norm(subs(Xasd1,a,tts(j+2)) - subs(Xasd1,a,tts(j+1)),'fro')/norm(subs(Xasd1,a,tts(j+1))-subs(Xasd1,a,tts(j)),'fro');
% end
% 
% figure(4)
% plot(lograte,'k')




%%%%%%%%%%%%%%%%%%%%%%%%%%%% rate of convergence plots (iterates)
N = length(tts);
ratesd1 = zeros(N-1,1);
ratesd2 = zeros(N-1,1);
ratesd3 = zeros(N-1,1);
for i=1:N-1
    ratesd1(i) = norm(X1(:,:,i+1) - Xsd1,'fro')/norm(X1(:,:,i) - Xsd1,'fro');
    ratesd2(i) = norm(X2(:,:,i+1) - Xsd2,'fro')/norm(X2(:,:,i) - Xsd2,'fro');
    ratesd3(i) = norm(X3(:,:,i+1) - Xsd3,'fro')/norm(X3(:,:,i) - Xsd3,'fro');
end

figure(4)
clf
plot(ratesd1,'k')
hold on
%plot(ratesd2,'r')
plot(ratesd3,'b')
legend('sd1','sd3','Location','best')
title('Convergence Rate of Empirical Iterates')
hold off

