function [Xrec , iterres , hiostep , erstep] = DR(Diagonal,Y_abs_noisy,imagetype,Index,itermax,ifer,X0,ifadaptive,ifdisplay)
%% Function of implementing DR followed by ER if random phase illumination
%% is exactly known, i.e. delta = 0

% Diagonal:       illumination
% Y_abs_noisy:    Noisy Fourier intensity measurements
% imagetype:      % 1 real positive
                  % 2 complex positive
                  % -2 general complex
% Index:          information about image sizes
% itermax:        maximal number of iterations for DR
% ifer:           maximal number of iterations for ER
% X0:             starting point for images
% ifadaptive:     adapt the termination of DR
% ifdisplay:      whether to display results at each iteration


%% parameters
i            = sqrt(-1);
N11          = Index(2,1);
N22          = Index(2,2);
N1start   = Index(3,1);
N1end     = Index(3,2);
N2start   = Index(4,1);
N2end     = Index(4,2);
P         = size(Diagonal,3);

%% initiation
Xk           = X0;
Imagek       = zeros(N11,N22);
Imagek(N1start:N1end,N2start:N2end) = projection(Xk(N1start:N1end,N2start:N2end),imagetype);
Y_abs_fro    = 0;
for k = 1 : P
    Y_abs_fro = Y_abs_fro + norm(Y_abs_noisy(:,:,k),'fro')^2;
end
Y_abs_fro    = sqrt(Y_abs_fro);

%% error
iterres      = [];        % relative residual at each iteration

%% to be adaptive
increase = 0;

%% ASR iteration
for iter = 1 : itermax
    %% Pre
    Xpre = Xk;
    %% Po
    XkO = Imagek;
    %% Pf
    XkF = Xk;
    for k = 1 : P
        FXk  = fft2(Diagonal(:,:,k).*XkF);
        FXkP = Y_abs_noisy(:,:,k).*exp(i*angle(FXk)); 
        XkF  = ifft2(FXkP)./Diagonal(:,:,k);
    end
    %% PoPf
    XkOF = zeros(N11,N22);
    XkOF(N1start:N1end,N2start:N2end) = projection(XkF(N1start:N1end,N2start:N2end),imagetype);
    %% updated Xk
    Xk = Xk - XkO - XkF + 2*XkOF;
    %% Imagek
    Imagek  = zeros(N11,N22);
    Imagek(N1start:N1end,N2start:N2end) = projection(Xk(N1start:N1end,N2start:N2end),imagetype);
    %% residual
    reres = 0;
    for k = 1 : P
        reres = reres + norm(abs(fft2(Diagonal(:,:,k).*Imagek))-Y_abs_noisy(:,:,k),'fro')^2;
    end
    reres = sqrt(reres)/Y_abs_fro;
    iterres(iter)  = reres;
    %% display the results
    if ifdisplay == 1
        fprintf('%6.0f res = %6.4f  \n',iter,reres)
    end
    %% Adapt DEER and terminate it when the residual keeps increasing in 5
    %% steps
    if ifadaptive == 1 
        if (imagetype == 1) || (imagetype == 2)
            %% stop if the residual doesn't drop in 5 steps
            if iterres(end) >= iterres(max(end-1,1))
                increase = increase + 1;
            else
                increase = 0;
            end
            if increase >= 5 && iter >= 10
                break;
            end
        end
    end
    %% The other adaptive rule that can be applied
% %     if ifadaptive == 1 && imagetype == 1
% %         if norm(Xk-Xpre,'fro')/norm(Xpre,'fro') < 0.01
% %             break;
% %         end
% %     end
end

hiostep  = length(iterres);

%% ER iteration
for iter = 1 : ifer
    %%
    Xpre = Xk;
    %% update X
    XkP   = Xk;
    for k = 1 : P
        FXk  = fft2(Diagonal(:,:,k).*XkP);
        FXkP = Y_abs_noisy(:,:,k).*exp(i*angle(FXk)); %*(FXk./abs(FXk));
        XkP  = ifft2(FXkP)./Diagonal(:,:,k);
    end
    %%%% update Xk to fit the the space domain constraint   
    Xk = zeros(N11,N22);
    Xk(N1start:N1end,N2start:N2end) = projection(XkP(N1start:N1end,N2start:N2end),imagetype);
    %% residual
    reres = 0;
    for k = 1 : P
        reres = reres + norm(abs(fft2(Diagonal(:,:,k).*Xk))-Y_abs_noisy(:,:,k),'fro')^2;
    end
    reres = sqrt(reres)/Y_abs_fro;
    iterres(hiostep+iter)  = reres;
    %% display the results
    if ifdisplay == 1
        fprintf('%6.0f res = %6.4f \n',iter,reres)
    end
    if ifadaptive  == 1
        if norm(Xk-Xpre,'fro')/norm(Xpre,'fro') < 0.0005
            break;
        end
    end
end

erstep  = length(iterres) - hiostep;
Xrec    = Xk;



