function [mus,cond_NUFFT,cond_SPIRiT,spfraction] = L1SPIRiT_ADMM_select_parameter(K_NUFFT,K_SPIRiT,spfraction_max,NUFFTgram_cs_spectrum,im_init,SPIRiTop,lambda,wavWeight,tvWeight,wavelet_args,tv_isotropic,mus)
%
% aims to find mus such that cond_NUFFT <= K_NUFFT and cond_SPIRiT <=
% K_SPIRiT and wavWeight/mus(2) and tvWeight/mus(2) <= threshold_max
%
% note that if lambda <= 0, data consistency is assumed, and mus(3) should
% be specified as an input
% 
% Copyright 2013, Daniel Weller, University of Michigan. All rights
% reserved. Email wellerds@umich.edu.

if K_NUFFT-1 < eps || K_SPIRiT-1 < eps, fail('target condition numbers should be > 1!'); end
if ~exist('mus','var'), mus = []; end

mus = col(mus);

% 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 lambda > 0
    cond_SPIRiT_fun = @(x) (SPIRiTmaxsv*lambda + x(1))/(SPIRiTminsv*lambda + x(1)); % depends only on mus(1)
    use_mu3 = false;
else
    cond_SPIRiT_fun = @(x) (SPIRiTmaxsv + x(1))/(SPIRiTminsv + x(1)); % depends only on mus(1)
    use_mu3 = true;
end

if length(mus) < 1 || mus(1) == 0 % initial guess for mus(1) based on cond_SPIRiT_fun
    mus(1) = exp(fzero(@(x) cond_SPIRiT_fun([exp(x),0]) - K_SPIRiT, 0));
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
        cond_NUFFT_fun = @(x) max(col(abs(mus(3).*NUFFTgram_cs_spectrum + x(1) + x(2).*NUFFT_mu2_spectrum)))/min(col(abs(mus(3).*NUFFTgram_cs_spectrum + x(1) + x(2).*NUFFT_mu2_spectrum)));
    else
        cond_NUFFT_fun = @(x) max(col(abs(NUFFTgram_cs_spectrum + x(1) + x(2).*NUFFT_mu2_spectrum)))/min(col(abs(NUFFTgram_cs_spectrum + x(1) + x(2).*NUFFT_mu2_spectrum)));
    end
    spfraction_fun = @(x) sum(absvals <= 1/x(2))/length(absvals);
else
    if use_mu3
        cond_NUFFT_fun = @(x) max(col(abs(mus(3).*NUFFTgram_cs_spectrum + x(1))))/min(col(abs(mus(3).*NUFFTgram_cs_spectrum + x(1))));
    else
        cond_NUFFT_fun = @(x) max(col(abs(NUFFTgram_cs_spectrum + x(1))))/min(col(abs(NUFFTgram_cs_spectrum + x(1))));
    end
    spfraction_fun = @(x) 0;
end

if use_mu2 
    if length(mus) < 2 || mus(2) == 0 % initial guess for mus(2) based on spfraction_fun
        mus(2) = exp(fzero(@(x) spfraction_fun([mus(1),exp(x)]) - spfraction_max, 0));
    end
else
    mus(2) = 0;
end

if use_mu3 && (length(mus) < 3 || mus(3) == 0)
    fail('need to specify mus(3) manually (there''s no bound on it)');
end

% so we have mus(1) that satisfies K_SPIRiT and mus(2) that satisfies
% spfraction_max; these mus form tight lower bounds on the optimal mus

function [c,ceq] = cond_funs(x)
    c = cond_NUFFT_fun(x) - K_NUFFT; % only remaining constraint we need to satisfy
    ceq = [];
end

% try to find smallest mu's that satisfy all constraints
options = optimset('Algorithm','interior-point','Display','notify-detailed','FinDiffType','central','MaxFunEvals',1000*length(mus));
[mus(1:2),~,exitflag] = fmincon(@(x) norm(x,2),mus(1:2),[],[],[],[],mus(1:2),[],@cond_funs,options);

% function [c,ceq] = cond_funs(x)
%     x = exp(x);
%     c(1) = cond_SPIRiT_fun(x) - K_SPIRiT;
%     c(2) = cond_NUFFT_fun(x) - K_NUFFT;
%     c(3) = spfraction_fun(x) - spfraction_max;
%     ceq = [];
% end

% try to find smallest mu's that satisfy this
% options = optimset('Algorithm','interior-point','Display','notify-detailed','FinDiffType','central','MaxFunEvals',1000*length(mus));
% [mus,~,exitflag] = fmincon(@(x) norm(exp(x),1),log(mus),[],[],[],[],[],[],@cond_funs,options);
% mus = exp(mus);

if exitflag < 0 % error
    warn('fmincon could not find suitable mu''s; maybe your K''s or sparse fraction are too small?');
end

cond_SPIRiT = cond_SPIRiT_fun(mus(1:2));
cond_NUFFT = cond_NUFFT_fun(mus(1:2));
spfraction = spfraction_fun(mus(1:2));

end
