function [obj, Awout, g, D] = loss_fun_and_grad(loss_func, grad_func, w, A, b, epsilon, Aw)
%Given z = (A*x), compute sum (individual_loss)
% 	if (~exist('epsilon'))
% 		epsilon = 1e-10;
% 	end
	n = length(w);
% 	m = length(b);
    idx = [];
	if (nargin < 7)
		if (nnz(w) < n / 4)
			idx = find(w);
		end

		if ~isempty(idx)
			Awout = (A(:,idx)*w(idx));
		else
			Awout = (A*w);
		end
	else
		Awout = Aw;
    end

    if (epsilon == 0)
        obj = loss_func(Awout, b);
    else
        obj = loss_func(Awout, b) + epsilon * w'*w / 2;
    end
    
	if (nargout >= 3)
		if (nargout == 4)
			[grad_coeff, D] = grad_func(Awout, b);
		else
			grad_coeff = grad_func(Awout, b);
		end
		g = A'* grad_coeff + epsilon * w;
	end
end
