function [iFR] = LNPm_sta_predict(SE,k,NLfunc,NLparams,nlin,xb)

% [iFR] = LNPm_sta_predict(SE,k,NLfunc,NLparams,nlin,xb)
%
% SE  - stimulus ensembles: dimensions x ensembes
% k   - filter kernel
% NLfunc   - @MyNonlinFunc --> parametric non-linearity
% NLparams - parameter for @MyNonlinFunc
% nlin - non-parametric non-linearity
% xb   - linear response output related to nlin (i.e., the x-axis)
% 
% Reference:
% Schwartz et al. (2006) Spike-triggered neural characterization. J Vison 6:484-507.
%
% Description: The script predicts the instantenous firing rate based on the
% spike-triggered average (STA) kernel. If all inputs are provided the
% non-linearity is estimated based on the parametric non-linearity. Use [] for 
% the third and fourth input to calculate the non-parametric non-linearity.
%
% ---------
%
%    Copyright (C) 2016, B. Herrmann
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
% ----------------------------------------------------------------------------
% B. Herrmann, Email: herrmann.b@gmail.com, 2016-04-11

% dot product / STA response / linear filter output
dprodSE  = k' * SE;

if nargin >= 4 && isa(NLfunc,'function_handle') && ~isempty(NLparams)
	% calculate the instantenous firing rate based on parametric non-linearity
	iFR = NLfunc(dprodSE,NLparams);
elseif nargin == 6 && ~isempty(nlin) && ~isempty(xb)
	% calculate the instantenous firing rate based on non-parametric non-linearity
	ixnan = isnan(nlin);
	xb(ixnan)   = [];
	nlin(ixnan) = [];

	[~,ix] = min(abs(repmat(xb',[1 length(dprodSE)])-repmat(dprodSE,[length(xb) 1])),[],1);
	iFR = nlin(ix);
else
	error('Error: Inputs are not correct!')
end
