  function [grad, curv] = de_pl_dercurv(data, sil, Gb, curvtype, coupled_material, iblock, nblock)
%|function [grad, curv] = de_pl_dercurv(data, sil, Gb, curvtype,  coupled_material, iblock, nblock)
%|
%| evaluate gradients and curvatures for negative log-likelihood function:
%| U(x) = \sum_i \sum_m -h_im(ybar_im)
%| h_im(ybar_im) = ybar_im - h_im *  log(ybar_im)
%| ybar_im = I_im * exp( - f_im(s_il) )
%|
%|
%| in
%|	data	{hi, I0, ftab}
%|				hi and I0 are [nb na]
%|
%|	sil	[nb na L]	material density line integrals
%|  Gb	[nd,np]		system matrix (see Gblock)
%|	curvtype		todo
%|	iblock, nblock		for OS type methods
%|  coupled_material   1 for PWSQS (pixel-wise separable, materials are coupled), 0 for SQS
%| out
%|	grad	[np]	\dot hi(l)
%|	curv	[np]	L-separable surrogate curvature for each U_i
%|
%|
%| Copyright 2010-10-25, Yong Long
%| 


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

[hi  I0 ftab] = deal(data{:});
[nb na M] = size(hi);

fit = ftab.fit;
xrs = ftab.xrs;
lin = ftab.lin;

% need modification 2011-10-11
if nargin == 7 % handle subset case
	ia = iblock:nblock:size(hi,2);
    Gb = Gb{iblock};
    if M == 1
        hi = hi(:, ia);
        I0 = I0(:, ia);
    else
        hi = hi(:, ia, :);
        I0 = I0(:, ia, :);
    end
%	warn 'not tested'
end
% precomputed curvature for SQS
if sum(sil(:)) == 0 && streq(curvtype, 'pc') && (coupled_material == 0) % precomputed curvature
    % g_i = sum_j g_ij.  caution: this requires real g_ij and g_ij >= 0
    gi = sum(Gb');
    gi = max(reale(gi), 0); % trick: just in case...
    gi = reshape(gi, [nb na]);
    
    % effective energies: xrs.eff [Ns M]
    % linear attenuation coef. at effective energies [Ns L M]
    for mm = 1 : M
        lac(:, :, mm) = lin.lac(xrs.eff(:, mm)); %[nb L M]
        lac_sum(:, mm) = sum(lac(:, :, mm), 2); %[nb M]
    end
    
    for ll = 1 : ftab.LL
        lac_eff = squeeze(lac(:, ll, :)) .* lac_sum; %[nb M]
        lac_eff = repmat(lac_eff, [1 1 na]); %[nb M na]
        lac_eff = permute(lac_eff, [1 3 2]);
        hi_lac = hi .* lac_eff; %[nb na M]
        hi_lac = sum(hi_lac, 3);%[nb na]
        curv(:, ll) = Gb' * col( gi .* hi_lac); %[np L]
    end
        

    grad = 0;%dummy
elseif sum(sil(:)) == 0 && streq(curvtype, 'pc') && (coupled_material == 1)
     % precomputed curvature for PWSQS  
     % g_i = sum_j g_ij.  caution: this requires real g_ij and g_ij >= 0
    gi = sum(Gb');
    gi = max(reale(gi), 0); % trick: just in case...
    gi = reshape(gi, [nb na]);
    
    if 1 
    % percomputed curvature matrix C^{PC}_{im}
    Cim = zeros(nb, ftab.LL, ftab.LL, M);
    for mm = 1 : M
        lac = fit.mac{mm}; %[Ne L]
        Ne = size(lac, 1); 
        lac_matrix = zeros(Ne, ftab.LL, ftab.LL); % mu*mu' LxL matrix at Ne energies
        for ee = 1 : Ne
            lac_matrix(ee, :, :) = lac(ee, :)' * lac(ee, :); % [Ne L L]
        end
        for ii = 1 : ftab.LL
            for jj = 1 : ftab.LL
                lac_matrix_element = squeeze(lac_matrix(:, ii, jj)); % [Ne 1]
                lac_matrix_element = repmat(lac_matrix_element, [1 nb]);% [Ne Ns]
                Cim(:, ii, jj, mm) = transpose(sum(fit.coef{mm} .* lac_matrix_element, 1));% [Ns 1]
            end
        end
    end
    end
    % about 10^15 magnitude larger than the old curvature
    % it has something to do with the Ide (not normalized)
    % maybe Pde is the right one to use
    % Yong Long, 2013-12-05
    if 0
    lac = lin.lac(xrs.en); % linear attenuation coefficients [Ne L]
    Ne = size(lac, 1);
    lac_matrix = zeros(Ne, ftab.LL, ftab.LL); % mu*mu' LxL matrix at Ne energies 
    for ee = 1 : Ne
        lac_matrix(ee, :, :) = lac(ee, :)' * lac(ee, :); % [Ne L L]
    end
    
    % percomputed curvature matrix C^{PC}_{im}
    Cim = zeros(nb, ftab.LL, ftab.LL, M);
    I = repmat(xrs.I, [1 1 Ne]); %[Ns M Ne]
    Pde = xrs.Ide ./ permute(I, [3 1 2]);%[Ne, Ns, M]
    % xrs.Ide [Ne, Ns, M]  I(E)dE
    % Ne = 140, but 80-140 = 0 for xrs.Ide of 80kVp
    for mm = 1 : M
        for ii = 1 : ftab.LL
            for jj = 1 : ftab.LL
                lac_matrix_element = squeeze(lac_matrix(:, ii, jj)); % [Ne 1]
                lac_matrix_element = repmat(lac_matrix_element, [1 nb]);% [Ne Ns]
                Cim(:, ii, jj, mm) = transpose(sum(Pde(:,:,mm) .* lac_matrix_element, 1));% [Ns 1]
            end
        end
    end
    end
    
    Cim = sum(Cim, 4); %[Ns L L]
    for ii = 1 : ftab.LL
        for jj = 1 : ftab.LL
            elementCim = squeeze(Cim(:, ii, jj)); %[Ns 1]
            elementCim = repmat(elementCim, [1 na]); %[nb na]
            curv(:, ii, jj) = Gb' * col( gi .* elementCim); %[np 1]
        end
    end
    grad = 0;%dummy
else
    if 1 == size(sil,2)
        sil = reshapee(sil, [], fit.LL);
    end
    fi = fit.fmfun(sil);
  

    [grad curv] = de_ftab_pl_grad_curv(fit, hi, I0, fi, sil, Gb, curvtype); % [*ns L]
   
    if any(isnan(grad(:))) || any(isnan(curv(:)))
        warn 'nan'
        keyboard
    end

end


%
% de_ftab_pl_grad_curv()
% gradient of likelihood-based cost function of the EF beam-hardening
% method
% in
%	fit	strum
%	hi	(ns)	measurements
%   I0  (ns)    intensity
%	fi   (ns)	f_i (log predictions from sil), possibly empty
%					fit.fmfun(sil)
%	sil	[(ns) L]	component material density integrals
%   Gb	[(ns), np]		system matrix (see Gblock)
%   curvtype         'pc' or 'fly'
% out
%	grad	[(ns) L]	gradient of likelihood-based cost function
%	curv	[(ns) L]	curvatures of likelihood-based cost function
%
function [grad curv] = de_ftab_pl_grad_curv(fit, hi, I0, fi, sil, Gb, curvtype)

L = fit.LL;
ns = size(sil);
M = fit.MM;
if L == 1
	if ns(end) == 1, ns = ns(1:end-1); end
else
	if fit.LL ~= ns(end), error 'sil size', end; ns = ns(1:end-1);
end

sil = reshapee(sil, [], L); % [*ns L]

if isempty(fi)
	fi = fit.fmfun(sil); % [*ns]
end

fgrad = fit.fgrad(sil); % [*ns L M] gradient of fi wrt si

if ~streq(fit.ctype, 'pre10')
	fail 'not done'
end

hgrad = reshapee(hi - I0 .* exp(-fi), size(fgrad, 1), M ); %[*ns] gradient of hi wrt fi

for ii = 1 : fit.LL
    sgrad = sum(squeeze(fgrad(:, ii, :) ).* hgrad, 2);
    grad(:, ii) = Gb' * sgrad; %[np L] gradient of hi wrt f(p)*p
end

if streq(curvtype, 'pc')
    curv = 0; %dummy, use precomouted curvature!
elseif streq(curvtype, 'fly') 
    
    Cim = zeros(prod(ns), L, L, M); %
    for mm = 1 : M
       % alf = fit.coef{mm}; % [Ne Ns]
        mac = fit.mac{mm}; % [Ne L]
        Ne = size(mac, 1); % # exponential terms in the fit
        mac_sil = sil * mac'; %[*Ns Ne]
        mac_sil(mac_sil < 0) = 0;
        A = exp(-mac_sil); % [*Ns Ne]
        % optimal curvature of f = exp(-x) + alpha * x
        cim_opt = 2 * (1 - A - mac_sil .* A) ./ (mac_sil .^ 2);
        cim_opt = max(0, cim_opt);
        cim_opt = min(cim_opt, 1);
        cim_opt(mac_sil == 0) = 1; % [*Ns Ne] 
        if size(ns, 2) == 1
            na = ns / size(fit.coef{mm}, 2);
            coef = reshape(repmat(fit.coef{mm}, [1 na]), [], ns);%[Ne Ns]->[Ne Ns Na]->[Ne *Ns] todo, Nt>1
        else
            coef = reshape(repmat(fit.coef{mm}, [1 ns(1) ns(2)]), [], prod(ns));%[Ne Ns]->[Ne Ns Na]->[Ne *Ns] todo, Nt>1
        end
        Q = cim_opt .* coef'; % [*Ns Ne]
        % curvature matrix
        mac_matrix = zeros(Ne, L, L); % mu*mu' LxL matrix at Ne energies
        for ee = 1 : Ne
            mac_matrix(ee, :, :) = mac(ee, :)' * mac(ee, :); % [Ne L L]
        end
        for ii = 1 : L
            for jj = 1 : L
                mac_matrix_element = squeeze(mac_matrix(:, ii, jj)); % [Ne 1]
                mac_matrix_element = repmat(mac_matrix_element, [1 prod(ns)]);% [Ne *Ns]
                Cim(:, ii, jj, mm) = transpose(sum(Q' .* mac_matrix_element, 1));% [*Ns 1]
            end
        end 
    end    
    
    Cim = sum(Cim, 4); %[Ns L L] 
    % g_i = sum_j g_ij.  caution: this requires real g_ij and g_ij >= 0
    gi = sum(Gb');
    gi = max(reale(gi), 0); % trick: just in case...
    
    for ii = 1 : L
        for jj = 1 : L
            elementCim = squeeze(Cim(:, ii, jj)); %[*Ns 1]
            curv(:, ii, jj) = Gb' * col( gi .* elementCim'); %[np 1]
        end
    end
    
else 
    error 'use precomputed curvature or on the fly curvature'
end




