function [WSURE WMSE MSE PSNR C u RGCV tE] = WSURE_IRLSMIL_CSMRI_Nobeta( y, beta, params )
% Only for Cartesian CSMRI Recon
% 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 Predicted MSE and Predicted SURE 
% 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

%% Load params
R = params.R;
mn = params.mn;
mx = params.mx;

switch params.Penalty
    case{'FP'}
        delta = params.delta;
        
    otherwise
        delta = 1;
end

Ninv = 1 / params.Npix;
Ndinv = 1 / params.Ndata;
psnrnum = max( abs( params.xtrue(:) ) .^ 2 );

%% Fatrices for some beta-dependent objects
Cresp = params.PSF2 + params.alpha; % Divide by beta for stability

Cinvresp = 1 ./ Cresp;
rho = max( params.RRr(:) .* Cinvresp(:) ) * params.rhofact; % Find max eigenvalue of R * C^-1 * R'
Cinv = Gblur_S( true( params.N ), 'psffft', Cinvresp, 'type', 'fft,same' ); % Corresponds to Cinv
CinvAA = Gblur_S( true( params.N ), 'psffft', Cinvresp .* params.PSF2, 'type', 'fft,same' ); % Corresponds to Cinv * A' * A
CinvAy = Cinv * params.Ay;

if params.dispinfo
    disp( [ 'Condition number of C = ' num2str( max( Cresp(:) ) / min( Cresp(:) ) ) '; rho = ' num2str( rho ) ] );
end

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

Jur = ( params.A )' * params.r;
Jucr = zeros( size( Jur ) );
CinvAr = Cinv * ( params.A )' * params.r;

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

Jbr = CinvAr + Jur - CinvAA * Jur;
Jbcr = Jucr - CinvAA * Jucr; 

RJbr = R * Jbr;
RJbcr = R * Jbcr;

v = Rb;
Jvr = RJbr;
Jvcr = RJbcr;

WSURE = zeros(params.niter + 1, 1); WMSE = WSURE; MSE = WMSE; PSNR = MSE; C = PSNR; div = C; RGCV = div; divtrue = div; tE = div;

%% Compute cost and display
switch params.Penalty
    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 = Ndinv * sum( abs( y(:) - Au(:) ) .^ 2 ); 
RGCVdeno = params.r .* ( params.A * Jur ); RGCVdeno = ( 1 - Ndinv * sum( real( RGCVdeno(:) ) ) ) ^ 2;

difs = params.xtrue - u;
Adifs = params.A * difs;
MSE(1) = Ninv * sum( abs( difs(:) ) .^ 2 );
WMSE(1) = Ndinv * sum( abs( Adifs(:) ) .^ 2 );
WSURE(1) = WMSE(1);
RGCV(1) = RGCVnume / RGCVdeno;
PSNR(1) = 10 * log10( psnrnum / MSE(1) );

if params.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 / ' int2str(params.niter) ' ; C=' num2str(cnew, '%0.5E') '; DC= ' num2str(dcost, '%0.5E') ...
        '; PSNR=' num2str(PSNR(1), '%0.4E') '; MSE=' num2str(MSE(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 = params.Ay;
    imagesc( abs( uini ), [mn mx] );
    colormap gray; colorbar; axis equal off; title( [ 'PSNR: ' num2str( PSNR(1) ) ] );
end

%% ========================= IRLS-MIL algo =========================
for ii = 1 : params.niter
    tS = tic;
    
    switch params.Penalty
        case('FP')
            Dinv = delta * ( delta + abs( Ru ) ) / beta;
            
        case('l1')
            Dinv = abs( Ru ) / beta;
            
        case{'TV'}
            error(' TV: Yet to be done ...');
    end
    
    Hinv = 1 ./ ( Dinv + rho );
    
    RJur = R * Jur;
    RJucr = R * Jucr;
    
    switch params.Penalty
        case{'FP'}
            RuD = Ru ./ abs( Ru );
            Jgammar = 0.5 / beta * delta * conj( RuD ) .* RJur + 0.5 / beta * delta * RuD .* conj( RJucr );
            Jgammacr = 0.5 / beta * delta * conj( RuD ) .* RJucr + 0.5 / beta * delta * RuD .* conj( RJur );
            
        case{'l1'}
            RuD = Ru ./ Dinv;
            Jgammar = 0.5 / beta ^ 2 * conj( RuD ) .* RJur + 0.5 / beta ^ 2 * RuD .* conj( RJucr );
            Jgammacr = 0.5 / beta ^ 2 * conj( RuD ) .* RJucr + 0.5 / beta ^ 2 * RuD .* conj( RJur );
            
        case{'TV'}
            % dtheta = -repmat( sum( Ru ./ Dinvz .* ( R * Jur ), 3 ), [1 1 dimsRu(3)] );
            error(' TV: Yet to be done ...');
    end
    JgH = Jgammar .* Hinv .^ 2;
    JgcH = Jgammacr .* Hinv .^ 2;
    
    for jj = 1 : params.niter_j
        for kk = 1 : params.niter_k 
            rhsv = Rb + rho * v - R * ( Cinv * ( R' * v ) );
            v = Hinv .* rhsv ;
            
            Jrhsvr = RJbr + rho * Jvr - R * ( Cinv * ( R' * Jvr ) );
            Jrhsvcr = RJbcr + rho * Jvcr - R * ( Cinv * ( R' * Jvcr ) );
            
            Jvr = Jrhsvr .* Hinv - rhsv .* JgH;
            Jvcr = Jrhsvcr .* Hinv - rhsv .* JgcH;
        end
        
        % Update u & Ju
        u = b - Cinv * ( R' * v );
        Jur = Jbr - Cinv * ( R' * Jvr );
        Jucr = Jbcr - Cinv * ( R' * Jvcr );
        
        b = CinvAy + u - CinvAA * u;
        Rb = R * b;
        
        Jbr = CinvAr + Jur - CinvAA * Jur; 
        Jbcr = Jucr - CinvAA * Jucr; % Divide by beta for stability
        
        RJbr = R * Jbr;
        RJbcr = R * Jbcr;
        
        Ru = R * u; % Needed for the algo
    end
    
    % Compute the divergence term
    divvr = params.r .* ( params.A * Jur ); % Predictive SURE
    divvr = sum( real( divvr(:) ) );
    div(ii + 1) = divvr;
    
    % Prepare to compute MSE, WMSE, WSURE, etc
    Au = params.A * u; % Needed for WMSE, WSURE and cost
    
    if params.dispinfo && ~params.realdata
        divtrue( ii + 1 ) = real( sum( conj( params.nois(:) ) .* Au(:) ) ) / params.sigma_sq;
    elseif params.realdata
        divtrue( ii + 1 ) = 1;
    end
        
    difs = params.xtrue - u;
    Adifs = params.A * difs;
    resd = Ndinv * sum( abs( y(:) - Au(:) ) .^ 2 );
    
    MSE(ii + 1) = Ninv * sum( abs( difs(:) ) .^ 2 );
    WMSE(ii + 1) = Ndinv * sum( abs( Adifs(:) ) .^ 2 );
    WSURE(ii + 1) = resd - params.sigma_sq + 2 * Ndinv * params.sigma_sq * div(ii + 1);
    RGCV(ii + 1) = resd / ( 1 - Ndinv * div(ii + 1) ) ^ 2;
    PSNR(ii + 1) = 10 * log10( psnrnum / MSE(ii + 1) );
    
    tE( ii + 1 ) = toc(tS);
    
    % Compute cost and display
    switch params.Penalty
        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 params.dispinfo
        dcost = cnew - cold; % Redundant, can be avoided if not monitoring cost
        disp(['IRLSMIL' int2str(ii) '/' int2str(params.niter) '; C=' num2str(cnew, '%0.5E') '; DC= ' num2str(dcost, '%0.5E') ...
            '; PSNR=' num2str(PSNR(ii), '%0.4E') '; MSE=' num2str(MSE(ii), '%0.4E') ...
            '; WMSE=' num2str(WMSE(ii), '%0.4E') '; WSURE=' num2str(WSURE(ii), '%0.4E') ...
            '; RGCV=' num2str(RGCV(ii + 1), '%0.4E') ...
            '; Time = ' num2str( tE( ii + 1 ) ) ...
            '; DIVest=' num2str(div(ii + 1), '%0.4E') '; DIVtrue=' num2str(divtrue(ii + 1), '%0.4E') ...
            '; ratio=' num2str(divtrue(ii + 1) / div(ii + 1), '%0.4E') ...
            ]);
        
        
        figure(2); clf; 
        imagesc( abs( 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, 4); cla; plot( log10( RGCV( 1 : ii + 1) ), 'k'); title('RGCV');
    end
end

if params.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, 4); cla; plot( log10( RGCV( 1 : ii + 1) ), 'k'); legend('RGCV in black');
    
    figure(2); clf;
    imagesc( abs( u ), [mn mx]); 
    colormap gray; colorbar; axis equal off; title(['PSNR: ' num2str(PSNR(ii))]);
end