swrtrig spect and some cleaning up of base cmw code

droumis · Oct 31, 2019 · bd1bbd6 · bd1bbd6
1 parent 47042e0
commit bd1bbd6
Show file tree

Hide file tree

Showing 12 changed files with 640 additions and 236 deletions.
diff --git a/DFAnalysis/dfa_reactivationPSTH.m b/DFAnalysis/dfa_reactivationPSTH.m
@@ -45,13 +45,17 @@
 
 timeBinEdges = epTimeRange(1):bin:epTimeRange(2);
 for s = 1:length(su(:,1))
-    iSpks = spikes{idx(1)}{idx(2)}{su(s,1)}{su(s,2)}.data;
+    iSpks = spikes{idx(1)}{idx(2)}{su(s,3)}{su(s,4)}.data;
     if length(iSpks) < numSpikesThresh
         fprintf('%s %d %d not enough spikes: %d < %d\n', animal, idx(1), idx(2), length(iSpks), numSpikesThresh);
         return
     end
+    outSpikes{s,1} = [repmat(s,length(iSpks),1) iSpks];
     binSpikes(s,:) = histc(iSpks,timeBinEdges);
 end
+out.su = su;
+out.binSpikes = binSpikes;
+out.spikes = cell2mat(outSpikes);
 %% Get Template and Match binned zSpikes
 templMask = logical(isExcluded(timeBinEdges, templateIntervals));
 fprintf('template pct eptime %.02f (%.03f s)\n',sum(templMask)/length(templMask)*100, ...
@@ -80,6 +84,7 @@
 %% Full Model version of Reactivation Strength
 if fullModel
     reactFull=diag(zSpikesMatch'*zCCtemplateNoDiag*zSpikesMatch)';
+%     reactFullNorm=(1/(2*numMatchBins))*sum(reactFull);
     out.reactFull = reactFull;
     %% all SWR psth
     d = setfiltertime(f, swrTimeFilter);
@@ -149,6 +154,9 @@
     burstIntervals = vec2list(~isExcluded(times, l.excludetime{1}{1}), times);
     allLicks = lick{idx(1)}{idx(2)}.starttime;
     lickBurstTime = allLicks(logical(isExcluded(allLicks, burstIntervals)));
+    % exclude licks within 1 second of a swr
+    [~, d] = knnsearch(swrBurstTime(:), lickBurstTime);
+    lickBurstTime = lickBurstTime(d > 1);
     lickBurstTime = lickBurstTime(lickBurstTime+min(idxWin)>0);
     lickBurstTime = lickBurstTime(lickBurstTime+max(idxWin)<length(reactFull));
     out.lickBurstTime = lickBurstTime;
@@ -277,6 +285,9 @@
 function out = init_out(idx, animal, varargin)
 out.idx = idx;
 out.animal = animal;
+out.su = [];
+out.binSpikes = [];
+out.spikes = [];
 out.reactFull = []; % full epoch 'match' intervals reactivation strength series
 out.reactTime = []; % full epoch 'match' time.. aka reactivation
 

diff --git a/DFFunctions/evaluatefilter.m b/DFFunctions/evaluatefilter.m
@@ -0,0 +1,46 @@
+function resultindex = evaluatefilter(cellvar, filterString)
+%
+% resultindex = evaluatefilter(cellvar, filterString)
+%
+% Return the indices of the cell structure that have fields that satisfy
+% the given constraints. Each row of the output is the index into one cell.
+%
+% CELLVAR - any cell structure, for example spikes{}{}{}{}.####
+% FILTERSTRING - A string containing an expression with the field variables
+%               of CELLVAR.  Each field variable must be preceded with a
+%               '$'.
+%
+% Example: index = evaluatefilter(cellinfo,'(($meanrate < 10) && isequal($area,''CA3''))')
+%          Returns all indices of cellinfo where the 'area' field is equal to 'CA3' and the 'meanrate' field is less than 10.
+
+
+if ~isequal(filterString, '')
+    newstruct = [];
+    [variables, expression] = parsefilterstring(filterString);
+    for i = 1:length(variables)
+        varfetch = cellfetch(cellvar, variables{i});
+        for j = 1:length(varfetch.values)
+
+                eval(['newstruct(',num2str(j),').',variables{i},' = varfetch.values{j};']);
+                %newstruct = setfield(newstruct,{j},variables{i},varfetch.values{j});
+
+
+        end
+    end
+
+    results = [];
+    for i = 1:length(newstruct)
+        structVar = newstruct(i);
+        eval(['tmpresult = ',expression,';']);
+        if isempty(tmpresult)
+            tmpresult = 0;
+        end
+        results(i) = tmpresult;
+    end
+    cellindex = cellfetch(cellvar, '');
+    resultindex = cellindex.index(find(results),:);
+else
+    cellindex = cellfetch(cellvar, '');
+    resultindex = cellindex.index;
+end
+
diff --git a/DFFunctions/load_filter_params.m b/DFFunctions/load_filter_params.m
@@ -119,6 +119,24 @@
             uselfptype = 'eeg';
         case 'unreferenced'
             uselfptype = 'eeggnd';
+
+        case 'lickbouts'
+            maxIntraBurstILI = 0.25;  % max burst ili threshold in seconds
+            minBoutLicks = 3; %filter out bouts with less than boutNum licks
+            % timefilt1: lick bouts, 
+            timefilter{end+1} = {'getLickBout', '($lickBout == 1)', ...
+                'maxIntraBurstILI', maxIntraBurstILI, ...
+                'minBoutLicks', minBoutLicks};
+
+        case 'nolickbouts'
+            % timefilt: immoble not lick bout
+            maxIntraBurstILI = 0.25;  % max burst ili threshold in seconds
+            minBoutLicks = 3; % filter out bouts with less than boutNum licks
+            minTimeFromLick = .5; % time duration from closest lick
+            timefilter{end+1} = {'getLickBout', ...
+                sprintf('($lickBout == 0) & ($timeFromLick > %d)',minTimeFromLick), ...
+                'maxIntraBurstILI', maxIntraBurstILI, 'minBoutLicks', minBoutLicks};
+
         case 'ripples'
             TF = 1;
             eventtype = 'rippleskons';
@@ -143,7 +161,7 @@
             %         Fp.welldist);
             %% filter function specific params
         case 'dfa_reactivationPSTH'
-            bin = 0.01; % seconds
+            bin = 0.025; % seconds
             win = [-2 2];
 
             consensus_numtets = 2;   % minimum # of tets for consensus event detection
@@ -174,7 +192,10 @@
 
         case 'wavelets4-300Hz'
             waveSet = '4-300Hz';
-            wp = getWaveParams(waveSet);
+%             wp = getWaveParams(waveSet);
+        case '4-300Hz_focusSWR'
+            waveSet = '4-300Hz_focusSWR';
+%             wp = getWaveParams(waveSet);
 
         case 'reactivationPLTH'
             iterator = 'singleepochanal';
@@ -190,24 +211,7 @@
             iterator = 'singlecellanal';
             datatypes = {'spikes'};
             options = {'savefigas', 'png', 'eventName',eventName};
-
-        case 'lickbouts'
-            maxIntraBurstILI = 0.25;  % max burst ili threshold in seconds
-            minBoutLicks = 3; %filter out bouts with less than boutNum licks
-            % timefilt1: lick bouts, 
-            timefilter{end+1} = {'getLickBout', '($lickBout == 1)', ...
-                'maxIntraBurstILI', maxIntraBurstILI, ...
-                'minBoutLicks', minBoutLicks};
-
-        case 'nolickbouts'
-            % timefilt: immoble not lick bout
-            maxIntraBurstILI = 0.25;  % max burst ili threshold in seconds
-            minBoutLicks = 3; % filter out bouts with less than boutNum licks
-            minTimeFromLick = .5; % time duration from closest lick
-            timefilter{end+1} = {'getLickBout', ...
-                sprintf('($lickBout == 0) & ($timeFromLick > %d)',minTimeFromLick), ...
-                'maxIntraBurstILI', maxIntraBurstILI, 'minBoutLicks', minBoutLicks};
-
+
         case 'dfa_lickphaseSUclustering'
             eventName = 'lick';
             filtfunction = 'dfa_lickphaseSUclustering';
@@ -319,7 +323,7 @@
             LFPtypes = {'eeggnd', 'eeg'};%, 'theta', 'ripple'}; %'lowgamma', 'fastgamma', 'ripple'};
             LFPrangesHz = {'1-400', '1-400'};%, '6-9', '140-250'}; %'20-50', '65-140',};
             srate = 1500;
-            win = [2 2];
+            win = [1.5 1.5];
             time = -win(1):1/srate:win(2);
             TF = 1;
             eventtype = 'rippleskons';

diff --git a/DFFunctions/load_plotting_params.m b/DFFunctions/load_plotting_params.m
@@ -34,13 +34,18 @@
             SupFontS = 12;
             tickSz = 8; % tick labels font size
             set(0,'defaultAxesFontSize',10)
+        case 'areaTFspect'
+%             area = {{'mec', 'sup'}, {'mec', 'deep'}, {'ca1', 'd'}};
+            position = [.1 .1 .8 .5];
+            win = [-.5 .5];
+            usecolormap = hot;
+            tickFsize = 8;
         case 'ReactivationStrength'
             position = [.1 .1 .5 .5];
-
+            winSE = [-2 2];
         case 'wetVDryILIphaseSWR'
             position = [.1 .1 .5 .5];
 
-
         case 'FeatureTracking'
             SpVt = 0.00;
             Msz = 30;
@@ -102,7 +107,7 @@
             stitFsize = 16;
             tickFsize = 8;
             sfTitFsize = 14;
-            pwin = [1 1];
+            pwin = [.5 .5];
             usecolormap = 'bone';
             contourRes = 40;
             position = [.1 .1 1 1];

diff --git a/DFFunctions/save_data.m b/DFFunctions/save_data.m
@@ -23,6 +23,7 @@ function save_data(F, filtOutputDirectory, filename, varargin)
     end
 end
 tic
+fprintf('saving...\n')
 if per_animal
     for an = 1:length(animals)
         if strcmp(filtOutputDirectory, 'filterframework')

diff --git a/DFFunctions/stack_riptriglfp.m b/DFFunctions/stack_riptriglfp.m
@@ -1,5 +1,6 @@
 function out = stack_riptriglfp(data, Fp, varargin)
 % make single data matrix per animal
+% exclude any events containing at least one nan in any tetrode
 
 % compile results from F(anim).output{dayep}.data{lfptype}{ripN: ntXsamp}
 % into out(ian).data{lfptype}(ntrode x sample x ripple)
@@ -11,69 +12,73 @@
 end
 
 out = struct;
-for ian = 1:length(data)
-    out(ian).animal = data(ian).animal{3};
+for ian = 1:length(data) % per animal
+    animal = data(ian).animal{3};
+    out(ian).animal = animal;
     out(ian).lfptypes = data(ian).datafilter_params.LFPtypes;
     try
         out(ian).LFPrangesHz = data(ian).datafilter_params.LFPrangesHz;
     catch
-    end    
+    end
     out(ian).data_dims = {'ntrode', 'sample', 'ripple'};
-    out(ian).dayeps = data(ian).epochs{1, 1};
-    out(ian).data = {};
-    out(ian).numrips_perep = {};
-    for ide = 1:length(out(ian).dayeps(:,1)) % per epoch
-        day = out(ian).dayeps(ide,1);
-        epoch = out(ian).dayeps(ide,2);
-        try
-        if isempty(data(ian).output{ide}.data)
-            fprintf('missing data %s %d %d\n', out(ian).animal, day, epoch);
-            continue
-        end
-        catch
-            fprintf('missing data %s %d %d\n', out(ian).animal, day, epoch);
-            continue
-        end
-        for t = 1:length(out(ian).lfptypes) % LFP type.. 
-            % TODO: need to seperate out the different lfp bands so i can load per type
-            % ntrode x sample X ripple
-            tmp = data(ian).output{1,ide}.data{t};
-            if isempty(tmp)
+    dayeps = data(ian).epochs{1, 1};
+    out(ian).dayeps = dayeps;
+    for t = 1:length(out(ian).lfptypes) % LFP type..
+        for ide = 1:size(dayeps,1) % per epoch
+            day = dayeps(ide,1);
+            epoch = dayeps(ide,2);
+            try
+                epData = data(ian).output{1,ide}.data{t};
+                if isempty(epData)
+                    fprintf('missing data %s %d %d\n', animal, day, epoch);
+                    continue
+                end
+            catch
+                fprintf('missing data %s %d %d\n', animal, day, epoch);
                 continue
-            else
-                out(ian).data{t}{1,ide} = cat(3,tmp{:}); % concat ripples into 3rd dim
             end
+            epDataTnsr = cat(3,epData{:}); % concat into [ntrode x sample x event]
+            % exclude any events with a nan
+            vEvents = find(squeeze(all(all(~isnan(epDataTnsr),1),2)));
+            numValid = numel(vEvents);
+            if max(vEvents) ~= numValid
+                fprintf('D:%d E:%d LFP:%s excluding bc nan: %d of %d events\n', ...
+                    day, epoch, out(ian).lfptypes{t}, max(vEvents)-length(vEvents), length(vEvents));
+            end
+            dataTnsr{1,ide} = epDataTnsr(:,:,vEvents);
+            evStartIdx{ide,1} = data(ian).output{ide}.eventStartIndices(vEvents);
+            evEndIdx{ide,1} = data(ian).output{ide}.eventEndIndices(vEvents);
+
+            % collect day/epoch level info
+            numEvPerEp{ide,1} = numValid;
+            evDay{ide,1} = day*ones(numValid, 1);
+            evEpoch{ide,1} = epoch*ones(numValid, 1);
+            evStart{ide,1} = data(ian).output{ide}.LFPtimes(evStartIdx{ide});
+            evEnd{ide,1} = data(ian).output{ide}.LFPtimes(evEndIdx{ide});
         end
-        % collect day/epoch level info
-        out(ian).numrips_perep{ide,1} = length(out(ian).data{t}{ide}(1,1,:));
-        out(ian).day{ide,1} = day*ones(length(out(ian).data{t}{ide}(1,1,:)), 1);
-        out(ian).epoch{ide,1} = epoch*ones(length(out(ian).data{t}{ide}(1,1,:)), 1);
-        out(ian).ripStartIdx{ide,1} = data(ian).output{ide}.eventStartIndices;
-        out(ian).ripEndIdx{ide,1} = data(ian).output{ide}.eventEndIndices;
-        out(ian).ripStartTime{ide,1} = data(ian).output{ide}.LFPtimes(...
-            data(ian).output{ide}.eventStartIndices);
-        out(ian).ripEndTime{ide,1} = data(ian).output{ide}.LFPtimes(...
-            data(ian).output{ide}.eventEndIndices);
-    end
-    % collapse data across collected epochs into one matrix
-    for t = 1:length(out(ian).lfptypes) % LFP type
-        out(ian).data{t} = cell2mat(permute(out(ian).data{t}, [1 3 2])); % stack all rips
+        out(ian).time{t} = data(ian).datafilter_params.time;
+        out(ian).day{t} = cell2mat(evDay);
+        out(ian).epoch{t} = cell2mat(evEpoch);
+        out(ian).ntrodes{t} = unique(cell2mat(cellfun(@(x) x(:,3), {data(ian).output{ide}.index}, ...
+            'un', 0)'), 'stable');
+
+        out(ian).data{t} = cell2mat(permute(dataTnsr, [1 3 2])); % stack all rips
+        out(ian).numEvPerEp{t} = cell2mat(numEvPerEp);
+        out(ian).evStartIdx{t} = cell2mat(evStartIdx);
+        out(ian).evEndIdx{t} = cell2mat(evEndIdx);
+        out(ian).evStart{t} = cell2mat(evStart);
+        out(ian).evEnd{t} = cell2mat(evEnd);
     end
-    out(ian).ntrodes = unique(cell2mat(cellfun(@(x) x(:,3), {data(ian).output{ide}.index}, ...
-        'un', 0)'), 'stable');
-    out(ian).time = data(ian).datafilter_params.time;
-    out(ian).numrips_perep = cell2mat(out(ian).numrips_perep);
-    out(ian).day = cell2mat(out(ian).day);
-    out(ian).epoch = cell2mat(out(ian).epoch);
-    out(ian).ripStartIdx = cell2mat(out(ian).ripStartIdx);
-    out(ian).ripEndIdx = cell2mat(out(ian).ripEndIdx);
-    out(ian).ripStartTime = cell2mat(out(ian).ripStartTime);
-    out(ian).ripEndTime = cell2mat(out(ian).ripEndTime);
+%     % collapse data across collected epochs into one matrix
+%     for t = 1:length(out(ian).lfptypes) % LFP type
+%         out(ian).data{t} = cell2mat(permute(out(ian).data{t}, [1 3 2])); % stack all rips
+%     end
 end
 if saveout
-    save_data(out, Fp.paths.resultsDirectory,['riptriglfpstack_',Fp.epochEnvironment]);
+    save_data(out, Fp.paths.resultsDirectory,['riptriglfpstack_',Fp.env]);
 end
 
 end
 
 
+