Description of dsPlotFR

0001 function handles = dsPlotFR(data,varargin)
0002 %PLOTFR - plot spike rates in various ways depending on what data was provided.
0003 %
0004 % Usage:
0005 %   dsPlotFR(data,'option',value)
0006 %
0007 % Inputs:
0008 %   - data: DynaSim data structure (see dsCheckData)
0009 %   - options: (same as dsCalcFR)
0010 %     'variable' : name of field containing data on which to calculate firing
0011 %                  rates (default: *_spikes or first variable in data.labels)
0012 %     'threshold': scalar threshold value for detecting events (default: 0)
0013 %     'bin_size' : size of temporal window over which to calculate rate [ms or
0014 %                  fraction of data set] (default: 5% of the data set)
0015 %     'bin_shift': how much to shift the bin before calculating rate again [ms
0016 %                  or fraction of data set] (default: 1% of the data set)
0017 %     'plot_type': options for sim study mode. Options include 'heatmap',
0018 %                  'heatmap_sorted' (default), 'meanFR,' 'meanFRdens', and
0019 %                  'summary'.
0020 %     'lock_gca' : Plots within currently active axis (gca); doesn't
0021 %                  open new figures or subplots.
0022 %
0023 % Examples:
0024 % dsPlotFR(data,'bin_size',30,'bin_shift',10);
0025 %
0026 % TODO: add rastergrams
0027 %
0028 % See also: dsCalcFR, dsSimulate, dsCheckData
0029 %
0030 % Author: Jason Sherfey, PhD <jssherfey@gmail.com>
0031 % Copyright (C) 2016 Jason Sherfey, Boston University, USA
0032 
0033 
0034 % Pull out lock_gca from varargin
0035 options=dsCheckOptions(varargin,{...
0036   'lock_gca',false,[true,false],...
0037   },false);
0038 lock_gca = options.lock_gca;
0039 
0040 data=dsCheckData(data, varargin{:});
0041 fields=fieldnames(data);
0042 handles=[];
0043 
0044 % calc firing rates if not already present in data
0045 if all(cellfun(@isempty,regexp(fields,'.*_FR$')))
0046   data=dsCalcFR(data,varargin{:}); % equivalent: data=dsAnalyzeStudy(data,@dsCalcFR);
0047   fields=fieldnames(data);
0048 end
0049 % get list of fields with firing rate data
0050 FR_fields=fields(~cellfun(@isempty,regexp(fields,'.*_FR$')));
0051 FR_fields=setdiff(FR_fields,'time_FR');
0052 % store time bins for firing rate data
0053 time=data.time_FR;
0054 nsets=length(FR_fields);
0055 
0056 % plot firing rates
0057 if numel(data)==1 || ~isfield(data,'varied')
0058   % plot separate firing rates vs time (one figure per element of data)
0059   for i=1:length(data)
0060     plotFR_SingleSim(i);
0061   end
0062 else
0063   % data contains results from a set of simulations varying something
0064   % plot average firing rates vs whatever was varied
0065   plotFR_SimStudy;
0066 end
0067 
0068   % NESTED FUNCTIONS
0069   function plotFR_SingleSim(i)
0070     % purpose: plot each data set data.(FR_fields{k})
0071     % plots for N populations (FR data sets) from this simulation
0072     ht=320; % height per subplot row (=per population or FR data set)
0073     if ~lock_gca; handles(end+1)=figure('position',[250 max(50,600-(nsets-1)*ht) 1400 min(ht*nsets,750)]); end
0074     for k=1:nsets % index of firing rate data field
0075       dat=data(i).(FR_fields{k});
0076       bins=0:1.05*max(dat(:));
0077       rlims=[min(bins) max(bins)];
0078       if rlims(1)==rlims(2), rlims=rlims(1)+[0 1e-5]; end
0079       tlims=[min(time) max(time)];
0080       % get population name from field (assumes: pop_*)
0081       popname=regexp(FR_fields{k},'^([a-zA-Z0-9]+)_','tokens','once');
0082       popname=popname{1};
0083       ncells=size(dat,2);
0084       if ncells==1
0085         nc=4; % number of plots for this FR data set
0086         % 1.0 plot firing rate trace
0087         subplot(nsets,nc,(1:nc-1)+(k-1)*nc); % plot(t,FR)
0088         plot(time,dat(:,1),'o-','linewidth',2); ylim(rlims); xlim(tlims);
0089         title([popname ': firing rate over time']);
0090         xlabel('time (ms)'); ylabel([popname ': firing rate (Hz)']);
0091         % 2.0 plot firing rate density
0092         subplot(nsets,nc,nc+(k-1)*nc); % hist(FR)
0093         if numel(bins)>1
0094           H=hist(dat(:,1),bins)/length(dat(:,1));
0095           h=bar(bins,H); hold on
0096           set(get(h,'children'),'FaceAlpha',0,'EdgeAlpha',.4,'linewidth',2);
0097           rn=ksr(bins,H,.75,length(bins));
0098           plot(bins,rn.f,'linewidth',2); title([popname ': firing rate density']);
0099           xlabel([popname ': firing rate (Hz)']); ylabel('fraction');
0100           xlim(rlims); ylim([0 1]);
0101         end
0102       else
0103         nc=4; % number of plots for this FR data set
0104         % 1.0 plot firing rate heat map
0105         subplot(nsets,nc,1+(k-1)*nc); % imagesc(t,cells,FR)
0106         imagesc(time,1:ncells,dat'); axis xy
0107         title([popname ': firing rates (Hz)']);
0108         xlabel('time (ms)'); ylabel([popname ' cell index']);
0109         caxis(rlims); xlim(tlims);
0110         if ncells<=10
0111           ytick=1:ncells;
0112         else
0113           ytick=round(linspace(1,ncells,5));
0114           if ~ismember(ncells,ytick)
0115             ytick=[ytick ncells];
0116           end
0117         end
0118         set(gca,'ytick',ytick,'yticklabel',ytick);
0119         % 2.0 plot sorted firing rate heat map
0120         subplot(nsets,nc,2+(k-1)*nc); % imagesc(t,cells_sorted,FR)
0121         tmp=sum(dat,1);
0122         [tmp,inds]=sort(tmp);
0123         imagesc(time,1:ncells,dat(:,inds)'); axis xy
0124         caxis(rlims); xlim(tlims);
0125         title([popname ': firing rates (Hz)']);
0126         xlabel('time (ms)'); ylabel([popname ' cell index (sorted by FR)']);
0127         if ncells<=10
0128           ytick=1:ncells;
0129           yticklabel=ytick(inds);
0130         else
0131           ytick=[];
0132           yticklabel=[];
0133         end
0134         set(gca,'ytick',ytick,'yticklabel',yticklabel);
0135         % 3.0 plot population-average firing rate trace
0136         subplot(nsets,nc,3+(k-1)*nc); % plot(t,<FR|pop>)
0137         plot(time,mean(dat,2),'o-','linewidth',2);
0138         title([popname ': population average firing rate']);
0139         xlabel('time (ms)'); ylabel([popname ': avg firing rate (Hz)']);
0140         ylim(rlims); xlim(tlims);
0141         % 4.0 plot firing rate density
0142         if numel(bins)>1
0143           subplot(nsets,nc,4+(k-1)*nc); % hist(FR)
0144           H=hist(dat(:),bins)/numel(dat);
0145           h=bar(bins,H); hold on
0146           set(get(h,'children'),'FaceAlpha',0,'EdgeAlpha',.4,'linewidth',2);
0147           rn=ksr(bins,H,.75,length(bins));
0148           plot(bins,rn.f,'linewidth',2); title([popname ': firing rate density']);
0149           xlabel([popname ': firing rate (Hz)']); ylabel('fraction');
0150           xlim(rlims); ylim([0 1]);
0151         end
0152       end
0153     end
0154   end
0155 
0156   function plotFR_SimStudy
0157     % calculate average firing rates across population and time
0158     varied=data(1).varied;
0159     % eliminate parameters that had the same value in all sims or were non-numeric
0160     keep=zeros(size(varied));
0161     for v=1:length(varied)
0162       if isnumeric(data(1).(varied{v})) && length(unique([data.(varied{v})]))>1
0163         keep(v)=1;
0164       end
0165     end
0166     varied=varied(keep==1);
0167     nvaried=length(varied); % number of model components varied across simulations
0168     nsims=length(data);
0169     FRmu=zeros(nsims,nsets);
0170     for i=1:nsims % simulations
0171       for k=1:nsets % populations
0172         FRmu(i,k)=mean(data(i).(FR_fields{k})(:));
0173       end
0174     end
0175     % collect info on parameters varied
0176     params=zeros(nsims,nvaried);
0177     for j=1:nvaried
0178       if isnumeric(data(1).(varied{j}))
0179         params(:,j)=[data.(varied{j})];
0180       else
0181         % todo: handle sims varying non-numeric model components
0182         % (eg, mechanisms)
0183       end
0184     end
0185     % plots for N populations and M varied elements
0186     % plot how avg firing rate for each pop varies with each parameter
0187     ht=320; % height per subplot row (=per population or FR data set)
0188     wt=500;
0189     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)]); end
0190     cnt=0;
0191     for k=1:nsets % populations
0192       popname=regexp(FR_fields{k},'^([a-zA-Z0-9]+)_','tokens','once');
0193       popname=popname{1};
0194       ncells=size(data(1).(FR_fields{k}),2);
0195       for j=1:nvaried % varied parameters
0196         % calculate mean avg FRs for each value of this varied parameter
0197         pvals=unique(params(:,j)); % unique values for this parameter j
0198         nv=length(pvals); % number of values used for this parameter j
0199         rvals=zeros(nv,1);
0200         for v=1:nv
0201           idx=(params(:,j)==pvals(v)); % sims with param j = value v
0202           % average across all sims with param j set to value v
0203           rvals(v)=mean(FRmu(idx,k)); % avg pop k FR given this value
0204         end
0205         npoints=length(find(idx)); % number of points in this average
0206         % plot avg FRs
0207         cnt=cnt+1; % subplot index
0208         subplot(nsets,nvaried,cnt);
0209         plot(pvals,rvals,'o-','linewidth',2);
0210         % todo: add confidence intervals / error bars
0211         axis tight; %ylim([0 max(FRmu(:))]);
0212         xlabel(strrep(varied{j},'_','\_'));
0213         ylabel([popname ': avg firing rate (Hz)']);
0214         title([popname ': firing rate averaged over pop, time, sims']);
0215         % add text: ncells, length(time), npoints
0216         xpos=min(xlim)+.7*diff(xlim);%.05*diff(xlim);
0217         ypos=min(ylim)+.2*diff(ylim);%.85*diff(ylim);
0218         text(xpos,ypos,sprintf('per point:\n #cells=%g\n #bins=%g\n #sims=%g',ncells,length(time),npoints));
0219       end
0220     end
0221     if length(varied)==2
0222       % plots for N populations and 2 varied elements
0223       % organize and imagesc FRmu(param 1, param 2)
0224       if ~lock_gca; handles(end+1)=figure('position',[1150 max(50,600-(nsets-1)*ht) 500 min(ht*nsets,750)]); end
0225       pvals1=unique(params(:,1)); nv1=length(pvals1);
0226       pvals2=unique(params(:,2)); nv2=length(pvals2);
0227       for k=1:nsets
0228         popname=regexp(FR_fields{k},'^([a-zA-Z0-9]+)_','tokens','once');
0229         popname=popname{1};
0230         FRmu2=zeros(nv1,nv2);
0231         for i=1:nv1
0232           for j=1:nv2
0233             idx=(params(:,1)==pvals1(i))&(params(:,2)==pvals2(j));
0234             FRmu2(i,j)=mean(FRmu(idx,k));
0235           end
0236         end
0237         % plot imagesc(param1,param2,FRmu)
0238         subplot(nsets,1,k);
0239         imagesc(pvals1,pvals2,FRmu2'); axis xy
0240         xlabel(strrep(varied{1},'_','\_')); ylabel(strrep(varied{2},'_','\_'));
0241         title([popname ': avg firing rate (Hz)']);
0242         caxis([0 max(FRmu(:))]); colorbar
0243       end
0244     end
0245   end
0246 end
dsPlotFR

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
