simulate_model_protos

 simulate_model_protos.m

 Creates a design matrix

 In combination with the get_exp_info* scripts, a number of different types of
 experiment designs can be simulated with relative ease, ranging from block
 designs to event designs, including cue-target type of designs.  Originally,
 the script was written to model a priming experiment using musical material,
 but that's just semantics. 

 In combination with run_post_process.m, this set of scripts can be used to
 examine properties of the design matrix and the repercussions that various
 design choices may have on the evaluation of data.

0001 % simulate_model_protos.m
0002 %
0003 % Creates a design matrix
0004 %
0005 % In combination with the get_exp_info* scripts, a number of different types of
0006 % experiment designs can be simulated with relative ease, ranging from block
0007 % designs to event designs, including cue-target type of designs.  Originally,
0008 % the script was written to model a priming experiment using musical material,
0009 % but that's just semantics.
0010 %
0011 % In combination with run_post_process.m, this set of scripts can be used to
0012 % examine properties of the design matrix and the repercussions that various
0013 % design choices may have on the evaluation of data.
0014 
0015 % 11/12/03 Petr Janata - adapted for simulating various design matrices
0016 
0017 global model_action
0018 
0019 clear tinfo
0020 
0021 nmodels = 1;
0022 nruns = 1;                % this will have to be set dynamically
0023 
0024 % Since we may want to try various types of models but not have to revise the
0025 % files containing the trial specifications each time we want to go back and
0026 % forth, we can refer to different models by different names and have the script
0027 % be smart about constructing the design matrix.  Here we create an array of
0028 % structure with experiment information. We can easily add to the array as we add
0029 % experiments we want to simulate.  We then select one of these "experiments"
0030 
0031 nexp = 0;
0032 
0033 % Create the first element in the exp_info structure array
0034 
0035 % increment the experiment list counter
0036 nexp = nexp+1;                
0037 
0038 % specify an arbitrary name for the model
0039 exp_info(nexp).name = 'block_1fact1level'; 
0040 
0041 % specify the member of the get_exp_info_* family of scripts we are going to call
0042 exp_info(nexp).script_stub = 'block'; % i.e., get_exp_info_block.m
0043 
0044 % List the conditions that will be in the model. Note, these names might refer
0045 % to the factors in the model, or factor/level combinations. Another portion of
0046 % this script will associate 'trial types' specified in the get_exp_info.m
0047 % script with each of the conditions in this list
0048 exp_info(nexp).cond_names = {'Tonic'};  % List of conditions we want to evaluate
0049 
0050 % Create the next element in the exp_info structure array
0051 nexp = nexp+1;
0052 exp_info(nexp).name = 'block_1fact2levels';
0053 exp_info(nexp).script_stub = 'block';
0054 exp_info(nexp).cond_names = {'Tonic','Subdominant'};  
0055 
0056 nexp = nexp+1;
0057 exp_info(nexp).name = 'block_2fact';
0058 exp_info(nexp).script_stub = 'block';
0059 exp_info(nexp).cond_names = {'Tonic','Subdominant','Consonant','Dissonant'};
0060 
0061 nexp = nexp+1;
0062 exp_info(nexp).name = 'block_2fact_bylevel';
0063 exp_info(nexp).script_stub = 'block';
0064 exp_info(nexp).cond_names = {'TC','TD','SC','SD'};
0065 
0066 nexp = nexp+1;
0067 exp_info(nexp).name = 'event_byfact';
0068 exp_info(nexp).script_stub = 'event';
0069 exp_info(nexp).cond_names = {'Cue silence','Cue music','Tonal','Atonal','Tonic','Subdominant','Random','Consonant','Dissonant'};
0070 
0071 nexp = nexp+1;
0072 exp_info(nexp).name = 'event_bylevel';
0073 exp_info(nexp).script_stub = 'event';
0074 exp_info(nexp).cond_names = {'Cue silence','Cue music','Tonal','Atonal','TC','TD','SC','SD','RC','RD'};
0075 
0076 nexp = nexp+1;
0077 exp_info(nexp).name = 'event_bylevel_w_silence';
0078 exp_info(nexp).script_stub = 'event';
0079 exp_info(nexp).cond_names = {'Cue silence','Cue music','Silence','Tonal','Atonal','TC','TD','SC','SD','RC','RD'};
0080 
0081 nexp = nexp+1;
0082 exp_info(nexp).name = 'event_bylevel_nosil'; % make sure to set number of
0083                                              % silent trials=0 in get_exp_info_
0084 exp_info(nexp).script_stub = 'event';
0085 exp_info(nexp).cond_names = {'Cue music','Tonal','Atonal','TC','TD','SC','SD','RC','RD'};
0086 
0087 % select which of the experiments we just specified is going to be modeled.
0088 use_exp = nexp;  
0089 
0090 % Run the script with all of the trial type and timing parameters, etc. for this experiment
0091 script_str = sprintf('get_exp_info_%s', exp_info(use_exp).script_stub);
0092 eval(script_str);
0093 
0094 % Build the appropriate condition list
0095 cond_names = exp_info(use_exp).cond_names;
0096 nc = length(cond_names);
0097 
0098 % Specify onset times and other info for each condition
0099 % Model everything as events with variable-length durations
0100 
0101 bf_ep_idx = zeros(1,nc);
0102 bf_ev_idx = ones(1,nc);
0103 
0104 %
0105 % This is the part of the code that specifies which trial types will be associated
0106 % with which conditions (members of the cond_names field in the exp_info
0107 % structure created above) in the design matrix.  All of the timing information will
0108 % get pulled from the structure and variables created in the get_exp_info
0109 % script. Note that some of the conditions model activity associated with the
0110 % prime part of each trial, whereas others reflect activity associated with the
0111 % targets.  The timing info is set accordingly.
0112 %
0113 % You can easily code novel trial_type combinations by adding a case statement
0114 % below.
0115 
0116 for ic = 1:nc
0117   switch cond_names{ic}
0118     %
0119     % Conditions that model the trial type cues
0120     %
0121     case 'Cue silence'
0122       % Find all of the desired trial types in the list of trials
0123       idxs = find(ismember(trial_list,find(ismember(tinfo.id,'Sil'))));
0124 
0125       % Get the event onset times for all of these trial types
0126       sot{ic} = onsets(idxs)/TR; 
0127 
0128       % Specify the events are events or epochs. Note: events with durations
0129       % are effectively epochs
0130       cond_types{ic} = 'events';
0131       
0132       % Specify the event durations in TR units
0133       durs{ic} = ones(size(idxs))*cue_dur/TR;  
0134 
0135     case 'Cue music'
0136       idxs = find(ismember(trial_list,find(~ismember(tinfo.id,'Sil'))));
0137       sot{ic} = onsets(idxs)/TR;
0138       cond_types{ic} = 'events';
0139       durs{ic} = ones(size(idxs))*cue_dur/TR;  % duration in msec
0140     
0141     %
0142     % Conditions that model the type of prime
0143     %
0144     case 'Tonal'  % Factor 1, Levels 1 and 2, irrespective of Factor 2 level
0145       idxs = find(ismember(trial_list, ...
0146       find(ismember(tinfo.id,{'F1L1_F2L1','F1L1_F2L2','F1L2_F2L1','F1L2_F2L2'})))); 
0147       sot{ic} = (onsets(idxs)+cue_stim_soa_list(idxs))/TR; % onset of prime
0148       cond_types{ic} = 'events';
0149       durs{ic} = ones(size(idxs))*prime_dur/TR; % duration of prime
0150 
0151     case 'Atonal'  % Factor 1, Levels 3 and 4, irrespective of Factor 2 level
0152       idxs = find(ismember(trial_list, ...
0153       find(ismember(tinfo.id,{'F1L3_F2L1','F1L3_F2L2','F1L4_F2L1','F1L4_F2L2'})))); 
0154       sot{ic} = (onsets(idxs)+cue_stim_soa_list(idxs))/TR;
0155       cond_types{ic} = 'events';
0156       durs{ic} = ones(size(idxs))*prime_dur/TR;
0157 
0158     case 'Silence'  % This could be thought of as Factor 1, Level 5
0159       idxs = find(ismember(trial_list, ...
0160       find(ismember(tinfo.id,{'Sil'})))); 
0161       sot{ic} = (onsets(idxs)+cue_stim_soa_list(idxs))/TR;
0162       cond_types{ic} = 'events';
0163       durs{ic} = ones(size(idxs))*prime_dur/TR;     
0164     
0165       %
0166       % Conditions that model the targets
0167       %
0168       
0169     case 'Tonic'
0170       idxs = find(ismember(trial_list, ...
0171       find(ismember(tinfo.id,{'F1L1_F2L1','F1L1_F2L2'}))));
0172       sot{ic} = (onsets(idxs)+cue_stim_soa_list(idxs) + prime_dur)/TR;
0173       cond_types{ic} = 'events';
0174       durs{ic} = ones(size(idxs))*ms_per_note/TR;
0175     
0176     case 'TC' % tonic, consonant
0177       idxs = find(ismember(trial_list, ...
0178       find(ismember(tinfo.id,{'F1L1_F2L1'}))));
0179       sot{ic} = (onsets(idxs)+cue_stim_soa_list(idxs) + prime_dur)/TR;
0180       cond_types{ic} = 'events';
0181       durs{ic} = ones(size(idxs))*ms_per_note/TR;
0182     
0183     case 'TD' % tonic, dissonant
0184       idxs = find(ismember(trial_list, ...
0185       find(ismember(tinfo.id,{'F1L1_F2L2'}))));
0186       sot{ic} = (onsets(idxs)+cue_stim_soa_list(idxs) + prime_dur)/TR;
0187       cond_types{ic} = 'events';
0188       durs{ic} = ones(size(idxs))*ms_per_note/TR;
0189 
0190     case 'Subdominant'
0191       idxs = find(ismember(trial_list, ...
0192       find(ismember(tinfo.id,{'F1L2_F2L1','F1L2_F2L2'}))));
0193       sot{ic} = (onsets(idxs)+cue_stim_soa_list(idxs) + prime_dur)/TR;
0194       cond_types{ic} = 'events';
0195       durs{ic} = ones(size(idxs))*ms_per_note/TR;
0196 
0197     case 'SC' % subdominant, consonant
0198       idxs = find(ismember(trial_list, ...
0199       find(ismember(tinfo.id,{'F1L2_F2L1'}))));
0200       sot{ic} = (onsets(idxs)+cue_stim_soa_list(idxs) + prime_dur)/TR;
0201       cond_types{ic} = 'events';
0202       durs{ic} = ones(size(idxs))*ms_per_note/TR;
0203 
0204     case 'SD' % subdominant, dissonant
0205       idxs = find(ismember(trial_list, ...
0206       find(ismember(tinfo.id,{'F1L2_F2L2'}))));
0207       sot{ic} = (onsets(idxs)+cue_stim_soa_list(idxs) + prime_dur)/TR;
0208       cond_types{ic} = 'events';
0209       durs{ic} = ones(size(idxs))*ms_per_note/TR;
0210 
0211     case 'Random'
0212       idxs = find(ismember(trial_list, ...
0213       find(ismember(tinfo.id,{'F1L3_F2L1','F1L3_F2L2','F1L4_F2L1','F1L4_F2L2'}))));
0214       sot{ic} = (onsets(idxs)+cue_stim_soa_list(idxs) + prime_dur)/TR;
0215       cond_types{ic} = 'events';
0216       durs{ic} = ones(size(idxs))*ms_per_note/TR;
0217       
0218     case 'RC'
0219       idxs = find(ismember(trial_list, ...
0220       find(ismember(tinfo.id,{'F1L3_F2L1','F1L4_F2L1'}))));
0221       sot{ic} = (onsets(idxs)+cue_stim_soa_list(idxs) + prime_dur)/TR;
0222       cond_types{ic} = 'events';
0223       durs{ic} = ones(size(idxs))*ms_per_note/TR;
0224     
0225     case 'RD'
0226       idxs = find(ismember(trial_list, ...
0227       find(ismember(tinfo.id,{'F1L3_F2L2','F1L4_F2L2'}))));
0228       sot{ic} = (onsets(idxs)+cue_stim_soa_list(idxs) + prime_dur)/TR;
0229       cond_types{ic} = 'events';
0230       durs{ic} = ones(size(idxs))*ms_per_note/TR;
0231       
0232     case 'Consonant'  % By Factor 2, Level 1, irrespective of prime
0233       idxs = find(ismember(trial_list, ...
0234       find(ismember(tinfo.id,{'F1L1_F2L1','F1L2_F2L1','F1L3_F2L1','F1L4_F2L1'}))));
0235       sot{ic} = (onsets(idxs)+cue_stim_soa_list(idxs) + prime_dur)/TR;
0236       cond_types{ic} = 'events';
0237       durs{ic} = ones(size(idxs))*ms_per_note/TR;
0238 
0239     case 'Dissonant'
0240       idxs = find(ismember(trial_list, ...
0241       find(ismember(tinfo.id,{'F1L1_F2L2','F1L2_F2L2','F1L3_F2L2','F1L4_F2L2'}))));
0242       sot{ic} = (onsets(idxs)+cue_stim_soa_list(idxs) + prime_dur)/TR;
0243       cond_types{ic} = 'events';
0244       durs{ic} = ones(size(idxs))*ms_per_note/TR;
0245 
0246   end % switch
0247 end % for ic=
0248 
0249 %
0250 % Fill in the generic model info
0251 %
0252 
0253 model_action_list = char({'specify','review','estimate','spec_and_estimate'});
0254 
0255 for imodel = 1:nmodels
0256   model(imodel) = struct( ...
0257       'types', strmatch(model_action,model_action_list), ... % 1=specify, 2=review, 3=estimate, 4=specify+estimate
0258       'nsess',          nruns, ...
0259       'RT',             TR/1000, ...
0260       'nscans',         [scans_per_iteration], ...    % this is actually specified in the calling script
0261       'replicated',     0, ...
0262       'same_time_param', 0, ...
0263       'conditions_nb',  ones(1,nruns) * nc(imodel), ...     
0264       'conditions',     1:nc, ...    
0265       'stochastics_flag', zeros(1,nruns), ...
0266       'stochastics',    [], ...
0267       'parametrics_type', {{'none','none'}}, ...
0268       'parametrics',    [], ...
0269       'regressors_nb',  zeros(1,nruns), ...
0270       'regressors',     zeros(1,nruns), ... 
0271       'global_effects', {'scaling'}, ...
0272       'burst_mode',     0, ...
0273       'HF_fil',         'none',  ...
0274       'HF_cut',         [], ...
0275       'LF_fil',         'none', ...
0276       'LF_cut',         [], ...
0277       'int_corr',       'none', ... 
0278       'trial_fcon',     0, ...
0279       'now_later',      1, ...        % 1 = now, 0 = later
0280       'stop_writing',   0, ...
0281       'files',          [] ...
0282       );
0283 
0284 %  model_prototype(imodel).nscans = model_prototype(imodel).nscans - SCAN_OFFSET;
0285 end
0286 
0287 %
0288 %  Conditions structures
0289 %
0290 
0291 for irun = 1:nruns
0292   conditions(irun) = struct( ...
0293       'volterra', 0, ...
0294       'variable_dur', 1, ...
0295       'names',   {cond_names}, ...
0296       'onsets',  {sot}, ...   
0297       'types',    {cond_types}, ... 
0298       'durations', {durs}, ...
0299       'bf_ev',   bf_ev_idx, ...
0300       'bf_ep',   bf_ep_idx ...
0301       );
0302 end
0303 
0304 %
0305 %  Basis function structures
0306 %
0307 
0308 bf_ev(1) = struct( ...
0309   'ev_type', 1, ...            % 1=hrf
0310   'win_len', [], ...
0311   'order', [], ...
0312   'gamma_idxs', [] ...
0313 );
0314