function [xs,info] = L1SPIRiT_ADMM_NonCartesian(kData, x0, SPIRiTop, NUFFTop, lambda, wavWeight, tvWeight, mus, varargin)

% solve the L1/TV regularized non-Cartesian SPIRiT reconstruction using
% ADMM:
%
% min lambda/2*||(G-I)x||^2 + 1/2*||NUFFTop*x-kData||^2
%  x    + wavWeight*||DWTop*x||_1,2 + tvWeight*TV(x)
%
% becomes
%
%     min     lambda/2*||(G-I)z||^2 + 1/2*||NUFFTop*x-kData||^2
% x,w,gx,gy,z   + wavWeight*||w||_1,2 + tvWeight*||[gx,gy]||_1,2
%
% with z = x, w = DWTop*x, gx = Dx*x, gy = Dy*x
%
% with ADMM, this becomes
% 
% x <- min ||NUFFTop*x-kData||^2 + mu(1)*||x - (z - u_z)||^2
%       x                             
%                    ||[DWTop]*x-[ w -  u_w]||^2
%            + mu(2)*||[  Dx ]   [gx - u_gx]||
%                    ||[  Dy ]   [gy - u_gy]||
%
% z <- min lambda*||(G-I)z||^2 + mu(1)*||z - (x + u_z)||^2
%       z
%
% w <- min wavWeight*||w||_1,2 + mu(2)/2*||w - (DWTop*x + u_w)||^2
%       w
%
% [gx,gy] <- min tvWeight*||[gx,gy]||_1,2 + mu(2)/2*||[gx] - [Dx*x + u_gx]||^2
%           gx,gy                                   ||[gy]   [Dy*x + u_gy]||
%
% u_z <- u_z + (x - z)
% u_w <- u_w + (DWTop*x - w)
% u_gx <- u_gx + (Dx*x - gx)
% u_gy <- u_gy + (Dy*x - gy)
%
% (we do the x-update last)
%
% Added code to interpret lambda = 0 as data-preserving case:
%
% min 1/2*||(G-I)x||^2
%  x    + wavWeight*||DWTop*x||_1,2 + tvWeight*TV(x)
%
% with kData = NUFFTop*x
%
% becomes
%
%     min     1/2*||(G-I)z||^2
% x,w,gx,gy,z   + wavWeight*||w||_1,2 + tvWeight*||[gx,gy]||_1,2
%
% with z = x, w = DWTop*x, gx = Dx*x, gy = Dy*x, kData = NUFFTop*x
%
% with ADMM, this becomes
% 
% x <- min mu(3)*||NUFFTop*x - (kData-u_d)||^2 + mu(1)*||x - (z - u_z)||^2
%       x                             
%                    ||[DWTop]*x-[ w -  u_w]||^2
%            + mu(2)*||[  Dx ]   [gx - u_gx]||
%                    ||[  Dy ]   [gy - u_gy]||
%
% z <- min ||(G-I)z||^2 + mu(1)*||z - (x + u_z)||^2
%       z
%
% w <- min wavWeight*||w||_1,2 + mu(2)/2*||w - (DWTop*x + u_w)||^2
%       w
%
% [gx,gy] <- min tvWeight*||[gx,gy]||_1,2 + mu(2)/2*||[gx] - [Dx*x + u_gx]||^2
%           gx,gy                                   ||[gy]   [Dy*x + u_gy]||
%
% u_d <- u_d + (NUFFTop*x - kData)
% u_z <- u_z + (x - z)
% u_w <- u_w + (DWTop*x - w)
% u_gx <- u_gx + (Dx*x - gx)
% u_gy <- u_gy + (Dy*x - gy)
% 
% Copyright 2013, Daniel Weller, University of Michigan. All rights
% reserved. Email wellerds@umich.edu.

opts.maxiters = 10;
opts.miniters = 1;
opts.itersave = []; % if not empty, iterations to save current output in "xs"
opts.tsave = []; % if not empty, times to save current output in "xs"
opts.abs_tol = 1e-4;
opts.rel_tol = 1e-4;
opts.xtol = 1e-4;
opts.update_x_solve = 'cg'; % also accepts 'exact' (for small # images)
opts.update_x_cg_opts = {}; % see options in pre_conj_grad.m
opts.update_x_precon = []; % (2-D) spectrum of circulant preconditioner to use for x-update
opts.wavelet_args = {'Daubechies',4,4}; % wavelet type, type spec, and # of levels
opts.tv_isotropic = true; % toggles isotropic/anisotropic TV
opts.update_z_solve = 'cg'; % also accepts 'exact' (for small # of coils)
opts.update_z_cg_opts = {}; % see options in pre_conj_grad.m
opts.update_z_precon = ''; % also accepts 'diag' and 'superdiag'
opts.output = 0; % 0 = quiet, 1 = summary + warnings only, 2 = iterations, 3 = display image
opts.compute_mus = ''; % if 'conds', uses L1SPIRiT_ADMM_select_parameter() to compute mus
opts.compute_mus_args = {}; % cell array of user arguments (e.g. condition numbers/threshold and NUFFTgram preconditioner spectrum for L1SPIRiT_ADMM_select_parameter())
opts.test_convergence = false; % set to true to compute r's and s's
opts.evalfun = []; % function out = evalfun(x) to call at each iteration

opts = vararg_pair(opts,varargin);

te = tic();
tinit = tic();

x = x0;
sz_images = size(x); sz_images(end+1:3) = 1;
sz_data = size(kData); sz_data(end+1:3) = 1;

u_d = 0;
if lambda <= 0
    NUFFT_x = NUFFTop*x;
end

z = x; u_z = 0;
if wavWeight > 0
    sz_impad = 2^nextpow2(max(sz_images(1:2)));
    sz_impad = [sz_impad,sz_impad,sz_images(3:end)];
    
    DWTop = Wavelet(opts.wavelet_args{:});
    DWT_x = DWTop*zpad(x,sz_impad);
    w = DWT_x;
    u_w = 0;
else
    w = []; u_w = []; DWTop = []; sz_impad = sz_images;
end
if tvWeight > 0
    TVop = TVOP(sz_images(1:2));
    TV_x = TVop*x;
    gxy = TV_x;
    u_gxy = 0;
else
    gxy = []; u_gxy = []; TVop = [];
end

res_z = []; res_x = []; % cg residuals

% compute condition number if needed
switch lower(opts.compute_mus)
    case '' % use specified mu
        info.mus = mus;
    case 'conds'
        if opts.output >= 2, fprintf(1,'Condition #-based ADMM tuning parameter selection: '); drawnow; end
        [mus,cond_NUFFT,cond_SPIRiT,spfraction] = L1SPIRiT_ADMM_select_parameter(opts.compute_mus_args{1:4},x0,SPIRiTop,lambda,wavWeight,tvWeight,opts.wavelet_args,opts.tv_isotropic,mus);
        if opts.output >= 2, fprintf(1,'K_x = %g, K_z = %g, spfraction = %g => mu_z = %g, mu_W = %g.\n',cond_NUFFT,cond_SPIRiT,spfraction,mus(1),mus(2)); drawnow; end
        info.mus = mus;
        info.cond_NUFFT = cond_NUFFT;
        info.cond_SPIRIT = cond_SPIRiT;
        info.spfraction = spfraction;
    otherwise
        warn('unknown mu computation method, ignoring...');
end

if (lambda > 0 && length(mus) < 2) || (lambda <= 0 && length(mus) < 3)
    fail('mu needs to be specified');
end

if lambda > 0, mus(3) = 1; end

% other initialization
if strcmpi(opts.update_z_solve,'exact')
    z_solve_helper = permute(getK(SPIRiTop),[4,3,1,2]); % [Nc Nc N1 N2]
    z_solve_helper = col(num2cell(z_solve_helper,[1,2]));
    z_solve_helper = cellfun(@(x) x-eye(sz_images(3)),z_solve_helper,'UniformOutput',false);
    z_solve_helper = cellfun(@(x) x'*x,z_solve_helper,'UniformOutput',false);
    if lambda > 0
        z_solve_helper = cellfun(@(x) lambda.*x + mus(1).*eye(sz_images(3)),z_solve_helper,'UniformOutput',false);
    else
        z_solve_helper = cellfun(@(x) x + mus(1).*eye(sz_images(3)),z_solve_helper,'UniformOutput',false);
    end
    z_solve_helper = cellfun(@(x) inv(x),z_solve_helper,'UniformOutput',false);
    z_solve_helper = reshape(permute(cat(3,z_solve_helper{:}),[3,2,1]),sz_images([1:3,3]));
else
    switch lower(opts.update_z_precon)
        case 'diag'
            % used for computing diagonal preconditioner
            z_solve_helper = permute(getK(SPIRiTop),[4,3,1,2]); % [Nc Nc N1 N2]
            z_solve_helper = col(num2cell(z_solve_helper,[1,2]));
            z_solve_helper = cellfun(@(x) x-eye(sz_images(3)),z_solve_helper,'UniformOutput',false);
            z_solve_helper = cellfun(@(x) sum(abs(x).^2,1),z_solve_helper,'UniformOutput',false); % [N1*N2 Nc]
            if lambda > 0
                z_solve_helper = lambda.*col(cat(1,z_solve_helper{:})) + mus(1); % [N1*N2*Nc]
            else
                z_solve_helper = col(cat(1,z_solve_helper{:})) + mus(1); % [N1*N2*Nc]
            end
            z_solve_helper = 1./z_solve_helper;
        case 'superdiag'
            % used for computing superdiagonal preconditioner
            z_solve_helper = permute(getK(SPIRiTop),[4,3,1,2]); % [Nc Nc N1 N2]
            z_solve_helper = col(num2cell(z_solve_helper,[1,2]));
            z_solve_helper = cellfun(@(x) x-eye(sz_images(3)),z_solve_helper,'UniformOutput',false);
            z_solve_helper = cellfun(@(x) x'*x,z_solve_helper,'UniformOutput',false);
            if lambda > 0
                z_solve_helper = cellfun(@(x) lambda.*x + mus(1).*eye(sz_images(3)),z_solve_helper,'UniformOutput',false);
            else
                z_solve_helper = cellfun(@(x) x + mus(1).*eye(sz_images(3)),z_solve_helper,'UniformOutput',false);
            end
            z_solve_helper = cellfun(@(x) diag(abs(x))./sum(abs(x).^2,2),z_solve_helper,'UniformOutput',false);
            z_solve_helper = col(cat(2,z_solve_helper{:}).');
    end
end

if strcmpi(opts.update_x_solve,'exact')
    x_solve_helper1 = 0; % spectrum of circulant W = mu*(I + DWT'*DWT + TV'*TV)
    if wavWeight > 0
        x_solve_helper1 = x_solve_helper1 + 1; % assume DWT is orthonormal
%         x_solve_helper1 = x_solve_helper1 + abs(getFFTs(DWTop)).^2;
    end
    if tvWeight > 0
        [DxFFT,DyFFT] = getFFTs(TVop);
        x_solve_helper1 = x_solve_helper1 + repmat(abs(DxFFT).^2,[1,sz_images(2)]) + repmat(abs(DyFFT).^2,[sz_images(1),1]);
    end
    x_solve_helper1 = mus(1) + mus(2).*x_solve_helper1;
    
    x_solve_helper2 = reshape(full(NUFFTop)',[sz_images(1),sz_images(2),prod(sz_data(1:2))]); % D' has dimensions [*N M]
    x_solve_helper2 = reshape(fft(fft(x_solve_helper2,[],1),[],2),[prod(sz_images(1:2)),prod(sz_data(1:2))]); % FFT(D') has dimensions [*N M]
    x_solve_helper2 = (prod(sz_images(1:2)).\x_solve_helper2')*(x_solve_helper2./repmat(x_solve_helper1(:),[1,prod(sz_data(1:2))])); % D*inv(W)*D' = (FFT(D'))'*(FFT(D')./W_spectrum) has dimensions [M M]
    x_solve_helper2 = x_solve_helper2 + eye(size(x_solve_helper2,1)); % [M M] matrix (I + D*W^-1*D')
else
    if ~isempty(opts.update_x_precon) % circulant
        x_solve_helper1 = mus(3).*opts.update_x_precon + mus(1);
        if wavWeight > 0
            x_solve_helper1 = x_solve_helper1 + mus(2); % assume orthonormal
%             x_solve_helper1 = x_solve_helper1 + mus(2).*abs(getFFTs(DWTop)).^2;
        end
        if tvWeight > 0
            [DxFFT,DyFFT] = getFFTs(TVop);
            x_solve_helper1 = x_solve_helper1 + mus(2).*(bsxfun(@plus,abs(DxFFT).^2,abs(DyFFT).^2));
        end
        x_solve_helper1 = 1./x_solve_helper1;
    end
end

tinit = toc(tinit);

% compute objective function
obj = ffun();

if opts.output == 3 && usejava('awt')
    h = figure;
    imshow(sos(x),[]); title('Initial guess'); drawnow;
else
    h = [];
end

% output storage
if opts.test_convergence
    r1s = zeros(opts.maxiters,1);
    s1s = zeros(opts.maxiters,1);
    r2s = zeros(opts.maxiters,1);
    s2s = zeros(opts.maxiters,1);
    r3s = zeros(opts.maxiters,1);
    converged = zeros(opts.maxiters,1);
end
normxdiffs = zeros(opts.maxiters,1);
objs = [obj;zeros(opts.maxiters,1)];
if strcmpi(opts.update_x_solve,'cg'), x_cg_info = cell(opts.maxiters,1); end
if strcmpi(opts.update_z_solve,'cg'), z_cg_info = cell(opts.maxiters,1); end
if ~isempty(opts.evalfun) && isa(opts.evalfun,'function_handle')
    evaloutputs = opts.evalfun(x);
end
ts = zeros(opts.maxiters,1);

xs = struct('x',x,'t',0,'iter',0);

iter = 1;
while iter <= opts.maxiters
    if opts.output == 1 % just update completion
        ticker(mfilename, iter, opts.maxiters);
    end
    
    titer = tic();
    
    % update z
    z_prev = z;
    switch lower(opts.update_z_solve)
        case 'cg'
            b_z = update_z_cg_AHbfun();
            if ~isempty(res_z), res_z = res_z - b_z; end
            if isempty(opts.update_z_precon)
                z_cg_opts = {'resvec',res_z};
            else
                z_cg_opts = {'M',@update_z_cg_Mfun,'resvec',res_z};
            end
            [z,z_cg_info{iter},res_z] = pre_conj_grad(col(z),@update_z_cg_AHAfun,b_z,opts.update_z_cg_opts{:},z_cg_opts{:});
            z = reshape(z,sz_images);
            res_z = res_z + b_z;
        case 'exact'
            for n=1:sz_images(3)
                z(:,:,n) = sum(z_solve_helper(:,:,:,n).*(mus(1).*(x + u_z)),3); 
            end
            
        otherwise
            fail('unknown z-solver!');
    end
    
    % update w
    if wavWeight > 0
        w_prev = w;
        w = DWT_x + u_w;
        wabs = sqrt(sum(abs(w).^2,3));
        wabs = 1 - (wavWeight/mus(2))./max(wabs,wavWeight/mus(2));
        wabs(~isfinite(wabs)) = 0; % safeguard for when wabs_n = wavWeight/mu = 0
        w = w.*repmat(wabs,[1,1,size(w,3)]);
    end
    
    % update gxy
    if tvWeight > 0
        gxy_prev = gxy;
        gxy = TV_x + u_gxy;
        if opts.tv_isotropic
            gxy_abs = sqrt(sum(sum(abs(gxy).^2,4),3));
        else
            gxy_abs = sqrt(sum(abs(gxy).^2,3));
        end
        gxy_abs = 1 - (tvWeight/mus(2))./max(gxy_abs,tvWeight/mus(2));
        gxy_abs(~isfinite(gxy_abs)) = 0; % safeguard for when wabs_n = tvWeight/mu = 0
        if opts.tv_isotropic
            gxy = gxy.*repmat(gxy_abs,[1,1,size(gxy,3),size(gxy,4)]);
        else
            gxy = gxy.*repmat(gxy_abs,[1,1,size(gxy,3),1]);
        end
    end
    
    % update u's, r's, s's, thresholds
    if opts.test_convergence
        A1 = x; B1 = -z; r1 = A1 + B1; u_z = u_z + r1;
        A1 = norm(col(A1)); B1 = norm(col(B1)); r1 = norm(col(r1));
        s1 = mus(1)*norm(col(z-z_prev));
        n_r1 = numel(u_z);
        Au1 = norm(col(u_z));
        
        if wavWeight > 0
            A_w = DWT_x; B_w = -w; r_w = A_w + B_w; u_w = u_w + r_w;
            A_w = norm(col(A_w)); B_w = norm(col(B_w)); r_w = norm(col(r_w));
            s_w = mus(2)*norm(col(crop(DWTop'*(w-w_prev),sz_images)));
            n_r_w = numel(u_w);
            Au_w = norm(col(crop(DWTop'*u_w,sz_images)));
            A2 = A_w; B2 = B_w; r2 = r_w; s2 = s_w; n_r2 = n_r_w; Au2 = Au_w; 
        else
            A2 = []; B2 = []; r2 = []; s2 = []; n_r2 = 0; Au2 = [];
        end
        
        if tvWeight > 0
            A_gxy = TV_x; B_gxy = -gxy; r_gxy = A_gxy + B_gxy; u_gxy = u_gxy + r_gxy;
            A_gxy = norm(col(A_gxy)); B_gxy = norm(col(B_gxy)); r_gxy = norm(col(r_gxy));
            s_gxy = mus(2)*norm(col(TVop'*(gxy-gxy_prev)));
            n_r_gxy = numel(u_gxy);
            Au_gxy = norm(col(TVop'*u_gxy));
            A2 = [A2,A_gxy]; B2 = [B2,B_gxy]; r2 = [r2,r_gxy]; s2 = [s2,s_gxy]; n_r2 = n_r2 + n_r_gxy; Au2 = [Au2,Au_gxy]; %#ok<AGROW>
        end
        
        if lambda > 0
            A3 = []; B3 = []; r3 = []; n_r3 = 0;
        else
            A3 = NUFFT_x; B3 = -kData; r3 = A3 + B3; u_d = u_d + r3;
            A3 = norm(col(A3)); B3 = norm(col(B3)); r3 = norm(col(r3));
            n_r3 = numel(u_d);
        end
        
        A2 = norm(A2); B2 = norm(B2); r2 = norm(r2); s2 = norm(s2); Au2 = norm(Au2);
        A3 = norm(A3); B3 = norm(B3); r3 = norm(r3);
        r1s(iter) = r1; s1s(iter) = s1; r2s(iter) = r2; s2s(iter) = s2; r3s(iter) = r3;
    else
        if lambda <= 0, u_d = u_d + (NUFFT_x - kData); end
        u_z = u_z + (x - z);
        if wavWeight > 0, u_w = u_w + (DWT_x - w); end
        if tvWeight > 0, u_gxy = u_gxy + (TV_x - gxy); end
    end
    
    % update x
    x_prev = x;
    switch lower(opts.update_x_solve)
        case 'cg'
            b_x = update_x_cg_AHbfun();
            if ~isempty(res_x), res_x = res_x - b_x; end
            if isempty(opts.update_x_precon)
                x_cg_opts = {'resvec',res_x};
            else
                x_cg_opts = {'M',@update_x_cg_Mfun,'resvec',res_x};
            end
            [x,x_cg_info{iter},res_x] = pre_conj_grad(col(x),@update_x_cg_AHAfun,b_x,opts.update_x_cg_opts{:},x_cg_opts{:});
            x = reshape(x,sz_images);
            res_x = res_x + b_x;
        case 'exact'
            x = mus(3).*(NUFFTop'*(kData - u_d)) + mus(1).*(z - u_z);
            if wavWeight > 0
                x = x + mus(2).*crop(DWTop'*(w - u_w),sz_images);
            end
            if tvWeight > 0
                x = x + mus(2).*(TVop'*(gxy - u_gxy));
            end
            x = mus(3).\ifft(ifft(repmat(x_solve_helper1,[1,1,sz_images(3:end)]).\fft(fft(x,[],1),[],2),[],2),[],1);
            x = x - mus(3).\ifft(ifft(repmat(x_solve_helper1,[1,1,sz_images(3:end)]).\fft(fft(...
                NUFFTop'*reshape(x_solve_helper2\reshape(NUFFTop*x,[prod(sz_data(1:2)),sz_images(3:end)]),sz_data)...
                ,[],1),[],2),[],2),[],1);
        otherwise
            fail('unknown x-solver!');
    end
    normxdiffs(iter) = norm(col(x-x_prev));
    
    % update intermediate variables
    if lambda <= 0
        NUFFT_x = NUFFTop*x;
    end
    if wavWeight > 0
        DWT_x = DWTop*zpad(x,sz_impad);
    end
    if tvWeight > 0
        TV_x = TVop*x;
    end

    if opts.test_convergence
        % residual tolerances from Boyd
        tol_r1 = sqrt(n_r1)*opts.abs_tol + max([A1,B1])*opts.rel_tol;
        tol_s1 = sqrt(prod(sz_images))*opts.abs_tol + mus(1)*Au1*opts.rel_tol;
        tol_r2 = sqrt(n_r2)*opts.abs_tol + max([A2,B2])*opts.rel_tol;
        tol_s2 = sqrt(prod(sz_images))*opts.abs_tol + mus(2)*Au2*opts.rel_tol;
        tol_r3 = sqrt(n_r3)*opts.abs_tol + max([A3,B3])*opts.rel_tol;

        if r1 <= tol_r1 && s1 <= tol_s1 && (~(wavWeight > 0 || tvWeight > 0) || (r2 <= tol_r2 && s2 <= tol_s2)) && (lambda > 0 || r3 <= tol_r3)
            converged(iter) = 1; % primal and dual residual convergence

            if normxdiffs(iter) <= opts.xtol*norm(col(x_prev))
                converged(iter) = 2; % x-convergence as well
            end
        end
    end
    
    ts(iter) = toc(titer);

    % compute updated objective
    obj = ffun();
    objs(iter+1) = obj;
    
    if ~isempty(opts.evalfun) && isa(opts.evalfun,'function_handle')
        evaloutputs(iter+1) = opts.evalfun(x);
    end
    
    if any(opts.itersave == iter) || any(opts.tsave > sum(ts(1:iter-1)) & opts.tsave <= sum(ts(1:iter)))
        xs(end+1) = struct('x',x,'t',sum(ts(1:iter)),'iter',iter); %#ok<AGROW>
    end
    
    % output
    if opts.output >= 2
        fprintf(1,'Iter #%d: ',iter);
        if opts.test_convergence
            if r1 <= tol_r1, str_r1 = '<='; else str_r1 = '>'; end
            if s1 <= tol_s1, str_s1 = '<='; else str_s1 = '>'; end
            fprintf(1,'r1 = %g %s %g, s1 = %g %s %g, ',r1,str_r1,tol_r1,s1,str_s1,tol_s1);
            if wavWeight > 0 || tvWeight > 0
                if r2 <= tol_r2, str_r2 = '<='; else str_r2 = '>'; end
                if s2 <= tol_s2, str_s2 = '<='; else str_s2 = '>'; end
                fprintf(1,'r2 = %g %s %g, s2 = %g %s %g, ',r2,str_r2,tol_r2,s2,str_s2,tol_s2);
            end
            if lambda <= 0
                if r3 <= tol_r3, str_r3 = '<='; else str_r3 = '>'; end
                fprintf(1,'r3 = %g %s %g, ',r3,str_r3,tol_r3);
            end
        end
        if normxdiffs(iter) <= opts.xtol*norm(col(x_prev)), str_x = '<='; else str_x = '>'; end
        fprintf(1,'dx = %g %s %g',normxdiffs(iter),str_x,opts.xtol*norm(col(x_prev)));

        if opts.test_convergence
            switch converged(iter)
                case 1
                    fprintf(1,' (residual convergence).\n');
                case 2
                    fprintf(1,' (x convergence).\n');
                otherwise
                    fprintf(1,'.\n');
            end
        else
            fprintf(1,'.\n');
        end
        
        if ~isempty(h) && ishandle(h)
            set(0,'CurrentFigure',h);
            imshow(sos(x),[]); title(sprintf('Iter #%d',iter));
        end
        drawnow;
    end
    
    if iter >= opts.miniters && opts.test_convergence && converged(iter) == 2
        break;
    end
    
    iter = iter + 1;
end

te = toc(te);
if iter > opts.maxiters, iter = opts.maxiters; end
info.te = te;
info.niters = iter;
if opts.test_convergence
    info.r1s = r1s(1:iter);
    info.s1s = s1s(1:iter);
    info.r2s = r2s(1:iter);
    info.s2s = s2s(1:iter);
    info.r3s = r3s(1:iter);
    info.converged = converged(1:iter);
end
info.objs = objs(1:iter+1);
info.normxdiffs = normxdiffs(1:iter);
if strcmpi(opts.update_x_solve,'cg'), info.x_cg_info = x_cg_info(1:iter); end
if strcmpi(opts.update_z_solve,'cg'), info.z_cg_info = z_cg_info(1:iter); end
if ~isempty(h) && ishandle(h), info.hfigure = h; end % handle of figure containing sos-image
if ~isempty(opts.evalfun) && isa(opts.evalfun,'function_handle')
    info.evaloutputs = evaloutputs(1:iter+1);
end
info.tinit = tinit;
info.times = ts(1:iter);
xs(1).x = x; xs(1).t = sum(ts(1:iter)); xs(1).iter = iter;

if opts.output == 1
    fprintf(1,'Iter #%d: ',iter);
    if opts.test_convergence
        if r1 <= tol_r1, str_r1 = '<='; else str_r1 = '>'; end
        if s1 <= tol_s1, str_s1 = '<='; else str_s1 = '>'; end
        fprintf(1,'r1 = %g %s %g, s1 = %g %s %g, ',r1,str_r1,tol_r1,s1,str_s1,tol_s1);
        if wavWeight > 0 || tvWeight > 0
            if r2 <= tol_r2, str_r2 = '<='; else str_r2 = '>'; end
            if s2 <= tol_s2, str_s2 = '<='; else str_s2 = '>'; end
            fprintf(1,'r2 = %g %s %g, s2 = %g %s %g, ',r2,str_r2,tol_r2,s2,str_s2,tol_s2);
        end
        if lambda <= 0
            if r3 <= tol_r3, str_r3 = '<='; else str_r3 = '>'; end
            fprintf(1,'r3 = %g %s %g, ',r3,str_r3,tol_r3);
        end
    end
    if normxdiffs(iter) <= opts.xtol*norm(col(x_prev)), str_x = '<='; else str_x = '>'; end
    fprintf(1,'dx = %g %s %g',normxdiffs(iter),str_x,opts.xtol*norm(col(x_prev)));

    if opts.test_convergence
        switch converged(iter)
            case 1
                fprintf(1,' (residual convergence).\n');
            case 2
                fprintf(1,' (x convergence).\n');
            otherwise
                fprintf(1,'.\n');
        end
    else
        fprintf(1,'.\n');
    end
    drawnow;
end

% nested helper functions (for CG and the like)
    function f = ffun()
        if lambda > 0
            f = 0.5*lambda*norm(col(SPIRiTop*x))^2 + 0.5*norm(col(NUFFTop*x - kData))^2;
        else
            f = 0.5*norm(col(SPIRiTop*x))^2;
        end
        if wavWeight > 0
            f = f + wavWeight*sum(sum(sqrt(sum(abs(DWTop*zpad(x,sz_impad)).^2,3)),2),1);
        end
        if tvWeight > 0
            if opts.tv_isotropic
                f = f + tvWeight*sum(sum(sqrt(sum(sum(abs(TVop*x).^2,4),3)),2),1);
            else
                f = f + tvWeight*sum(sum(sum(sqrt(sum(abs(TVop*x).^2,3)),4),2),1);
            end
        end
    end

    function out = update_z_cg_AHAfun(in)
        % out <- (lambda*(G-I)'*(G-I) + mu(1)*I)*in
        
        in = reshape(in,sz_images);
        if lambda > 0
            out = lambda.*(SPIRiTop'*(SPIRiTop*in)) + mus(1).*in;
        else
            out = SPIRiTop'*(SPIRiTop*in) + mus(1).*in;
        end
        out = out(:);
        drawnow;
        
    end

    function out = update_z_cg_AHbfun()
        % out <- mu(1)*(x + u_z)
        
        out = mus(1).*col(x + u_z);
        
    end

    function out = update_z_cg_Mfun(in)
        % out <- M.\in
        
        if ~isempty(z_solve_helper), out = in.*z_solve_helper; else out = in; end
    end

    function out = update_x_cg_AHAfun(in)
        % out <- (mu(3)*NUFFTop'*NUFFTop + mu(1)*(I + DWT'*DWT + TV'*TV))*in
        
        in = reshape(in,sz_images);
        out = mus(3).*(NUFFTop'*(NUFFTop*in)) + mus(1).*in;
        if wavWeight > 0
            out = out + mus(2).*in; % assume DWT is orthonormal
%             out = out + mu.*crop(DWTop'*(DWTop*zpad(x,sz_impad)),sz_images);
        end
        if tvWeight > 0
            out = out + mus(2).*(TVop'*(TVop*in));
        end
        out = out(:);
        drawnow;
        
    end

    function out = update_x_cg_AHbfun()
        % out <- mu(3)*NUFFTop'*(kData-u_d) + mu(1)*(z-u_z + DWT'*(w-u_w) + TV'*(gxy-u_gxy))
        
        out = mus(3).*(NUFFTop'*(kData-u_d)) + mus(1).*(z - u_z);
        if wavWeight > 0
            out = out + mus(2).*crop(DWTop'*(w - u_w),sz_images);
        end
        if tvWeight > 0
            out = out + mus(2).*(TVop'*(gxy-u_gxy));
        end
        out = out(:);
        
    end

    function out = update_x_cg_Mfun(in)
        % out <- ifft(precon_spectrum.\fft(in))
        % where precon_spectrum = spectrum of (mu(3)*D'D + mu(1)*(I + DWT'*DWT + TV'*TV)
        
        in = reshape(in,sz_images);
        if iscell(opts.update_x_precon)
            out = sum(bsxfun(@times,x_solve_helper1,fft(fft(bsxfun(@times,x_solve_helper2,in),[],1),[],2)),4);
            out = sum(bsxfun(@times,x_solve_helper2,ifft(ifft(bsxfun(@times,x_solve_helper1,out),[],2),[],1)),4);
        else
            out = ifft(ifft(bsxfun(@times,x_solve_helper1,fft(fft(in,[],1),[],2)),[],2),[],1);
        end

        out = out(:);
        drawnow;
        
    end

end
