itload('../merger_debug.it');

ndat = 1000; % check if true!

XL = [-1 4];
YL= XL;

iters = 0:9

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,:),w_is_' si ')']);
    if it==iters(1),
    ylabel('Importance density'); 
    end
    set(gca,'XLim',XL);
    set(gca,'YLim',YL);
    noi = noi+1;
end

for it=0:25
    si = num2str(it);
    eval(['std_w(it+1) = std(w_is_' si '*200);']);
    eval(['prop_mean(1:2,it+1) = Mpred_mean' si ';']);
end

figure(2)
subplot(1,3,1)
hold off
plot(std_w)
xlabel('iteration');
ylabel('standard deviation');
title('Non-normalized importance weights')

subplot(1,3,2);
hist (w_is_15*200);
title('Histogram of importance weights');

subplot(1,3,3);
hold off
plot(prop_mean');
hold on
plot(ones(size(prop_mean,2),1)*[3 2],'--')
xlabel('iteration');
title('Mean value of the merged density');

% 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