function [WSURE,Y_L1SPIRiT] = L1SPIRiT_WSURE(k_acq,Y_SPIRiT,b,Y_SPIRiT_b,SPIRiTop,nSPIRiTiters,wavWeight,epsilon,Y_L1SPIRiT)

sz_full = size(k_acq); sz_full(end+1:3) = 1;
K = abs(k_acq) > 0; % mask of acquired data

Y_SPIRiT_mask = Y_SPIRiT;
Y_SPIRiT_mask(K) = 0;

% perform L1SPIRiT with just acquired data
if ~exist('Y_L1SPIRiT','var') || isempty(Y_L1SPIRiT)
    Y_L1SPIRiT = pocsSPIRiT(k_acq,SPIRiTop,nSPIRiTiters,k_acq,wavWeight,0);
end
Y_L1SPIRiT_mask = Y_L1SPIRiT;
Y_L1SPIRiT_mask(K) = 0;

% WSURE = ||B*f_gamma(y)||^2 - 2*real(y'*A'*B*f_gamma(y))
WSURE = abs(Y_L1SPIRiT_mask(:)'*Y_L1SPIRiT_mask(:)) - 2*real(Y_SPIRiT_mask(:)'*Y_L1SPIRiT_mask(:));

% do the same with data + complex Bernoulli noise
Y_L1SPIRiT_yb = pocsSPIRiT(k_acq+epsilon.*b,SPIRiTop,nSPIRiTiters,k_acq+epsilon.*b,wavWeight,0);

varrho = (Y_L1SPIRiT_yb - Y_L1SPIRiT)./epsilon;
varrho(K) = 0;

% mask the output (if preserving data)
Y_SPIRiT_b_mask = Y_SPIRiT_b;
Y_SPIRiT_b_mask(K) = 0;

WSURE = WSURE + 2*real(Y_SPIRiT_b_mask(:)'*varrho(:));

end
