  function [cost, fit, reg] = de_pl_cost(xs, G, data, R,  mask)
%|function [cost, fit, reg] = de_pl_cost(xs, G, data, R, mask)
%|
%| compute penalized-likelihood cost function for each column of x
%| in
%|	xs	[np,niter]	iterates
%|	G	[nd,np]		system matrix
%|	data	{cell}		whatever data is needed for the likelihood
%|				data{1} must have size "[nb,na]", so for 3D
%|				one must use reshaper(yi, '2d').
%|	R	penalty objects for material images(see Robject.m)
%|				
%| mask
%|  ftab  beam hardening function
%|
%| out
%|	cost	[niter,1]	cost
%|
%|
%| Yong Long, 2011-10-25

if nargin < 5, help(mfilename), error(mfilename), end

if ~isvar('R'), R = []; end

[yi  I0 ftab] = deal(data{:});
[nb na] = size(yi);
LL = ftab.LL + 1;

if isvar('mask') & ~isempty(mask)
	xs = reshapee(xs, prod(size(mask)), LL, []);	% [(*N), LL, niter]
	xs = xs(mask(:), : ,:);				% [np,LL,niter]
    np = size(xs, 1);
end

niter = size(xs,3);
reg = zeros(niter,1);

if 1 == size(yi,2)
	yi = reshapee(yi, [], ftab.fit.MM);
end

if 1 == size(I0,2)
	I0 = reshapee(I0, [], ftab.fit.MM);
end


if isempty(R)
	warning 'empty R means no penalty'
elseif issparse(R) | isa(R, 'Fatrix')
	C = R; % trick!
	if size(C,2) == size(C,1)
		warning 'square C is quite unusual!?'
	end
	for ii=1:niter
		reg(ii) = sum(abs(C * xs(:,ii)).^2)/2;
	end

elseif isstruct(R) | isa(R, 'strum')
	for ii = 1 : niter
        reg(ii) = 0;
        for ll = 1 : LL 
            reg(ii) = reg(ii) + R.penal(R, xs(:,ll,ii));
        end
    end
elseif iscell(R)
    for ii = 1 : niter
        reg(ii) = 0;
        for ll = 1 : LL
            reg(ii) = reg(ii) + R{ll}.penal(R{ll}, xs(:,ll,ii));
        end
    end   
else
	keyboard
	error 'bad R'
end

fit = zeros(niter,1);
for ii=1:niter
    for ll = 1 : LL - 1
        sil(:, ll) = G * xs(:,ll,ii);
    end
    fim = ftab.fit.fmfun(sil);
    fim = reshape(fim, size(I0));          
    yibar = I0 .* exp( - fim);
    fit(ii) = sum(col(yibar - yi .* log(yibar)));
end
       
cost = fit + reg;
cost = reale(cost, 'warn'); % trick: x'*x is not always real for complex values
