function  [params,RMSE] = PLE_curvefit(fun,ydata,params0,xdata,lb,ub,optns)

%PLE_CURVEFIT -- Call LSQCURVEFIT on batches of PLExp1 data.
%
%  [params,RMSE] = PLE_curvefit(fun,ydata,params0,[xdata,lb,ub,optns])
%
%  FUN is either 'PLE_exponent' or 'PLE_2exp'.
%  Each column in YDATA is treated as a separate data set, [32xN].
%  PARAMS0 is the initial parameter vector. See PLE_EXPONENT and PLE_2EXP.
%  XDATA is a column vector suitable passed to FUN. Default=FUN's default.
%  PARAMS is the resulting matrix of N optimal parameter row-vectors.
%  RMSE are the corresponding root mean square errors, [Nx1].
%
%  Example (Switch-cost params for each of the 13 subjects):
%    ydata=zeros(32,13);for k=1:13 ydata(:,k)=mean(S(k).dprime,2);end
%    [p13,RMSE] = PLE_curvefit('PLE_exponent',ydata,[2 1 .1 -.3 .5])
%
%  Example (Model2 params for the three contrast levels):
%    d=HUMAN_ZCONGR';d=[d(:,1)+d(:,2),d(:,3)+d(:,4),d(:,5)+d(:,6)];
%    [p3,RMSE] = PLE_curvefit('PLE_2exp',d,[2 1 .1 +.5 .5])
%
%  See also PLE_EXPONENT, PLE_2EXP, LSQCURVEFIT.

%  Original coding by Alexander Petrov, UC Irvine.
%  $Revision: 1.1 $  $Date: 2003/04/07 18:52 $

%-- Supply XDATA defaults depending on FUN
if (nargin < 4)
  if (strcmp(fun,'PLE_exponent'))
    xdata = [1:32 ; Inf,1:8,1:8,1:8,1:6,1]' ;   % PLE_EXPONENT's defaults
  elseif (strcmp(fun,'PLE_2exp'))
    if (length(params0)==3)                     % PLE_2EXP: Model 1
      xdata = [1,1:8,2:9,9:16,10:15,17]' ;
    else
      xdata = [1:32 ; 1,1:8,2:9,9:16,10:15,17]' ;         % Model 2
    end
  else
    error(sprintf('I don''t know XDATA defaults for FUN=''%s''',fun)) ;
  end
end

if (nargin < 6)  lb = [] ; ub = [] ;  end

if (nargin < 7)  % default options for the LSQCURVEFIT optimizer
   optns=optimset('lsqcurvefit') ; 
   optns.MaxFunEvals = 1000 ;
   %optns.Display = 'off' ;
end
  
%-- Main loop
[M,N] = size(ydata) ;           % N data sets of M elements each
params = zeros(N,length(params0)) ;
RMSE = zeros(N,1) ;

for k = 1:N
  edpr = ydata(:,k) ;           % empirical d', [32x1]
  p = lsqcurvefit(fun,params0,xdata,edpr,lb,ub,optns) ;   % pass 1
  [p,SSE] = lsqcurvefit(fun,p,xdata,edpr,lb,ub,optns) ;   % pass 2
  params(k,:) = p ;
  RMSE(k) = sqrt(SSE/M) ;
end

%--- Return PARAMS and RMSE
%%%%% End of PLE_CURVEFIT .M
