itload('iter_cond_debug.it');

ndat = 1000; % check if true!

XL = [-3 3];
YL= XL;

iters =[];% [0 2 39]

figure(1);
niters = length(iters);
noi = 1;
for it=iters
%    figure(it+1)
    subplot(5,niters,noi);
    si = num2str(it);
    eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),exp(Mpdf' si '''))']);
    if it==iters(1),
    ylabel('Proposal density ');
    end
    set(gca,'XLim',XL);
    set(gca,'YLim',YL);
        
    subplot(5,niters,noi+niters);
    eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),exp(lW' si '(1,:)))']);
    if it==iters(1),
    ylabel('First source'); 
    end
    set(gca,'XLim',XL);
    set(gca,'YLim',YL);
    
    subplot(5,niters,noi+2*niters);
    eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),exp(lW' si '(2,:)))']);
    if it==iters(1),
    ylabel('Second source');
    end
    set(gca,'XLim',XL);
    set(gca,'YLim',YL);
    
    
    subplot(5,niters,noi+3*niters);
    hold off
    eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),w' si ''')']);
    if it==iters(1),
    ylabel('Merged density'); 
    end
    set(gca,'XLim',XL);
    set(gca,'YLim',YL);
    
    subplot(5,niters,noi+4*niters);
    hold off
    eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),exp(lw_m' si '''))']);
    if it==iters(1),
    ylabel('Importance density'); 
    end
    set(gca,'XLim',XL);
    set(gca,'YLim',YL);
    noi = noi+1;
end

prop_mean=zeros(2,40);
prop_var=zeros(2,40);
for it=0:39
    si = num2str(it);
    eval(['eff_w(it+1) = 1/(w' si '''*w' si ');']);
    eval(['prop_mean(1:2,it+1) = Mpred_mean' si ';']);
    eval(['prop_var(1:2,it+1) = Mpred_var' si '(1:2);']);
end

figure(2)
subplot(1,2,1)
hold off
plot(eff_w)
xlabel('iteration');
ylabel('Effective sample size');
%title('Non-normalized importance weights')

subplot(1,2,2);
hold off
plot(prop_mean');
hold on
plot(prop_mean'+2*sqrt(prop_var'),'.');
plot(prop_mean'-2*sqrt(prop_var'),'.');

%plot(ones(size(prop_mean,2),1)*[3 2],'--')
xlabel('iteration');
title('Mean value of the joint merger $\tilde{f}(\bm{x})$');

% itload('../merger_iter_test.it');
% XG = reshape(Grid(1,:),Npoints,Npoints);
% YG = reshape(Grid(2,:),Npoints,Npoints);
% 
% figure(2);
% M1 = reshape(exp(Res1),Npoints,Npoints);
% contour(XG,YG,M1,7);
% 
% figure(3);
% hold off
% mm = max(max(Res2));
% for i=1:size(Res2,1)
%     M2 = reshape(Res2(i,:),Npoints,Npoints);
%     contour(XG,YG,M2,[0:mm/7:mm]);
%     hold on
% end