clear, clc, close all
startup
load('satellite')
A = psfMatrix(PSF);
xex = x_true(:);
b = A*xex;
nl = 10^-1;
ber = addRELnoise(b,nl);
% setting the initial guess
b0 = ber;
b0(b0<0)=0;
m_in = 20; m_out = 20;
display('beginning of MFCGLS - first test...')
stopC_1.kind = 'discr';
stopC_1.noise = nl;
[Xiter_RS1, RelErr_RS1, RelRes_RS1, Alpha_RS1, tot_it1, cycle_it1, stopC_it1] = fcgls_rp_nn_ReSt (A, ber, m_in, m_out, xex, b0, 1, stopC_1);
display('beginning of MFCGLS - second test...')
stopC_2.kind = 'stab';
stopC_2.tau = 1e-10;
[Xiter_RS2, RelErr_RS2, RelRes_RS2, Alpha_RS2, tot_it2, cycle_it2, stopC_it2] = fcgls_rp_nn_ReSt (A, ber, m_in, m_out, xex, b0, 1, stopC_2);
display('beginning of MRNSD...')
options = IRset('x_true', xex, 'MaxIter', 200, 'x0', b0);
[X_MRNSD, IterInfo_MRNSD] = IRmrnsd(A, ber, options);
options = IRset(options, 'StepLength', -1);
display('beginning of NN SD...')
[X_SDNN, IterInfo_SDNN] = IRgdnn(A, ber, options);
%
% plotting the results (Relative Errors History),
% highlighting the different stopping iterations for MFCGLS
nzst1 = find(stopC_it1~=0,1); % discrepancy principle
if nzst1 > 1
    totit_sc1 = sum(cycle_it1(1:nzst1-1))+stopC_it1(nzst1);
else
    totit_sc1 = stopC_it1(nzst1);
end
nzst2 = find(stopC_it2~=0,1); % stabilization of the residual
if nzst2 > 1
    totit_sc2 = sum(cycle_it2(1:nzst2-1))+stopC_it2(nzst2);
else
    totit_sc2 = stopC_it2(nzst2);
end
figure, semilogy(RelErr_RS1)
hold on
semilogy(IterInfo_MRNSD.Enrm, '-g')
semilogy(IterInfo_SDNN.Enrm, '-r')
title('relative error history, satellite test problem, Gaussian noise')
legend('MFCGLS', 'MRNSD', 'SDNN')
semilogy(totit_sc1, RelErr_RS1(totit_sc1), 'om')
semilogy(totit_sc2, RelErr_RS2(totit_sc2), 'hm')
%
% showing the best restored images
[minv, indv] = min(RelErr_RS1);
figure, imshow(reshape(xex,256,256),[]), title('exact image')
figure, imshow(reshape(ber,256,256),[]), title('corrupted image')
figure, imshow(reshape(Xiter_RS1(:,indv),256,256),[]), title('MFCGLS')
figure, imshow(reshape(X_MRNSD,256,256),[]), title('MRNSD')
figure, imshow(reshape(X_SDNN,256,256),[]), title('SDNN')
