function [WSURE WMSE MSE PSNR C u RGCV PJSURE PJMSE] = WSURE_IRLSMIL_Nobeta(y, beta, params)
% Only for Deblur
% Computes RGCV (a Monte-Carlo version of GCV for nonlinear algos) for comparison; RGCV = (1 / Nsamp) * || y - A * f(y) ||^2 / ( 1 - (1 / Nsamp) trace{ J_f(y) } )^2
% Computes Projected MSE and Projected SURE where Projection is computed as A' * (AA')^-1 A and a tol is used in the pseudo inverse to avoid division by numerically small values
% Jaocbian matrix estimation not implemented for the complex case. Derivatives w.r.t. y and conj(y) need to be computed for the complex case
%
% Sathish Ramani, University of Michigan, 2011

%% Load params
xtrue = params.xtrue;
dispinfo = params.dispinfo;
trW = params.trW;
WA = params.WA;
VA = params.VA;
V = params.V;
A = params.A;
R = params.R;
Ay = params.Ay;
sigma_sq = params.sigma_sq;
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;
NPJSURE = params.ndata_PJSURE;
nx1 = N(1);
nx2 = N(2);
mn = params.mn;
mx = params.mx;

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

Ninv = 1 / Npix;
NinvPJSURE = 1 / NPJSURE;
NinvPJMSE = 1 / NPJSURE;
isnrnum = Ninv * sum( abs( y(:) - xtrue(:) ) .^ 2 );

%% Fatrices for some beta-dependent objects
Cresp = ( AAresp + alpha ); % ( + alpha * I ) works well for deblurring
Cresp = Cresp; % Divide by beta for stability

Cinvresp = 1 ./ Cresp;
CinvAAresp = Cinvresp .* AAresp;
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
CinvAA = Gblur_S(true([nx1 nx2]), 'psffft', CinvAAresp, 'type', 'fft,same'); % Corresponds to Cinv * A' * A

CinvAresp = Cinvresp .* conj( PSF );
CinvA = Gblur_S(true([nx1 nx2]), 'psffft', CinvAresp, 'type', 'fft,same'); % Corresponds to Cinv * A'
CinvAy = real( CinvA * y ); % Divide by beta for stability

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, WSURE and cost
Ru = R * u; % Needed for the algo
Jur = A' * r;

CinvAr = ( CinvA * r ); % Divide by beta for stability

b = CinvAy + u - ( CinvAA * u ); % Divide by beta for stability
Rb = R * b;

Jb = CinvAr + Jur - ( CinvAA * Jur ); % Divide by beta for stability
RJb = R * Jb;

v = Rb;
Jvr = RJb;

WSURE = zeros(niter + 1, 1); PJSURE = WSURE; WMSE = WSURE; PJMSE = WMSE; MSE = WMSE; PSNR = MSE; C = PSNR; div = C; PJdiv = C; RGCV = div;
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;

RGCVnume = Ninv * sum( abs( y(:) - Au(:) ) .^ 2 ); 
RGCVdeno = r .* ( A * Jur ); RGCVdeno = ( 1 - Ninv * sum( real( RGCVdeno(:) ) ) ) ^ 2;

difs = xtrue - u;
Adifs = A * difs;
VAdifs = VA * difs;
MSE(1) = Ninv * sum( abs( difs(:) ) .^ 2 );
WMSE(1) = Ninv * sum( abs( Adifs(:) ) .^ 2 );
PJMSE(1) = NinvPJMSE * sum( abs( VAdifs(:) ) .^ 2 );
WSURE(1) = WMSE(1);
PJSURE(1) = PJMSE(1);
RGCV(1) = RGCVnume / RGCVdeno;
PSNR(1) = 10 * log10( isnrnum / MSE(1) );

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') ...
        '; PSNR=' num2str(PSNR(1), '%0.4E') '; MSE=' num2str(MSE(1), '%0.4E') ...
        '; PJMSE=' num2str(PJMSE(1), '%0.4E') '; PJSURE=' num2str(PJSURE(1), '%0.4E') ...
        '; WMSE=' num2str(WMSE(1), '%0.4E') '; WSURE=' num2str(WSURE(1), '%0.4E') ...
        '; RGCV=' num2str(RGCV(1), '%0.4E') ...
        ]);
    figure(3); clf; 
    subplot(2, 2, 1); cla;
    uini = Ay;
    imagesc( uini, [mn mx]);
    colormap gray; colorbar; axis equal off; title(['PSNR: ' num2str(PSNR(1))]);
end

%% ========================= IRLS-MIL algo =========================
for ii = 1 : niter
    switch(penalt)
        case('FP')
            Dinv = delta * ( delta + abs( Ru ) ) / beta;
            
        case('l1')
            Dinv = abs( Ru ) / beta;
            
        case{'TV'}
            Dinv = repmat( sqrt( sum( abs( Ru ) .^ 2, 3 ) ), [1 1 dimsRu(3)] ) / beta;
            Dinvz = Dinv; Dinvz(Dinvz == 0) = eps; % Makesure 1/Dinv does not blow later in dtheta calc
    end
    
    Hinv = 1 ./ ( Dinv + lambda_max );
    % disp(['MinHinv = ' num2str(min(Hinv(:))) '; MaxHinv = ' num2str(max(Hinv(:))) '; CondHinv = ' num2str(max(Hinv(:)) / min(Hinv(:)))]);
    
    switch(penalt)
        case{'FP'}
            dtheta = -delta / beta * sign( Ru ) .* ( R * Jur );
            
        case{'l1'}
            dtheta = -sign( Ru ) .* ( R * Jur ) / beta;
            
        case{'TV'}
            dtheta = -repmat( sum( Ru ./ Dinvz .* ( R * Jur ) / beta, 3 ), [1 1 dimsRu(3)] );
    end
    
    for kk = 1 : niter_in
        rhsv = Rb + lambda_max * v - R * ( Cinv * ( R' * v ) );
        Jrhsv = ( RJb + lambda_max * Jvr - R * ( Cinv * ( R' * Jvr ) ) );
        v = Hinv .* rhsv ;
        Jvr = rhsv .* dtheta .* Hinv .^ 2 + Jrhsv .* Hinv;
    end
    
    % Update u & Ju
    u = b - Cinv * ( R' * v );
    Jur = Jb - Cinv * ( R' * Jvr );

    b = CinvAy + u - ( CinvAA * u ); % Divide by beta for stability
    Rb = R * b;
    
    Jb = CinvAr + Jur - ( CinvAA * Jur ); % Divide by beta for stability
    RJb = R * Jb;
    
    Ru = R * u; % Needed for the algo
    
    % Compute the divergence term
    divvr = r .* ( A * Jur ); % Predictive SURE
    divvr = sum( divvr(:) );
    div(ii + 1) = divvr;
    
    divvr = r .* ( WA * Jur ); % Projected SURE
    divvr = sum( divvr(:) );
    PJdiv(ii + 1) = divvr;

    % Prepare to compute MSE, WMSE, WSURE, etc
    Au = A * u; % Needed for WMSE, WSURE and cost

    difs = xtrue - u;
    Adifs = A * difs;
    VAdifs = VA * difs;
    Vdifs = V * (y - Au);
    
    MSE(ii + 1) = Ninv * sum( abs( difs(:) ) .^ 2 );
    WMSE(ii + 1) = Ninv * sum( abs( Adifs(:) ) .^ 2 );
    PJMSE(ii + 1) = NinvPJMSE * sum( abs( VAdifs(:) ) .^ 2 );
    WSURE(ii + 1) = Ninv * sum( abs( y(:) - Au(:) ) .^ 2 ) - sigma_sq + 2 * Ninv * sigma_sq * div(ii + 1);
    PJSURE(ii + 1) = NinvPJSURE * sum( abs( Vdifs(:) ) .^ 2 ) - sigma_sq * trW * NinvPJSURE + 2 * NinvPJSURE * sigma_sq * PJdiv(ii + 1);
    RGCV(ii + 1) = Ninv * sum( abs( y(:) - Au(:) ) .^ 2 ) / ( 1 - Ninv * div(ii + 1) ) ^ 2;
    PSNR(ii + 1) = 10 * log10( isnrnum / MSE(ii + 1) );
    
    % 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') ...
            '; PSNR=' num2str(PSNR(ii), '%0.4E') '; MSE=' num2str(MSE(ii), '%0.4E') ...
            '; PJMSE=' num2str(PJMSE(ii), '%0.4E') '; PJSURE=' num2str(PJSURE(ii), '%0.4E') ...
            '; WMSE=' num2str(WMSE(ii), '%0.4E') '; WSURE=' num2str(WSURE(ii), '%0.4E') ...
            '; RGCV=' num2str(RGCV(ii + 1), '%0.4E') ...
            ]);
        
        
        figure(2); clf; 
        imagesc( u, [mn mx]);
        colormap gray; colorbar; axis equal off; title(['PSNR: ' num2str(PSNR(ii))]);

        figure(3);
        subplot(2, 2, 1); cla; plot( log10( MSE( 1 : ii + 1) ) ); title('MSE');
        subplot(2, 2, 2); cla; plot( log10( WMSE( 1 : ii + 1) ), 'r'); hold on; plot( log10( abs(WSURE( 1 : ii + 1 ) ) ), 'k'); title('WMSE');
        subplot(2, 2, 3); cla; plot( log10( PJMSE( 1 : ii + 1) ), 'r'); hold on; plot( log10( abs(PJSURE( 1 : ii + 1 ) ) ), 'k'); title('PJMSE');
        subplot(2, 2, 4); cla; plot( log10( RGCV( 1 : ii + 1) ), 'k'); title('RGCV');
    end
end

if(dispinfo)
    figure(3);
    subplot(2, 2, 1); cla; plot( log10( MSE ) );
    subplot(2, 2, 2); cla; plot( log10( WMSE ), 'r'); hold on; plot( log10( abs( WSURE ) ), 'k'); legend('WMSE in red', 'WSURE in black');
    subplot(2, 2, 3); cla; plot( log10( PJMSE ), 'r'); hold on; plot( log10( abs( PJSURE ) ), 'k'); legend('PJMSE in red', 'PJSURE in black');
    subplot(2, 2, 4); cla; plot( log10( RGCV( 1 : ii + 1) ), 'k'); legend('RGCV in black');
    
    figure(2); clf;
    imagesc( u, [mn mx]);
    colormap gray; colorbar; axis equal off; title(['PSNR: ' num2str(PSNR(ii))]);
end