From 9d22944623b2a3080b0a1fdd8ccdd90ed444f219 Mon Sep 17 00:00:00 2001 From: raltodo Date: Mon, 25 Nov 2024 15:50:40 +0100 Subject: [PATCH] small changes --- SharpWaveRipples/DetectSWR.m | 3 +- lfp/FilterLFP.m | 2 +- lfp/FindThetaCycles.m | 14 ++++++- lfp/WaveletSpectrogram.m | 4 +- plotting/PlotIntervals.m | 53 +++++++++++++------------ plotting/RasterPlot.m | 2 +- plotting/semplot.m | 22 ++++++----- reactivation/ActivityTemplates.m | 2 +- spikes/DeconvolvePETH.m | 36 ++++++----------- spikes/FindBursts.m | 40 +++++++++---------- spikes/PETH.m | 55 +++++++++++--------------- spikes/PPC.m | 47 ++++++++++------------ spikes/PairOrder.m | 28 ++++++------- statistics/CorrByColumn.m | 28 ++++++++----- utilities/Bins.m | 4 +- utilities/Dist.m | 25 +++++++----- utilities/GLMgain.m | 34 ++++++++-------- utilities/Group.m | 13 +++--- utilities/RelativeTimes.m | 12 +++--- utilities/Shrink.m | 12 +++--- utilities/Unfind.m | 2 +- utilities/Unshift.m | 2 +- utilities/bz_Diff.m | 2 +- utilities/intervalsC+/FindInterval.m | 10 ++--- utilities/intervalsC+/SplitIntervals.m | 51 ++++++++++++++---------- 25 files changed, 255 insertions(+), 248 deletions(-) diff --git a/SharpWaveRipples/DetectSWR.m b/SharpWaveRipples/DetectSWR.m index b9bfa09a..0535ab17 100644 --- a/SharpWaveRipples/DetectSWR.m +++ b/SharpWaveRipples/DetectSWR.m @@ -1,4 +1,4 @@ -function [ripples] = DetectSWR(Channels, varargin) +function ripples = DetectSWR(Channels, varargin) %DetectSWR - Detect Sharp Wave Ripple in epochs of continuous data. % @@ -1483,6 +1483,7 @@ end +% ------------------------------- Helper functions ------------------------------- function logAnalysisFile(AnalysisFileName, writePath) % The purpose of this utility is to provide a record of analyses run by diff --git a/lfp/FilterLFP.m b/lfp/FilterLFP.m index e77b6fcc..d4262583 100644 --- a/lfp/FilterLFP.m +++ b/lfp/FilterLFP.m @@ -31,7 +31,7 @@ % gamma [30 80] % ripples [100 250] -% Copyright (C) 2004-2021 by Michaël Zugaro +% Copyright (C) 2004-2021 by Michaël Zugaro (C) 2021 by Ralitsa Todorova % % This program is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by diff --git a/lfp/FindThetaCycles.m b/lfp/FindThetaCycles.m index 492c9d00..2fc9f2af 100644 --- a/lfp/FindThetaCycles.m +++ b/lfp/FindThetaCycles.m @@ -118,4 +118,16 @@ troughs = troughs(ok); peaktopeak = peaktopeak(ok,:); amplitude = amplitude(ok,2); -end \ No newline at end of file + + +% ------------------------------- Helper functions ------------------------------- +function peaks = FindLocalMaxima(signal) +%FindLocalMaxima - find local peaks +if isdmatrix(signal,'@2'), + iPeaks = FindLocalMaxima(signal(:,2)); + peaks = signal(iPeaks,1); + return +end + +d = [nan;diff(signal(:))>0]; +peaks = strfind(d',[1 0])'; diff --git a/lfp/WaveletSpectrogram.m b/lfp/WaveletSpectrogram.m index b1c09d28..d966fc3c 100644 --- a/lfp/WaveletSpectrogram.m +++ b/lfp/WaveletSpectrogram.m @@ -16,7 +16,7 @@ % ------------------------------------------------------------------------- % 'range' frequency range (in Hz) (default = all) % 'resolution' desired frequency resolution (sub-octaves per octave, default = 0.05) -% 'step' desired step between successive windows (default = window/2): +% 'step' desired step between successive windows (default = window/2): % output will be downsampled to this value % 'interpolate' interpolate to create equidistant frequency bins instead % of bins in power of 2 (default = false) @@ -167,7 +167,7 @@ % loop through all scales and compute transform for a1 = 1:J1+1 - [daughter,fourier_factor,coi,dofmin]=wave_bases(mother,k,scale(a1),param); + [daughter,fourier_factor,coi,dofmin]=wave_bases(mother,k,scale(a1),param); wave(a1,:) = ifft(f.*daughter); % wavelet transform[Eqn(4)] end diff --git a/plotting/PlotIntervals.m b/plotting/PlotIntervals.m index a892fbc8..740ab48b 100644 --- a/plotting/PlotIntervals.m +++ b/plotting/PlotIntervals.m @@ -25,7 +25,8 @@ % ========================================================================= % -% Copyright (C) 2008-2013 by Gabrielle Girardeau & Michaël Zugaro +% Copyright (C) 2008-2013 by Gabrielle Girardeau & Michaël Zugaro, +% (C) 2023 by Ralitsa Todorova (graphics optimization for large numbers of intervals) % % This program is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by @@ -39,57 +40,57 @@ direction = 'v'; yLim = ylim; -if nargin < 1, +if nargin < 1 error('Incorrect number of parameters (type ''help PlotIntervals'' for details).'); end -if size(intervals,2) ~= 2, +if size(intervals,2) ~= 2 error('Incorrect list of intervals (type ''help PlotIntervals'' for details).'); end % Backward compatibility: previous versions used the syntax PlotIntervals(intervals,style,direction) parsed = false; -if (nargin == 2 || nargin == 3) && isstring_FMAT(lower(varargin{1}),'rectangles','bars'), +if (nargin == 2 || nargin == 3) && isastring(lower(varargin{1}),'rectangles','bars') style = lower(varargin{1}); parsed = true; end -if nargin == 3 && isstring_FMAT(lower(varargin{2}),'h','v'), +if nargin == 3 && isastring(lower(varargin{2}),'h','v') direction = lower(varargin{2}); parsed = true; end % Parse parameter list -if ~parsed, - for i = 1:2:length(varargin), - if ~ischar(varargin{i}), +if ~parsed + for i = 1:2:length(varargin) + if ~ischar(varargin{i}) error(['Parameter ' num2str(i+1) ' is not a property (type ''help PlotIntervals'' for details).']); end - switch(lower(varargin{i})), - case 'style', + switch(lower(varargin{i})) + case 'style' style = lower(varargin{i+1}); - if ~isstring_FMAT(style,'bars','rectangles'), + if ~isastring(style,'bars','rectangles') error('Incorrect value for property ''style'' (type ''help PlotIntervals'' for details).'); end - case 'direction', + case 'direction' direction = lower(varargin{i+1}); - if ~isstring_FMAT(direction,'h','v'), + if ~isastring(direction,'h','v') error('Incorrect value for property ''direction'' (type ''help PlotIntervals'' for details).'); end - case 'color', + case 'color' color = lower(varargin{i+1}); - if ~isstring_FMAT(color,'r','g','b','c','m','y','k','w') && ~isdvector(color,'#3','>=0','<=1'), + if ~isastring(color,'r','g','b','c','m','y','k','w') && ~isdvector(color,'#3','>=0','<=1') error('Incorrect value for property ''direction'' (type ''help PlotIntervals'' for details).'); end - case 'alpha', + case 'alpha' alphaValue = varargin{i+1}; - if ~isdscalar(alphaValue,'>=0','<=1'), + if ~isdscalar(alphaValue,'>=0','<=1') error('Incorrect value for property ''alpha'' (type ''help PlotIntervals'' for details).'); end - case 'ylim', + case 'ylim' yLim = varargin{i+1}; - if ~isdvector(yLim,'<'), + if ~isdvector(yLim,'<') error('Incorrect value for property ''yLim'' (type ''help PlotIntervals'' for details).'); end - otherwise, + otherwise error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help PlotIntervals'' for details).']); end end @@ -97,9 +98,9 @@ hold on; xLim = xlim; -if strcmp(style,'bars'), - for i = 1:size(intervals,1), - if strcmp(direction,'v'), +if strcmp(style,'bars') + for i = 1:size(intervals,1) + if strcmp(direction,'v') plot([intervals(i,1) intervals(i,1)],yLim,'Color',[0 0.75 0]); plot([intervals(i,2) intervals(i,2)],yLim,'Color',[0.9 0 0]); else @@ -108,8 +109,8 @@ end end else - for i=1:size(intervals,1), - if strcmp(direction,'v'), + for i=1:size(intervals,1) + if strcmp(direction,'v') dx = intervals(i,2)-intervals(i,1); dy = yLim(2)-yLim(1); rec(i) = patch(intervals(i,1)+[0 0 dx dx],yLim(1)+[0 dy dy 0],color,'LineStyle','none'); @@ -127,7 +128,7 @@ end end -if nargout>0, +if nargout>0 varargout{1} = rec; end diff --git a/plotting/RasterPlot.m b/plotting/RasterPlot.m index ae812739..937a8f39 100644 --- a/plotting/RasterPlot.m +++ b/plotting/RasterPlot.m @@ -13,7 +13,7 @@ % the Free Software Foundation; either version 3 of the License, or % (at your option) any later version. -if nargin==1, +if nargin==1 height = 1; end diff --git a/plotting/semplot.m b/plotting/semplot.m index 06a4b2db..56709f35 100644 --- a/plotting/semplot.m +++ b/plotting/semplot.m @@ -1,4 +1,4 @@ -function varargout = semplot(x,y,color,smooth) +function varargout = semplot(x,y,color,smooth,solid) %semplot - plot mean (line) +/- s.e.m. (shaded area) of a matrix "y" % semplot(x,y,color,smooth) @@ -12,32 +12,34 @@ % the Free Software Foundation; either version 3 of the License, or % (at your option) any later version. -solid = false; % one day I will make this into an option +if ~exist('solid','var') + solid = false; +end -if nargin<2 || (nargin==2 && ischar(y)), - if nargin==2, +if nargin<2 || (nargin==2 && ischar(y)) + if nargin==2 color = y; end y = x; x = 1:size(y,2); end -if ~exist('color','var'), +if ~exist('color','var') color = [0 0 0]; end -if ~exist('smooth','var'), +if ~exist('smooth','var') smooth = 0; end -if size(y,2)~=length(x), +if size(y,2)~=length(x) y = y'; - if size(y,2)~=length(x), + if size(y,2)~=length(x) error('Y should have one column for each element in X'); end end -if isvector(y), +if isvector(y) handles = plot(x,Smooth(y,smooth),'color',color); if nargout>0, varargout{1} = handles; end hold on; @@ -54,7 +56,7 @@ y = Smooth(nanmean(y),smooth); handles = fill(xx,yy,color); -if solid, +if solid set(handles,'FaceColor',mean([color;1 1 1]),'edgeAlpha',0); else % transparent set(handles,'FaceAlpha',0.5,'edgeAlpha',0); diff --git a/reactivation/ActivityTemplates.m b/reactivation/ActivityTemplates.m index b19fa105..70dce87c 100644 --- a/reactivation/ActivityTemplates.m +++ b/reactivation/ActivityTemplates.m @@ -45,7 +45,7 @@ % % See also ReactivationStrength. -% Copyright (C) 2016-2022 by Michaël Zugaro, Ralitsa Todorova +% Copyright (C) 2016-2024 by Michaël Zugaro, Ralitsa Todorova % % This program is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by diff --git a/spikes/DeconvolvePETH.m b/spikes/DeconvolvePETH.m index adfe74d0..21c17127 100644 --- a/spikes/DeconvolvePETH.m +++ b/spikes/DeconvolvePETH.m @@ -1,12 +1,12 @@ function [deconvolved,t] = DeconvolvePETH(signal,events,varargin) -% [DeconvolvePETH] - [Compute a deconvolved version of PETH which removes the -% smoothing effects of the events' autocorrelogram] +%DeconvolvePETH - Compute a deconvolved version of PETH which removes the +% smoothing effects of the events' autocorrelogram % % INPUTS -% [signal] [signal to find events for] -% [events] [events for PETH] -% optional list of property-value pairs (see table below) +% signal signal to find events for +% events events for PETH +% optional list of property-value pairs (see table below) % ========================================================================= % Properties Values % ------------------------------------------------------------------------- @@ -30,15 +30,15 @@ % ========================================================================= % % OUTPUTS -% [deconvolved] [deconvolved PETH] -% [t] [times of PETH] +% deconvolved deconvolved PETH +% t times of PETH % % SEE ALSO % -% [PETH. This is using similar script as PETH with deconvultion] +% PETH % -% [Ralitsa Todorova] [2021-2022] +% Copyright (C) 2021-2024 by Ralitsa Todorova % % This program is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by @@ -47,7 +47,7 @@ % default values: durations = [-1 1]; -nBins = [100]; +nBins = 100; argsToPassOn = {}; for i = 1:2:length(varargin) @@ -85,7 +85,7 @@ [~,t] = PETH(signal,events,'durations',durations,'nBins',nBins,argsToPassOn{:}); [autoResponse,~] = PETH(events,events,'durations',durations*2,'nBins',nBins*2-1,argsToPassOn{:}); % make sure the binsize stays the same while you double the duration -[rawResponse,~] = PETH(signal,events,'durations',durations*2,'nBins',nBins*2-1,argsToPassOn{:}); +[rawResponse,~] = PETH(signal,events,'durations',durations*2,'nBins',nBins*2-1,argsToPassOn{:}); autocorrelogram = sum(autoResponse); rawPeth = sum(rawResponse); @@ -102,17 +102,3 @@ % add the baseline to the final PETH deconvolved = T \ rawPeth(round(nBins/2 + 0):round(nBins/2*3 - 1))' + const/length(events); end - -%% Old code: -% function deconvolved = DeconvolvePETH(rawPeth,autocorrelogram) -% -% Note: "rawPETH" should be of double the duration of the autocorrelogram for this to work/ - -% Example: -% [autoResponse,t]= PETH(events,events,'durations',durations,'nBins',nBins); -% [rawResponse,~]= PETH(signal,events,,'durations',durations*2,'nBins',nBins*2-1); % make sure the binsize stays the same while you double the duration -% deconvolved = DeconvolvePETH(sum(rawResponse),sum(autoResponse)); -% plot(t,deconvolved); -% -% T = toeplitz([autocorrelogram(:); zeros(numel(rawPeth)-numel(autocorrelogram), 1)], [autocorrelogram(1), zeros(1, length(autocorrelogram)-1)]); -% deconvolved = T \ rawPeth(:); \ No newline at end of file diff --git a/spikes/FindBursts.m b/spikes/FindBursts.m index 1045bfb9..9029b999 100644 --- a/spikes/FindBursts.m +++ b/spikes/FindBursts.m @@ -34,7 +34,7 @@ % be discarded (default = [0 Inf]). % ========================================================================= % -% Copyright (C) 2016-2023 by Ralitsa Todorova, Michaël Zugaro +% Copyright (C) 2016-2024 by Ralitsa Todorova, Michaël Zugaro % % This program is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by @@ -53,34 +53,34 @@ basepath = []; save_mat = false; -for i = 1:2:length(varargin), - if ~ischar(varargin{i}), +for i = 1:2:length(varargin) + if ~ischar(varargin{i}) error(['Parameter ' num2str(i+2) ' is not a property (type ''help FindBursts'' for details).']); end - switch(lower(varargin{i})), - case 'thresholds', + switch(lower(varargin{i})) + case 'thresholds' thresholds = varargin{i+1}; - if ~isdvector(thresholds,'#2','<'), + if ~isdvector(thresholds,'#2','<') error('Incorrect value for property ''thresholds'' (type ''help FindBursts'' for details).'); end case 'binsize' binSize = varargin{i+1}; - if ~isdvector(binSize,'#1'), + if ~isdvector(binSize,'#1') error('Incorrect value for property ''binSize'' (type ''help FindBursts'' for details).'); end - case 'smooth', + case 'smooth' smooth = varargin{i+1}; - if ~isdvector(smooth,'#1','>0'), + if ~isdvector(smooth,'#1','>0') error('Incorrect value for property ''smooth'' (type ''help FindBursts'' for details).'); end - case 'intervals', + case 'intervals' intervals = varargin{i+1}; - if ~isdmatrix(intervals) || size(intervals,2) ~= 2, + if ~isdmatrix(intervals) || size(intervals,2) ~= 2 error('Incorrect value for property ''intervals'' (type ''help FindBursts'' for details).'); end - case 'durations', + case 'durations' durations = varargin{i+1}; - if ~isdmatrix(durations) || size(durations,2) ~= 2, + if ~isdmatrix(durations) || size(durations,2) ~= 2 error('Incorrect value for property ''durations'' (type ''help FindBursts'' for details).'); end case 'basepath' @@ -93,7 +93,7 @@ if ~islogical(save_mat) error('Incorrect value for property ''save_mat'' (type ''help FindBursts'' for details).'); end - otherwise, + otherwise error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help FindBursts'' for details).']); end end @@ -202,13 +202,13 @@ function save_mat_file(bursts,basepath,intervals,binSize,durations,thresholds,sm % dN (output binned spike counts as a matrix defined on bins starting with the % earliest spike across all channels and ending with the latest spike) % t (lower limit of each bin) -if nargin < 2; error('Need at least two input arguments'); end; +if nargin < 2; error('Need at least two input arguments'); end binSize=1/Fs; binSizemp=''; -if isstruct(data); +if isstruct(data) C=length(data); fnames=fieldnames(data); - if nargin <3 || isempty(t); + if nargin <3 || isempty(t) mintime=zeros(1,C); maxtime=zeros(1,C); for ch=1:C @@ -232,14 +232,14 @@ function save_mat_file(bursts,basepath,intervals,binSize,durations,thresholds,sm % if maxtimech > max(t); t=[t maxtimech+binSize]; end; end t=linspace(mintime,maxtime,1+(maxtime-mintime)/binSize); - for ch=1:C; + for ch=1:C eval(['binSizemp=data(ch).' fnames{1} ';']) x=histc(binSizemp,t); dN(:,ch)=x(:); end else binSizemp=data; - if nargin < 3; + if nargin < 3 mintime=min(binSizemp); maxtime=max(binSizemp); else @@ -253,7 +253,7 @@ function save_mat_file(bursts,basepath,intervals,binSize,durations,thresholds,sm end end -function timestamps = Unshift(timestamps, intervals), +function timestamps = Unshift(timestamps, intervals) % The opposite of the option 'shift' in Restrict % diff --git a/spikes/PETH.m b/spikes/PETH.m index 6a75912a..55f5c89e 100644 --- a/spikes/PETH.m +++ b/spikes/PETH.m @@ -73,7 +73,7 @@ % % See also Sync, SyncHist, SyncMap, PlotSync, PETHTransition. % -% Copyright (C) 2018-2022 by Ralitsa Todorova & Michaël Zugaro +% Copyright (C) 2018-2024 by Ralitsa Todorova & Michaël Zugaro % % This program is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by @@ -94,57 +94,48 @@ end mode = 'l'; -for i = 1:2:length(varargin), - switch(lower(varargin{i})), - case 'durations', +for i = 1:2:length(varargin) + switch(lower(varargin{i})) + case 'durations' duration = varargin{i+1}; - if ~isvector(duration) || length(duration) ~= 2, + if ~isvector(duration) || length(duration) ~= 2 error('Incorrect value for property ''durations'' (type ''help PETH'' for details).'); end - case 'duration', + case 'duration' duration = varargin{i+1}; - if ~isvector(duration) || length(duration) ~= 2, + if ~isvector(duration) || length(duration) ~= 2 error('Incorrect value for property ''durations'' (type ''help PETH'' for details).'); end - case 'nbins', + case 'nbins' nBins = varargin{i+1}; - if ~isvector(nBins) || length(nBins) ~= 1, + if ~isvector(nBins) || length(nBins) ~= 1 error('Incorrect value for property ''nBins'' (type ''help PETH'' for details).'); end - case 'show', + case 'show' show = varargin{i+1}; - if ~isastring(show,'on','off'), + if ~isastring(show,'on','off') error('Incorrect value for property ''show'' (type ''help PETH'' for details).'); end - case 'mode', + case 'mode' mode = varargin{i+1}; - if ~isastring(mode,'l','c'), + if ~isastring(mode,'l','c') error('Incorrect value for property ''mode'' (type ''help PETH'' for details).'); end - case 'title', + case 'title' namestring = varargin{i+1}; - if ~isastring(namestring), + if ~isastring(namestring) error('Incorrect value for property ''title'' (type ''help PETH'' for details).'); end - case 'smooth', + case 'smooth' smooth = varargin{i+1}; - if ~isvector(smooth) || length(smooth) ~= 1, + if ~isvector(smooth) || length(smooth) ~= 1 error('Incorrect value for property ''smooth'' (type ''help PETH'' for details).'); end - otherwise, + otherwise pictureoptions = varargin(i:end); break end end -seconds_around_event = diff(duration)/2; -% if ~exist('nBins', 'var'), -% if seconds_around_event < 1, -% nBins = seconds_around_event*300; -% else -% nBins = seconds_around_event*100; -% end -% end - if size(samples,2)==2 % if the provided data is a signal rather than events t = linspace(duration(:,1),duration(2),nBins); mat_t = bsxfun(@plus,events,t); @@ -163,7 +154,7 @@ return else % the samples are a point process [sync, j] = Sync(samples, events, 'durations', duration); - if nargout>0, + if nargout>0 s = Bin(sync(:,1),duration,nBins); mat = zeros(size(events,1),nBins); if isempty(sync) @@ -176,13 +167,13 @@ varargout{1} = mat; varargout{2} = t; end - if strcmpi(show,'on'), % compute 'm' that we will plot + if strcmpi(show,'on') % compute 'm' that we will plot [m, ~, t] = SyncHist(sync, j, 'nBins', nBins, 'smooth', smooth, 'mode', 'mean', 'durations', duration); end end -if strcmpi(show,'on'), - if isempty(pictureoptions), +if strcmpi(show,'on') + if isempty(pictureoptions) if exist('linetype','var'), PlotXY(t', m, linetype); else, PlotXY(t', m); end title([namestring ', ' num2str(numel(j)) ' x ' num2str(numel(unique(j))) ' instances']); else @@ -190,7 +181,7 @@ end end -if nargout>0, +if nargout>0 varargout{1} = mat; varargout{2} = t; if nargout>2 diff --git a/spikes/PPC.m b/spikes/PPC.m index 3253eebf..65449dbb 100644 --- a/spikes/PPC.m +++ b/spikes/PPC.m @@ -1,32 +1,29 @@ -function ppc = PPC(phases,optionalInput) +function ppc = PPC(phases,input2) -% [PPC - compute Pairwise Phase Consistency across spikes] +%PPC - compute Pairwise Phase Consistency across spikes % -% [ppc = PPC(phases) computes the Pairwise Phase Consistency as defined by +% ppc = PPC(phases) computes the Pairwise Phase Consistency as defined by % Vinck, Wingerden, Womelsdorf, Fries, Pennartz (2010, NeuroImage) % ppc = PPC(phases,trialID) computes P(hat)2 Pairwise Phase Consistency % as defined by Vinck, Battaglia, Womelsdorf, Pennartz (2012, J Comp Neurosci) % ppc = PPC(phases1,phases2) computes the Pairwise Phase Consistency across units -% as defined by Vinck, Wingerden, Womelsdorf, Fries, Pennartz (2010, -% NeuroImage)] +% as defined by Vinck, Wingerden, Womelsdorf, Fries, Pennartz (2010, NeuroImage) % % % INPUTS -% [phases] [phases from spikes, lfp] -% -% optional list of property-value pairs (see table below) -% ========================================================================= -% Properties Values -% ------------------------------------------------------------------------- -% ['trial'] [if 'optionalInput' is vector, trials are assumed] -% ========================================================================= +% phases a vector containing phases (in radians) +% input2 it can either be a +% (1) vector of phases (in radians) emitted by another unit +% to compute PPC between the two cells, or +% (2) a vector of trial IDs (one for each phase). The PPC +% will be computed independently for each trial. % % OUTPUTS -% [ppc] [pairwise phase consistency score] +% ppc pairwise phase consistency score % % % SEE ALSO % -% [Ralitsa Todorova] [2016-2022] +% Copyright (C) 2016-2024 by Ralitsa Todorova % % This program is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by @@ -46,28 +43,26 @@ % i.e. taking the mean is equilalent to the sum, normalised by n*(n-1) as in Vinck et al (2010) ppc = nanmean(cos(difference(:))); % as cos(a-b) = cos(a)*cos(b)+sin(a)*sin(b), Vinck et al's definition of PPC else - if isivector(optionalInput) % if it's a vector of integers, assume it's trial ID - trial = optionalInput; - phases = phases(:); trial = (:); - if numel(trial)~=numel(phases), + if isivector(input2) % if it's a vector of integers, assume it's trial ID + phases = phases(:); trialID = input2(:); + if numel(trialID)~=numel(phases) error('Please provide a trial number for each phase'); end - nTrials = max(trial); + nTrials = max(trialID); ppc = nan(nTrials); - for m=1:nTrials, - for l=1:nTrials, - if l~=m, - ppc(m,l) = PPC2(phases(trial==m),phases(trial==l)); + for m=1:nTrials + for l=1:nTrials + if l~=m + ppc(m,l) = PPC2(phases(trialID==m),phases(trialID==l)); end end end - % as the diagonal is now composed of nans, there are n*(n-1) non-nan elements in the matrix 'ppc' % i.e. taking the mean is equilalent to the sum, normalised by n*(n-1) as in Vinck et al (2012) ppc = nanmean(ppc(:)); else % assume it's phases of another unit and compute cross-unit PPC - phases1 = phases; phases2 = optionalInput; + phases1 = phases; phases2 = input2; difference = bsxfun(@minus,phases1,phases2'); ppc = nanmean(cos(difference(:))); % as cos(a-b) = cos(a)*cos(b)+sin(a)*sin(b), Vinck et al's definition of PPC end diff --git a/spikes/PairOrder.m b/spikes/PairOrder.m index 7a11e32c..d9b0cd53 100644 --- a/spikes/PairOrder.m +++ b/spikes/PairOrder.m @@ -1,41 +1,41 @@ function [score] = PairOrder(targetSequence, sequenceList, varargin) -% [Method used by Gupta et al (Redish) 2010] -% [sequenceList can be a matrix with NaNs, for example: +%Apply the method used by Gupta et al (Redish) 2010 +% sequenceList can be a matrix with NaNs, for example: % % 1 3 4 2 1 6 NaN % NaN 2 NaN 1 2 3 5 % 4 1 1 5 2 6 6 % where each row is a sequence. 'score' will have a score (row) for each of -% the rows in sequenceList] +% the rows in sequenceList % % EXAMPLE: -% [[~,peakBin] = max(firingMaps,[],2);] -% [[~,correctOrder] = sort(peakBin); % the order of place fields (as -% ordered by the location of their peak)] -% [sequenceList = GetSequenceList(spikes, thetaCycles)] -% [thetaScores = PairOrder(correctOrder,sequenceList)] +% [~,peakBin] = max(firingMaps,[],2); +% [~,correctOrder] = sort(peakBin); % the order of place fields (as +% ordered by the location of their peak) +% sequenceList = GetSequenceList(spikes, thetaCycles) +% thetaScores = PairOrder(correctOrder,sequenceList) % % % INPUTS -% [targetSequence] [sequence with intended order] -% [sequenceList] [sequence list actually recorder. See above] +% targetSequence sequence with intended order +% sequenceList sequence list actually recorder. See above % optional list of property-value pairs (see table below) % ========================================================================= % Properties Values % ------------------------------------------------------------------------- -% ['circular'] [if data is circular 'on.' If not 'off' (Default)] -% ['normalize'] [normalize score] +% 'circular' if data is circular 'on.' If not 'off' (Default) +% 'normalize' normalize score % ========================================================================= % % OUTPUTS -% [score] [Description of input 1] +% score Description of input 1 % % % SEE ALSO % -% [Ralitsa Todorova] [2016-2022] +% Copyright (C) 2016-2022 by Ralitsa Todorova % % This program is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by diff --git a/statistics/CorrByColumn.m b/statistics/CorrByColumn.m index 9bdffb7f..648d3e87 100644 --- a/statistics/CorrByColumn.m +++ b/statistics/CorrByColumn.m @@ -1,20 +1,21 @@ -function [c,slope,rSquared,SSres,startstop] = CorrByColumn(A,B) +function [c,slope,rSquared,SSres] = CorrByColumn(A,B,parametric) -% CorrByColumn - correlate matrices A and B column by column +% CorrByColumn - correlate each column of matrix A to each column of matrix B % -% i.e. c(i) = corr(A(:,i),B(:,i)) +% This is a faster implementation of the loop: +% for i=1:size(A,1), c(i) = corr(A(:,i),B(:,i)); end % % Input: -% A i,j matrix -% B i,j matrix +% A matrix +% B matrix +% parametric true/false toggle; parametric = true (default) will compute +% the Pearson coefficient, whereas parametric = +% false will compute the Spearman coefficient. % Output: % c correlation (r) % slope slope of regression % rSquared coefficient of variation (R2) % SSres sum of squares residual -% startstop ??? -% -% TODO: determine what startstop is % % Copyright (C) 2016 by Ralitsa Todorova % @@ -28,6 +29,13 @@ B(bad) = nan; nonconclusive = sum(double(~bad))'<3; +if ~exist('parametric','var'), parametric = true; end + +if ~parametric % transform the values to ranks + A = tiedrank(A); + B = tiedrank(B); +end + An = bsxfun(@minus,A,nanmean(A,1)); Bn = bsxfun(@minus,B,nanmean(B,1)); An = bsxfun(@times,An,1./sqrt(nansum(An.^2,1))); @@ -44,6 +52,6 @@ rSquared = 1 - sum((An-bsxfun(@times,Bn,c')).^2)'; SSres = (1-rSquared).*(sum(~bad)' - 1).*var(A)'; -Afit = bsxfun(@plus,bsxfun(@times,bsxfun(@minus,B,nanmean(B,1)),slope'),nanmean(A,1)); -startstop = Afit([1 end],:)'; +% Afit = bsxfun(@plus,bsxfun(@times,bsxfun(@minus,B,nanmean(B,1)),slope'),nanmean(A,1)); +% startstop = Afit([1 end],:)'; end \ No newline at end of file diff --git a/utilities/Bins.m b/utilities/Bins.m index 15d9a062..0282fda8 100644 --- a/utilities/Bins.m +++ b/utilities/Bins.m @@ -11,8 +11,8 @@ % the Free Software Foundation; either version 3 of the License, or % (at your option) any later version. -if nargin<4, - if nargin<3, +if nargin<4 + if nargin<3 window = 1; end step = window; diff --git a/utilities/Dist.m b/utilities/Dist.m index 8a989a59..99bd4e78 100644 --- a/utilities/Dist.m +++ b/utilities/Dist.m @@ -5,6 +5,9 @@ % where 'd' is the normalised distribution of the data for points 'x'. % Dist(resolution, data) without output arguments plots said distributions. % +% This is very similar to "hist" but it normalizes the distribitions and +% accepts a variety of ways to input the data. +% % USAGE % % FOR VECTORS: @@ -27,28 +30,30 @@ barplot = 'off'; -if strcmp(varargin{end}, 'bar'), +if strcmp(varargin{end}, 'bar') barplot = 'on'; varargin = varargin(1:end-1); end -if length(varargin)==1 && ~isvector(varargin{1}), - for i=2:size(varargin{1},2), +if length(varargin)==1 && ~isvector(varargin{1}) + for i=2:size(varargin{1},2) varargin{i} = varargin{1}(:,i); end varargin{1} = varargin{1}(:,1); elseif length(varargin)>1 && (strcmp(varargin{2}, 'group') || strcmp(varargin{2}, 'groups') || strcmp(varargin{2}, 'grouped')), - for i=max(varargin{1}(:,end)):-1:1 %flipped so 1 comes last and I don't override the data + for i=max(varargin{1}(:,end)):-1:1 %flipped so 1 comes last and we don't override the data varargin{i} = varargin{1}(varargin{1}(:,end)==i,1); end if (strcmp(varargin{2}, 'group') || strcmp(varargin{2}, 'groups') || strcmp(varargin{2}, 'grouped')), varargin(2) = []; end +else % treat each input as a separate variable. Make sure they are vertical vectors: + for i=1:length(varargin), varargin{i} = varargin{i}(:); end end n = length(varargin); -for i=1:n, - if isempty(varargin{i}), +for i=1:n + if isempty(varargin{i}) minmaxima(i,1:2) = nan; else minmaxima(i,1) = min(varargin{i}(:)); @@ -60,15 +65,15 @@ limits = [limits(1)-diff(limits)/5 limits(2)+diff(limits)/5]; % broaden by 40% [~, t] = hist(limits, resolution); -for i=1:n, +for i=1:n [h(:,i), ~] = hist(varargin{i}, t); end h = h./repmat(sum(h), size(h,1), 1); -if nargout==0, +if nargout==0 - if strcmp(barplot, 'on'), + if strcmp(barplot, 'on') handle = bar(t,h); else handle = plot(t, h); @@ -76,7 +81,7 @@ legend(num2str([1:n]')); end -if nargout>0, +if nargout>0 varargout{1} = h; varargout{2} = t; end \ No newline at end of file diff --git a/utilities/GLMgain.m b/utilities/GLMgain.m index b581433d..2a281f90 100644 --- a/utilities/GLMgain.m +++ b/utilities/GLMgain.m @@ -77,36 +77,36 @@ mode = 'mean'; fixedSets = false; -for i = 1:2:length(varargin), - if ~ischar(varargin{i}), +for i = 1:2:length(varargin) + if ~ischar(varargin{i}) error(['Parameter ' num2str(i+2) ' is not a property (type ''help GLMgain'' for details).']); end - switch(lower(varargin{i})), - case 'link', + switch(lower(varargin{i})) + case 'link' link = varargin{i+1}; - case 'dist', + case 'dist' dist = varargin{i+1}; - case 'mode', + case 'mode' mode = varargin{i+1}; - case 'setid', + case 'setid' setID = varargin{i+1}; nSets = max(setID); fixedSets = true; - case 'nsets', + case 'nsets' nSets = varargin{i+1}; if ~isscalar(nSets) || mod(nSets,1)>0 error('Incorrect value for property ''nSets'' (type ''help GLMgain'' for details).'); end - case 'niterations', + case 'niterations' nIterations = varargin{i+1}; if ~isscalar(nIterations) || mod(nIterations,1)>0 error('Incorrect value for property ''nIterations'' (type ''help GLMgain'' for details).'); end - case 'minevents', + case 'minevents' minEvents = varargin{i+1}; if ~isscalar(minEvents) || mod(minEvents,1)>0 error('Incorrect value for property ''minEvents'' (type ''help GLMgain'' for details).'); end - otherwise, + otherwise error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help GLMgain'' for details).']); end end @@ -129,7 +129,7 @@ sh = nan(nUnits,nIterations); w = nan(nUnits,size(source,2)+1,nSets); predictions = nan(size(target)); -for unit = 1:nUnits, +for unit = 1:nUnits %% Split data into nSets balanced sets if ~fixedSets, setID = nan(nEvents,1); end zero = target(:,unit)==0; @@ -146,7 +146,7 @@ end shuffled = nan(nSets,nIterations); errors = nan(nSets,1); - for set = 1:nSets, + for set = 1:nSets lastwarn(''); weights = glmfit(source(setID~=set,:),target(setID~=set,unit),dist,'link',link); if strcmp(lastwarn,'Iteration limit reached.') % if the algorhithm failed to converge @@ -156,7 +156,7 @@ prediction = linkFunction(source(setID==set,:),weights); predictions(setID==set,unit) = prediction; errors(set,1) = average(abs(target(setID==set,unit) - prediction)); - for iteration = 1:nIterations, + for iteration = 1:nIterations % get randomly ordered numbers [~,indices] = sort(rand(sum(setID==set),1)); shuffled(set,iteration) = average(abs(target(setID==set,unit) - prediction(indices))); @@ -169,7 +169,7 @@ end warning(state); -if sum(~notempty)>0, +if sum(~notempty)>0 gain0 = gain; shGain0 = shGain; er0 = er; @@ -197,11 +197,11 @@ % the Free Software Foundation; either version 3 of the License, or % (at your option) any later version. -if size(x,1) == 1, %if a row vector is provided +if size(x,1) == 1 %if a row vector is provided x = x(:); %change it to a column vector end -if ~exist('n', 'var'), +if ~exist('n', 'var') n = size(x,1); end diff --git a/utilities/Group.m b/utilities/Group.m index d95b8882..ed804565 100644 --- a/utilities/Group.m +++ b/utilities/Group.m @@ -4,11 +4,9 @@ % in the same format as spikes are grouped ([values id]) % vector1 1 % vector2 2 -% ... and so on. Useful for creating grouped events. -% ATTENTION: rows are automatically sorted, unless you call -% Group(vector1, vector2, ..., vectorN, 'sort', false); +% ... and so on. This is useful for creating grouped events. % -% Copyright (C) 2018-2022 by Ralitsa Todorova +% Copyright (C) 2018-2024 by Ralitsa Todorova % % This program is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by @@ -17,10 +15,10 @@ grouped = []; doSort = false; % Are we dealing with vectors? -if all(cellfun(@(x) min(size(x)), varargin)<=1), +if all(cellfun(@(x) min(size(x)), varargin)<=1) vectors = true; else vectors = false; end -for i=1:length(varargin), +for i=1:length(varargin) if ischar(varargin{i}) && strcmpi(varargin{i},'sort'), doSort = varargin{i+1}; break; end @@ -35,7 +33,6 @@ end end - -if doSort, +if doSort grouped = sortrows(grouped); end \ No newline at end of file diff --git a/utilities/RelativeTimes.m b/utilities/RelativeTimes.m index c193009c..1ea467e6 100644 --- a/utilities/RelativeTimes.m +++ b/utilities/RelativeTimes.m @@ -46,11 +46,11 @@ % (at your option) any later version. -if nargin<3, +if nargin<3 values = 0:size(intervals,2)-1; end -rt = nan(length(t),size(intervals,2)-1); +[rt,intervalID] = deal(nan(length(t),size(intervals,2)-1)); for i=1:size(intervals,2)-1 [rt(:,i),intervalID(:,i)] = RelativeTime(t,intervals(:,[i i+1])); end @@ -62,7 +62,7 @@ [~,columnsToKeep] = max(abs(0.5-rt),[],2); indicesToKeep = sub2ind(size(rt),(1:length(t))',columnsToKeep); -for i=1:size(intervals,2)-1, +for i=1:size(intervals,2)-1 rt(:,i) = rt(:,i)*(values(i+1)-values(i)) + values(i); end @@ -84,16 +84,16 @@ rt = nan(size(t)); -if size(intervals,2)==3, +if size(intervals,2)==3 intervalID = zeros(length(t), size(intervals,2)-1); - for i=1:size(intervals,2)-1, + for i=1:size(intervals,2)-1 theseIntervals = intervals(:,i:i+1); [inIntervals,intervalID(:,i)] = InIntervals(t, theseIntervals); rt(inIntervals) = RelativeTime(t(inIntervals), theseIntervals)+i-1; end return end -[inIntervals intervalID] = InIntervals(t, intervals); +[inIntervals, intervalID] = InIntervals(t, intervals); rt(inIntervals,1) = (t(inIntervals) - intervals(intervalID(inIntervals),1))./... (intervals(intervalID(inIntervals),end) - intervals(intervalID(inIntervals),1)); diff --git a/utilities/Shrink.m b/utilities/Shrink.m index 88b26b36..4b162824 100644 --- a/utilities/Shrink.m +++ b/utilities/Shrink.m @@ -19,12 +19,12 @@ n = ceil(n0./[rowBinSize columnBinSize]).*[rowBinSize columnBinSize]; % This is a heuristic to make the code run faster: -if columnBinSize>rowBinSize, +if columnBinSize>rowBinSize shrunk = Shrink(matrix',columnBinSize,rowBinSize)'; return; end -if n0(1)~=n(1), +if n0(1)~=n(1) if n0(1)==1, matrix = [matrix; nan(n(1)-1,n0(2))]; else nanRows = round(linspace(1,n(1),2+n(1)-n0(1))); nanRows([1 end]) = []; % Equally spaced rows of nans matrix0 = matrix; @@ -33,7 +33,7 @@ end end -if n0(2)~=n(2), +if n0(2)~=n(2) if n0(2)==1, matrix = [matrix nan(n(1),n(2)-1)]; else nanColumns = round(linspace(1,n(2),2+n(2)-n0(2))); nanColumns([1 end]) = []; % Equally spaced columns of nans matrix0 = matrix; @@ -43,14 +43,14 @@ end -if columnBinSize == 1 && rowBinSize == 1, +if columnBinSize == 1 && rowBinSize == 1 shrunk = matrix; return end order = []; -for i=1:columnBinSize, - if istops(end), +if isempty(stops) || starts(end)>stops(end) stops = [stops length(logical)]; end -if starts(1)>stops(1), +if starts(1)>stops(1) starts = [1 starts]; end diff --git a/utilities/intervalsC+/SplitIntervals.m b/utilities/intervalsC+/SplitIntervals.m index d5d84d1b..d66e8c6a 100644 --- a/utilities/intervalsC+/SplitIntervals.m +++ b/utilities/intervalsC+/SplitIntervals.m @@ -37,14 +37,21 @@ % ========================================================================= % Hint: to get overlapping windows, simply use % sortrows([SplitIntervals(intervals, window); SplitIntervals(intervals+window/2, window)]) -% EXAMPLE: -% SplitIntervals([0 2.1],'pieceSize',1,'mode','keep') generates [0 1; 1 2; 2 2.1] +% +% EXAMPLE 1 (with nPieces): +% SplitIntervals([2 5; 6 10], 'nPieces', 2) (i.e. split each of these intervals into 2) +% gives intervals = [2 3.5; 3.5 5; 6 8; 8 0] and ids = [1; 1; 2; 2] +% +% EXAMPLE 2 (with pieceSize): +% SplitIntervals([2 5; 6 10], 'pieceSize', 2) (i.e. split these intervals into 2-second pieces) +% gives intervals = [2 4; 6 8; 8 10] and ids = [1; 2; 2] % SplitIntervals([0 2.1],'pieceSize',1,'mode','discard') generates [0 1; 1 2], discarding the -% additional bin [2 2.1] which is too short. +% bin [2 2.1] which is too short (desired duration is 1 second). +% SplitIntervals([0 2.1],'pieceSize',1,'mode','keep') generates [0 1; 1 2; 2 2.1] % SplitIntervals([0 2.1],'pieceSize',1,'mode','keep','extend',true) generates [0 1; 1 2; 2 3] % extending the last piece so that it is also 1 second long. % -% Copyright (C) 2019-2023 Ralitsa Todorova +% Copyright (C) 2019-2023 Ralitsa Todorova & (C) 2023 by Federica Lareno Faccini % % This program is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by @@ -56,29 +63,29 @@ extend = false; nPieces = []; -for i = 1:2:length(varargin), - if ~ischar(varargin{i}), +for i = 1:2:length(varargin) + if ~ischar(varargin{i}) error(['Parameter ' num2str(i+2) ' is not a property (type ''help SplitIntervals'' for details).']); end - switch(lower(varargin{i})), - case 'piecesize', + switch(lower(varargin{i})) + case 'piecesize' pieceSize = varargin{i+1}; - if ~isvector(pieceSize) || length(pieceSize) ~= 1, + if ~isvector(pieceSize) || length(pieceSize) ~= 1 error('Incorrect value for property ''pieceSize'' (type ''help SplitIntervals'' for details).'); end - case 'mode', + case 'mode' mode = lower(varargin{i+1}); - if ~isastring(lower(mode),'discard','keep','round'), + if ~isastring(lower(mode),'discard','keep','round') error('Incorrect value for property ''mode'' (type ''help SplitIntervals'' for details).'); end - case 'npieces', + case 'npieces' nPieces = varargin{i+1}; - if ~isvector(nPieces) || length(nPieces) ~= 1, + if ~isvector(nPieces) || length(nPieces) ~= 1 error('Incorrect value for property ''nPieces'' (type ''help SplitIntervals'' for details).'); end - case 'extend', + case 'extend' extend = varargin{i+1}; - if isastring(lower(extend),'on','off'), + if isastring(lower(extend),'on','off') extend = strcmpi(extend,'on'); % transform to logical end if length(extend)==1 @@ -86,16 +93,17 @@ else error('Incorrect value for property ''discard'' (type ''help SplitIntervals'' for details).'); end - otherwise, + otherwise error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help SplitIntervals'' for details).']); end end -%% SPLIT EACH INTERVAL INTO EQUAL PIECES (nPieces) -if ~isempty(nPieces), +%% Split each interval into we equal pieces (nPieces) + +if ~isempty(nPieces) matrix = nan(size(intervals,1),nPieces+1); matrix(:,[1 end]) = intervals; - for i=2:(nPieces), + for i=2:(nPieces) matrix(:,i) = matrix(:,1)+(i-1)*(matrix(:,end)-matrix(:,1))/nPieces; end pieces = [reshape(matrix(:,1:end-1)',[],1) reshape(matrix(:,2:end)',[],1)]; @@ -103,7 +111,8 @@ return end -%% SPLIT INTERVALS INTO PIECES OF EQUAL ABSULUTE LENGTHS (pieceSize) +%% Split intervals into pieces of equal durations (pieceSize) + d = diff(intervals,[],2); switch mode case 'discard' @@ -134,7 +143,7 @@ try ids = indicesNotEmpty(ids); catch - warning('There is an issue with the indices'); + warning('There is an issue with the indices. Check results manually'); end