function  F = PLE_exponent(params,xdata)

%PLE_EXPONENT -- Perceptual learning: exponential law with switch costs.
%
% F = PLE_exponent(params,xdata)
% F = PLE_exponent(params)
%
% This function is designed to be minimized via PLE_CURVEFIT or LSQCURVEFIT.
% It implements an additive mixture of two exponential curves -- 
% "main carrier" and "switch cost". The basic equation is:
%   F(T,R) = D.*(1 - g*exp(-(T-1)/tau) - s*exp(-(R-1)/tau_s)
%
% XDATA is a two-column matrix, typically with 32 rows (one per block).
% XDATA(:,1) is the block number T since the beginning of the experiment.
% XDATA(:,2) is the number of blocks R since the last switch.
% If not supplied, XDATA defaults to [1:32 ; Inf,1:8,1:8,1:8,1:6,1]'.
%
% PARAMS is a vector. There are several "modes" depending on its length L:
%  L=3 (baseline): params=[D,g,tau] -- ignores switch costs altogether
%  L=5 (standard): params=[D,g,tau, s,tau_s] -- see the equation above
%      - D is the asymptotic F level (when T=Inf)
%      - g is the relative magnitude of the learning effect: g=(D-F(1))/D
%      - tau is the time constant of the main exponent
%      - s is the (relative) switch cost (when R=1)
%      - tau_s is the time constant of the switch-cost exponential decay
%  L=7 (3-curves, short): params=[D1,D2,D3,g,tau, s,tau_s] -- fits
%    three curves in parallel, each with its own D but w/ equal tau, etc.
%    This is used to fit the different stimulus contrast levels together.
%  L=9 (3-curves, long):  params[D1,D2,D3,g1,g2,g3,tau, s,tau_s] --
%    three parallel curves, each with its own D and g; common tau and s.
%
% F is a column vector with the same number of rows as XDATA (default 32)
% when L=length(PARAMS) is 3 or 5. It has thrice as many rows when L=7,9.
%
% See also PLE_2EXP, PLE_CURVEFIT, PLE_POWER, LSQCURVEFIT.

% Original coding by Alexander Petrov, UC Irvine.
% $Revision: 2.0 $  $Date: 2003/04/03 10:32 $

%-- Decipher arguments
if (nargin < 2)    % XDATA not supplied?
   t = [1:32]' ;                   % time since beginning of experiment
   r = [Inf,1:8,1:8,1:8,1:6,1]' ;  % time since last switch
else
   t = xdata(:,1) ;
   r = xdata(:,2) ;
end

%-- Dispatch based on the number of PARAMS
L = length(params) ;
if     (L==3)
   D=params(1) ; g=params(2) ; tau=params(3) ;
   F = D.*(1 - g.*exp((1-t)./tau)) ;
%elseif (L==4)
%   D=params(1) ; g=params(2) ; tau=params(3) ; s=params(4) ; 
%   tau_s = 1.18 ;  % constant for fitting RT data; based on d' fits
%   F = D.*(1 - g.*exp((1-t)./tau) - s.*exp((1-r)./tau_s)) ;
elseif (L==5)
   D=params(1) ; g=params(2) ; tau=params(3) ; 
   s=params(4) ; tau_s=params(5) ;
   F = D.*(1 - g.*exp((1-t)./tau) - s.*exp((1-r)./tau_s)) ;
elseif (L==7)
   D1=params(1) ; D2=params(2) ; D3=params(3) ;
   g=params(4) ; tau=params(5) ; s=params(6) ; tau_s=params(7) ;
   f = (1 - g.*exp((1-t)./tau) - s.*exp((1-r)./tau_s)) ;  % [32x1]
   F = [D1.*f ; D2.*f ; D3.*f] ;                          % [96x1]
elseif (L==9)
   D1=params(1) ; D2=params(2) ; D3=params(3) ;
   g1=params(4) ; g2=params(5) ; g3=params(6) ;
   tau=params(7) ; s=params(8) ; tau_s=params(9) ;
   fs = (1 -s.*exp((1-r)./tau_s)) ;
   fg = exp((1-t)./tau) ;
   F = [D1.*(fs - g1.*fg) ; D2.*(fs - g2.*fg) ; D3.*(fs - g3.*fg)] ;
else
   error('PARAMS should have 3, 5, 7 or 9 elements.') ;
end

%--- Return F
%%%%% End of PLE_EXPONENT.M
