dsPlot

DSPlot - plot data in various ways depending on what data was provided and what options are defined.

 This function is wrapped by dsPlotWaveforms, PlotPower, etc. to provide a single
 function for organizing and displaying data.

 Usage:
   handles=dsPlot(data,'option',value)

 Inputs:
   - data: DynaSim data structure (see dsCheckData)
   - options:
     'plot_type'       : what to plot {'waveform' (default),'rastergram','rates','power'}
     'variable'        : name of field containing data to plot (default: all
                         pops with state variable of variable in data.labels)
     'time_limits'     : in units of data.time {[beg,end]}
     'max_num_overlaid': maximum # of waveforms to overlay per plot
     'max_num_rows'    : maximum # of subplot rows per figure
     'xlim'            : x-axis limits {[XMIN XMAX]} (default: all data)
     'yscale'          : whether to plot linear or log scale {'linear','log','log10'}
     'visible'         : {'on','off'}
     'lock_gca'        : Plots within currently active axis (gca); doesn't
                         open new figures or subplots.
     - NOTE: analysis options are available depending on plot_type
       - see see dsCalcFR options for plot_type 'rastergram' or 'rates'
       - see dsCalcPower options for plot_type 'power'

 Outputs:
   - handles: graphic handles to figures

 Notes:
   if Nsims>1: one sim per row
   elseif Npops>1: one pop per row
   else: one cell per row

 Examples for specifying 'variable' option:
   []      : plot all data.labels with same variable name as first element of
             data.labels (eg, 'pop1_v' and 'pop2_v')
   '*'     : plot all data.labels
   '*_v'   : plot all data.labels ending in _v (i.e., all state variables 'v'
             for all populations)
   'pop1_*': plot all data.labels starting with pop1_ (i.e., all variables for
             population 'pop1')
   'pop1_v': plot only variable 'pop1_v'
   'v'     : look for all data.labels ending in _v then starting with v_ (eg,
             all populations with variable 'v')
   'pop1'  : look for all data.labels ending in _pop1 then starting with pop1_
             (eg, all variables for population 'pop1')
   '*_iNa_*': plot all data.labels for the 'iNa' mechanism (for all populations)

 Examples:
   - Example 1: One cell:
       data=dsSimulate('dv/dt=@current+10; {iNa,iK}','tspan',[0 500]);
       dsPlot(data); % plot first state variable ('v')
       dsPlot(data,'variable','*'); % plot all state variables
       % plot all variables and time 30-60ms
       dsPlot(data,'variable','*','time_limits',[30 60]);
       % plot power spectrum
       dsPlot(data,'variable','v','plot_type','power');
       dsPlot(data,'variable','*','plot_type','power');

   - Example 2: One population with noisy input:
       data=dsSimulate('dv[5]/dt=@current+10*(1+randn(1,Npop)); {iNa,iK}','tspan',[0 250]);
       dsPlot(data);
       dsPlot(data,'variable','*'); % plot all state variables (all cells)
       dsPlot(data,'variable','m'); % plot state variable 'm' (all cells)
       % plot all variables and time 30-60ms
       dsPlot(data,'variable','*','time_limits',[30 60]);
       % plot power spectrum
       dsPlot(data,'variable','v','plot_type','power');
       dsPlot(data,'variable','*','plot_type','power');
       % plot rastergram
       dsPlot(data,'variable','v','plot_type','rastergram');
       dsPlot(data,'variable','*','plot_type','rastergram');

   - Example 3: One population varying one parameter (input amplitude):
       eqns='dv[5]/dt=@current+amp*(1+randn(1,Npop)); {iNa,iK}';
       vary={'','amp',[0 10 20]};
       data=dsSimulate(eqns,'vary',vary,'tspan',[0 200]);
       dsPlot(data);
       dsPlot(data,'variable','m');
       dsPlot(data,'variable','*');
       % plot power spectrum
       dsPlot(data,'variable','v','plot_type','power');
       % plot rastergram
       dsPlot(data,'variable','v','plot_type','rastergram');

   - Example 4: One population varying two parameters (input amplitude and
                membrane capacitance):
       eqns='dv[5]/dt=@current/Cm+amp*(1+randn(1,Npop)); {iNa,iK}';
       vary={'','Cm',[1 2]; '','amp',[0 10 20]};
       data=dsSimulate(eqns,'vary',vary,'tspan',[0 200]);
       dsPlot(data);
       dsPlot(data,'variable','*');
       % plot power spectrum
       dsPlot(data,'variable','v','plot_type','power');
       % plot rastergram
       dsPlot(data,'variable','v','plot_type','rastergram');

   - Example 5: Two populations: noisy input to E and excitatory connection from E to I
       spec=[];
       spec.populations(1).name='E1';
       spec.populations(1).equations='dv[5]/dt=@current+amp*(1+randn(1,Npop)); amp=10; {iNa,iK}';
       spec.populations(2).name='E2';
       spec.populations(2).equations='dv[2]/dt=@current; {iNa,iK}';
       spec.connections(1).direction='E1->E2';
       spec.connections(1).mechanism_list='iAMPA';
       data=dsSimulate(spec,'tspan',[0 200]);
       dsPlot(data); % plot first state variable
       dsPlot(data,'variable','*');
       % plot monitored synaptic current with post-synaptic voltages:
       dsPlot(data,'variable',{'E2_v','ISYN'});
       % plot monitored synaptic current with pre- and post-synaptic voltages:
       dsPlot(data,'variable',{'v','ISYN'});
       % plot power spectrum
       dsPlot(data,'variable','v','plot_type','power');
       dsPlot(data,'variable',{'E2_v','ISYN'},'plot_type','power');
       % plot rastergram
       dsPlot(data,'variable','v','plot_type','rastergram');

   - Example 6: Two populations varying one parameter (input amplitude):
       vary={'E1','amp',[0 10 20]};
       data=dsSimulate(spec,'vary',vary,'tspan',[0 200]);
       dsPlot(data);
       dsPlot(data,'variable','*');
       dsPlot(data,'variable','*_iNa_*');
       % plot power spectrum
       dsPlot(data,'variable','v','plot_type','power');
       % plot rastergram
       dsPlot(data,'variable','v','plot_type','rastergram');

   - Example 7: Two populations varying two parameters (input amplitude and
                synaptic conductance):
       vary={'E1','amp',[0 10 20]; 'E1->E2','gSYN',[0 .05 .1]};
       data=dsSimulate(spec,'vary',vary,'tspan',[0 200]);
       % plot voltage waveforms
       dsPlot(data,'variable','v','plot_type','power');
       % plot voltage power spectrum
       dsPlot(data,'variable','v','plot_type','waveform');
       % plot voltage-derived rastergram
       dsPlot(data,'variable','v','plot_type','rastergram');
       % more plots
       dsPlot(data,'variable','ISYN');
       dsPlot(data,'variable','E1_v');
       dsPlot(data,'variable','*');

 See also: dsCalcFR, dsCalcPower, dsPlotWaveforms, dsCheckData

 Author: Jason Sherfey, PhD <jssherfey@gmail.com>
 Copyright (C) 2016 Jason Sherfey, Boston University, USA

0001 function handles = dsPlot(data,varargin)
0002 %DSPlot - plot data in various ways depending on what data was provided and what options are defined.
0003 %
0004 % This function is wrapped by dsPlotWaveforms, PlotPower, etc. to provide a single
0005 % function for organizing and displaying data.
0006 %
0007 % Usage:
0008 %   handles=dsPlot(data,'option',value)
0009 %
0010 % Inputs:
0011 %   - data: DynaSim data structure (see dsCheckData)
0012 %   - options:
0013 %     'plot_type'       : what to plot {'waveform' (default),'rastergram','rates','power'}
0014 %     'variable'        : name of field containing data to plot (default: all
0015 %                         pops with state variable of variable in data.labels)
0016 %     'time_limits'     : in units of data.time {[beg,end]}
0017 %     'max_num_overlaid': maximum # of waveforms to overlay per plot
0018 %     'max_num_rows'    : maximum # of subplot rows per figure
0019 %     'xlim'            : x-axis limits {[XMIN XMAX]} (default: all data)
0020 %     'yscale'          : whether to plot linear or log scale {'linear','log','log10'}
0021 %     'visible'         : {'on','off'}
0022 %     'lock_gca'        : Plots within currently active axis (gca); doesn't
0023 %                         open new figures or subplots.
0024 %     - NOTE: analysis options are available depending on plot_type
0025 %       - see see dsCalcFR options for plot_type 'rastergram' or 'rates'
0026 %       - see dsCalcPower options for plot_type 'power'
0027 %
0028 % Outputs:
0029 %   - handles: graphic handles to figures
0030 %
0031 % Notes:
0032 %   if Nsims>1: one sim per row
0033 %   elseif Npops>1: one pop per row
0034 %   else: one cell per row
0035 %
0036 % Examples for specifying 'variable' option:
0037 %   []      : plot all data.labels with same variable name as first element of
0038 %             data.labels (eg, 'pop1_v' and 'pop2_v')
0039 %   '*'     : plot all data.labels
0040 %   '*_v'   : plot all data.labels ending in _v (i.e., all state variables 'v'
0041 %             for all populations)
0042 %   'pop1_*': plot all data.labels starting with pop1_ (i.e., all variables for
0043 %             population 'pop1')
0044 %   'pop1_v': plot only variable 'pop1_v'
0045 %   'v'     : look for all data.labels ending in _v then starting with v_ (eg,
0046 %             all populations with variable 'v')
0047 %   'pop1'  : look for all data.labels ending in _pop1 then starting with pop1_
0048 %             (eg, all variables for population 'pop1')
0049 %   '*_iNa_*': plot all data.labels for the 'iNa' mechanism (for all populations)
0050 %
0051 % Examples:
0052 %   - Example 1: One cell:
0053 %       data=dsSimulate('dv/dt=@current+10; {iNa,iK}','tspan',[0 500]);
0054 %       dsPlot(data); % plot first state variable ('v')
0055 %       dsPlot(data,'variable','*'); % plot all state variables
0056 %       % plot all variables and time 30-60ms
0057 %       dsPlot(data,'variable','*','time_limits',[30 60]);
0058 %       % plot power spectrum
0059 %       dsPlot(data,'variable','v','plot_type','power');
0060 %       dsPlot(data,'variable','*','plot_type','power');
0061 %
0062 %   - Example 2: One population with noisy input:
0063 %       data=dsSimulate('dv[5]/dt=@current+10*(1+randn(1,Npop)); {iNa,iK}','tspan',[0 250]);
0064 %       dsPlot(data);
0065 %       dsPlot(data,'variable','*'); % plot all state variables (all cells)
0066 %       dsPlot(data,'variable','m'); % plot state variable 'm' (all cells)
0067 %       % plot all variables and time 30-60ms
0068 %       dsPlot(data,'variable','*','time_limits',[30 60]);
0069 %       % plot power spectrum
0070 %       dsPlot(data,'variable','v','plot_type','power');
0071 %       dsPlot(data,'variable','*','plot_type','power');
0072 %       % plot rastergram
0073 %       dsPlot(data,'variable','v','plot_type','rastergram');
0074 %       dsPlot(data,'variable','*','plot_type','rastergram');
0075 %
0076 %   - Example 3: One population varying one parameter (input amplitude):
0077 %       eqns='dv[5]/dt=@current+amp*(1+randn(1,Npop)); {iNa,iK}';
0078 %       vary={'','amp',[0 10 20]};
0079 %       data=dsSimulate(eqns,'vary',vary,'tspan',[0 200]);
0080 %       dsPlot(data);
0081 %       dsPlot(data,'variable','m');
0082 %       dsPlot(data,'variable','*');
0083 %       % plot power spectrum
0084 %       dsPlot(data,'variable','v','plot_type','power');
0085 %       % plot rastergram
0086 %       dsPlot(data,'variable','v','plot_type','rastergram');
0087 %
0088 %   - Example 4: One population varying two parameters (input amplitude and
0089 %                membrane capacitance):
0090 %       eqns='dv[5]/dt=@current/Cm+amp*(1+randn(1,Npop)); {iNa,iK}';
0091 %       vary={'','Cm',[1 2]; '','amp',[0 10 20]};
0092 %       data=dsSimulate(eqns,'vary',vary,'tspan',[0 200]);
0093 %       dsPlot(data);
0094 %       dsPlot(data,'variable','*');
0095 %       % plot power spectrum
0096 %       dsPlot(data,'variable','v','plot_type','power');
0097 %       % plot rastergram
0098 %       dsPlot(data,'variable','v','plot_type','rastergram');
0099 %
0100 %   - Example 5: Two populations: noisy input to E and excitatory connection from E to I
0101 %       spec=[];
0102 %       spec.populations(1).name='E1';
0103 %       spec.populations(1).equations='dv[5]/dt=@current+amp*(1+randn(1,Npop)); amp=10; {iNa,iK}';
0104 %       spec.populations(2).name='E2';
0105 %       spec.populations(2).equations='dv[2]/dt=@current; {iNa,iK}';
0106 %       spec.connections(1).direction='E1->E2';
0107 %       spec.connections(1).mechanism_list='iAMPA';
0108 %       data=dsSimulate(spec,'tspan',[0 200]);
0109 %       dsPlot(data); % plot first state variable
0110 %       dsPlot(data,'variable','*');
0111 %       % plot monitored synaptic current with post-synaptic voltages:
0112 %       dsPlot(data,'variable',{'E2_v','ISYN'});
0113 %       % plot monitored synaptic current with pre- and post-synaptic voltages:
0114 %       dsPlot(data,'variable',{'v','ISYN'});
0115 %       % plot power spectrum
0116 %       dsPlot(data,'variable','v','plot_type','power');
0117 %       dsPlot(data,'variable',{'E2_v','ISYN'},'plot_type','power');
0118 %       % plot rastergram
0119 %       dsPlot(data,'variable','v','plot_type','rastergram');
0120 %
0121 %   - Example 6: Two populations varying one parameter (input amplitude):
0122 %       vary={'E1','amp',[0 10 20]};
0123 %       data=dsSimulate(spec,'vary',vary,'tspan',[0 200]);
0124 %       dsPlot(data);
0125 %       dsPlot(data,'variable','*');
0126 %       dsPlot(data,'variable','*_iNa_*');
0127 %       % plot power spectrum
0128 %       dsPlot(data,'variable','v','plot_type','power');
0129 %       % plot rastergram
0130 %       dsPlot(data,'variable','v','plot_type','rastergram');
0131 %
0132 %   - Example 7: Two populations varying two parameters (input amplitude and
0133 %                synaptic conductance):
0134 %       vary={'E1','amp',[0 10 20]; 'E1->E2','gSYN',[0 .05 .1]};
0135 %       data=dsSimulate(spec,'vary',vary,'tspan',[0 200]);
0136 %       % plot voltage waveforms
0137 %       dsPlot(data,'variable','v','plot_type','power');
0138 %       % plot voltage power spectrum
0139 %       dsPlot(data,'variable','v','plot_type','waveform');
0140 %       % plot voltage-derived rastergram
0141 %       dsPlot(data,'variable','v','plot_type','rastergram');
0142 %       % more plots
0143 %       dsPlot(data,'variable','ISYN');
0144 %       dsPlot(data,'variable','E1_v');
0145 %       dsPlot(data,'variable','*');
0146 %
0147 % See also: dsCalcFR, dsCalcPower, dsPlotWaveforms, dsCheckData
0148 %
0149 % Author: Jason Sherfey, PhD <jssherfey@gmail.com>
0150 % Copyright (C) 2016 Jason Sherfey, Boston University, USA
0151 
0152 % Check inputs
0153 data=dsCheckData(data, varargin{:});
0154   % NOTE: calling dsCheckData() at beginning enables analysis/plotting functions to
0155   %       accept data matrix [time x cells] in addition to DynaSim data structure.
0156 
0157 % get options
0158 options=dsCheckOptions(varargin,{...
0159   'time_limits',[-inf inf],[],...
0160   'variable',[],[],...
0161   'max_num_overlaid',50,[],...
0162   'max_num_rows',20,[],...
0163   'plot_mode','trace',{'trace','image'},...
0164   'plot_type','waveform',{'waveform','rastergram','raster','power','rates'},...
0165   'xlim',[],[],...
0166   'ylim',[],[],...
0167   'figwidth',[1],[],...
0168   'figheight',[1],[],...
0169   'lock_gca',[false],[false, true],...
0170   'yscale','linear',{'linear','log','log10','log2'},...
0171   'visible','on',{'on','off'},...
0172   'auto_gen_test_data_flag',0,{0,1},...
0173   },false);
0174 
0175 %% auto_gen_test_data_flag argin
0176 if options.auto_gen_test_data_flag
0177   varargs = varargin;
0178   varargs{find(strcmp(varargs, 'auto_gen_test_data_flag'))+1} = 0;
0179   varargs(end+1:end+2) = {'unit_test_flag',1};
0180   argin = [{data}, varargs]; % specific to this function
0181 end
0182 
0183 
0184 %% varied fields
0185 
0186 % Make sure there is no empty
0187 labels = data(1).labels;
0188 inds = arrayfun(@(s) ~isempty(s.(labels{1})),data);
0189 data = data(inds);
0190 
0191 
0192 fields=fieldnames(data);
0193 if any(strcmp(fields, 'varied'))
0194   % get varied labels
0195   vary_labels = data(1).varied; % data(1).simulator_options.vary;
0196   no_vary_labels = length(vary_labels);
0197   vary_params = nan(length(data), no_vary_labels);
0198   vary_vectors = cell(no_vary_labels, 1);
0199   vary_lengths = nan(no_vary_labels, 1);
0200 
0201   % get varied params
0202   for v = 1:no_vary_labels
0203     vary_params(:, v) = [data.(vary_labels{v})];
0204     vary_vectors{v} = unique(vary_params(:, v));
0205     vary_lengths(v) = length(vary_vectors{v});
0206   end
0207 
0208   [effective_vary_indices, ~] = dsCheckCovary(vary_lengths, vary_params, varargin{:});
0209 
0210   if prod(vary_lengths(effective_vary_indices)) == length(data)
0211 
0212       dimensions_varied = sum(effective_vary_indices);
0213       vary_params = vary_params(:, effective_vary_indices);
0214       vary_vectors = vary_vectors(effective_vary_indices);
0215       vary_lengths = vary_lengths(effective_vary_indices);
0216 
0217   else
0218 
0219       warning('unable to determine which parameters are covaried. Data will be plotted as a lattice.')
0220 
0221       dimensions_varied = 1;
0222 
0223   end
0224 
0225   if dimensions_varied > 2
0226     no_figures = prod(vary_lengths(3:end));
0227 
0228     figure_params = nan(dimensions_varied - 2, 1);
0229 
0230     vary_lengths_cp = cumprod(vary_lengths);
0231 
0232     for f = 1:no_figures
0233       figure_params(1) = vary_vectors{3}(mod(f - 1, vary_lengths(3)) + 1);
0234 
0235       for v = 4:dimensions_varied
0236         figure_params(v - 2) = vary_vectors{v}(ceil(f/vary_lengths_cp(v - 3)));
0237       end
0238 
0239       figure_data_index = ones(size(vary_params, 1), 1);
0240       for v = 3:dimensions_varied
0241         figure_data_index = figure_data_index & vary_params(:, v) == figure_params(v - 2);
0242       end
0243 
0244       vary_title = '';
0245 
0246       for v = 1:(dimensions_varied - 2)
0247         vary_title = [vary_title, sprintf('%s = %f ', vary_labels{v + 2}, figure_params(v))];
0248       end
0249 
0250       handles = dsPlot(data(figure_data_index), varargin{:});
0251       for h = 1:length(handles)
0252         % figure(handles(h))
0253         mtit(handles(h), vary_title, 'FontSize', 14, 'yoff', .2)
0254       end
0255 
0256     end % no_figures
0257     return
0258   end % dimensions_varied
0259 end
0260 
0261 data=dsCheckData(data, varargin{:});
0262 handles=[];
0263 
0264 lock_gca = options.lock_gca;
0265 
0266 % TODO: add option 'plot_mode' {'trace','image'}
0267 
0268 % variables to plot
0269 var_fields=dsSelectVariables(data(1).labels,options.variable, varargin{:});
0270 tmp=regexp(var_fields,'_(.+)$','tokens','once');
0271 variables=unique([tmp{:}]);
0272 
0273 % populations to plot
0274 pop_names={data(1).model.specification.populations.name}; % list of populations
0275 % restrict to populations with variables to plot
0276 pop_indices=[];     % indices of populations to plot
0277 pop_var_indices={}; % indices of var_fields to plot per population
0278   % pop_var_indices{pop}(variable): index into var_fields
0279 for i=1:length(pop_names)
0280   % do any variables start with this population name?
0281   var_inds=find(~cellfun(@isempty,regexp(var_fields,['^' pop_names{i}])));
0282   if any(var_inds)
0283     inds=cellfun(@(x)find(~cellfun(@isempty,regexp(var_fields(var_inds),['_' x '$'],'once'))),variables,'uni',0);
0284     varsel=~cellfun(@isempty,inds);
0285     fldind=unique([inds{:}]);
0286     pop_indices(end+1)=i;
0287     pop_var_indices{end+1}=nan(1,length(variables));
0288     pop_var_indices{end}(varsel)=var_inds(fldind);
0289   end
0290 end
0291 pop_names=pop_names(pop_indices);
0292 
0293 % data set info
0294 time=data(1).time; % time vector
0295 pop_sizes=[data(1).model.specification.populations(pop_indices).size];
0296 num_pops=length(pop_names); % number of populations to plot
0297 num_sims=length(data); % number of simulations
0298 num_vars=length(variables);
0299 num_labels=length(var_fields); % number of labels to plot
0300 num_times=length(time);
0301 % set x-axis limits
0302 if isempty(options.time_limits)
0303   options.time_limits=[min(time) max(time)];
0304 end
0305 
0306 % do any analysis if necessary and set x-data
0307 switch options.plot_type
0308   case 'waveform'   % plot VARIABLE
0309     xdata=time;
0310     xlab='time (ms)'; % x-axis label
0311   case 'power'      % plot VARIABLE_Power_SUA.Pxx
0312     if any(cellfun(@isempty,regexp(var_fields,'.*_Power_SUA$')))
0313       data=dsCalcPower(data,varargin{:});
0314     end
0315     xdata=data(1).([var_fields{1} '_Power_SUA']).frequency;
0316     xlab='frequency (Hz)'; % x-axis label
0317     % set default x-limits for power spectrum
0318     if isempty(options.xlim)
0319       options.xlim=[0 200]; % Hz
0320     end
0321   case {'rastergram','raster'} % raster VARIABLE_spike_times
0322     if any(cellfun(@isempty,regexp(var_fields,'.*_spike_times$')))
0323       spike_fields=cellfun(@(x)[x '_spikes'],var_fields,'uni',0);
0324       idx=cellfun(@(x)isfield(data,x),spike_fields);
0325       if any(idx)
0326         % get spike times from binary spike matrix
0327         inds=find(idx);
0328         for s=1:length(data) % sims
0329           for i=1:length(inds) % pops
0330             spike_fld=[var_fields{inds(i)} '_spikes'];
0331             spike_time_fld=[var_fields{inds(i)} '_spike_times'];
0332             for j=1:size(data(s).(spike_fld),2) % cells
0333               data(s).(spike_time_fld){j}=data(s).time((1==data(s).(spike_fld)(:,j)));
0334             end % cells
0335           end % pops
0336         end % sims
0337 %       rmfields=cellfun(@(x)[x '_spikes'],var_fields,'uni',0);
0338 %       idx=cellfun(@(x)isfield(data,x),rmfields);
0339 %       if any(idx)
0340 %          data=rmfield(data,rmfields);
0341 %          data.labels=setdiff(data.labels,rmfields,'stable');
0342 %       end
0343       else
0344         % find spikes from threshold crossing
0345         data=dsCalcFR(data,varargin{:});
0346       end
0347     end
0348     xdata=time;
0349     xlab='time (ms)'; % x-axis label
0350   case 'rates'      % plot VARIABLE_FR
0351     if any(cellfun(@isempty,regexp(var_fields,'.*_FR$')))
0352       data=dsCalcFR(data,varargin{:});
0353     end
0354     xdata=data.time_FR;
0355     xlab='time (ms, bins)'; % x-axis label
0356 end
0357 if isempty(options.xlim)
0358   options.xlim=[min(xdata) max(xdata)];
0359 end
0360 
0361 MRPF = options.max_num_rows; % max rows per fig
0362 MTPP = options.max_num_overlaid; % max traces per plot
0363 
0364 % how many plots:
0365 if num_sims==1 && num_pops==1 && num_vars==1 && ~lock_gca
0366   num_fig_sets=1; num_figs=ceil(pop_sizes/MRPF); num_rows=min(pop_sizes,MRPF);
0367 elseif num_sims==1 && num_pops==1 && num_vars==1 && lock_gca
0368   num_fig_sets=1; num_figs=1; num_rows=1;
0369 elseif num_sims==1 && num_pops==1 && num_vars>1
0370   num_fig_sets=1; num_figs=ceil(num_vars/MRPF); num_rows=min(num_vars,MRPF);
0371 elseif num_sims==1 && num_pops>1 && num_vars==1 && ~lock_gca
0372   num_fig_sets=1; num_figs=ceil(num_pops/MRPF); num_rows=min(num_pops,MRPF);
0373 elseif num_sims==1 && num_pops>1 && num_vars==1 && lock_gca
0374   num_fig_sets=1; num_figs=ceil(num_pops/MRPF); num_rows=1;
0375 elseif num_sims==1 && num_pops>1 && num_vars>1
0376   num_fig_sets=num_vars; num_figs=ceil(num_pops/MRPF); num_rows=min(num_pops,MRPF);
0377 elseif num_sims>1 && num_pops==1 && num_vars==1
0378   num_fig_sets=1; num_figs=ceil(num_sims/MRPF); num_rows=min(num_sims,MRPF);
0379 elseif num_sims>1 && num_pops==1 && num_vars>1
0380   num_fig_sets=num_vars; num_figs=ceil(num_sims/MRPF); num_rows=min(num_sims,MRPF);
0381 elseif num_sims>1 && num_pops>1 && num_vars==1
0382   num_fig_sets=1; num_figs=ceil(num_sims/MRPF); num_rows=min(num_sims,MRPF);
0383 elseif num_sims>1 && num_pops>1 && num_vars>1
0384   num_fig_sets=num_vars; num_figs=ceil(num_sims/MRPF); num_rows=min(num_sims,MRPF);
0385 else
0386   error('unrecognized dimensions');
0387 end
0388 
0389 % If are doing rastergram, pops can be greater than 1 when doing lock_gca
0390 if lock_gca && (num_sims>1 || num_vars>1)
0391     error('Option lock_gca cannot with more than one simulation or variable');
0392 end
0393 
0394 % make subplot adjustments for varied parameters
0395 if num_sims>1 && isfield(data,'varied')
0396   % collect info on parameters varied
0397   varied=data(1).varied;
0398   num_varied=length(varied); % number of model components varied across simulations
0399   num_sims=length(data); % number of data sets (one per simulation)
0400   % collect info on parameters varied
0401   param_mat=zeros(num_sims,num_varied); % values for each simulation
0402   param_cell=cell(1,num_varied); % unique values for each parameter
0403   % loop over varied components and collect values
0404   for j=1:num_varied
0405     if isnumeric(data(1).(varied{j}))
0406       param_mat(:,j)=[data.(varied{j})]; % values for each simulation
0407       param_cell{j}=unique([data.(varied{j})]); % unique values for each parameter
0408     else
0409       % TODO: handle sims varying non-numeric model components
0410       % (eg, mechanisms) (also in dsPlotFR and dsSelect)
0411     end
0412   end
0413   param_size=cellfun(@length,param_cell); % number of unique values for each parameter
0414   % varied parameter with most elements goes along the rows (everything else goes along columns)
0415   row_param_index=find(param_size==max(param_size),1,'first');
0416   row_param_name=varied{row_param_index};
0417   row_param_values=param_cell{row_param_index};
0418   num_rows=length(row_param_values);
0419   %num_cols=num_sims/num_rows;
0420   num_figs=ceil(num_rows/MRPF);
0421   % collect sims for each value of the row parameter
0422   indices={};
0423   for row=1:num_rows
0424     indices{row}=find(param_mat(:,row_param_index)==row_param_values(row));
0425   end
0426   num_per_row=cellfun(@length,indices);
0427   num_cols=max(num_per_row);
0428   sim_indices=nan(num_cols,num_rows);
0429   % arrange sim indices for each row in a matrix
0430   for row=1:num_rows
0431     sim_indices(1:num_per_row(row),row)=indices{row};
0432   end
0433 %   sim_indices=[];
0434 %   for row=1:num_rows
0435 %     sim_indices=[sim_indices find(param_mat(:,row_param_index)==row_param_values(row))];
0436 %   end
0437 else
0438   sim_indices=ones(1,num_rows); % index into data array
0439   num_cols=1;
0440 end
0441 
0442 max_legend_entries=10;
0443 for figset=1:num_fig_sets
0444   for fig=1:num_figs
0445     ylims=[nan nan];
0446     % create figure
0447     if ~lock_gca
0448         handles(end+1)=figure('units','normalized','outerposition',[0 0 options.figwidth, options.figheight],'visible',options.visible);
0449         % position axes
0450         haxes=tight_subplot(num_rows,num_cols,[.01 .03],[.05 .01],[.03 .01]);
0451     else
0452         handles = gcf;
0453         haxes = gca;
0454     end
0455 
0456     axis_counter=0;
0457     AuxData=[];
0458     vlines=[];
0459     allspikes={};
0460     text_string=''; % string to add to subplot (set below)
0461     legend_strings=''; % legend for subplot (set below)
0462     shared_ylims_flag=1;
0463     % draw plots
0464     for row=1:num_rows
0465       for col=1:num_cols
0466         dat=[];
0467         sim_index=sim_indices(col,row); % index into data array for this subplot
0468         axis_counter=axis_counter+1; % number subplot axis we're on
0469         if isnan(sim_index)
0470           continue;
0471         end
0472         % #################################################################
0473         % what to plot
0474         % -----------------------------------------------------------------
0475 
0476         if num_sims==1 && num_pops==1 && num_vars==1 && ~lock_gca
0477         % -----------------------------------------------------------------
0478           % one cell per row: dat = data(s=1).(var)(:,c=r) where var=vars{v=1}
0479           var=var_fields{1};
0480           switch options.plot_type
0481             case 'waveform'
0482               dat=data(sim_index).(var)(:,row);
0483             case 'power'
0484               AuxData=data(sim_index).([var '_Power_MUA']).Pxx;
0485               vlines=data(sim_index).([var '_Power_MUA']).PeakFreq;
0486               AuxDataName={'MUA Power'};
0487               var=[var '_Power_SUA'];
0488               dat=data(sim_index).(var).Pxx(:,row);
0489               legend_strings={'SUA','MUA'};
0490             case {'rastergram','raster'}
0491               set_name=regexp(var,'^([a-zA-Z0-9]+)_','tokens','once');
0492               allspikes{1}{1}=data(sim_index).([var '_spike_times']){row};
0493                 % one pop, cell array of spike times for each cell in population
0494           end
0495           if num_rows>1
0496             text_string{row,col}=sprintf('cell %g',row);
0497           end
0498 
0499         elseif num_sims==1 && num_pops==1 && num_vars==1 && lock_gca
0500           % one population per row: dat = data(s=1).(var)(:,1:MTPP) where var=vars{v=r}
0501           var=var_fields{1};
0502           switch options.plot_type
0503             case 'waveform'
0504               dat=data(sim_index).(var);
0505             case 'power'
0506               AuxData=data(sim_index).([var '_Power_MUA']).Pxx;
0507               vlines=data(sim_index).([var '_Power_MUA']).PeakFreq;
0508               AuxDataName={'MUA Power'};
0509               var=[var '_Power_SUA'];
0510               dat=data(sim_index).(var).Pxx;
0511             case {'rastergram','raster'}
0512               set_name=regexp(var,'^([a-zA-Z0-9]+)_','tokens','once');
0513               allspikes{1}=data(sim_index).([var '_spike_times']);
0514           end
0515         % -----------------------------------------------------------------
0516         elseif num_sims==1 && num_pops==1 && num_vars>1
0517         % -----------------------------------------------------------------
0518           % one variable per row: dat = data(s=1).(var)(:,1:MTPP) where var=vars{v=r}
0519           var=var_fields{row};
0520           switch options.plot_type
0521             case 'waveform'
0522               dat=data(sim_index).(var);
0523             case 'power'
0524               AuxData=data(sim_index).([var '_Power_MUA']).Pxx;
0525               vlines=data(sim_index).([var '_Power_MUA']).PeakFreq;
0526               AuxDataName={'MUA Power'};
0527               var=[var '_Power_SUA'];
0528               dat=data(sim_index).(var).Pxx;
0529             case {'rastergram','raster'}
0530               set_name=regexp(var,'^([a-zA-Z0-9]+)_','tokens','once');
0531               allspikes{1}=data(sim_index).([var '_spike_times']);
0532           end
0533           shared_ylims_flag=0;
0534         % -----------------------------------------------------------------
0535         elseif num_sims==1 && num_pops>1 && num_vars==1 && ~lock_gca
0536         % -----------------------------------------------------------------
0537           % one population per row: dat = data(s=1).(var)(:,1:MTPP) where var=vars{v=r}
0538           var=var_fields{row};
0539           switch options.plot_type
0540             case 'waveform'
0541               dat=data(sim_index).(var);
0542             case 'power'
0543               AuxData=data(sim_index).([var '_Power_MUA']).Pxx;
0544               vlines=data(sim_index).([var '_Power_MUA']).PeakFreq;
0545               AuxDataName={'MUA Power'};
0546               var=[var '_Power_SUA'];
0547               dat=data(sim_index).(var).Pxx;
0548             case {'rastergram','raster'}
0549               set_name=regexp(var,'^([a-zA-Z0-9]+)_','tokens','once');
0550               allspikes{1}=data(sim_index).([var '_spike_times']);
0551           end
0552         % -----------------------------------------------------------------
0553         elseif num_sims==1 && num_pops>1 && num_vars==1 && lock_gca
0554         % -----------------------------------------------------------------
0555           % one simulation per row, overlay pops: dat = <data(s=r).(var)(:,1:MTPP),2|vars>
0556           switch options.plot_type
0557             case 'waveform'
0558               % calculate averages across populations
0559               dat=nan(num_times,num_pops);
0560               if ~strcmp(reportUI,'matlab') && exist('nanmean') ~= 2 % 'nanmean is not in Octave's path
0561                 try
0562                   pkg load statistics; % trying to load octave forge 'statistics' package before using nanmean function
0563                 catch
0564                   error('nanmean function is needed, please install the statistics package from Octave Forge');
0565                 end
0566               end
0567               for k=1:num_pops
0568                 dat(:,k)=nanmean(data(sim_index).(var_fields{k}),2);
0569               end
0570               var=['<' variables{1} '>'];
0571             case 'power'
0572               dat=nan(length(xdata),num_pops);
0573               AuxData=nan(length(xdata),num_pops);
0574               AuxDataName={}; vlines=[];
0575               if ~strcmp(reportUI,'matlab') && exist('nanmean') ~= 2 % 'nanmean is not in Octave's path
0576                 try
0577                   pkg load statistics; % trying to load octave forge 'statistics' package before using nanmean function
0578                 catch
0579                   error('nanmean function is needed, please install the statistics package from Octave Forge');
0580                 end
0581               end
0582               for k=1:num_pops
0583                 dat(:,k)=nanmean(data(sim_index).([var_fields{k} '_Power_SUA']).Pxx,2);
0584                 AuxData(:,k)=data(sim_index).([var_fields{k} '_Power_MUA']).Pxx;
0585                 AuxDataName{end+1}=strrep([var_fields{k} '_Power_MUA'],'_','\_');
0586                 vlines(end+1)=data(sim_index).([var_fields{k} '_Power_MUA']).PeakFreq;
0587               end
0588               var=['<' variables{1} '_Power_SUA>'];
0589             case {'rastergram','raster'}
0590               set_name={};
0591               for k=1:num_pops
0592                 tmp=regexp(var_fields{k},'^([a-zA-Z0-9]+)_','tokens','once');
0593                 set_name{k}=tmp{1};
0594                 allspikes{k}=data(sim_index).([var_fields{k} '_spike_times']);
0595               end
0596               var=['<' variables{1} '>'];
0597 
0598           end
0599         % -----------------------------------------------------------------
0600         elseif num_sims==1 && num_pops>1 && num_vars>1
0601         % -----------------------------------------------------------------
0602           % one population per row: dat = data(s=1).(var)(:,1:MTPP) where var=vars{these(p=r)}
0603           if isnan(pop_var_indices{row}(figset))
0604             continue;
0605           end
0606           var=var_fields{pop_var_indices{row}(figset)};
0607           switch options.plot_type
0608             case 'waveform'
0609               dat=data(sim_index).(var);
0610             case 'power'
0611               AuxData=data(sim_index).([var '_Power_MUA']).Pxx;
0612               vlines=data(sim_index).([var '_Power_MUA']).PeakFreq;
0613               AuxDataName={'MUA Power'};
0614               var=[var '_Power_SUA'];
0615               dat=data(sim_index).(var).Pxx;
0616             case {'rastergram','raster'}
0617               set_name=regexp(var,'^([a-zA-Z0-9]+)_','tokens','once');
0618               allspikes{1}=data(sim_index).([var '_spike_times']);
0619           end
0620         % -----------------------------------------------------------------
0621         elseif num_sims>1 && num_pops==1 && num_vars==1
0622         % -----------------------------------------------------------------
0623           % one simulation per row: dat = data(s=r).(var)(:,1:MTPP) where var=vars{v=1}
0624           var=var_fields{1};
0625           switch options.plot_type
0626             case 'waveform'
0627               dat=data(sim_index).(var);
0628             case 'power'
0629               AuxData=data(sim_index).([var '_Power_MUA']).Pxx;
0630               vlines=data(sim_index).([var '_Power_MUA']).PeakFreq;
0631               AuxDataName={'MUA Power'};
0632               var=[var '_Power_SUA'];
0633               dat=data(sim_index).(var).Pxx;
0634             case {'rastergram','raster'}
0635               set_name=regexp(var,'^([a-zA-Z0-9]+)_','tokens','once');
0636               allspikes{1}=data(sim_index).([var '_spike_times']);
0637           end
0638         % -----------------------------------------------------------------
0639         elseif num_sims>1 && num_pops==1 && num_vars>1
0640         % -----------------------------------------------------------------
0641           % one simulation per row: dat = data(s=r).(var)(:,1:MTPP) where var=vars{v++}
0642           if isnan(pop_var_indices{1}(figset))
0643             continue;
0644           end
0645           var=var_fields{pop_var_indices{1}(figset)};
0646           switch options.plot_type
0647             case 'waveform'
0648               dat=data(sim_index).(var);
0649             case 'power'
0650               AuxData=data(sim_index).([var '_Power_MUA']).Pxx;
0651               vlines=data(sim_index).([var '_Power_MUA']).PeakFreq;
0652               AuxDataName={'MUA Power'};
0653               var=[var '_Power_SUA'];
0654               dat=data(sim_index).(var).Pxx;
0655             case {'rastergram','raster'}
0656               set_name=regexp(var,'^([a-zA-Z0-9]+)_','tokens','once');
0657               allspikes{1}=data(sim_index).([var '_spike_times']);
0658           end
0659         % -----------------------------------------------------------------
0660         elseif num_sims>1 && num_pops>1 && num_vars==1
0661         % -----------------------------------------------------------------
0662           % one simulation per row, overlay pops: dat = <data(s=r).(var)(:,1:MTPP),2|vars>
0663           switch options.plot_type
0664             case 'waveform'
0665               % calculate averages across populations
0666               dat=nan(num_times,num_pops);
0667               if ~strcmp(reportUI,'matlab') && exist('nanmean') ~= 2 % 'nanmean is not in Octave's path
0668                 try
0669                   pkg load statistics; % trying to load octave forge 'statistics' package before using nanmean function
0670                 catch
0671                   error('nanmean function is needed, please install the statistics package from Octave Forge');
0672                 end
0673               end
0674               for k=1:num_pops
0675                 dat(:,k)=nanmean(data(sim_index).(var_fields{k}),2);
0676               end
0677               var=['<' variables{1} '>'];
0678             case 'power'
0679               dat=nan(length(xdata),num_pops);
0680               AuxData=nan(length(xdata),num_pops);
0681               AuxDataName={}; vlines=[];
0682               if ~strcmp(reportUI,'matlab') && exist('nanmean') ~= 2 % 'nanmean is not in Octave's path
0683                 try
0684                   pkg load statistics; % trying to load octave forge 'statistics' package before using nanmean function
0685                 catch
0686                   error('nanmean function is needed, please install the statistics package from Octave Forge');
0687                 end
0688               end
0689               for k=1:num_pops
0690                 dat(:,k)=nanmean(data(sim_index).([var_fields{k} '_Power_SUA']).Pxx,2);
0691                 AuxData(:,k)=data(sim_index).([var_fields{k} '_Power_MUA']).Pxx;
0692                 AuxDataName{end+1}=strrep([var_fields{k} '_Power_MUA'],'_','\_');
0693                 vlines(end+1)=data(sim_index).([var_fields{k} '_Power_MUA']).PeakFreq;
0694               end
0695               var=['<' variables{1} '_Power_SUA>'];
0696             case {'rastergram','raster'}
0697               set_name={};
0698               for k=1:num_pops
0699                 tmp=regexp(var_fields{k},'^([a-zA-Z0-9]+)_','tokens','once');
0700                 set_name{k}=tmp{1};
0701                 allspikes{k}=data(sim_index).([var_fields{k} '_spike_times']);
0702               end
0703               var=['<' variables{1} '>'];
0704 
0705           end
0706         % -----------------------------------------------------------------
0707         elseif num_sims>1 && num_pops>1 && num_vars>1
0708         % -----------------------------------------------------------------
0709           % one simulation per row, overlay pops: dat = <data(s=r).(var)(:,1:MTPP),2|vars(these)>
0710           switch options.plot_type
0711             case 'waveform'
0712               % calculate averages across populations
0713               dat=nan(num_times,num_pops);
0714               if ~strcmp(reportUI,'matlab') && exist('nanmean') ~= 2 % 'nanmean is not in Octave's path
0715                 try
0716                   pkg load statistics; % trying to load octave forge 'statistics' package before using nanmean function
0717                 catch
0718                   error('nanmean function is needed, please install the statistics package from Octave Forge');
0719                 end
0720               end
0721               for k=1:num_pops
0722                 if isnan(pop_var_indices{k}(figset))
0723                   continue;
0724                 end
0725                 var=var_fields{pop_var_indices{k}(figset)};
0726                 dat(:,k)=nanmean(data(sim_index).(var),2);
0727               end
0728               var=['<' variables{figset} '>'];
0729             case 'power'
0730               dat=nan(length(xdata),num_pops);
0731               AuxData=nan(length(xdata),num_pops);
0732               AuxDataName={}; vlines=[];
0733               if ~strcmp(reportUI,'matlab') && exist('nanmean') ~= 2 % 'nanmean is not in Octave's path
0734                 try
0735                   pkg load statistics; % trying to load octave forge 'statistics' package before using nanmean function
0736                 catch
0737                   error('nanmean function is needed, please install the statistics package from Octave Forge');
0738                 end
0739               end
0740               for k=1:num_pops
0741                 if isnan(pop_var_indices{k}(figset))
0742                   continue;
0743                 end
0744                 var=var_fields{pop_var_indices{k}(figset)};
0745                 dat(:,k)=nanmean(data(sim_index).([var '_Power_SUA']).Pxx,2);
0746                 AuxData(:,k)=data(sim_index).([var '_Power_MUA']).Pxx;
0747                 AuxDataName{end+1}=strrep([var '_Power_MUA'],'_','\_');
0748                 vlines(end+1)=data(sim_index).([var '_Power_MUA']).PeakFreq;
0749               end
0750               var=['<' variables{figset} '_Power_SUA>'];
0751             case {'rastergram','raster'}
0752               set_name={};
0753               for k=1:num_pops
0754                 if isnan(pop_var_indices{k}(figset))
0755                   continue;
0756                 end
0757                 var=var_fields{pop_var_indices{k}(figset)};
0758                 tmp=regexp(var,'^([a-zA-Z0-9]+)_','tokens','once');
0759                 set_name{k}=tmp{1};
0760                 allspikes{k}=data(sim_index).([var '_spike_times']);
0761               end
0762               var=['<' variables{figset} '>'];
0763           end
0764         end
0765         % #################################################################
0766         if size(dat,2)>1
0767           legend_strings=cellfun(@(x)['cell ' num2str(x)],num2cell(1:min(size(dat,2),max_legend_entries)),'uni',0);
0768         end
0769 
0770         if isfield(data,'varied')
0771           if num_sims>1
0772             % list the parameter varied along the rows first
0773             str=[row_param_name '=' num2str(row_param_values(row)) ': '];
0774             for k=1:num_varied
0775               fld=data(sim_index).varied{k};
0776               if ~strcmp(fld,row_param_name)
0777                 val=data(sim_index).(fld);
0778                 str=[str fld '=' num2str(val) ', '];
0779               end
0780             end
0781             if num_pops>1
0782               legend_strings=cellfun(@(x)[x ' (mean)'],pop_names,'uni',0);
0783             end
0784           else
0785             str='';
0786             for k=1:length(data.varied)
0787               fld=data(sim_index).varied{k};
0788               str=[str fld '=' num2str(data(sim_index).(fld)) ', '];
0789             end
0790           end
0791           text_string{row,col}=['(' strrep(str(1:end-2),'_','\_') ')'];
0792         end
0793         if ~isempty(AuxData) && length(legend_strings)<=max_legend_entries
0794           legend_strings=cat(2,legend_strings,AuxDataName);
0795         end
0796         % plot data
0797         %axes(haxes(axis_counter));
0798         set(gcf,'CurrentAxes',haxes(axis_counter));
0799         switch options.plot_type
0800           case {'waveform','power'}
0801             % finish preparing data
0802             if ~strcmp(options.yscale,'linear')
0803               dat=feval(options.yscale,dat); % log or log10
0804               % alternative approach: use semilogy for log10
0805             end
0806             if length(options.xlim)==2
0807               sel=(xdata>=options.xlim(1)&xdata<=options.xlim(2));
0808             else
0809               sel=1:length(xdata);
0810             end
0811             % plot traces
0812             if strcmp(options.plot_mode,'trace')
0813               % select max subset allowed
0814               dat=dat(:,1:min(size(dat,2),MTPP)); % select max subset to plot
0815               plot(xdata(sel),dat(sel,:));
0816               set(gca,'ticklength',get(gca,'ticklength')/2) %make ticks shorter
0817             else
0818               imagesc(dat);
0819             end
0820           case {'rastergram','raster'}
0821             % draw spikes
0822             ypos=0; % y-axis position tracker
0823             yticks=[]; % where to position population names
0824             yticklabels={}; % population names
0825             for p=1:length(allspikes) % loop over populations
0826               spikes=allspikes{p}; % spikes for one population
0827               for c=1:length(spikes) % loop over cells in population p
0828                 spks=spikes{c}; % spikes for one cell
0829                 for k=1:length(spks) % loop over spikes for cell c
0830                   spk=spks(k); % time of spike k in cell c of population p
0831                   line([spk spk],[c+ypos-.5 c+ypos+.5],'color','k'); hold on
0832                 end
0833               end
0834               % record position for population tick name
0835               yticks(end+1)=ypos+c/2+.5;
0836               yticklabels{end+1}=set_name{p};
0837               % draw line separating populations
0838               if length(allspikes)>1
0839                 pos=c+ypos+.5;
0840                 line([min(time) max(time)],[pos pos],'color','k','linewidth',3);
0841                 if p<length(allspikes)
0842                   % increment y-position for next population
0843                   ypos=ypos+c;
0844                 end
0845               end
0846             end
0847             % artificially set "dat" to get correct ylims below
0848             dat=[.5 ypos+c+.5];
0849             shared_ylims_flag=0;
0850             legend_strings='';
0851             % set y-ticks to population names
0852             set(gca,'ytick',yticks,'yticklabel',yticklabels);
0853             % set x-ticks
0854             plot([min(time) max(time)],[.5 .5],'w');
0855             nticks=length(get(gca,'xtick'));
0856             xticks=linspace(options.xlim(1),options.xlim(2),nticks);
0857             set(gca,'xtick',xticks,'xticklabel',xticks);
0858             %set(gca,'xticklabel',get(gca,'ytick'));
0859             set(gca,'ticklength',get(gca,'ticklength')/2) %make ticks shorter
0860         end % end switch options.plot_type
0861         % plot auxiliary data
0862         if ~isempty(AuxData) %strcmp(options.plot_type,'power')
0863           hold on
0864           plot(xdata(sel),AuxData(sel,:),'-','linewidth',3);%,'o-','linewidth',3);%'--.');
0865         end
0866         % format axes
0867         if row==num_rows
0868           xlabel(xlab);
0869         else
0870           set(haxes(axis_counter),'XTickLabel','');
0871           %set(haxes(row),'YTickLabel','');
0872         end
0873         xlim(options.xlim);
0874         if ~strcmp(options.plot_type,'rastergram')
0875           ylabel(strrep(var,'_','\_'));
0876         end
0877         if ~isempty(options.ylim)
0878           ylims=options.ylim;
0879         elseif shared_ylims_flag
0880           % update max/min
0881           ylims(1)=min(ylims(1),min(dat(:)));
0882           ylims(2)=max(ylims(2),max(dat(:)));
0883         else
0884           % set ylim to max/min of this data set
0885           ylims=[min(dat(:)) max(dat(:))];
0886           if ylims(1)~=ylims(2)
0887             ylim(ylims);
0888           end
0889           % add text
0890           if ~isempty(text_string)
0891             xmin=min(xlim); xmax=max(xlim);
0892             ymin=min(dat(:)); ymax=max(dat(:));
0893             text_xpos=xmin+.05*(xmax-xmin);
0894             text_ypos=ymin+.9*(ymax-ymin);
0895             try
0896               if any(strcmp(options.plot_type, {'rastergram','raster'}))
0897                 xlims = double(get(gca,'xlim'));
0898                 ylims = double(get(gca,'ylim'));
0899                 text(0.05*xlims(end),0.9*ylims(end),text_string{row,col});
0900               else
0901                 text(text_xpos,text_ypos,text_string{row,col});
0902               end
0903             end
0904           end
0905         end
0906         % plot lines and text (used for power)
0907         if ~isempty(vlines)
0908           for k=1:length(vlines)
0909             if ~isnan(vlines(k))
0910               line([vlines(k) vlines(k)],ylim,'color','k','linestyle','--');
0911               ymax=max(ylim);
0912               text(double(vlines(k) + 0.1*range(xlim)), 0.9*ymax, sprintf('MUA Sxx Peak F: %.f', vlines(k)))
0913             end
0914           end
0915         end
0916         % add legend
0917         if ~isempty(legend_strings) && axis_counter==1
0918           legend(legend_strings);
0919         end
0920       end % end loop over subplot columns
0921     end % end loop over subplot rows
0922     % set y-limits to max/min over data in this figure
0923     if shared_ylims_flag || ~isempty(options.ylim)
0924       if ylims(1)~=ylims(2)
0925         set(haxes,'ylim',ylims);
0926       end
0927       if ~isempty(text_string)
0928         axis_counter=0;
0929         for row=1:num_rows
0930           for col=1:num_cols
0931             if ~ischar(text_string{row,col})
0932               continue;
0933             end
0934             axis_counter=axis_counter+1;
0935             %axes(haxes(axis_counter));
0936             set(gcf,'CurrentAxes',haxes(axis_counter));
0937             xmin=min(xlim); xmax=max(xlim);
0938             ymin=min(ylim); ymax=max(ylim);
0939             text_xpos=double(xmin+.05*(xmax-xmin));
0940             text_ypos=ymin+.9*(ymax-ymin);
0941             text(text_xpos,text_ypos,text_string{row,col});
0942           end
0943         end
0944       end
0945     end
0946 
0947     %link x axes
0948     if numel(haxes) > 1
0949       linkaxes(haxes, 'x')
0950     end
0951 
0952   end % end loop over figures in this set
0953 end % end loop over figure sets
0954 
0955 %% auto_gen_test_data_flag argout
0956 if options.auto_gen_test_data_flag
0957   argout = {handles}; % specific to this function
0958 
0959   dsUnitSaveAutoGenTestDir(argin, argout);
0960 end
0961 
0962 % 1 sim, 1 pop, 1 var (X)
0963 %     N=1        one fig, one row (plot var X)
0964 %     N>1        one row per cell (plot var X),             enough figs for all cells
0965 %     (nsims=1, num_pops=1, num_vars=1, pop_sizes>=1): num_fig_sets=1, num_figs=ceil(pop_sizes/MRPF), num_rows=min(pop_sizes,MRPF): row r: dat = data(s=1).(var)(:,c=r) where var=vars{v=1}
0966 %
0967 % 1 sim, 1 pop, >1 vars (X,Y,...)
0968 %     N=1        one row per var (plot cell 1),             enough figs for all vars
0969 %     N>1        one row per var (overlay cells),         enough figs for all vars
0970 %     (nsims=1, num_pops=1, num_vars>=1, pop_sizes>=1): num_fig_sets=1, num_figs=ceil(num_vars/MRPF), num_rows=min(num_vars,MRPF): row r: dat = data(s=1).(var)(:,1:MTPP) where var=vars{v=r}
0971 %
0972 % 1 sim, >1 pops, 1 var (X)
0973 %     all N=1        one row per pop (plot var X, cell 1),        enough figs for all pops
0974 %     any N>1        one row per pop (plot var X, overlay cells),     enough figs for all pops
0975 %     (nsims=1, num_pops>=1, num_vars=1, pop_sizes>=1): num_fig_sets=1, num_figs=ceil(num_pops/MRPF), num_rows=min(num_pops,MRPF): row r: dat = data(s=1).(var)(:,1:MTPP) where var=vars{v=r}
0976 %
0977 % 1 sim, >1 pops, >1 vars (X,Y,...)
0978 %     all N=1        one row per pop (plot var X, cell 1),         enough figs for all pops, separate figs for each var
0979 %     any N>1        one row per pop (plot var X, overlay cells),     enough figs for all pops, separate figs for each var
0980 %     (nsims=1, num_pops>=1, num_vars>=1, pop_sizes>=1): num_fig_sets=num_vars, num_figs=ceil(num_pops/MRPF), num_rows=min(num_pops,MRPF): row r: dat = data(s=1).(var)(:,1:MTPP) where var=vars{these(p=r)}
0981 %
0982 % >1 sim, 1 pop, 1 var (X)
0983 %     N=1        one row per sim (plot var X, cell 1),         enough figs for all sims
0984 %     N>1        one row per sim (plot var X, overlay cells),     enough figs for all sims
0985 %     (nsims>1, num_pops=1, num_vars=1, pop_sizes>=1): num_fig_sets=1, num_figs=ceil(nsims/MRPF), num_rows=min(nsims,MRPF): row r: dat = data(s=r).(var)(:,1:MTPP) where var=vars{v=1}
0986 %
0987 % >1 sim, 1 pop, >1 vars (X,Y,...)
0988 %     N=1        one row per sim (plot var X, cell 1),         enough figs for all sims, separate figs for each var
0989 %     N>1        one row per sim (plot var X, overlay cells),     enough figs for all sims, separate figs for each var
0990 %     (nsims>1, num_pops=1, num_vars=1, pop_sizes>=1): num_fig_sets=num_vars, num_figs=ceil(nsims/MRPF), num_rows=min(nsims,MRPF): row r: dat = data(s=r).(var)(:,1:MTPP) where var=vars{v++}
0991 %
0992 % >1 sim, >1 pops, 1 var (X)
0993 %     all N=1        one row per sim (plot var X, overlay pops),    enough figs for all sims
0994 %     any N>1        one row per sim (plot var <X>, overlay pops),    enough figs for all sims
0995 %     (nsims>1, num_pops=1, num_vars=1, pop_sizes>=1): num_fig_sets=1, num_figs=ceil(nsims/MRPF), num_rows=min(nsims,MRPF): row r: dat = <data(s=r).(var)(:,1:MTPP),2|vars>
0996 %
0997 % >1 sim, >1 pops, >1 vars (X,Y,...)
0998 %     all N=1        one row per sim (plot var X, overlay pops),    enough figs for all sims, separate figs for each var
0999 %     any N>1        one row per sim (plot var <X>, overlay pops),    enough figs for all sims, separate figs for each var
1000 %     (nsims>1, num_pops=1, num_vars=1, pop_sizes>=1): num_fig_sets=num_vars, num_figs=ceil(nsims/MRPF), num_rows=min(nsims,MRPF): row r: dat = <data(s=r).(var)(:,1:MTPP),2|vars(these)>
1001 
1002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5