function thresh=ZEST(response,init)
% thresh=ZEST(response,init)
%
% This routine generalizes the ZEST procedure. The first invocation should
% have a second parameter "init" (actually a structure) that specifies the
% various parameters to initialize the ZEST routine. In particular, init
% should have the following fields:
%   zestA   PDF scale factor
%   zestB   falling PDF slope
%   zestC   rising PDF slope
%   zestmaxrange    the maximum value for possible estimates (usually in
%                   log units)
%   zestminrange    the minimum value for possible estimates (usually in
%                   log units)
%   zestfa  False alarm rate estimate (0<= fa <=1)
%   zestmiss    Miss rate estimate (0<= miss <=)
%   zestbeta    Response function slope
%   zesteta     "sweat factor" or response criterion parameter
% If there is a passed structure init, then the value of "response" (though
% necessary) is ignored.
%
% Otherwise, just response (=0 or 1) is passed to ZEST. Based on teh
% current set of parameters, the routine returns the most probable mean
% vlaue of the PDF that can either be used for the next trial OR as a final
% threshold estimate.
%
% To be done: return the variance of the current PDF so one could stop at a
% particular confidence interval rather than after so many trials.
%
% Note that this code is currently written with default values for a
% modulation detection task. 
% 
% For details about the parameters, see Marvit, et al. (2003), JASA
% 113(6):3348-3361.
%
% In the main code, we specify modulation depth as percentage, but it is
% better expressed and manipulated in a dB scale. So: Convert modulation
% depth, m, to delta L in dB 
%   DLdB=20*log10((1+m)./(1-m))

% CHeck Default values. Currently for 2AFC
persistent T numbersteps A B C maxrange minrange fa miss beta eta q meanpdf convert

% Initialize
if nargin > 1
	% Parameters for Initial PDF for ZEST
	if isfield(init,'zestA') A=init.zestA;, else, 
		A=1;    % Scale factor to start with unity pdf area
	end
	if isfield(init,'zestB') B=init.zestB;, else, 
		B=3.0;  % Falling slope
	end
	if isfield(init,'zestC') C=init.zestC;, else, 
		C=1.5;  % Rising slope
	end	
	if isfield(init,'zestmaxrange') maxrange=init.zestmaxrange;, else, 
		maxrange=-2.5;
	end
	if isfield(init,'zestminrange') minrange=init.zestminrange;, else, 
		minrange=2.5;
	end
	
	% Parameters for Response function (Weibull function) (P(yes) given stimulus)
	if isfield(init,'zestfa') fa=init.zestfa;, else, 
		fa=0.50;   %gamma in the text, false alarm rate (guess rate for 2AFC)
	end
	if isfield(init,'zestmiss') miss=init.zestmiss;, else, 
		miss=0.01; %delta in the text, miss rate (1/2 inattention rate for 2AFC)
	end
	if isfield(init,'zestbeta') beta=init.zestbeta;, else, 
		beta=.6;    %beta in the text, slope of response function
	end
	if isfield(init,'zesteta') eta=init.zesteta;, else, 
		eta=0;     %eta in the text, "sweat factor" or response criterion parameter
	end
	
	% Starting params
	if isfield(init,'zestinit_AM') 
        m=init.zestinit_AM; % Comes in as a percentage
    else 
		m=.5; % Initially, start with 50% AMdepth
	end
    % But convert it to log(dB)
    init_t=log10(20*log10((1+m)/(1-m))); % Use initial "guess" as midpoint of PDF
	
	% Finally, flag to know how to return the threshold
	% POssibilities include 'm' (AMdepth->(0,100)%), 'mdB', 'logmdB' 
	if isfield(init,'zestconvert') convert=init.convert, else,
		convert='m';
	end
		
	% Create a discrete vector of the extent/range of possible DLs
	T=linspace(minrange,maxrange,10000); % Linear DL's can vary from .01 dB to 20 dB in .01 dB increments
	numbersteps=length(T);
	
	% Calculate the initial PDF
	q=A./(B*exp(-C*(T-init_t)) + C*exp(B*(T-init_t)));
	meanpdf=init_t;
else % if nargin>1
	% If we just have a response, calculate the next thresh. The prior
	% thresh estimate (stimulus value) was meanpdf.

	%   Psychometric function (Weibull):p is  model prob of resp given log_lev_diff if true log_DL is T
	p=1-miss-((1-fa-miss)*exp(-10.^(beta*(meanpdf-T+eta))));
	if response==0
		p=1-p;
	end

	% Compute the next q (the next pdf)
	q=p.*q;
	meanpdf=sum(T.*q)/sum(q); % Calculate new midpoint
	%thresh=meanpdf;
end %if nargin>1

% FInal conversions, if any.
if strcmp(convert,'m')
    % De-logify
    antimean=10^meanpdf;
    % Now, de-dBify
	thresh=((10^(antimean/20))-1)/((10^(antimean/20))+1);
elseif strcmp(convert,'mdB')
    % Just de-logify
	thresh=10^meanpdf;
else
	thresh=meanpdf;
end