% creating the test problem
clear, clc, close all
startup
load('star_cluster')
xex=x_true(:);
N = size(xex,1); 
n = sqrt(N);
A = psfMatrix(PSF25);
A.imsize = [n,n];
b = A*xex;
sigma = 10^-2;
ber = addRELnoise(b,sigma); 
%
% defining some additional inputs for MFCGLS
b0 = ber; b0(b0<0)=0;
x0 = b0; % initial guess for MFCGLS
m_in = 20; % maximum number of inner iterations
m_out = 25; % maximum number of outer iterations
% considering a stopping criterion based on the discrepancy principle
display('beginning of MFCGLS - first test...')
stopC_1.kind = 'discr';
stopC_1.noise = sigma;
[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);
% considering a stopping criterion based on the stabilization of the residual
display('beginning of MFCGLS - second test...')
stopC_2.kind = 'stab';
stopC_2.tau = 1e-5;
[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);
%
% defining some additional inputs for gradient-descent-type methods
display('beginning of MRNSD...')
options = IRset('x_true', xex, 'MaxIter', 400, 'x0', b0);
[X_MRNSD, IterInfo_MRNSD] = IRmrnsd(A, ber, options);
display('beginning of NN SD...')
options = IRset(options, 'StepLength',0);
[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, star\_cluster test problem')
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(sqrt(abs(reshape(xex,256,256))),[0 20]), title('exact image')
figure, imshow(sqrt(abs(reshape(ber,256,256))),[0 20]), title('corrupted image')
figure, imshow(sqrt(abs(reshape(Xiter_RS1(:,indv),256,256))),[0 20]), title('MFCGLS')
figure, imshow(sqrt(abs(reshape(X_MRNSD,256,256))),[0 20]), title('MRNSD')
figure, imshow(sqrt(abs(reshape(X_SDNN,256,256))),[0 20]), title('SDNN')
