diff --git a/SharpWaveRipples/DetectSWR.m b/SharpWaveRipples/DetectSWR.m
index b9bfa09..0535ab1 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 e77b6fc..d426258 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 492c9d0..2fc9f2a 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 b1c09d2..d966fc3 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 a892fbc..740ab48 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 ae81273..937a8f3 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 06a4b2d..56709f3 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 b19fa10..70dce87 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 adfe74d..21c1712 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 1045bfb..9029b99 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 6a75912..55f5c89 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 3253eeb..65449db 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 7a11e32..d9b0cd5 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 9bdffb7..648d3e8 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 15d9a06..0282fda 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 8a989a5..99bd4e7 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 b581433..2a281f9 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 d95b888..ed80456 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 c193009..1ea467e 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 88b26b3..4b16282 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 d5d84d1..d66e8c6 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