diff --git a/DFAnalysis/dfa_reactivationPSTH.m b/DFAnalysis/dfa_reactivationPSTH.m index 711a6b57..552c72b3 100644 --- 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) %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 index 27533d24..5ccfb335 100644 --- 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 index 0c67d7a1..7c123851 100644 --- 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 index 06735fb4..ec051ca4 100644 --- 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 + diff --git a/Functions/computeAnalyticSignal.m b/Functions/computeAnalyticSignal.m index c36b79df..2c37f447 100644 --- a/Functions/computeAnalyticSignal.m +++ b/Functions/computeAnalyticSignal.m @@ -10,23 +10,24 @@ % Demetris Roumis 2019 saveOutput = 1; -uselfptype = 'eeg'; +lfptype = 'eeg'; waveSet = '4-300Hz'; -epochEnvironment = 'ep'; +env = 'ep'; if ~isempty(varargin) assign(varargin{:}); end for ian = 1:length(lfpstack) animal = lfpstack(ian).animal; - fprintf('computing Analytic Signal for %s %s\n', animal, uselfptype); + fprintf('computing Analytic Signal for %s %s\n', animal, lfptype); fprintf('waveSet: %s\n', waveSet); - wp = getWaveParams(waveSet); - itypeidx = find(strcmp(lfpstack(ian).lfptypes, uselfptype)); - dsampdata = lfpstack(ian).data{itypeidx}(:,1:wp.dsamp:end,:); + wp = getWaveParams(waveSet); + t = find(strcmp(lfpstack(ian).lfptypes, lfptype)); + + dsampdata = lfpstack(ian).data{t}(:,1:wp.dsamp:end,:); + nNtrodes = size(dsampdata,1); nevents = size(dsampdata,3); nsamps = size(dsampdata,2); - nNtrodes = size(dsampdata,1); nData = nsamps*nevents*nNtrodes; timeWin = wp.win(1):1/(wp.srate/wp.dsamp):wp.win(2); baseind(1,1) = dsearchn(timeWin',wp.basewin(1)); @@ -53,47 +54,47 @@ % Phase fprintf('getting phase \n'); - phaseout.ph = single(angle(as)); % convert to single to save diskspace - phaseout.wp = wp; - phaseout.animal = animal; - phaseout.day = lfpstack(ian).day; - phaseout.epoch = lfpstack(ian).epoch; - phaseout.ripStartTime = lfpstack(ian).ripStartTime; - phaseout.ripEndTime = lfpstack(ian).ripEndTime; - phaseout.ntrode = lfpstack(ian).ntrodes; - phaseout.frequency = wp.frex; - phaseout.lfptype = uselfptype; - phaseout.dsampsrate = wp.win(1):1/(wp.srate/wp.dsamp):wp.win(2); - phaseout.dsamp = wp.dsamp; - phaseout.time = timeWin; - phaseout.dims = {'ntrode', 'sample', 'ripple', 'frequency'}; + phaseout(ian).ph = single(angle(as)); % convert to single to save diskspace + phaseout(ian).wp = wp; + phaseout(ian).animal = animal; + phaseout(ian).day = lfpstack(ian).day{t}; + phaseout(ian).epoch = lfpstack(ian).epoch{t}; + phaseout(ian).evStart = lfpstack(ian).evStart{t}; + phaseout(ian).evEnd = lfpstack(ian).evEnd{t}; + phaseout(ian).ntrode = lfpstack(ian).ntrodes{t}; + phaseout(ian).frequency = wp.frex; + phaseout(ian).lfptype = lfptype; + phaseout(ian).dsampsrate = wp.win(1):1/(wp.srate/wp.dsamp):wp.win(2); + phaseout(ian).dsamp = wp.dsamp; + phaseout(ian).time = timeWin; + phaseout(ian).dims = {'ntrode', 'sample', 'event', 'frequency'}; % Power fprintf('getting power \n'); - powerout.pwr = single(abs(as).^2); - powerout.wp = wp; - powerout.animal = animal; - powerout.day = lfpstack(ian).day; - powerout.epoch = lfpstack(ian).epoch; - powerout.ripStartTime = lfpstack(ian).ripStartTime; - powerout.ripEndTime = lfpstack(ian).ripEndTime; - powerout.ntrode = lfpstack(ian).ntrodes; - powerout.frequency = wp.frex; - powerout.lfptype = uselfptype; - powerout.dsampsrate = wp.win(1):1/(wp.srate/wp.dsamp):wp.win(2); - powerout.dsamp = wp.dsamp; - powerout.time = timeWin; - powerout.dims = {'ntrode', 'sample', 'ripple', 'frequency'}; - fprintf('%.02f seconds to compute pwr,ph AS\n', toc); + powerout(ian).pwr = single(abs(as).^2); + powerout(ian).wp = wp; + powerout(ian).animal = animal; + powerout(ian).day = lfpstack(ian).day{t}; + powerout(ian).epoch = lfpstack(ian).epoch{t}; + powerout(ian).evStart = lfpstack(ian).evStart{t}; + powerout(ian).evEnd = lfpstack(ian).evEnd{t}; + powerout(ian).ntrode = lfpstack(ian).ntrodes{t}; + powerout(ian).frequency = wp.frex; + powerout(ian).lfptype = lfptype; + powerout(ian).dsampsrate = wp.win(1):1/(wp.srate/wp.dsamp):wp.win(2); + powerout(ian).dsamp = wp.dsamp; + powerout(ian).time = timeWin; + powerout(ian).dims = {'ntrode', 'sample', 'event', 'frequency'}; + fprintf('%.02f seconds to compute AS, pwr, ph \n', toc); %% ---------------- Save Output --------------------------------------------------- if saveOutput == 1 fprintf('saving AS\n') andef = animaldef(lower('Demetris')); - save_data(phaseout, sprintf('%s/analyticSignal/', andef{2}), ... - sprintf('AS_waveSet-%s_%s_%s_phase', waveSet, uselfptype, epochEnvironment)) - save_data(powerout, sprintf('%s/analyticSignal/', andef{2}), ... - sprintf('AS_waveSet-%s_%s_%s_power', waveSet, uselfptype, epochEnvironment)) + save_data(phaseout(ian), sprintf('%s/analyticSignal/', andef{2}), ... + sprintf('AS_waveSet-%s_%s_%s_phase', waveSet, lfptype, env)) + save_data(powerout(ian), sprintf('%s/analyticSignal/', andef{2}), ... + sprintf('AS_waveSet-%s_%s_%s_power', waveSet, lfptype, env)) end end end @@ -110,7 +111,10 @@ gaus_win = exp(-timeWin.^2./(2*bn^2)); % gaussian wavelet = sine_wave .* gaus_win; % gaus_fwhm = exp(-(4*log(2)*timeWin).^2/fwhm(fi).^2); - wavefft = fft(wavelet,nConv2pow); % fft of the wavelet, pad to next power of 2 +% i think this will ensure that the wavelett and data fft have same + % freq resolution + wavefft = fft(wavelet,nConv2pow); % fft of the wavelet, pad to next power of 2 of data + % normalize wavelet to a maximum of 1 to ensure convolution units are same as data waveletFFT(:,fi) = (wavefft ./ max(wavefft))'; fprintf('creating wavelet for freq %d of %d \n',fi,wp.numfrex); diff --git a/Functions/getPower.m b/Functions/getPower.m index c472a539..0c0ccd5c 100644 --- a/Functions/getPower.m +++ b/Functions/getPower.m @@ -9,7 +9,7 @@ me = animaldef('Demetris'); noiseEvents = struct; lfptype = 'eeg'; -expvars = {'onlywdays'}; %,'rewarded', 'unrewarded', 'inbound' , 'outbound', 'proximalWell', ... +expvars = {}; %,'rewarded', 'unrewarded', 'inbound' , 'outbound', 'proximalWell', ... % 'distalWell'}; saveout = 1; run_perm = 1; @@ -17,7 +17,9 @@ if ~isempty(varargin) assign(varargin{:}); end - +if isempty(expvars) + expvars = expvarCat(1).expvars; +end for ian = 1:length(rawpwr) wp = getWaveParams(Fp.waveSet); animal = Fp.animals{ian}; @@ -58,7 +60,7 @@ out(ian).meandbpower = meandbpower; if saveout - save_data(out, [pconf.andef{2},outdir,'/'], [outdir, '_', Fp.epochEnvironment]) + save_data(out, [pconf.andef{2},outdir,'/'], [outdir, '_', Fp.env]) end end end diff --git a/Functions/getWaveParams.m b/Functions/getWaveParams.m index 9f7b9f37..9746d6fd 100644 --- a/Functions/getWaveParams.m +++ b/Functions/getWaveParams.m @@ -64,6 +64,30 @@ % frexres = 2; % numfrex = floor((max_freq - min_freq)/frexres); % frex = logspace(log10(min_freq),log10(max_freq),numfrex); + +case '4-300Hz_focusSWR' + srate = 1500; + win = [-1.5 1.5]; + dsamp = 2; + basewin = [-.5 -.2]; % in seconds, period before event start to use as baseline + pval = 0.05; + voxel_pval = 0.01; + mcc_voxel_pval = 0.05; + mcc_cluster_pval = 0.05; + zval = abs(norminv(pval)); % convert p-value to Z value + n_permutes = 100; + + min_freq = 4; + max_freq = 300; + numfrex = 25; + range_cycles = [3 14]; % more decreases temp, increases freq precision + + frex = logspace(log10(min_freq),log10(max_freq),numfrex); + nWavecycles = logspace(log10(range_cycles(1)),log10(range_cycles(end)),numfrex); %number of wavelet cycles per freq +% hws = (length(timeWin)-1)/2; % half_wave_size + + + case '4-30Hz' srate = 1500; win = [-2 2]; diff --git a/Functions/makeExpvarCatDesignMat.m b/Functions/makeExpvarCatDesignMat.m index 48f18adc..27ec3f5b 100644 --- a/Functions/makeExpvarCatDesignMat.m +++ b/Functions/makeExpvarCatDesignMat.m @@ -4,75 +4,82 @@ % return design matrix of experiment variables [rip var] % Demetris Roumis 2019 pconf = paramconfig; - expvars = {'onlywdays', 'rewarded', 'unrewarded', 'inbound' , 'outbound', ... + expvars = {'all', 'rewarded', 'unrewarded', 'inbound' , 'outbound', ... 'proximalWell', 'distalWell', 'rewarded_outbound', 'rewarded_inbound'}; saveout = 1; outdir = 'expvarCat'; defaults = {'wtrackdays', 'excludeNoise','excludePriorFirstWell'}; + lfptype = 'eeg'; if ~isempty(varargin) assign(varargin{:}) end - + fprintf('defaults: %s\n', defaults{:}) outpath = [pconf.andef{2},outdir,'/']; for ian = 1:length(lfpstack) + t = find(strcmp(lfpstack(ian).lfptypes, lfptype)); animal = lfpstack(ian).animal; out(ian).animal = animal; out(ian).dims = {'ripple', 'expvar'}; - out(ian).ripStartTime = lfpstack(ian).ripStartTime; - out(ian).ripEndTime = lfpstack(ian).ripEndTime; - out(ian).dayeps = [lfpstack(ian).day lfpstack(ian).epoch]; + out(ian).evStart = lfpstack(ian).evStart{t}; + out(ian).evEnd = lfpstack(ian).evEnd{t}; + out(ian).dayeps = [lfpstack(ian).day{t} lfpstack(ian).epoch{t}]; out(ian).expvars = expvars; - out(ian).dm = zeros(length(out(ian).ripStartTime), length(expvars)); + out(ian).dm = zeros(length(out(ian).evStart), length(expvars)); Fp = struct; for ss = 1:length(expvars) - Fp.add_params = defaults; + Fp.params = defaults; switch expvars{ss} -% case 'all' -% Fp.add_params = {'excludeNoise','excludePriorFirstWell'}; - case 'onlywdays' + case 'all' % all defaults + case 'lickbouts' + Fp.params{end+1} = 'lickbouts'; + case 'nolickbouts' + Fp.params{end+1} = 'nolickbouts'; case 'rewarded' - Fp.add_params{end+1} = 'correcttrials'; + Fp.params{end+1} = 'correcttrials'; case 'unrewarded' - Fp.add_params{end+1} = 'errortrials'; + Fp.params{end+1} = 'errortrials'; case 'outbound' - Fp.add_params{end+1} = 'outbound'; + Fp.params{end+1} = 'outbound'; case 'inbound' - Fp.add_params{end+1} = 'inbound'; + Fp.params{end+1} = 'inbound'; case 'rewarded_outbound' - Fp.add_params{end+1} = 'correcttrials'; - Fp.add_params{end+1} = 'outbound'; + Fp.params{end+1} = 'correcttrials'; + Fp.params{end+1} = 'outbound'; case 'unrewarded_outbound' - Fp.add_params{end+1} = 'errortrials'; - Fp.add_params{end+1} = 'outbound'; + Fp.params{end+1} = 'errortrials'; + Fp.params{end+1} = 'outbound'; case 'rewarded_inbound' - Fp.add_params{end+1} = 'correcttrials'; - Fp.add_params{end+1} = 'inbound'; + Fp.params{end+1} = 'correcttrials'; + Fp.params{end+1} = 'inbound'; case 'unrewarded_inbound' - Fp.add_params{end+1} = 'errortrials'; - Fp.add_params{end+1} = 'inbound'; + Fp.params{end+1} = 'errortrials'; + Fp.params{end+1} = 'inbound'; case 'proximalWell' - Fp.add_params{end+1} = 'proximalWell'; + Fp.params{end+1} = 'proximalWell'; case 'distalWell' - Fp.add_params{end+1} = 'distalWell'; + Fp.params{end+1} = 'distalWell'; otherwise error('condition string %s unrecognized', expvars{ss}); end - Fp = load_filter_params(Fp, 'add_params', Fp.add_params); + fprintf('============ %s =============\n', expvars{ss}); + Fp = load_filter_params(Fp); F = createfilter('animal', {out(ian).animal}, 'epochs', ... Fp.epochfilter, 'excludetime', Fp.timefilter); for de = 1:length(F.epochs{1}(:,1)) day = F.epochs{1}(de,1); epoch = F.epochs{1}(de,2); + % events in current epoch ieprips = ismember(out(ian).dayeps, [day epoch], 'rows'); + % events in current timefilter out(ian).dm(ieprips,ss) = ... - ~isExcluded(out(ian).ripStartTime(ieprips), ... + ~isExcluded(out(ian).evStart(ieprips), ... F.excludetime{1}{de}); end end end if saveout - save_data(out, outpath, [outdir,'_',Fp.epochEnvironment]); + save_data(out, outpath, [outdir,'_',Fp.env]); end end diff --git a/notebooks/reactivationBurstSwr_20191026.m b/notebooks/reactivationBurstSwr_20191026.m index ac54d86e..93a42478 100644 --- a/notebooks/reactivationBurstSwr_20191026.m +++ b/notebooks/reactivationBurstSwr_20191026.m @@ -6,27 +6,31 @@ swr's are real by comparing the neural activity during the burst swr to non-swr times in the lick burst.. and further, that there is some structure during burst swr's that is different than the activity during nonswr bursts -at the same ILI phase.. +at the same ILI phase.. -- compare intra-burst SU PSwrTH to PLickTH (basically just a SU modulation).. +- compare intra-burst SU PSwrTH to PLickTH (basically just a SU modulation).. - compare intra-burst Reactivation PSwrTH to PLickTH - i think the easier first step is to do the psth relative to time rather -than ILIphase.. +than ILIphase.. %} -create_filter = 1; -run_ff = 1; +pconf = paramconfig; +create_filter = 0; +run_ff = 0; load_ffdata = 0; -createSummaryData = 0; - -plotfigs = 0; -pausefigs = 1; -savefigs = 0; +createSummaryData = 1; + +plotfigs = 1; +plotETA = 1; +plotTrace = 0; +plotPCdemo = 0; +pausefigs = 0; +savefigs = 1; %% FF -Fp.animals = {'D10', 'D12', 'JZ1', 'JZ2', 'JZ4'}; %, 'JZ2', 'JZ4'}; +Fp.animals = {'D10', 'D13', 'JZ1', 'JZ2', 'JZ4'}; Fp.filtfunction = 'dfa_reactivationPSTH'; Fp.params = {'>4cm/s', 'ca1SU', 'wtrackdays', 'excludePriorFirstWell', ... - Fp.filtfunction}; %'exemplar_wepochs', + Fp.filtfunction}; %'exemplar_wepochs', Fp = load_filter_params(Fp); %% if create_filter @@ -50,10 +54,11 @@ %% create summary data if createSummaryData + D = struct; for i = 1:length(F) data = [F(i).output{:}]; D(i).animal = F(i).animal; - D(i).etTime = F(i).output{1}.reactPSTHtime; + D(i).etTime = F(i).output{end}.reactPSTHtime; D(i).swrReactETAfull = cell2mat(cellfun(@(x) x', {data.swrReactETAfull}', 'un', 0)')'; D(i).swrReactETAfullMean = nanmean(D(i).swrReactETAfull); D(i).swrReactETAfullShufs = cell2mat(cellfun(@(x) x', {data.swrReactETAfullShufs}', 'un', 0)')'; @@ -73,82 +78,120 @@ %% plot if plotfigs - %% plot Full ETA w sem errorbars - figname = 'ReactivationStrength'; - for i = 1:length(D) - Pp=load_plotting_params({'defaults',figname}, 'pausefigs', pausefigs, ... - 'savefigs', savefigs); - % plot swr ETA with ebars, shuff - time = D(i).etTime; - subplot(3,1,1) - plot(time, D(i).swrReactETAfullMean,'k','linewidth',3); - hold on; - plot(time, D(i).swrReactETAfullShufMean,'r','linewidth',3); - errorbar(time, D(i).swrReactETAfullMean, nanstd(D(i).swrReactETAfull)/... - sqrt(size(D(i).swrReactETAfull,1)), 'k','linewidth',0.5) - errorbar(time, D(i).swrReactETAfullShufMean, nanstd(D(i).swrReactETAfullShufs)/... - sqrt(size(D(i).swrReactETAfullShufs,1)), 'r','linewidth',0.5) - ylabel('rxn strength') - xticks([]); -% xlabel('Time (s)') - title('SWR') - hold off; - - % plot swr burst ETA with ebars, shuff - subplot(3,1,2) - plot(time, D(i).swrBurstReactETAfullMean,'k','linewidth',3); - hold on; - plot(time, D(i).swrBurstReactETAfullShufMean,'r','linewidth',3); - errorbar(time, D(i).swrBurstReactETAfullMean,nanstd(D(i).swrBurstReactETAfull)/... - sqrt(size(D(i). swrBurstReactETAfull,1)), 'k','linewidth',0.5) - errorbar(time, D(i).swrBurstReactETAfullShufMean, ... - nanstd(D(i).swrBurstReactETAfullShufs)/... - sqrt(size(D(i).swrBurstReactETAfullShufs,1)),'r','linewidth',0.5) - ylabel('rxn strength') -% xlabel('Time (s)') - xticks([]); - title('SWR in Burst') - hold off; - - % plot lick ETA with ebars, shuff - subplot(3,1,3) - plot(time, D(i).lickReactETAfullMean,'k','linewidth',3); - hold on; - plot(time,D(i).lickReactETAfullShufMean,'r','linewidth',3); - errorbar(time,D(i).lickReactETAfullMean,nanstd(D(i).lickReactETAfull)/... - sqrt(size(D(i).lickReactETAfull,1)),'k','linewidth',0.5) - errorbar(time,D(i).lickReactETAfullShufMean,nanstd(D(i).lickReactETAfullShufs)/... - sqrt(size(D(i).lickReactETAfullShufs,1)),'r','linewidth',0.5) - ylabel('rxn strength') - xlabel('Time (s)') - title('lick in Burst') - hold off; - %% - stit = [figname D(i).animal{3} Fp.env]; - setSuperAxTitle(stit); - if pausefigs - pause; - end - if savefigs - save_figure([pconf.andef{4} '/' figname '/'], stit); + if plotETA + %% plot Full ETA w sem errorbars + figname = 'ReactivationStrength'; + for i = 1:length(D) + Pp=load_plotting_params({'defaults',figname}, 'pausefigs', pausefigs, ... + 'savefigs', savefigs); +% plot swr ETA with ebars, shuff + time = D(i).etTime'; + subplot(1,1,1) + winidx = knnsearch(time, Pp.winSE'); + time = time(winidx(1):winidx(2)); + Rs = D(i).swrReactETAfull(:,winidx(1):winidx(2)); + R = D(i).swrReactETAfullMean(:,winidx(1):winidx(2)); + shs = D(i).swrReactETAfullShufs(:,winidx(1):winidx(2)); + sh = D(i).swrReactETAfullShufMean(:,winidx(1):winidx(2)); + Rsem = nanstd(Rs)/sqrt(size(Rs,1)); + shssem = nanstd(shs)/sqrt(size(shs,1)); + plot(time, R, 'g', 'linewidth', 1.5) +% errorbar(time, R, Rsem, 'b','linewidth',0.5,'CapSize',0) + hold on; + plot(time, sh, 'k', 'linewidth', 1.5) +% errorbar(time, sh, shssem,'k', 'linewidth',0.5, 'CapSize',0) + fill([time; flipud(time)],[R'+Rsem'; flipud(R'-Rsem')],'g','linestyle','none', ... + 'facealpha', .1) + fill([time; flipud(time)]',[sh-shssem flipud(sh+shssem)]','k','linestyle','none', ... + 'facealpha', .1) +% axis tight +% line([0 0], ylim, 'color', 'k', 'linestyle', '--', 'linewidth', .5) + ylabel('rxn strength') +% xticks([]); +% title('SWR') +% hold off; + + %% plot swr burst ETA with ebars, shuff + subplot(1,1,1) + Rs = D(i).swrBurstReactETAfull(:,winidx(1):winidx(2)); + R = D(i).swrBurstReactETAfullMean(:,winidx(1):winidx(2)); + shs = D(i).swrBurstReactETAfullShufs(:,winidx(1):winidx(2)); + sh = D(i).swrBurstReactETAfullShufMean(:,winidx(1):winidx(2)); + Rsem = nanstd(Rs)/sqrt(size(Rs,1)); + shssem = nanstd(shs)/sqrt(size(shs,1)); + plot(time, R, 'b', 'linewidth', 1.5) +% errorbar(time, R, Rsem, 'b','linewidth',0.5,'CapSize',0) + hold on; + plot(time, sh, 'k', 'linewidth', 1.5) +% errorbar(time, sh, shssem,'k', 'linewidth',0.5, 'CapSize',0) + fill([time; flipud(time)],[R'+Rsem'; flipud(R'-Rsem')],'b','linestyle','none', ... + 'facealpha', .1) + fill([time; flipud(time)]',[sh-shssem flipud(sh+shssem)]','k','linestyle','none', ... + 'facealpha', .1) + ylabel('rxn strength') +% xticks([]); + axis tight + line([0 0], ylim, 'color', 'k', 'linestyle', '--', 'linewidth', .5) + ylabel('rxn strength') + xlabel('Time from SWR (s)') + title('rxn SWR ETA') +% hold off; + %% plot lick ETA with ebars, shuff + subplot(1,1,1) + Rs = D(i).lickReactETAfull(:,winidx(1):winidx(2)); + R = D(i).lickReactETAfullMean(:,winidx(1):winidx(2)); + shs = D(i).lickReactETAfullShufs(:,winidx(1):winidx(2)); + sh = D(i).lickReactETAfullShufMean(:,winidx(1):winidx(2)); + Rsem = nanstd(Rs)/sqrt(size(Rs,1)); + shssem = nanstd(shs)/sqrt(size(shs,1)); + plot(time, R, 'm', 'linewidth', 1.5) +% errorbar(time, R, Rsem, 'b','linewidth',0.5,'CapSize',0) + hold on; + plot(time, sh, 'k', 'linewidth', 1.5) +% errorbar(time, sh, shssem,'k', 'linewidth',0.5, 'CapSize',0) + fill([time; flipud(time)],[R'+Rsem'; flipud(R'-Rsem')],'m','linestyle','none', ... + 'facealpha', .1) + fill([time; flipud(time)]',[sh-shssem flipud(sh+shssem)]','k','linestyle','none', ... + 'facealpha', .1) + axis tight + line([0 0], ylim, 'color', 'k', 'linestyle', '--', 'linewidth', .5) + ylabel('rxn strength') + xlabel('Time from lickDin (s)') + title('rxn Lick ETA') +% legend({'allSWR', 'allSWRshuff', 'burstSWR', 'burstSWRshuff', 'lick', 'lickshuf'}); + hold off; + %% + allAxesInFigure = findall(gcf,'type','axes'); + linkaxes(allAxesInFigure, 'y'); + try + stit = [figname D(i).animal{3} Fp.env]; + catch + stit = [figname D(i).animal Fp.env]; + end + setSuperAxTitle(stit); + if pausefigs + pause; + end + if savefigs + save_figure([pconf.andef{4} '/' figname], stit); + end end end - - %% plot time snippet of spikes + zbinSpikes traces + reactivation result (full,pc) - exWin = []; - % plot spikes in time window - - % plot zbinspikes traces in time window - % plot like lfp - - % plot reactivation result in time window - % reactPerPC (1 trace per sig pc) - % reactFull - + if plotTrace + %% plot time snippet of spikes + zbinSpikes traces + reactivation result (full,pc) + exWin = []; + % plot spikes + + % plot full reactivation trace + + % plot lick times + + % plot swr times + end %% plot per PC method Demo % plot the steps and examples of the main parts % templateCC Eigvec, val, sig PC's, match z spikes, per pc react - if 0 + if plotPCdemo figure subplot(2,3,1); imagesc(zCCtemplateNoDiag); ax = gca; ax.YDir = 'normal'; title('zCCtemplateNoDiag') ylabel('cell') @@ -178,9 +221,3 @@ - - - - - - diff --git a/notebooks/swrtrigspect_20191030.m b/notebooks/swrtrigspect_20191030.m new file mode 100644 index 00000000..cde8e47a --- /dev/null +++ b/notebooks/swrtrigspect_20191030.m @@ -0,0 +1,258 @@ + + +%{ +just get a per-animal, all-swr-trig spect.. focusing on the higher +frequencies + +%} +% get data +pconf = paramconfig; +create_filter = 0; +run_ff = 0; +load_ffdata = 0; +stack_swrLFP = 0; +load_swrLFPstack = 0; +% run cmwc +make_ASig = 0; +load_rawpwr = 0; +load_phase = 0; +% create condition design mat +make_expvarCat = 0; +load_expvarCat = 0; +makeNoiseEvents = 0; +loadNoiseEvents = 0; +% compute per condition +make_expvarCatMeanPwr = 0; +load_expvarCatMeanPwr = 0; +% combine per area +combineArea = 0; +% plot +plot_expvarCatMeanPwr = 0; +plot_ByArea = 1; +pausefigs = 0; +savefigs = 1; + +%% +Fp.animals = {'D13'}; %{'D10','D13', 'JZ1', 'JZ2', 'JZ3', 'JZ4'}; %, 'JZ4'}; %;{'D10'}; % +Fp.filtfunction = 'dfa_riptriglfp'; +% Fp.add_params = {'sleepwtrackdays', 'excludeNoise', '<4cm/s', 'wavelets4-300Hz'}; +Fp.params = {'referenced','wtrackdays', 'excludePriorFirstWell','<4cm/s', ... + '4-300Hz_focusSWR','excludeAfterLastWell', Fp.filtfunction}; +Fp = load_filter_params(Fp); +wp = getWaveParams(Fp.waveSet); +rs = ''; % ripstate +area = {{'ca1', 'd'}, {'mec', 'deep'}, {'mec', 'supf'}}; +%% +if create_filter + F = createfilter('animal', Fp.animals, 'epochs', Fp.epochfilter, 'eegtetrodes', ... + Fp.tetfilter, 'excludetime', Fp.timefilter, 'iterator', Fp.iterator); + F = setfilterfunction(F, Fp.filtfunction, Fp.datatypes, Fp.options{:}); +end +if run_ff + F = arrayfun(@(x) setfield(F(x),'datafilter_params',Fp),1:length(F), 'un', 1); + F = runfilter(F); + save_data(F, Fp.paths.filtOutputDirectory, Fp.paths.filenamesave, ... + 'filetail', ['_' Fp.env]); +end +if load_ffdata + F = load_data(Fp.paths.filtOutputDirectory, Fp.paths.filenamesave, Fp.animals, ... + 'filetail', ['_' Fp.env]); +end +%% vectorize swr-trig lfp +if stack_swrLFP + lfpstack = stack_riptriglfp(F, Fp); +end +if load_swrLFPstack + lfpstack = load_data(Fp.paths.resultsDirectory, ... + ['riptriglfpstack_',Fp.env], Fp.animals); +end +%% make rawpwr (all trials) [ntrode time rip freq] +if make_ASig + [rawpwr, ~] = computeAnalyticSignal(lfpstack, 'waveSet', Fp.waveSet, 'saveOutput',1, ... + 'lfptype', Fp.uselfptype, 'env', Fp.env); +end +if load_rawpwr + rawpwr = load_data(sprintf('%s/analyticSignal/', pconf.andef{2}), ... + sprintf('AS_waveSet-%s_%s_%s_power', wp.waveSet, Fp.uselfptype, Fp.env), ... + Fp.animals); +end +%% make design mat to slice the rawpwr trials +if make_expvarCat + outdir = 'expvarCatBurst'; + expvarCat = makeExpvarCatDesignMat(lfpstack, 'outdir', outdir, 'expvars', {'all', ... + 'lickbouts', 'nolickbouts'}, 'lfptype', Fp.uselfptype); +end +if load_expvarCat + outdir = 'expvarCatBurst'; + outpath = [pconf.andef{2},outdir,'/']; + expvarCat = load_data(outpath, [outdir,'_',Fp.env], Fp.animals); +end +%% exclude noise +if makeNoiseEvents + noiseEvents = make_noiseEvents(lfpstack); +end +if loadNoiseEvents + noiseEvents = load_data('filterframework','noiseEvents', Fp.animals, 'animpos', 0); +end +%% get mean power per condition +if make_expvarCatMeanPwr % :expvarCat @mean /ntTF $time + expvarCatMeanPwr = getPower(expvarCat, rawpwr, Fp, 'run_perm', 0, 'noiseEvents',... + noiseEvents); +end +if load_expvarCatMeanPwr + outdir = 'expvarCatMeanPwr'; + outpath = [pconf.andef{2},outdir,'/']; + expvarCatMeanPwr = load_data(outpath, [outdir,'_', Fp.env] ,Fp.animals); +end +%% mean per area per condition +if combineArea + for ani = 1:length(expvarCatMeanPwr) % for each animal + animal = expvarCatMeanPwr(ani).animal; + aninfo = animaldef(animal); + ntinfo = loaddatastruct(aninfo{2}, animal, 'tetinfo'); + ntrodes = evaluatefilter(ntinfo, 'strcmp($valid, ''yes'')'); + ntrodes = unique(ntrodes(:,3)); + for ia = 1:length(area) + ntsInArea = evaluatefilter(ntinfo,... + sprintf('isequal($area,''%s'') && isequal($subarea,''%s'')', area{ia}{1}, ... + area{ia}{2})); + ntsA = unique(ntsInArea(:,3)); + for iv = 1:length(expvarCatMeanPwr(ani).expvars) + areaData = expvarCatMeanPwr(ani).meandbpower{iv}.pwr_mean_db(ntsA,:,:); + meanPwrArea{ia}{iv} = squeeze(nanmean(areaData,1))'; + end + end + end +end + +%% swr example lfp trace over ca1 area all heatmap + + +%% plot per area +if plot_ByArea + figname = 'expvarCatMeanPwrByArea'; + for ani = 1:length(expvarCatMeanPwr) % for each animal + animal = expvarCatMeanPwr(ani).animal; + + for iv = 1:length(expvarCatMeanPwr(ani).expvars) + + Pp=load_plotting_params({'defaults','areaTFspect'}, 'pausefigs', pausefigs, ... + 'savefigs', savefigs); + f = 0; + for ia = 1:length(area) + f = f+1; + sf{ia} = subaxis(1,length(area),f); + frequency = round(expvarCatMeanPwr(ani).wp.frex); + time = wp.win(1):1/(wp.srate/wp.dsamp):wp.win(2); + trim = knnsearch(time',Pp.win(1)):knnsearch(time', Pp.win(2)); + ttime = time(trim); + tdata = meanPwrArea{ia}{iv}(:,trim); + contourf(ttime,frequency,tdata,40,'linecolor','none') + set(gca,'ydir','normal','yscale','log'); + colormap(Pp.usecolormap) + caxis(sf{ia}, 'auto') + cax = caxis; + caxis([0 cax(2)]) + ytickskip = 2:4:wp.numfrex; + set(gca,'ytick', round(wp.frex(ytickskip)), 'FontSize', Pp.tickFsize) + hold on + title(sprintf('%s %s',area{ia}{1},area{ia}{2})) + yl = ylim; + xl = xlim; + line([0 0], yl, 'Color', [0.8 0.8 0.8],'LineStyle','--', 'LineWidth', 1); + colorbar + end + ylabel(sf{1},'frequency Hz') + xlabel(sf{1}, 'time s') + %% + stit = sprintf('%s %s %s %s %s',figname, expvarCatMeanPwr(ani).expvars{iv}, animal, ... + Fp.env, Fp.waveSet); + setSuperAxTitle(stit); + if pausefigs + pause + end + if savefigs + save_figure([pconf.andef{4} figname], stit); + end + end + + end + end + +%% plot MeanPower varCat TFzmap /nt +if plot_expvarCatMeanPwr + figname = 'expvarCatMeanPwr'; + for ani = 1:length(expvarCatMeanPwr) % for each animal + Pp=load_plotting_params({'defaults','powerTFmap'}, 'pausefigs', pausefigs, ... + 'savefigs', savefigs); + animal = expvarCatMeanPwr(ani).animal; + aninfo = animaldef(animal); + ntinfo = loaddatastruct(aninfo{2}, animal, 'tetinfo'); + ntrodes = evaluatefilter(ntinfo, 'strcmp($valid, ''yes'')'); + ntrodes = unique(ntrodes(:,3)); + den = cellfetch(ntinfo, 'area'); + matidx = unique(den.index(:,3)); + matidx = circshift(matidx,-1); + anidx = find(strcmp({expvarCatMeanPwr.animal}, animal)); + evanidx = find(strcmp({expvarCatMeanPwr.animal}, animal)); + % exclude invalid tets + invalidtets = evaluatefilter(ntinfo, 'isequal($valid, ''no'')'); + invalidtets = unique(invalidtets(:,3)); + if isempty(rs) + rs = expvarCatMeanPwr(evanidx).expvars; + end + for i = 1:length(rs); + iv = find(strcmp(rs{i}, expvarCatMeanPwr(evanidx).expvars)); + numcols = 8; + numrows = 4; %ceil(length(ntrodes) / numcols); + for nti = 1:length(ntrodes) + nt = ntrodes(nti); + if ismember(nt, invalidtets) + continue + end + area = ntinfo{1}{1}{nt}.area; + subarea = ntinfo{1}{1}{nt}.subarea; + cann = ntinfo{1}{1}{nt}.cannula; + ntxy = ntinfo{1}{1}{nt}.ntxy; + if cann == 'ca1' + ntxy(1) = ntxy(1)+4; % offset ca1 to the right + end + ntp = (ntxy(2)-1)*8+ntxy(1); + sf = subaxis(numrows,numcols,ntp); + if isnumeric(subarea) + subarea = num2str(subarea); + end + ntidx = find(matidx == nt); + idata2plot = squeeze(... + expvarCatMeanPwr(anidx).meandbpower{iv}.pwr_mean_db(ntidx,:,:))'; + idata2plot = trim2win(idata2plot, Fp.srate, Pp.pwin, ... + 'dsamp', expvarCatMeanPwr(anidx).wp.dsamp); + time = linspace(-Pp.pwin(1), Pp.pwin(2), length(idata2plot(1,:))); + contourf(sf, time, wp.frex, idata2plot, Pp.contourRes, ... + 'linecolor','none'); + set(gca,'ydir','normal','yscale','log'); + colormap(hot); %Pp.usecolormap) + caxis(sf, 'auto') + hold on; + ytickskip = 2:4:wp.numfrex; + set(gca,'ytick', round(wp.frex(ytickskip)), 'FontSize', 8) + title(sprintf('%s%s nt%d',area,subarea,nt), 'FontSize',14,... + 'FontWeight',Pp.FontW, 'FontName', Pp.FontNm) + yl = ylim; + line([0 0], yl, 'Color', [0.8 0.8 0.8],'LineStyle','--', 'LineWidth', 1); + end + %% + % allAxesInFigure = findall(gcf,'type','axes'); + % linkaxes(allAxesInFigure, 'xy'); + stit = sprintf('%s %s %s %s %s',figname, expvarCatMeanPwr(anidx).expvars{iv}, animal, ... + Fp.env, Fp.waveSet); + setSuperAxTitle(stit); + if pausefigs + pause + end + if savefigs + save_figure([pconf.andef{4} figname], stit); + end + end + end +end