function [AGCV C u] = AGCV_IRLSMIL(beta, params)
% Only for Deblur
% Computes AGCV - a Monte-Carlo approximate GCV suggested by Giryes Elad Eldar et al for nonlinear algos applied to inverse problems such as deblurring

%% Load params
dispinfo = params.dispinfo;
A = params.A;
R = params.R;
penalt = params.Penalty;
niter = params.niter;
niter_in = params.niter_in;
r = params.r;
PSF = params.PSF;
mnAA = params.mnAA;
mxAA = params.mxAA;
AAresp = params.AAr;
RRresp = params.RRr;
alpha = params.alpha;
N = params.N;
Npix = params.Npix;
nx1 = N(1);
nx2 = N(2);

switch(penalt)
    case{'FP'}
        delta = params.delta;
        
    otherwise
        delta = 1;
end

Ninv = 1 / Npix;

y = r;
Ay = A' * y;

%% Fatrices for some beta-dependent objects
Cresp = ( AAresp + alpha ) / beta;
Cinvresp = 1 ./ Cresp;
CAAresp = Cinvresp .* AAresp;
CAresp = Cinvresp .* conj(PSF);
lambda_max = max(RRresp(:) .* Cinvresp(:)); % Find max eigenvalue of R * C^-1 * R'

Cinv = Gblur_S(true([nx1 nx2]), 'psffft', Cinvresp, 'type', 'fft,same'); % Corresponds to Cinv * A' * A
CAA = Gblur_S(true([nx1 nx2]), 'psffft', CAAresp, 'type', 'fft,same'); % Corresponds to Cinv * A' * A
CA = Gblur_S(true([nx1 nx2]), 'psffft', CAresp, 'type', 'fft,same'); % Corresponds to Cinv * A'
CAy = real( CA * y ) / beta;

if(dispinfo)
    disp(['Min AA = ' num2str(mnAA) '; Max AA = ' num2str(mxAA) ...
    '; Condition number of AA = ' num2str(mxAA / mnAA) ...
    '; Condition number of C = ' num2str(max(Cresp(:)) / min(Cresp(:))) ...
    '; alpha = ' num2str(alpha) ...
    '; Lambda_max = ' num2str(lambda_max) ...
    ]);
end

%% Initial estimates for AL algo
u = Ay;
Au = A * u; % Needed for WMSE, AGCV and cost
Ru = R * u; % Needed for the algo

b = CAy + u - CAA * ( u / beta );
Rb = R * b;
v = Rb;

AGCV = zeros(niter + 1, 1); C = AGCV; tr = AGCV;
dimsRu = size(Ru);

%% Compute cost and display
switch(penalt)
    case{'FP'}
        regterm = sum( abs( Ru(:) ) ) / delta - sum( log( 1 + abs( Ru(:) ) / delta ) ); % Redundant, can be avoided if not monitoring cost
        
    case{'l1'}
        regterm = sum( abs( Ru(:) ) ) ; % Redundant, can be avoided if not monitoring cost
        
    case{'TV'}
        regterm = sqrt( sum( abs( Ru ) .^ 2, 3 ) ); regterm = sum( regterm(:) );
end
cnew = 0.5 * sum( abs( y(:) - Au(:) ) .^ 2 ) + beta * regterm;
C(1) = cnew;

AGCVnume = Ninv * sum( abs( y(:) - Au(:) ) .^ 2 ); 
AGCVdeno = r(:) .* ( Au(:) ); AGCVdeno = ( 1 - Ninv * sum( AGCVdeno(:) ) ) ^ 2;
AGCV(1) = AGCVnume / AGCVdeno;

if(dispinfo)
    cold = 2 * cnew; % Redundant, can be avoided if not monitoring cost
    dcost = cnew - cold; % Redundant, can be avoided if not monitoring cost
    disp(['IRLSMIL0; C=' num2str(cnew, '%0.5E') '; DC= ' num2str(dcost, '%0.5E') ...
        '; AGCV=' num2str(AGCV(1), '%0.4E') ...
        ]);
    figure(6); clf; imagesc( u ); colormap gray; colorbar; axis equal off;
end

%% ========================= IRLS-MIL algo =========================
for ii = 1 : niter
    switch(penalt)
        case('FP')
            Dinv = delta * ( delta + abs( Ru ) );
            
        case('l1')
            Dinv = abs( Ru );
            
        case{'TV'}
            Dinv = repmat( sqrt( sum( abs( Ru ) .^ 2, 3 ) ), [1 1 dimsRu(3)] );
    end
    
    b = CAy + u - CAA * ( u / beta ); % Update b
    Rb = R * b;
    Hinv = 1 ./ ( Dinv + lambda_max );
    
    for kk = 1 : niter_in
        v = Hinv .* ( Rb + (lambda_max * v - R * ( Cinv * ( R' * v ) ) ) );
    end
    
    % Update u & Ju
    u = b - Cinv * ( R' * v );

    % Prepare to compute MSE, WMSE, AGCV, etc
    Ru = R * u; % Needed for the algo
    Au = A * u; % Needed for WMSE, AGCV and cost

    % Compute the divergence term
    tr(ii + 1) = sum( r(:) .* Au(:) );
    AGCV(ii + 1) = Ninv * sum( abs( y(:) - Au(:) ) .^ 2 ) / ( 1 - Ninv * tr(ii + 1) ) ^ 2;
    
    % Compute cost and display
    switch(penalt)
        case{'FP'}
            regterm = sum(abs(Ru(:))) / delta - sum( log(1 + abs(Ru(:)) / delta) ); % Redundant, can be avoided if not monitoring cost
            
        case{'l1'}
            regterm = sum(abs(Ru(:))) ; % Redundant, can be avoided if not monitoring cost
            
        case{'TV'}
            regterm = sqrt( sum( abs( Ru ) .^ 2, 3 ) ); regterm = sum( regterm( : ) );
    end
    cold = cnew; % Redundant, can be avoided if not monitoring cost
    cnew = 0.5 * sum( abs( y(:) - Au(:) ) .^ 2 ) + beta * regterm;
    C(ii + 1) = cnew;
    if(dispinfo)
        dcost = cnew - cold; % Redundant, can be avoided if not monitoring cost
        disp(['IRLSMIL' int2str(ii) '; C=' num2str(cnew, '%0.5E') '; DC= ' num2str(dcost, '%0.5E') ...
            '; AGCV=' num2str(AGCV(ii + 1), '%0.4E') ...
            ]);
        
        figure(6); clf; imagesc( u ); colormap gray; colorbar; axis equal off;
        figure(7); clf; plot(AGCV(1 : ii + 1));
    end
end