% Mahmood Amintoosi, IUST 2009
% The previous version only runs for an specified snr indicated in run_homo_MAT_04
% This version is adapted for multi snrs

    alg.list = {'homo_fa_MAT'  'homo_fa_lm_MAT' 'homo_fa_lm_MAT_SSIM'};
    alg.name = {    'LK'    'LK-LM'   'LK-SSIM-LM'};
SNRs = [20,30,50,70];%10:20:100;

desc = 'Tokyo';%'Afsaran';%'Closing';
nAlgorithms = length(alg.name);
maxRunNo = 60;
n_iters =10;


for snrNo = 1:length(SNRs)
    Converging{snrNo} = zeros(maxRunNo,nAlgorithms);
    Times{snrNo} = zeros(maxRunNo,nAlgorithms);
    for algNo=1:nAlgorithms
      rms{snrNo,algNo} = zeros(maxRunNo,n_iters);
      ssim{snrNo,algNo} = rms{algNo};
      ConvergedIter{snrNo,algNo} = zeros(maxRunNo,n_iters);
    end
end

% ConvergedIter{algoNo}{runNo,:) is a vector with size n_iters, which only
% one element of it (i-th) may be equal 1 indicating that the algorithm algNo in
% run number runNo, converged on i-th iteration.
disp('Starting Time')
datestr(now)

% for snrNo = 1:length(SNRs)
%     snr = SNRs(snrNo);
%     for runNo = 1:maxRunNo
%         [snr runNo]
%         %LR image name for tokyo sequence, for simple two images, comment
%         %it and delete the correspondent parameter from function
%         lrImName = sprintf('%04d',(runNo-1) * 5); 
%         [converged,convergedIter,rms_error,ssim_error,times] = run_homo_MAT_10_LM(n_iters,alg,snr,lrImName);
%         Converging{snrNo}(runNo,:) = converged;
%         for algNo=1:nAlgorithms
%           rms{snrNo,algNo}(runNo,:) = rms_error{algNo};
%           ssim{snrNo,algNo}(runNo,:) = ssim_error{algNo};
%           ConvergedIter{snrNo,algNo}(runNo,:) = convergedIter(algNo,:);
%           Times{snrNo}(runNo,algNo) = times(algNo);
%         end
%     end
% end

matlabpool('open','local',4)

%Using parfor for parallel executing on 4 cores
parfor snrNo = 1:length(SNRs)
    snr = SNRs(snrNo);
    for runNo = 1:maxRunNo
        [snr runNo]
        %LR image name for tokyo sequence, for simple two images, comment
        %it and delete the correspondent parameter from function
        lrImName = sprintf('%04d',(runNo-1) * 5); 
        [converged,convergedIter,rms_error,ssim_error,times] = run_homo_MAT_10_LM(n_iters,alg,snr,lrImName);
        Converging{snrNo}(runNo,:) = converged;
        for algNo=1:nAlgorithms
%             switch algNo
%                 case 1
                    rms_par{snrNo}(runNo,algNo,:) = rms_error{algNo};
                    ssim_par{snrNo}(runNo,algNo,:) = ssim_error{algNo};
                    ConvergedIter_par{snrNo}(runNo,algNo,:) = convergedIter(algNo,:);
% %                 case 2
%                     rms{snrNo,2}(runNo,:) = rms_error{algNo};
%                     ssim{snrNo,2}(runNo,:) = ssim_error{algNo};
%                     ConvergedIter{snrNo,2}(runNo,:) = convergedIter(algNo,:);
%             end
          Times{snrNo}(runNo,algNo) = times(algNo);
        end
    end
end
matlabpool('close')


disp('Ending Time')
datestr(now)

for snrNo = 1:length(SNRs)
    for runNo = 1:maxRunNo
        for algNo=1:nAlgorithms
            rms{snrNo,algNo}(runNo,:) = rms_par{snrNo}(runNo,algNo,:);
            ssim{snrNo,algNo}(runNo,:) = ssim_par{snrNo}(runNo,algNo,:);
            ConvergedIter{snrNo,algNo}(runNo,:) = ConvergedIter_par{snrNo}(runNo,algNo,:);
        end
    end
end

save result_Tokyo60.mat rms ssim Converging ConvergedIter Times

%% Plotting
% sum(Converging)

% markers = {'rp:','g+-.','bs--','k.-','+:','s:','d:'};
markers = {'rd--','bp-','ks--','m.-','+:','s:','d:'};

% In the following plot we compare the convergence frequency
h=figure(30); clf;
for snrNo = 1:length(SNRs)
    convFreq{snrNo} = zeros(1,nAlgorithms);
    for algNo=1:nAlgorithms
        convFreq{snrNo}(algNo)=(sum(sum(ConvergedIter{snrNo,algNo}))/maxRunNo)*100;
    end
end
temp = reshape([convFreq{:}],[nAlgorithms,length(SNRs)])';
bar(temp,'group');
% for algNo=1:nAlgorithms
%     text(algNo-.05,convFreq(algNo)-3,sprintf('%4.1f%%',convFreq(algNo)));
% end
ylabel('Convergence Frequency'); 
xlabel('Registration Method');
legend(alg.name);
for snrNo = 1:length(SNRs)
  snrText{snrNo} = num2str(SNRs(snrNo));
end
set(gca,'XTickLabel',snrText);
fileName = sprintf('%s_convergence_frequency.eps',desc);
print_EPS(h,fileName);

%% ------ Which algorithm converged sooner
% h=figure(31); clf;
% bar(sum(Converging)./size(Converging,1)*100,'g');
% ylabel('Percentage '); 
% xlabel('Registration Method');
% set(gca,'XTickLabel', alg.name);
% fileName = sprintf('%s_convergence_faster.eps',desc);
% print_EPS(h,fileName);

% In the following plot we compare the mean convergence iteration number,
% when all of the algorithms converged.
% a=sum(sum((ConvergedIter{algNo}))) is the number of runs which algorithm
% algNo is converged.
% sum(sum(ConvergedIter{algNo}).*(1:n_ites))/a is the average iteration no
% which algorithm algNo converged, but it is not a usefull criterion,
% because the rate of convergence is different, hence we compare the
% mentioned criterion when all of the methods converged.
clear allConverged meanConvIter
for snrNo = 1:length(SNRs)
    idx = false(maxRunNo,nAlgorithms);
    for algNo=1:nAlgorithms
        %Indicating whether algorithm algNo converged in runNo or not
        idx(:,algNo) = sum(ConvergedIter{algNo},2)>0; 
    end
    allConverged{snrNo} = true(maxRunNo,1);
    for algNo=1:nAlgorithms
        allConverged{snrNo} = allConverged{snrNo} & idx(:,algNo);
    end
end
% ------ Average Cycles needed for convergence
% http://jasonlphd.blogspot.com/2008/01/on-importing-matlab-figures-into-pa
% pers.html
h=figure(32); clf;
for snrNo = 1:length(SNRs)
    meanConvIter{snrNo} = zeros(1,nAlgorithms);
    for algNo=1:nAlgorithms
        tmp=sum(sum((ConvergedIter{snrNo,algNo}(allConverged{snrNo},:))));
        meanConvIter{snrNo}(algNo)=(sum(sum(ConvergedIter{snrNo,algNo}(allConverged{snrNo},:)).*(1:n_iters))/tmp);
    end
end

temp = reshape([meanConvIter{:}],[nAlgorithms,length(SNRs)])';
bar(temp,'group');
% bar(meanConvIter,'c');
% for algNo=1:nAlgorithms
%     text(algNo-.05,meanConvIter(algNo)-.3,sprintf('%4.1f',meanConvIter(algNo)));
% end
ylabel('Average Iteration No.'); 
xlabel('SNR (dB)');
set(gca,'XTickLabel',snrText);
legend(alg.name);
% set(gca,'XTickLabel', alg.name);
% set(gca,'FontName', 'Arial');

fileName = sprintf('%s_mean_conv_iter',desc);
print_EPS(h,fileName);

% print('-dpdf','-loose',sprintf('%s_mean_conv_iter.pdf',desc))

%% -------------------- RMS and SSIM Error
set(h,'DefaultLineLineWidth',2);
% Diagram of specified SNRs
SNRlist = [1:length(SNRs)];
% plotNo = 1;
for snrNo = SNRlist
    h=figure(33); clf;
%     subplot(length(SNRlist),1,plotNo); plotNo = plotNo+1;
    hold on;
    for algNo=1:nAlgorithms
        plot(mean(rms{snrNo,algNo}(allConverged{snrNo},:)),markers{algNo});
    end
    hold off
    legend(alg.name);
    ylabel('RMS Error'); 
    xlabel('Iteration No.');
    title(sprintf('SNR=%2d dB',SNRs(snrNo)));
    fileName = sprintf('%s_rms_error_%2d_dB.eps',desc,SNRs(snrNo));
    print_EPS(h,fileName,[4.3 3]);
end
% title('RMS-Error'); 

%% 
h=figure(34); clf;
set(h,'DefaultLineLineWidth',2);
hold on;
snrNo = length(SNRs);
for algNo=1:nAlgorithms
    plot(mean(rms{snrNo,algNo}),markers{algNo});
end
hold off
% title('RMS-Error'); 
legend(alg.name);
ylabel('RMS Error'); 
xlabel('Iteration No.');
fileName = sprintf('%s_biased_rms_error.eps',desc);
print_EPS(h,fileName,[4.3 3]);

% h=figure(34); clf;
% set(h,'DefaultLineLineWidth',2);
% % set(gcf,'name',sprintf('%s, snr=%d',desc,snr));
% % subplot(211);
% hold on;
% for algNo=1:nAlgorithms
%     plot(mean(ssim{algNo}(allConverged,:)),markers{algNo});
% end
% hold off
% title('SSIM-Error'); 
% legend(alg.name);
% print('-depsc',sprintf('%s_ssim_error.eps',desc))

% figure(34); clf;
% bar(mean(Times),'y');
% title('Times'); 
% set(gca,'XTickLabel', alg.name);

% maxRunNo alg.name alg.list

