  function fit = de_ftab_fit_yl(xrs, mac, varargin)
%|function fit = de_ftab_fit_yl(xrs, mac, [options])
%|
%| Fit to DE table fm(), suitable for subsequent interpolation / extrapolation.
%| Uses either classic polynomial basis functions,
%| or an exponential model: -log(sum_k p_k exp(-m_k . s))
%|
%| in
%|  xrs	strum	X-ray spectra; see xray_read_spectra.m
%|  mac    [Ne L]	mass attenuation coefficient, at spectrum energies
%|
%| option
%|
%|	'show'	1|0		plot? (default false)
%|	'ntype' ''		options for negative sl values (default '')
%|				'abs' - experimental symmetrizing option
%|	'macbar' [M L]		if nonempty, use it to match derivative at s=0.
%|
%| out
%|	fit	strum		strum object for fitted f_m(s_1,...,s_L)
%|				fit.coef is [nbasis M] for 'poly'
%|	methods:
%|				(sll denotes stacked array [(Ns) L] or [*Ns L])
%|	.fmfun(sll [() L])	fm function evaluation [() L] -> [() M]
%|	.fgrad(sll)		fm gradient evaluation [() L] -> [() L M]
%|	.fhess(sll)		fm hessian evaluation [() L] -> [() L L M]
%|	.show_sp([en, sp])	plot true spectrum vs fitted spectrum
%|	.show_fm(sl, [fm])	mesh plot of fm and its fit
%|	.show_err(sl, fm)	mesh plot of fit error
%|	.mac_eff		effective mass atten coef based on fit
%|				(valid for 'exp' only)
%|
%| Copyright 2006-3-3, Jeff Fessler, University of Michigan
%|
%|
%| Yong Long, 2010-04-06,  no approximation to f_im

if nargin == 1 && streq(sl, 'test'), de_ftab_fit_test, return, end
if nargin < 2, help(mfilename), error(mfilename), end

arg.show = false;
arg.ntype = ''; % how to handle negatives
arg.macbar = [];
arg = vararg_pair(arg, varargin);
 
MM = size(xrs.sp, 2); % # of spectra
LL = ncol(mac);
Ne = size(xrs.en, 1);
kev = xrs.en; 
I = repmat(xrs.I, [Ne 1]); %[Ne M]
Pde = xrs.Ide ./ I;%[Ne M]

% no approximation to beam-hardening function
fit = de_ftab_fit_no(MM, mac, kev, Pde);

fmfun = @de_ftab_fit_no_eval; 
fgrad = @de_ftab_fit_no_grad; %todo
fhess = @de_ftab_fit_exp_hess; %todo

fit.mac_eff = zeros(MM, LL); % [M L] mac effective
for mm=1:MM
	fit.mac_eff(mm,:) = fit.mac{mm}' * fit.coef{mm};
end

     
fit.LL = LL;
fit.MM = MM;
fit.ntype = arg.ntype;

meth = {'fmfun', fmfun, '(sll [() L]) -> [() M]'; ...
    'fgrad', fgrad, '(sll [() L]) -> [() L M]'; ... 
	'fhess', fhess, '(sll [() L]) -> [() L L M]'; ...%todo
	'show_err', @de_ftab_fit_show_err, '(sl, fm)'; ...%todo
	'show_fm', @de_ftab_fit_show_fm, '(sl, [fm])'; ...%todo
	'show_sp', @de_ftab_fit_show_sp, '([en, sp])'; ...%todo
	};
fit = strum(fit, meth);

if arg.show && im
	fit.show_fm(sl, fm);
end

end % de_ftab_fit()

%
% de_ftab_fit_no()
% no approximation to beam-hardening function
%
function fit = de_ftab_fit_no(MM, mac, kev, Pde)


if isempty(mac)
	if isempty(mtype), error 'mac or mtype required', end
	mac = xray_read_atten(mtype, kev); % [Ne L]
else
	if size(mac,1) ~= length(kev), fail 'size mismatch', end
end

fit.kev = cell(1,MM); % [Ne M]
fit.mac = cell(1,MM);  % [Ne L M] 
fit.coef = cell(1,MM); %[Ne M]
for mm=1:MM
    % find key energies. intensities at low energies (<10Kev) are very
    % small, but mass attenuation coefficients are big (>1000 for <2Kev),
    % which may result in inf values when compute exp(-beta * s)
    ie = mean(Pde(:, mm),2) > 1e-6; 
    %1e-5,there was an unit error. now use cm, hope 1e-6 or 1e-7 works
    
	fit.kev{mm} = kev(ie); % [Ne 1] (now Ne may be smaller than before)
	fit.mac{mm} = mac(ie,:); % [Ne L] 
    
    % final coefficients at enerties of energy vector:
    fit.coef{mm} = Pde(ie, mm); %[Ne, 1]
end

end % de_ftab_fit_no()


%
% de_ftab_fit_no_eval()
% evaluate 
% in
%	sll	[(Ns) L]	stackup of s1,s2,...,s_L
% out
%	f	[(Ns) M]	stackup of f1,f2,...,f_M
%
function f = de_ftab_fit_no_eval(fit, sll)
LL = fit.LL;
if LL == 1
	Ns = size(sll);
	if Ns(end) == 1, Ns = Ns(1:end-1); end
else
	Ns = size(sll); if LL ~= Ns(end), fail 'bug', end; Ns = Ns(1:end-1);
end
sll = reshapee(sll, [], LL); % [*Ns L]
MM = fit.MM;

persistent warned
if ~isvar('warned') || isempty(warned)
	warned = 0;
end

f = zeros(prod(Ns),MM);
for mm=1:MM

	switch fit.ntype
	case '' % standard
		
        A = exp(-sll * fit.mac{mm}'); % [*Ns Ne]
		tmp = -log(A * fit.coef{mm}); % [*Ns 1]
		f(:,mm) = tmp;
        
	case 'abs' % trick for negatives
		if LL == 1
			A = exp(-abs(sll) * fit.mac{mm}'); % [*Ns Ne]
			tmp = -log(A * fit.coef{mm}); % [*Ns 1]
			f(:,mm) = -sign(sll) .* log(A * fit.coef{mm}); % [*Ns 1]
%			bad = sll < 0;
%			f(bad,mm) = sll(bad) * fit.mac_eff(mm,1);
		else
			if ~warned
				warned = 1;
				warn 'negative trick not done for L>1'
			end
		end

	otherwise
		fail('bad ntype: %s', fit.ntype)
	end

end
f = reshape(f, [Ns MM]);

end % de_ftab_fit_no_eval()


%
% de_ftab_fit_no_grad()
% evaluate gradient of f for each of the given s vectors.
% in
%	sll	[(Ns) L]	stackup of s1,s2,...,s_L
% out
%	g	[(Ns) L M]	stackup of gradients of f(s)
%
function g = de_ftab_fit_no_grad(fit, sll)
LL = fit.LL;
Ns = size(sll);
if LL == 1
	if Ns(end) == 1, Ns = Ns(1:end-1); end
else
	if LL ~= Ns(end); error 'bug', end
	Ns = Ns(1:end-1);
end
sll = reshape(sll, [], LL); % [*Ns L]
MM = fit.MM;

persistent warned
if ~isvar('warned') || isempty(warned)
	warned = 0;
end

g = zeros(prod(Ns), LL, MM);
for mm=1:fit.MM
	alf = fit.coef{mm}; % [Ne 1]
	mac = fit.mac{mm}; % [Ne L]
	Ne = length(alf); % # exponential terms in the fit

	switch fit.ntype
%	case '' % standard
	case 'abs'
		if any(sll(:) < 0) % trick for negatives
			if LL == 1
				sll = abs(sll);
%				bad = sll < 0;
%				g(bad,1,mm) = fit.mac_eff(mm,1);
			else
				if ~warned
					warned = 1;
					warn 'negative trick not done for L>1'
				end
			end
		end
	end

	A = exp(-sll * mac'); % [*Ns Ne]
   	Q = A .* repmat(alf', [prod(Ns) 1]); % [*Ns Ne]
    denom = sum(Q,2); % [*Ns 1]
    Q = Q ./ repmat(denom, [1 Ne]);
	g(:,:,mm) = Q * mac;

end
g = reshape(g, [Ns LL MM]);

end % de_ftab_fit_exp_grad()


%
% de_ftab_fit_exp_hess()
% evaluate hessian of f for each of the given s vectors.
% in
%	sll	[(Ns) L]	stackup of s1,s2,...,s_L
% out
%	h	[(Ns) L L M]	stackup of hessians of f(s)
%
function h = de_ftab_fit_exp_hess(fit, sll)
Ns = size(sll); LL = Ns(end); Ns = Ns(1:end-1);
sll = reshape(sll, [], LL); % [*Ns L]
MM = fit.MM;

persistent warned
if ~isvar('warned') || isempty(warned)
	warned = 0;
end

h = zeros(prod(Ns), LL, LL, MM);
for mm=1:fit.MM
	mac = fit.mac{mm}; % [Ne L]
	alf = fit.coef{mm}; % [Ne 1]
	Ne = length(alf); % # exponential terms in the fit

	switch fit.ntype
%	case '' % standard
	case 'abs'
		if any(sll(:) < 0) % trick for negatives
			if LL == 1
				sll_sign = sign(sll);
				sll = abs(sll);
%				bad = sll < 0;
%				hm(bad,1,1) = 0;
			else
				if ~warned
					warned = 1;
					warn 'negative trick not done for L>1'
				end
			end
		end
	end

	A = exp(-sll * fit.mac{mm}'); % [*Ns Ne]

	Q = A .* repmat(alf', [prod(Ns) 1]); % [*Ns Ne]
	denom = sum(Q,2); % [*Ns 1]
	Q = Q ./ repmat(denom, [1 Ne]);
	g = Q * mac; % [*Ns L] gradient

	% loop over LL because it is smaller than *Ns or Ne
	hm = zeros(prod(Ns), LL, LL);
	for l1=1:LL
		for l2=1:LL
			hm(:,l1,l2) = g(:,l1) .* g(:,l2) ...
				- Q * (mac(:,l1) .* mac(:,l2));
		end
	end

	switch fit.ntype
%	case '' % standard
	case 'abs'
		if any(sll(:) < 0) % trick for negatives
			if LL == 1
				hm(:,1,1) = hm(:,1,1) .* sll_sign;
%				bad = sll < 0;
%				hm(bad,1,1) = 0;
			end
		end
	end

	h(:,:,:,mm) = hm;

end
h = reshape(h, [Ns LL LL MM]);

end % de_ftab_fit_exp_hess()


%
% de_ftab_fit_show_sp()
% compare true spectra to fitted spectra
%
function out = de_ftab_fit_show_sp(fit, en, sp)
if nargin < 1, fail 'de_ftab_fit_show_sp(fit, en, sp)', end

if ~streq(fit.type, 'exp'), printm 'show_sp only done for exp', return, end
if im
	clf, pl = (fit.MM+1)*100 + 10 + 1;
	if nargin > 1
		subplot(pl)
		plot(en, sp * diag(1 ./ max(sp)))
		tmp = num2cell(char([1:fit.MM] + '0'));
		legend(tmp{:})
		axisx(minmax(en))
	end

	if isfield(fit, 'kev')
		for mm=1:fit.MM
			subplot(pl+mm)
			bar(fit.kev{mm}, fit.coef{mm})
			xtick(fit.kev{mm}(fit.coef{mm} > 0))
			axis tight
			if nargin > 1, axisx(minmax(en)), end
		end
	else
		warn 'kev unknown'
	end
end

if nargout, out = []; end

end % de_ftab_fit_show_sp()


%
% de_ftab_fit_show_fm()
% compare fit to sampled fm
%
function out = de_ftab_fit_show_fm(fit, sl, fm, sll)
if nargin < 2, error 'de_ftab_fit_show_fm(fit, sl, fm)', end

do_fm = nargin > 2;
if ~isvar('sll') || isempty(sll)
	sll = ndgrid_jf('mat', sl{:});
end

fh = fit.fmfun(sll);

if ~iscell(sl), fail 'sl must be cell', end
switch length(sl) % LL
case 1
	s1 = sl{1};
	im clf
	arg = {};
	leg = {};
	order = 'bgrcm';
	if do_fm
		for mm=1:fit.MM
			arg = {arg{:}, s1, fm(:,mm), [order(mm) '.']};
			leg = {leg{:}, sprintf('true m=%d', mm)};
		end
	end
	for mm=1:fit.MM
		arg = {arg{:}, sll, fh(:,mm), [order(mm) '-']};
		leg = {leg{:}, sprintf('fit m=%d', mm)};
	end
	plot(arg{:})
	xlabel 's1'
	ylabel 'fm(s1)'
	legend(leg{:}, 4)

case 2
	s1 = sl{1};
	s2 = sl{2};
	smax(1) = max(s1);
	smax(2) = max(s2);
	if do_fm
		fmax = max(fm(:));
	else
		fmax = max(fh(:));
	end

	ax = [0 smax(1) 0 smax(2) 0 fmax];
	cax = [0 fmax];

	im clf
	im('pl', 2, fit.MM)
	for mm=1:fit.MM
		if do_fm
		show(mm, s1, s2, fm(:,:,mm), ax, cax, sprintf('f_%d(s)', mm))
		end
		m0 = mm + fit.MM;
		show(m0, s1, s2, fh(:,:,mm), ax, cax, sprintf('f_%d fit', mm))
	end
%	% text(-60, 10, '[cm^2/g]')

case 3
	s1 = sl{1};
	s2 = sl{2};
	s3 = sl{3};
	smax(1) = max(s1);
	smax(2) = max(s2);
	smax(3) = max(s3);
	if do_fm
		fmax = max(fm(:));
	else
		fmax = max(fh(:));
	end

	ax = [0 smax(2) 0 smax(3) 0 fmax];
	cax = [0 fmax];

	im clf
	im('pl', 2, fit.MM)
	for mm=1:fit.MM
		if do_fm
		show(mm, s2, s3, fm(1,:,:,mm), ax, cax, sprintf('f_%d(s)', mm))
		end
		m0 = mm + fit.MM;
		show(m0, s2, s3, fh(1,:,:,mm), ax, cax, sprintf('f_%d fit', mm))
	end
%	% text(-60, 10, '[cm^2/g]')

otherwise
	fail 'not done'
end

if nargout, out = []; end

end % de_ftab_fit_show_fm()


%
% de_ftab_fit_show_err()
% show fit errors, and relate to HU
% f = mac s, mac(70kev,H2O) = 0.2 cm^2 / g = 1000 HU
% mh = mac * 1000 HU / (0.2 g/cm^2)
% so Df/Dmh = Df/Dmac * Dmac/Dmh = (50 g/cm^2) (0.2 cm^2/g) / 1000 HU = 1/100 HU
% so Dmh = 100 HU * Df for 50 cm of H2O.
%
function out = de_ftab_fit_show_err(fit, sl, fm)
if nargin < 3, error 'de_ftab_fit_show_err(fit, sl, fm)', end

sll = ndgrid_jf('mat', sl{:});

fh = fit.fmfun(sll);
err = fh - fm;
printm('worst model error = %g of %g', max(abs(err(:))), max(fm(:)))
printm('worst model error = %g HU over 50cm H2O', 100*max(abs(err(:))))
for mm=1:fit.MM
	ee = stackpick(err,mm);
	printm('worst error (mm=%d): %g', mm, max(abs(ee(:))))
end

if max(abs(err(:))) == 0
	if nargout, out = []; end
return
end

if 0
	e1 = err(:,:,1);
	e2 = err(:,:,2);
	disp([minmax(e1); minmax(e2)]')
	printm('worst error1 %g', max(col(abs(err(:,1,:)))))
	printm('worst error2 %g', max(col(abs(err(:,2,:)))))
end
%err = abs(err);

switch length(sl) % LL
case 1
	s1 = sl{1};
	arg = {};
	leg = {};
	order = 'bgrcm';
	for mm=1:fit.MM
		arg = {arg{:}, s1, err(:,mm), [order(mm) '-']};
		leg = {leg{:}, sprintf('m=%d', mm)};
	end
	plot(arg{:})
	xlabel 's1'
	ylabel 'error'
	legend(leg{:}, 4)

case 2
	s1 = sl{1};
	s2 = sl{2};

	elim = minmax(err(:))';
	%elim = [-1 1] * 0.01; % +/- 1 HU
	ax = [0 max(s1) 0 max(s2) elim];
	im clf, im('pl', 1, fit.MM)
	for mm=1:fit.MM
		show(mm, s1, s2, err(:,:,mm), ax, elim, sprintf('f_%d error', mm))
	end

case 3
	s1 = sl{1};
	s2 = sl{2};
	s3 = sl{3};

	elim = minmax(err(:))';
	ax = [0 max(s2) 0 max(s3) elim];
	im clf, im('pl', 1, fit.MM)
	for mm=1:fit.MM
		show(mm, s2, s3, err(1,:,:,mm), ax, elim, sprintf('f_%d error', mm))
	end

otherwise
	fail 'not done'
end

if nargout, out = []; end

end % de_ftab_fit_show_err()


%
% show()
%
function show(pl, x, y, f, ax, cax, ti)
f = squeeze(f); % for LL=3 case
if ~im, return, end
im('subplot', pl)
if 1
	mesh(x,y,f')
	colormap hsv, caxis(cax), cbar
	axis(ax)
	xtick, ytick, ztick, zwhite
else
	im(x,y,f), cbar
	xtick, ytick
end
xlabel 's_1', ylabel 's_2', title(ti)

end % show()


%
% de_ftab_fit_test()
%
function de_ftab_fit_test

%stype = 'mono,70,100';
%stype = 'ps1';
%stype = 'poly1,60';
stype = 'poly1,80,100,160'; % stress test
xrs = xray_read_spectra(stype);

list.sl{1} = linspace(0, 50, 26); % coarse for fast fitting
list.sl{2} = linspace(0, 30, 31);
list.sl{3} = linspace(0, 2, 11);
%list.slfine{1} = linspace(0, 5, 261); % fine for evaluating
%list.slfine{2} = linspace(0, 3, 511);
%list.slfine{3} = linspace(0, 2, 201);
list.mtype = {'water', 'bone', 'iodine'};
list.wls_simplex_arg = {{'reg', 1e-16}, {'reg', 1e-9}, {'reg', 1e-4}}; % more reg for LL>1 case

for LL=1:3
	sl = {list.sl{1:LL}};
	mtype = {list.mtype{1:LL}};
	warg = list.wls_simplex_arg{LL};

	mas = xray_read_mac(mtype);
	mac = mas.mac(xrs.en);
	if im
		clf, semilogy(xrs.en, mac), legend(mtype{:})
	end
	sll = ndgrid_jf('mat', sl{:});
	fm = de_ftab_fm(sll, mac, xrs.Ide);
%	fit = de_ftab_fit(sl, fm, 'show', 1, 'type', 'poly')
	fit = de_ftab_fit(sl, fm, 'show', 1, 'type', 'exp', 'mtype', mtype, ...
		'wls_simplex_arg', warg); 

	fh = fit.fmfun(sll);
	g = fit.fgrad(sll);
	h = fit.fhess(sll);

	if 0 && LL == 2 % test gradient
		slfine = {list.slfine{1:LL}};
		slf = ndgrid_jf('mat', slfine{:});
		fh = fit.fmfun(slf);
		g = fit.fgrad(slf);
		i1 = 5;
		i2 = 9;
		d1 = slfine{1}(i1+1) - slfine{1}(i1);
		d2 = slfine{2}(i2+1) - slfine{2}(i2);
		gh(:,1) = (fh(i1+1,i2,:) - fh(i1,i2,:)) / d1;
		gh(:,2) = (fh(i1,i2+1,:) - fh(i1,i2,:)) / d2;
		pr transpose(gh)
		gg = squeeze(g(i1,i2,:,:));
		pr gg
		pr transpose(fit.mac_eff) % at origin
%		gg - fit.mac_eff'
	end

	if im
		fit.show_sp(xrs.en, xrs.sp);
		prompt
		fit.show_fm(sl, fm);
		prompt
		fit.show_err(sl, fm);
		prompt
	end

end

end % de_ftab_fit_test()
