dsPlotFR2

PLOTFR2 - plot spike rates in various ways depending on what data was provided.

 As with dsPlotFR, but has additional options for controlling output of SimStudy mode.

 Usage:
   dsPlotFR2(data,'option',value)

0001 function handles = dsPlotFR2(data,varargin)
0002 %PLOTFR2 - plot spike rates in various ways depending on what data was provided.
0003 %
0004 % As with dsPlotFR, but has additional options for controlling output of SimStudy mode.
0005 %
0006 % Usage:
0007 %   dsPlotFR2(data,'option',value)
0008 
0009 %
0010 % Inputs:
0011 %   - data: DynaSim data structure (see dsCheckData)
0012 %   - options: (same as dsCalcFR)
0013 %     'variable' : name of field containing data on which to calculate firing
0014 %                  rates (default: *_spikes or first variable in data.labels)
0015 %     'threshold': scalar threshold value for detecting events (default: 0)
0016 %     'bin_size' : size of temporal window over which to calculate rate [ms or
0017 %                  fraction of data set] (default: 5% of the data set)
0018 %     'bin_shift': how much to shift the bin before calculating rate again [ms
0019 %                  or fraction of data set] (default: 1% of the data set)
0020 %     'plot_type': options for sim study mode. Options include 'heatmap',
0021 %                  'heatmap_sorted' (default), 'meanFR,' 'meanFRdens', and
0022 %                  'summary'.
0023 %     'lock_gca' : Plots within currently active axis (gca); doesn't
0024 %                  open new figures or subplots.
0025 %
0026 %
0027 % Examples:
0028 % dsPlotFR(data,'bin_size',30,'bin_shift',10);
0029 %
0030 % TODO: add rastergrams
0031 %
0032 % See also: dsCalcFR, dsSimulate, dsCheckData
0033 
0034 data=dsCheckData(data, varargin{:});
0035 fields=fieldnames(data);
0036 handles=[];
0037 
0038 options=dsCheckOptions(varargin,{...
0039   'plot_type','heatmap_sorted',{'heatmap','heatmap_sorted','meanFR','meanFRdens','summary'},...
0040   'variable',[],[],...
0041   'threshold',1e-5,[],... % slightly above zero in case variable is point process *_spikes {0,1}
0042   'bin_size',.05,[],...  % 30
0043   'bin_shift',.01,[],... % 10
0044   'exclude_data_flag',0,{0,1},...
0045   'lock_gca',false,[true,false],...
0046   'visible','on',[],...
0047   },false);
0048 
0049 lock_gca = options.lock_gca;
0050 keyvals=dsOptions2Keyval(rmfield(options,{'plot_type'}));
0051 
0052 
0053 % calc firing rates if not already present in data
0054 if all(cellfun(@isempty,regexp(fields,'.*_FR$')))
0055   data=dsCalcFR(data,keyvals{:}); % equivalent: data=dsAnalyzeStudy(data,@dsCalcFR);
0056   fields=fieldnames(data);
0057 end
0058 % get list of fields with firing rate data
0059 FR_fields=fields(~cellfun(@isempty,regexp(fields,'.*_FR$')));
0060 FR_fields=setdiff(FR_fields,'time_FR');
0061 % store time bins for firing rate data
0062 time=data.time_FR;
0063 nsets=length(FR_fields);
0064 
0065 % plot firing rates
0066 if (numel(data)==1 || ~isfield(data,'varied')) && ~lock_gca
0067   % plot separate firing rates vs time (one figure per element of data)
0068   for i=1:length(data)
0069     plotFR_SingleSim(i);
0070   end  
0071 else
0072   % data contains results from a set of simulations varying something
0073   % plot average firing rates vs whatever was varied
0074   if strcmp(options.plot_type,'summary') && ~lock_gca
0075       plotFR_SimStudy;
0076   else
0077       if numel(data) > 1 && lock_gca
0078           error('Cannot lock gca if number of elements in data is greater than 1');
0079       end
0080       plotFR_SimStudy_specialized;
0081       
0082   end
0083 end
0084 
0085   % NESTED FUNCTIONS
0086   function plotFR_SingleSim(i)
0087     % purpose: plot each data set data.(FR_fields{k})
0088     % plots for N populations (FR data sets) from this simulation
0089     ht=320; % height per subplot row (=per population or FR data set)
0090     if ~lock_gca; handles(end+1)=figure('position',[250 max(50,600-(nsets-1)*ht) 1400 min(ht*nsets,750)],'visible',options.visible); end
0091     for k=1:nsets % index of firing rate data field
0092       dat=data(i).(FR_fields{k});
0093       bins=0:1.05*max(dat(:));
0094       rlims=[min(bins) max(bins)];
0095       if rlims(1)==rlims(2), rlims=rlims(1)+[0 1e-5]; end
0096       tlims=[min(time) max(time)];
0097       % get population name from field (assumes: pop_*)
0098       popname=regexp(FR_fields{k},'^([a-zA-Z0-9]+)_','tokens','once');
0099       popname=popname{1};
0100       ncells=size(dat,2);
0101       if ncells==1
0102         nc=4; % number of plots for this FR data set
0103         % 1.0 plot firing rate trace
0104         subplot(nsets,nc,(1:nc-1)+(k-1)*nc); % plot(t,FR)
0105         plot(time,dat(:,1),'o-','linewidth',2); ylim(rlims); xlim(tlims);
0106         title([popname ': firing rate over time']);
0107         xlabel('time (ms)'); ylabel([popname ': firing rate (Hz)']);
0108         % 2.0 plot firing rate density
0109         subplot(nsets,nc,nc+(k-1)*nc); % hist(FR)
0110         if numel(bins)>1
0111           H=hist(dat(:,1),bins)/length(dat(:,1));
0112           h=bar(bins,H); hold on
0113           set(get(h,'children'),'FaceAlpha',0,'EdgeAlpha',.4,'linewidth',2);
0114           rn=ksr(bins,H,.75,length(bins));
0115           plot(bins,rn.f,'linewidth',2); title([popname ': firing rate density']);
0116           xlabel([popname ': firing rate (Hz)']); ylabel('fraction');
0117           xlim(rlims); ylim([0 1]);
0118         end
0119       else
0120         nc=4; % number of plots for this FR data set
0121         % 1.0 plot firing rate heat map
0122         subplot(nsets,nc,1+(k-1)*nc); % imagesc(t,cells,FR)
0123         imagesc(time,1:ncells,dat'); axis xy
0124         title([popname ': firing rates (Hz)']);
0125         xlabel('time (ms)'); ylabel([popname ' cell index']);
0126         caxis(rlims); xlim(tlims);
0127         if ncells<=10
0128           ytick=1:ncells;
0129         else
0130           ytick=round(linspace(1,ncells,5));
0131           if ~ismember(ncells,ytick)
0132             ytick=[ytick ncells];
0133           end
0134         end
0135         set(gca,'ytick',ytick,'yticklabel',ytick);
0136         % 2.0 plot sorted firing rate heat map
0137         subplot(nsets,nc,2+(k-1)*nc); % imagesc(t,cells_sorted,FR)
0138         tmp=sum(dat,1);
0139         [tmp,inds]=sort(tmp);
0140         imagesc(time,1:ncells,dat(:,inds)'); axis xy
0141         caxis(rlims); xlim(tlims);
0142         title([popname ': firing rates (Hz)']);
0143         xlabel('time (ms)'); ylabel([popname ' cell index (sorted by FR)']);
0144         if ncells<=10
0145           ytick=1:ncells;
0146           yticklabel=ytick(inds);
0147         else
0148           ytick=[];
0149           yticklabel=[];
0150         end
0151         set(gca,'ytick',ytick,'yticklabel',yticklabel);
0152         % 3.0 plot population-average firing rate trace
0153         subplot(nsets,nc,3+(k-1)*nc); % plot(t,<FR|pop>)
0154         plot(time,mean(dat,2),'o-','linewidth',2);
0155         title([popname ': population average firing rate']);
0156         xlabel('time (ms)'); ylabel([popname ': avg firing rate (Hz)']);
0157         ylim(rlims); xlim(tlims);
0158         % 4.0 plot firing rate density
0159         if numel(bins)>1
0160           subplot(nsets,nc,4+(k-1)*nc); % hist(FR)
0161           H=hist(dat(:),bins)/numel(dat);
0162           h=bar(bins,H); hold on
0163           set(get(h,'children'),'FaceAlpha',0,'EdgeAlpha',.4,'linewidth',2);
0164           rn=ksr(bins,H,.75,length(bins));
0165           plot(bins,rn.f,'linewidth',2); title([popname ': firing rate density']);
0166           xlabel([popname ': firing rate (Hz)']); ylabel('fraction');
0167           xlim(rlims); ylim([0 1]);
0168         end
0169       end
0170     end
0171   end
0172 
0173   function plotFR_SimStudy
0174     % calculate average firing rates across population and time
0175     varied=data(1).varied;
0176     % eliminate parameters that had the same value in all sims or were non-numeric
0177     keep=zeros(size(varied));
0178     for v=1:length(varied)
0179       if isnumeric(data(1).(varied{v})) && length(unique([data.(varied{v})]))>1
0180         keep(v)=1;
0181       end
0182     end
0183     varied=varied(keep==1);
0184     nvaried=length(varied); % number of model components varied across simulations
0185     nsims=length(data);
0186     FRmu=zeros(nsims,nsets);
0187     for i=1:nsims % simulations
0188       for k=1:nsets % populations
0189         FRmu(i,k)=mean(data(i).(FR_fields{k})(:));
0190       end
0191     end
0192     % collect info on parameters varied
0193     params=zeros(nsims,nvaried);
0194     for j=1:nvaried
0195       if isnumeric(data(1).(varied{j}))
0196         params(:,j)=[data.(varied{j})];
0197       else
0198         % todo: handle sims varying non-numeric model components
0199         % (eg, mechanisms)
0200       end
0201     end
0202     % plots for N populations and M varied elements
0203     % plot how avg firing rate for each pop varies with each parameter
0204     ht=320; % height per subplot row (=per population or FR data set)
0205     wt=500;
0206     if ~lock_gca; handles(end+1)=figure('position',[250 max(50,600-(nsets-1)*ht) min(1500,500+(nvaried-1)*wt) min(ht*nsets,750)],'visible',options.visible); end
0207     cnt=0;
0208     for k=1:nsets % populations
0209       popname=regexp(FR_fields{k},'^([a-zA-Z0-9]+)_','tokens','once');
0210       popname=popname{1};
0211       ncells=size(data(1).(FR_fields{k}),2);
0212       for j=1:nvaried % varied parameters
0213         % calculate mean avg FRs for each value of this varied parameter
0214         pvals=unique(params(:,j)); % unique values for this parameter j
0215         nv=length(pvals); % number of values used for this parameter j
0216         rvals=zeros(nv,1);
0217         for v=1:nv
0218           idx=(params(:,j)==pvals(v)); % sims with param j = value v
0219           % average across all sims with param j set to value v
0220           rvals(v)=mean(FRmu(idx,k)); % avg pop k FR given this value
0221         end
0222         npoints=length(find(idx)); % number of points in this average
0223         % plot avg FRs
0224         cnt=cnt+1; % subplot index
0225         subplot(nsets,nvaried,cnt);
0226         plot(pvals,rvals,'o-','linewidth',2);
0227         % todo: add confidence intervals / error bars
0228         axis tight; %ylim([0 max(FRmu(:))]);
0229         xlabel(strrep(varied{j},'_','\_'));
0230         ylabel([popname ': avg firing rate (Hz)']);
0231         title([popname ': firing rate averaged over pop, time, sims']);
0232         % add text: ncells, length(time), npoints
0233         xpos=min(xlim)+.7*diff(xlim);%.05*diff(xlim);
0234         ypos=min(ylim)+.2*diff(ylim);%.85*diff(ylim);
0235         text(xpos,ypos,sprintf('per point:\n #cells=%g\n #bins=%g\n #sims=%g',ncells,length(time),npoints));
0236       end
0237     end
0238     if length(varied)==2
0239       % plots for N populations and 2 varied elements
0240       % organize and imagesc FRmu(param 1, param 2)
0241       if ~lock_gca; handles(end+1)=figure('position',[1150 max(50,600-(nsets-1)*ht) 500 min(ht*nsets,750)],'visible',options.visible);end
0242       pvals1=unique(params(:,1)); nv1=length(pvals1);
0243       pvals2=unique(params(:,2)); nv2=length(pvals2);
0244       for k=1:nsets
0245         popname=regexp(FR_fields{k},'^([a-zA-Z0-9]+)_','tokens','once');
0246         popname=popname{1};
0247         FRmu2=zeros(nv1,nv2);
0248         for i=1:nv1
0249           for j=1:nv2
0250             idx=(params(:,1)==pvals1(i))&(params(:,2)==pvals2(j));
0251             FRmu2(i,j)=mean(FRmu(idx,k));
0252           end
0253         end
0254         % plot imagesc(param1,param2,FRmu)
0255         subplot(nsets,1,k);
0256         imagesc(pvals1,pvals2,FRmu2'); axis xy
0257         xlabel(strrep(varied{1},'_','\_')); ylabel(strrep(varied{2},'_','\_'));
0258         title([popname ': avg firing rate (Hz)']);
0259         caxis([0 max(FRmu(:))]); colorbar
0260       end
0261     end
0262   end
0263     function plotFR_SimStudy_specialized    % Specifc plot for all sims, based on options.plot_type
0264         %%
0265         
0266         % New code (imported from dsPlot)
0267         num_sims=length(data); % number of simulations
0268         
0269         
0270         % make subplot adjustments for varied parameters
0271         if num_sims>1 && isfield(data,'varied')
0272           % collect info on parameters varied
0273           varied=data(1).varied;
0274           num_varied=length(varied); % number of model components varied across simulations
0275           num_sims=length(data); % number of data sets (one per simulation)
0276           % collect info on parameters varied
0277           param_mat=zeros(num_sims,num_varied); % values for each simulation
0278           param_cell=cell(1,num_varied); % unique values for each parameter
0279           % loop over varied components and collect values
0280           for j=1:num_varied
0281             if isnumeric(data(1).(varied{j}))
0282               param_mat(:,j)=[data.(varied{j})]; % values for each simulation
0283               param_cell{j}=unique([data.(varied{j})]); % unique values for each parameter
0284             else
0285               % todo: handle sims varying non-numeric model components
0286               % (eg, mechanisms) (also in dsPlotFR and dsSelect)
0287             end
0288           end
0289           param_size=cellfun(@length,param_cell); % number of unique values for each parameter
0290           % varied parameter with most elements goes along the rows (everything else goes along columns)
0291           row_param_index=find(param_size==max(param_size),1,'first');
0292           row_param_name=varied{row_param_index};
0293           row_param_values=param_cell{row_param_index};
0294           num_rows=length(row_param_values);
0295           %num_cols=num_sims/num_rows;
0296           num_figs=1;
0297           % collect sims for each value of the row parameter
0298           indices={};
0299           for row=1:num_rows
0300             indices{row}=find(param_mat(:,row_param_index)==row_param_values(row));
0301           end
0302           num_per_row=cellfun(@length,indices);
0303           num_cols=max(num_per_row);
0304           sim_indices=nan(num_cols,num_rows);
0305           % arrange sim indices for each row in a matrix
0306           for row=1:num_rows
0307             sim_indices(1:num_per_row(row),row)=indices{row};
0308           end
0309         %   sim_indices=[];
0310         %   for row=1:num_rows
0311         %     sim_indices=[sim_indices find(param_mat(:,row_param_index)==row_param_values(row))];
0312         %   end
0313         else
0314             num_rows = 1;
0315             sim_indices=ones(1,num_rows); % index into data array
0316             num_cols=1;
0317         end
0318         
0319         ht=320; % height per subplot row (=per population or FR data set)
0320         if ~lock_gca; handles(1) = figure('units','normalized','position',[0,1-min(.33*num_rows,1),min(.25*num_cols,1) min(.33*num_rows,1)],'visible',options.visible); end
0321         if ~lock_gca; hsp = subplot_grid(num_rows,num_cols);  end
0322         
0323         axis_counter = 0;
0324         for row=1:num_rows
0325             for col=1:num_cols
0326                 
0327                 sim_index=sim_indices(col,row); % index into data array for this subplot
0328                 axis_counter=axis_counter+1; % number subplot axis we're on
0329                 if isnan(sim_index)
0330                   continue;
0331                 end
0332                 
0333                 if ~lock_gca; hsp.set_gca(axis_counter); end
0334                 
0335                 num_pops = 1;
0336                 if isfield(data,'varied')
0337                   if num_sims>1
0338                     % list the parameter varied along the rows first
0339                     str=[row_param_name '=' num2str(row_param_values(row)) ': '];
0340                     for k=1:num_varied
0341                       fld=data(sim_index).varied{k};
0342                       if ~strcmp(fld,row_param_name)
0343                         val=data(sim_index).(fld);
0344                         str=[str fld '=' num2str(val) ', '];
0345                       end
0346                     end
0347                     if num_pops>1
0348                       legend_strings=cellfun(@(x)[x ' (mean)'],pop_names,'uni',0);
0349                     end
0350                   else
0351                     str='';
0352                     for k=1:length(data.varied)
0353                       fld=data(sim_index).varied{k};
0354                       str=[str fld '=' num2str(data(sim_index).(fld)) ', '];
0355                     end
0356                   end
0357                   text_string{row,col}=['(' strrep(str(1:end-2),'_','\_') ')'];
0358                 end
0359         
0360 
0361                 k=1;
0362                 i=sim_index;
0363                 dat=data(i).(FR_fields{k});
0364                 bins=0:1.05*max(dat(:));
0365                 rlims=[min(bins) max(bins)];
0366                 if rlims(1)==rlims(2), rlims=rlims(1)+[0 1e-5]; end
0367                 tlims=[min(time) max(time)];
0368                 % get population name from field (assumes: pop_*)
0369                 popname=regexp(FR_fields{k},'^([a-zA-Z0-9]+)_','tokens','once');
0370                 popname=popname{1};
0371                 ncells=size(dat,2);
0372                 nc=4; % number of plots for this FR data set
0373                 % 1.0 plot firing rate heat map
0374                 if strcmp(options.plot_type,'heatmap')
0375                     %subplot(nsets,nc,1+(k-1)*nc); % imagesc(t,cells,FR)
0376                     imagesc(time,1:ncells,dat'); axis xy
0377                     if ~lock_gca; hsp.figtitle([popname ': firing rates (Hz) ']); title(text_string{row,col}); end
0378                     if row == num_rows; xlabel('time (ms)'); end; ylabel([popname ' cell index']);
0379                     caxis(rlims); xlim(tlims);
0380                     if ncells<=10
0381                       ytick=1:ncells;
0382                     else
0383                       ytick=round(linspace(1,ncells,5));
0384                       if ~ismember(ncells,ytick)
0385                         ytick=[ytick ncells];
0386                       end
0387                     end
0388                     set(gca,'ytick',ytick,'yticklabel',ytick);
0389                 end
0390                 % 2.0 plot sorted firing rate heat map
0391                 if strcmp(options.plot_type,'heatmap_sorted')
0392                     %subplot(nsets,nc,2+(k-1)*nc); % imagesc(t,cells_sorted,FR)
0393                     tmp=sum(dat,1);
0394                     [tmp,inds]=sort(tmp);
0395                     imagesc(time,1:ncells,dat(:,inds)'); axis xy
0396                     caxis(rlims); xlim(tlims);
0397                     if ~lock_gca; hsp.figtitle([popname ': firing rates (Hz) ']); title(text_string{row,col}); end;
0398                     
0399                     if row == num_rows; xlabel('time (ms)'); end; ylabel([popname ' cell index (sorted by FR)']);
0400                     if ncells<=10
0401                       ytick=1:ncells;
0402                       yticklabel=ytick(inds);
0403                     else
0404                       ytick=[];
0405                       yticklabel=[];
0406                     end
0407                     set(gca,'ytick',ytick,'yticklabel',yticklabel);
0408                 end
0409                 % 3.0 plot population-average firing rate trace
0410                 if strcmp(options.plot_type,'meanFR')
0411                     %subplot(nsets,nc,3+(k-1)*nc); % plot(t,<FR|pop>)
0412                     plot(time,mean(dat,2),'o-','linewidth',2);
0413                     if ~lock_gca; hsp.figtitle([popname ': pop. avg FR']); title(text_string{row,col}); end
0414                     if row == num_rows; xlabel('time (ms)'); end; ylabel([popname ': avg firing rate (Hz)']);
0415                     ylim(rlims); xlim(tlims);
0416                 end 
0417                 % 4.0 plot firing rate density
0418                 if strcmp(options.plot_type,'meanFRdens')
0419                     if numel(bins)>1
0420                       %subplot(nsets,nc,4+(k-1)*nc); % hist(FR)
0421                       H=hist(dat(:),bins)/numel(dat);
0422                       h=bar(bins,H); hold on
0423                       set(get(h,'children'),'FaceAlpha',0,'EdgeAlpha',.4,'linewidth',2);
0424                       rn=ksr(bins,H,.75,length(bins));
0425                       plot(bins,rn.f,'linewidth',2); if ~lock_gca; hsp.figtitle([popname ': FR density']); end; title(text_string{row,col});
0426                       if row == num_rows; xlabel([popname ': firing rate (Hz)']); end; ylabel('fraction');
0427                       xlim(rlims); ylim([0 1]);
0428                     end
0429                 end
0430             end
0431         end
0432     end
0433 end