% ts1_model.m

global dataroot spm_path subject_ids proc_subs use_model_proto sinfo
global SCAN_OFFSET
global regressors2add

global created_model_specs		% trying to make things more efficient
                                        % by only constructing structures once
model_spec_fname = 'model_specs.mat';

if created_model_specs
  disp(sprintf('Loading model from %s', model_spec_fname));
  load(model_spec_fname)
else
  nsub = length(subject_ids);
  nsub_proc = length(proc_subs);

  subj_nsess = length(sinfo(1).cond_order);

  for isub = 1:nsub_proc
    sub_idx = proc_subs(isub);
    disp(sprintf('Mapping files for %s', sinfo(sub_idx).id))
    for isess = 1:subj_nsess
      sessions(sub_idx).images{isess} = spm_get('Files', ... 
	  sprintf('%s/%s/epi/run%d/', dataroot, char(subject_ids(sub_idx)), isess), ...
	  sprintf('sn%s*.img', char(subject_ids(sub_idx))) ...
	  );   
    end  % for isess=
  end  % for isub=

  ts1_model_protos;
  
  model_init = model_prototype(use_model_proto);
  conditions = conditions_prototype{use_model_proto};

  nregress_idx = 0;			% Overall number or entries in
  % regressor specification
  % structure. Since most regressors are
  % likely to differ between subjects,
  % every subject has a unique set of
  % regressor structure entries.  The
  % number of entries/subject equals the
  % number of runs
  for isub = 1:nsub_proc
    sub_idx = proc_subs(isub);
    
    disp(sprintf('Specifying model parameters for %s', sinfo(sub_idx).id));
    
    model(sub_idx) = model_init;
    model(sub_idx).files = sessions(sub_idx).images;
    
    % Set the correct run order for this subject
    model(sub_idx).conditions = sinfo(sub_idx).cond_order;
    model(sub_idx).nscans = sinfo(sub_idx).nvol-SCAN_OFFSET;
    
    %
    % Add regressors if requested
    %
    
    nr2add = length(regressors2add);
    
    nruns = length(sinfo(sub_idx).cond_order);
    for irun = 1:nruns
      nregress_idx = nregress_idx+1;
      regress_nb(irun) = 0;		% The number of regressors in this run
      sub_regressor_idxs(irun) = nregress_idx; % regressor indices specific to this subject
      
      for ir = 1:nr2add
	switch regressors2add{ir}
	  case 'linear'
	    run_regress_start_idx = regress_nb(irun)+1;
	    
	    regressors(nregress_idx).names(run_regress_start_idx) = {'linear'};

	    vals = 0:sinfo(sub_idx).nvol(irun)-SCAN_OFFSET-1;
	    vals = vals'/max(vals);
	    regressors(nregress_idx).values(1:sinfo(sub_idx).nvol(irun) - ...
		SCAN_OFFSET, run_regress_start_idx) = vals(:); 
	    
	    % Global tallies
	    regress_nb(irun) = regress_nb(irun) + 1;
	    
	  case 'motion'
	    run_regress_start_idx = regress_nb(irun)+1;

	    motion_names = {'x','y','z','pitch','roll','yaw'};
	    nmotion_regress = length(motion_names);
	    
	    % Entries specific to this regressor
	    regressors(nregress_idx).names(run_regress_start_idx: ...
		run_regress_start_idx + nmotion_regress-1) = motion_names;
	    
	    % Load regressors
	    fname = spm_get('Files',sprintf('%s/%s/epi/run%d/', dataroot, char(subject_ids(sub_idx)), irun),'realignment_params*.txt');
	    
	    vals = load(fname);
	    ncol = size(vals,2);
	    if ncol ~= nmotion_regress
	      error('Problem getting expected number of motion regressors')
	    else
	      regressors(nregress_idx).values(1:(sinfo(sub_idx).nvol(irun) - ...
		  SCAN_OFFSET), regress_nb(irun)+1:regress_nb(irun)+ncol) = vals;   
	    end
	    
	    % Global tallies
	    regress_nb(irun) = regress_nb(irun) + nmotion_regress;
	    
	  otherwise
	    disp(sprintf('Unknown regressor type: %s', regressors2add{ir}))
	end % switch regressor type
      end % end adding of regressors
    end % run loop
    
    model(sub_idx).regressors = sub_regressor_idxs;
    model(sub_idx).regressors_nb = regress_nb;
  end % end of subject loop

  % Save model information to a specs file to save time
  if ~exist('regressors'), regressors = []; end
  if ~exist('stochastics'), stochastics = []; end
  if ~exist('parametrics'), parametrics = []; end
  
  disp(sprintf('Saving model specifications to %s', model_spec_fname))
  save(model_spec_fname,'model','conditions','regressors','stochastics','parametrics','bf_ep','bf_ev')
  created_model_specs = 1;
end % if created_model_specs
