function results = MCarlo_psearch(model_data,search_params,prev_results)

% MCARLO_PSEARCH -- Monte Carlo PARAMSEARCH over many synthetic data sets.
%
% results = MCarlo_psearch(model_data,search_params)
% results = MCarlo_psearch(model_data,search_params,prev_results)
%
% MODEL_DATA is a structure array having one struct per "subject".
% In addition to the other fields, it must have a field 'eparams' keeping
% the "empirical parameters" that PARAMSEARCH has produced from the
% observed responses for this subject. Each MODEL_DATA struct is used
% in two ways: within the generative model, and within the model-tracer.
%
% SEARCH_PARAMS is the same structure as in BIG_PARAMSEARCH (that is, with
% fields  'model_name', 'params', 'p2v_templ', 'v2p_templ', 'bounds,
% 'clip_prob', 'optns', 'deadline', 'fname', and 'verbose'). 
% In addition it must have the following Monte-Carlo specific fields:
%  11. gener_name -- M-file implementing the generative model (e.g. 'anchor2')
%        [Note: Do not confuse with the model-tracer (e.g. 'anchor2mt'), ]
%        [      whose name is specified in SEARCH_PARAMS.model_name.     ]
%  12. MC_iter -- number of Monte Carlo iterations per subject
%  13. MC_order -- when set to 'iter' (the default), the algorithm does
%        one iteration for all subjects, then another iteration, etc;
%        when set to 'sbj', the algorithm finishes all iterations for the
%        first subject, then moves on to the second subject, etc.
%  14. resetp -- when 1 (the default), the function resets the seed of
%        MATLAB's random number generator   % random('seed',sum(100*clock))
%
% The function produces a structure RESULTS with the following fields:
%  1. descriptor -- a string, e.g. '04-Jun-2001 18:55:25' (see DATESTR(NOW))
%  2. search_params -- a copy of the input argument
%  3. opt_X -- [NxpxI] matrix, where N is the number of "subjects" (that is,
%        the size of the MODEL_DATA array), p is the number of search
%        dimensions (as specified in SEARCH_PARAMS.pv2_templ), and I is
%        the number of Monte Carlo iterations
%  4. L -- [NxI] vector of L values at the optimum (see LOGL, PARAMSEARCH)
%  5. curr_state = [curr_sbj curr_iter] -- where to resume after interrupts
%
% When the search has been interrupted (by the user or by deadline), it
% can later be resumed on the basis of the partially completed structure
% saved on the disk.  To that end, pass this structure as PREV_RESULTS.
% MCARLO_PSEARCH will read the 'curr_state' field and continue from there.
%
% Reference: WH Press, SA Teukolsky, WT Vettering, BP Flannery (1992)
%   "Numerical Recipes in C" (2nd Ed, Section 15.6, p. 691, Fig. 15.6.2)
%
% See also BIG_PSEARCH, PARAMSEARCH, LOGL, VECTOR2PARAMS.

% Original coding by Alex Petrov, Ohio State University
% $Revision: 1.0 $  $Date: 2001/08/11 19:30 $
%
% Part of the utils toolbox version 1.1 for MATLAB version 5 and up.
% http://alexpetrov.com/softw/utils/
% Copyright (c) Alexander Petrov 1999-2006, http://alexpetrov.com
% Please read the LICENSE and NO WARRANTY statement in ../utils_license.m

%--  Handle the input arguments
N_sbj = length(model_data) ;                % number of subjects
N_dim = length(search_params.p2v_templ) ;   % number of active parameters
N_iter = search_params.MC_iter ;            % Monte Carlo iterations per sbj

if (nargin < 3)              % PREV_RESULTS not supplied
  prev_results = [] ;        % build up from scratch
end
if (isempty(prev_results))
  results.descriptor = datestr(now) ;
  results.search_params = search_params ;
  results.opt_X = zeros(N_sbj,N_dim,N_iter) ;
  results.L = zeros(N_sbj,N_iter) ;
  results.curr_state = [0 0] ;      % [curr_sbj curr_iter]
else
  results = prev_results ;
  results.search_params = search_params ;   % in case they've changed
end

%--  Supply the optional fields of SEARCH_PARAMS if necessary
if (isfield(search_params,'deadline'))
  deadline = search_params.deadline ;
else
  deadline = +Inf ;      % no deadline, default
end
if (isfield(search_params,'fname'))
  fname = search_params.fname ;
else
  fname = ['MCarlo-' date] ;   % e.g. 'MCarlo-04-Jun-2001'
end
if (isfield(search_params,'verbose'))
  verbose = search_params.verbose ;
else
  verbose = 1 ;
end
if (isfield(search_params,'MC_order'))
  MC_order = lower(search_params.MC_order) ;
else
  MC_order = 'iter' ;
end
if (isfield(search_params,'resetp'))
  resetp = search_params.resetp ;
else
  resetp = 1 ;
end

%--  Reset the seed of MATLAB's random number generator
%  MATLAB starts each session with the same seed, so if you run two
%  parallel Monte Carlo searches on two machines, you risk visiting
%  the same points twice.
if (resetp)
  seed = sum(100*clock) ;
  fprintf('The random-generator seed set to %g.\n',seed) ;
  rand('state',seed) ;
end

%--  Main loop
curr_state = next_state(results.curr_state,MC_order,[N_sbj N_iter]) ;
while (~isempty(curr_state))   % finished?
  curr_sbj = curr_state(1) ;
  curr_iter = curr_state(2) ;
  if (now > deadline)
    if verbose  
      fprintf('%s: Time out, cannot proceed with subject %d iter %d.\n\n',...
              datestr(now),curr_sbj,curr_iter) ;
    end
    break
  else
    %- do the real work: call DO-ONE-ITER for this subject
    [X,L,mresp] = do_one_iter(model_data(curr_sbj),search_params) ;
    %- record the results for this subject, ignoring MRESP
    results.descriptor = datestr(now) ;     % update the time stamp
    results.opt_X(curr_sbj,:,curr_iter) = X ;
    results.L(curr_sbj,curr_iter) = L ;
    results.curr_state = curr_state ;
    %- save partial RESULTS to disk and report to the console
    if (isempty(fname))        % saving option switched off?
      if verbose
        fprintf('%s: Done iter %3d, sbj %2d, L=%g\n',...
                results.descriptor,curr_iter,curr_sbj,L);
      end
    else
      save(fname,'results') ;  % save partial RESULTS to disk
      if verbose
        fprintf('%s: Done iter %3d, sbj %2d, L=%g, saved to %s.mat\n',...
                results.descriptor,curr_iter,curr_sbj,L,fname);
      end
    end %if isempty(fname)
  end %if (now > deadline)  
curr_state = next_state(results.curr_state,MC_order,[N_sbj N_iter]) ;
end  % while (~isempty(curr_state))

%--- Return RESULTS
%%%% End of the main function

%--- Subfunctions
function  new_state = next_state(old_state,order,max_state)

switch order
case 'iter'          % SBJ is the inner loop, ITER the outer
  new_sbj = old_state(1) + 1 ;
  new_iter = old_state(2) ;
  if (new_sbj > max_state(1))
    new_sbj = 1 ;
    new_iter = new_iter + 1 ;
  end
case 'sbj'           % ITER is the inner loop, SBJ the outer
  new_sbj = old_state(1) ;
  new_iter = old_state(2) + 1 ;
  if (new_iter > max_state(2))
    new_iter = 1 ;
    new_sbj = new_sbj + 1 ;
  end
otherwise
  error('serac_params.MC_order must be either ''iter'' or ''sbj''.') ;
end

if ((new_sbj > max_state(1)) | (new_iter > max_state(2)))
  new_state = [] ;    % end signal
else
  new_state = max([1 1],[new_sbj new_iter]) ;
end
%--  Return NEW_STATE
%%%  End of subfunction NEXT_STATE

function  [X,L,mresp] = do_one_iter(mdata,search_params)

%-- Phase 1: Generate synthetic responses
eparams = mdata.eparams ;       % parameters derived from the empirical data
mdata.params = eparams ;
mresp = feval(search_params.gener_name,mdata) ;  % call the generative model

%-- Phase 2: Run PARAMSEARCH over these synthetic responses
mdata.params = search_params.params ;   % just in case
mdata.resp = mresp ;
[opt_params,L,X] = paramsearch(mdata,search_params) ;  % <-- heavy load

%--  Return X, L, and MRESP
%%%  End of subfunction DO_ONE_ITER

%%% End of file MCARLO_PSEARCH.M
