function [cond_NUFFT,cond_SPIRiT,spfraction] = L1SPIRiT_ADMM_examine_parameters(NUFFTgram_cs_spectrum,im_init,SPIRiTop,lambda,wavWeight,tvWeight,wavelet_args,tv_isotropic,mu1s,mu2s,mu3s)
% 
% Copyright 2013, Daniel Weller, University of Michigan. All rights
% reserved. Email wellerds@umich.edu.

% get output size
sz_mu1s = size(mu1s);
sz_mu2s = size(mu2s);
sz_mu3s = size(mu3s);

use_mu2 = wavWeight > 0 || tvWeight > 0;
use_mu3 = lambda <= 0;

if use_mu3 && ~isscalar(mu3s)
    sz_out = sz_mu3s;
elseif use_mu2 && ~isscalar(mu2s)
    sz_out = sz_mu2s;
else
    sz_out = sz_mu1s;
end
if use_mu3 && isscalar(mu3s)
    mu3s = repmat(mu3s,sz_out);
end
if use_mu2 && isscalar(mu2s)
    mu2s = repmat(mu2s,sz_out);
end
if isscalar(mu1s)
    mu1s = repmat(mu1s,sz_out);
end

cond_SPIRiT = zeros(sz_out);
cond_NUFFT = zeros(sz_out);
spfraction = zeros(sz_out);

% condition number for SPIRiT equation (lambda*(G-I)'*(G-I) + mu(1)*I)
SPIRiTkernels = permute(getK(SPIRiTop),[4,3,1,2]); % [Nc Nc N1 N2]
nCoils = size(SPIRiTkernels,1);
SPIRiTkernels = col(num2cell(SPIRiTkernels,[1,2]));
SPIRiTkernels = cellfun(@(x) x-eye(nCoils),SPIRiTkernels,'UniformOutput',false);
SPIRiTkernels = cellfun(@(x) x'*x,SPIRiTkernels,'UniformOutput',false);
SPIRiTsvs = cellfun(@(x) svd(x),SPIRiTkernels,'UniformOutput',false);
SPIRiTmaxsv = max(col(cellfun(@(x) max(col(abs(x))),SPIRiTsvs)));
SPIRiTminsv = min(col(cellfun(@(x) min(col(abs(x))),SPIRiTsvs)));
if use_mu3
    for imu = 1:prod(sz_out)
        cond_SPIRiT(imu) = (SPIRiTmaxsv + mu1s(imu))/(SPIRiTminsv + mu1s(imu)); % depends only on mus(1)
    end
else
    for imu = 1:prod(sz_out)
        cond_SPIRiT(imu) = (SPIRiTmaxsv*lambda + mu1s(imu))/(SPIRiTminsv*lambda + mu1s(imu)); % depends only on mus(1)
    end
end

% condition number for NUFFT equation => mus(1:2)
%
% singular values = NUFFTgram_cs_spectrum + mu(1) + mu(2)*(TVfft + DWTffts)
%
sz_image = size(NUFFTgram_cs_spectrum);
sz_images = size(im_init); 

% assemble singular values "matrix"
if wavWeight > 0
    sz_impad = 2^nextpow2(max(sz_images(1:2)));
    sz_impad = [sz_impad,sz_impad,sz_images(3:end)];
    DWTop = Wavelet(wavelet_args{:});
    
    NUFFT_mu2_spectrum = 1; % assume orthogonal DWT
    DWT_abs = DWTop*zpad(im_init,sz_impad);
    DWT_abs = sqrt(sum(abs(DWT_abs).^2,3));
    absvals = col(DWT_abs)./wavWeight;
    absvals = sort(absvals,'descend');
    absvals = absvals(1:size(im_init,1)*size(im_init,2));
    use_mu2 = true;
else
    NUFFT_mu2_spectrum = 0;
    absvals = [];
    use_mu2 = false;
end
if tvWeight > 0
    TVop = TVOP(sz_images(1:2));
    [DxFFT,DyFFT] = getFFTs(TVop);
    NUFFT_mu2_spectrum = NUFFT_mu2_spectrum + (repmat(abs(DxFFT).^2,[1,sz_image(2)]) + repmat(abs(DyFFT).^2,[sz_image(1),1]));
    TV_abs = TVop*im_init;
    if tv_isotropic
        TV_abs = sqrt(sum(sum(abs(TV_abs).^2,4),3));
    else
        TV_abs = sqrt(sum(abs(TV_abs).^2,3));
    end
    absvals = [absvals;col(TV_abs)./tvWeight];
    use_mu2 = true;
end

if use_mu2
    if use_mu3
        for imu = 1:prod(sz_out)
            cond_NUFFT(imu) = max(col(abs(mu3s(imu).*NUFFTgram_cs_spectrum + mu1s(imu) + mu2s(imu).*NUFFT_mu2_spectrum)))/min(col(abs(mu3s(imu).*NUFFTgram_cs_spectrum + mu1s(imu) + mu2s(imu).*NUFFT_mu2_spectrum)));
        end
    else
        for imu = 1:prod(sz_out)
            cond_NUFFT(imu) = max(col(abs(NUFFTgram_cs_spectrum + mu1s(imu) + mu2s(imu).*NUFFT_mu2_spectrum)))/min(col(abs(NUFFTgram_cs_spectrum + mu1s(imu) + mu2s(imu).*NUFFT_mu2_spectrum)));
        end
    end
    for imu = 1:prod(sz_out)
        spfraction(imu) = sum(absvals <= 1/mu2s(imu))/length(absvals);
    end
else
    if use_mu3
        for imu = 1:prod(sz_out)
            cond_NUFFT(imu) = max(col(abs(mu3s(imu).*NUFFTgram_cs_spectrum + mu1s(imu))))/min(col(abs(mu3s(imu).*NUFFTgram_cs_spectrum + mu1s(imu))));
        end
    else
        for imu = 1:prod(sz_out)
            cond_NUFFT(imu) = max(col(abs(NUFFTgram_cs_spectrum + mu1s(imu))))/min(col(abs(NUFFTgram_cs_spectrum + mu1s(imu))));
        end
    end
    for imu = 1:prod(sz_out)
        spfraction(imu) = 0;
    end
end

end
