diff --git a/DFAnalysis/dfa_reactivationPSTH.m b/DFAnalysis/dfa_reactivationPSTH.m index 552c72b3..3ed39125 100644 --- a/DFAnalysis/dfa_reactivationPSTH.m +++ b/DFAnalysis/dfa_reactivationPSTH.m @@ -5,7 +5,7 @@ % uses the lick data to generate lickBurstIntervals -function out = dfa_reactivationPSTH(idx, templateIntervals, varargin) +function out = dfa_reactivationPSTH(idx, nonSwrIntervals, varargin) % check for required data in varargin reqData = {'spikes', 'lick', 'ca1rippleskons', 'cellinfo'}; @@ -18,7 +18,8 @@ % assign params numCellsThresh = 5; numSpikesThresh = 50; -numShuffs = 100; +minLickTimeFromSwr = 1; +numShufs = 100; bin = 0.01; % seconds fullModel = 1; perPC = 0; @@ -28,26 +29,25 @@ % init output out = init_out(idx, animal); +day = idx(1); +ep = idx(2); -%% Filter and bin All su spikes +%% Filter for SU su = evaluatefilter(cellinfo, cellfilter); numCells=size(su,1); if numCells0); +swrTimeIdx = swrTimeIdx(swrTimeIdx+max(idxWin)0); +lickBurstTime = lickBurstTime(lickBurstTime+max(idxWin)minLickTimeFromSwr,:); + +% exclude events, intervals too close to ep edges +[~, swrBurstTimeIdx] = histc(swrBurstTime(:,1), matchTimeEdges); +swrBurstTimeIdx = swrBurstTimeIdx(swrBurstTimeIdx+min(idxWin)>0); +swrBurstTimeIdx = swrBurstTimeIdx(swrBurstTimeIdx+max(idxWin)0); +lickBurstTimeIdx = lickBurstTimeIdx(lickBurstTimeIdx+max(idxWin)0); - swrTimeIdx = swrTimeIdx(swrTimeIdx+max(idxWin)0); - swrBurstTimeIdx = swrBurstTimeIdx(swrBurstTimeIdx+max(idxWin) 1); - lickBurstTime = lickBurstTime(lickBurstTime+min(idxWin)>0); - lickBurstTime = lickBurstTime(lickBurstTime+max(idxWin)0); - lickBurstTimeIdx = lickBurstTimeIdx(lickBurstTimeIdx+max(idxWin)eigValSigThresh); - eigVecSig = eigVecSort(:,eigValSort>eigValSigThresh); - numEigValSig = length(eigValSig); + eigValSortSig = eigValSort(eigValSort>eigValSigThresh); + eigVecSortSig = eigVecSort(:,eigValSort>eigValSigThresh); + numEigValSig = length(eigValSortSig); for iEig = 1:numEigValSig - PC = eigVecSig(:,iEig)*eigVecSig(:,iEig)'*eigValSig(iEig); + PC = eigVecSortSig(:,iEig)*eigVecSortSig(:,iEig)'*eigValSortSig(iEig); PCsNoDiag(iEig,:,:) = PC - diag(diag(PC)); end - reactPerPC = []; - for it=1:numMatchBins - % outer product of the instantaneous ensemble firing - % in each i,j is the product of firing of cell i and cell j (in this time bin). - ensOutProd=zSpikesMatch(:,it)*zSpikesMatch(:,it)'; - for jEig=1:numEigValSig + out.eigVecSortSig = eigVecSortSig; + out.eigValSortSig = eigValSortSig; + %% compute reactivation trace per sig PC + for jEig=1:numEigValSig + for it=1:numMatchBins + % outer product of the instantaneous ensemble firing + % in each i,j is the product of firing of cell i and cell j (in this time bin). + ensOutProd=zSpikesMatch(:,it)*zSpikesMatch(:,it)'; % dot-product ensembleOuterProd with corr matrix (whose diagonal has been 0'd) - reactPerPC(jEig,it)=sum(sum(ensOutProd.*squeeze(PCsNoDiag(jEig,:,:)))); - end - end - out.reactPerPC = reactPerPC; % per pc full epoch 'match' intervals reactivation strength series - %% per PC shuffle and test psth - - reactSignPCs = []; - for iEig=1:numEigValSig - % comparing reactivation strength to per PC circshift match zbinspikes - reactPerPCShufmean=[]; - reactPerPCShuf=[]; - for qq=1:numShuffs - zSpikesMatchShuf=[]; - for tr=1:numCells - shiftBy=round(rand(1)*numMatchBins); - zSpikesMatchShuf(tr,:)=circshift(zSpikesMatch(tr,:),shiftBy,2); - end - - for i2=1:numMatchBins - ensOutProdShuf=zSpikesMatchShuf(:,i2)*zSpikesMatchShuf(:,i2)'; - reactPerPCShuf(i2)=sum(sum(ensOutProdShuf.*squeeze(PCsNoDiag(iEig,:,:)))); - end - reactPerPCShufmean =[reactPerPCShufmean mean(reactPerPCShuf)]; - end - reactSignPC=mean(mean(reactPerPC(iEig,:))swrReactETAperPCShufsTrim)>0.95 - sigxcorr=1; - else - sigxcorr=0; + %% SWR rxn psth per PC + [out.swrReactPSTHPerPC{jEig}, out.swrReactETAPerPC{jEig}, out.swrReactETAPerPCShufs{jEig}] = ... + stackPSTHwShuf(out.reactPerPC{jEig}, swrTimeIdx, idxWin, numShufs); + %% SWR-Burst rxn psth + if ~isempty(swrBurstTimeIdx) + [out.swrBurstReactPSTHPerPC{jEig}, out.swrBurstReactETAPerPC{jEig}, out.swrBurstReactETAPerPCShufs{jEig}] = ... + stackPSTHwShuf(out.reactPerPC{jEig}, swrBurstTimeIdx, idxWin, numShufs); end - - % need to add burst swr and lick to this psth,shuffle - % swrBurstReactPSTHperPC = - % swrBurstReactETAperPC = - % swrBurstReactPSTHperPCshuff = - % swrBurstReactETAperPCshuff = - % - % lickReactPSTHperPC = - % lickReactETAperPC = - % lickReactPSTHperPCshuff = - % lickReactETAperPCshuff = - + %% lick rxn psth + [out.lickReactPSTHPerPC{jEig}, out.lickReactETAPerPC{jEig}, out.lickReactETAPerPCShufs{jEig}] = ... + stackPSTHwShuf(out.reactPerPC{jEig}, lickBurstTimeIdx, idxWin, numShufs); end end @@ -288,24 +198,61 @@ out.su = []; out.binSpikes = []; out.spikes = []; +out.swrTime = []; +out.swrBurstTime = []; +out.etaTime = []; + out.reactFull = []; % full epoch 'match' intervals reactivation strength series out.reactTime = []; % full epoch 'match' time.. aka reactivation -out.swrTime = []; out.swrReactPSTHfull = []; out.swrReactETAfull = []; -out.reactPSTHtime = []; out.swrReactETAfullShufs = []; -out.swrBurstTime = []; out.swrBurstReactPSTHfull = []; out.swrBurstReactETAfull = []; out.swrBurstReactETAfullShufs = []; +out.burstIntervals = []; out.lickBurstTime = []; out.lickReactPSTHfull = []; out.lickReactETAfull = []; out.lickReactETAfullShufs = []; -end +% per pc +out.eigVecSortSig = []; +out.eigValSortSig = []; +out.reactPerPC = []; + +out.swrReactPSTHPerPC = []; +out.swrReactETAPerPC = []; +out.swrReactETAPerPCShufs = []; +out.swrBurstReactPSTHPerPC = []; +out.swrBurstReactETAPerPC = []; +out.swrBurstReactETAPerPCShufs = []; + +out.lickReactPSTHPerPC = []; +out.lickReactETAPerPC = []; +out.lickReactETAPerPCShufs = []; +end + +function [swrReactPSTHfull, swrReactETAfull, swrReactETAfullShufs] = ... + stackPSTHwShuf(reactFull, swrTimeIdx, idxWin, numShufs) +% stack to psth +swrReactPSTHfull = cell2mat(arrayfun(@(r) reactFull(r+idxWin), swrTimeIdx, 'un', 0)); +swrReactETAfull = nanmean(swrReactPSTHfull); +% time shuffle psth +swrReactETAfullShufs = []; +for sh=1:numShufs + swrReactPSTHfullShuf=[]; + for qq=1:length(swrTimeIdx) + shiftBy=round(rand(1)*size(swrReactPSTHfull,2)); + swrReactPSTHfullShuf(qq,:)=circshift(swrReactPSTHfull(qq,:),shiftBy,2); + end + swrReactETAfullShufs(sh,:) = nanmean(swrReactPSTHfullShuf); +end +swrReactETAfullShufsMean = nanmean(swrReactETAfullShufs); +swrReactETAfullShufsSEM = nanstd(swrReactETAfullShufs)/... + sqrt(size(swrReactETAfullShufs,1)); +end \ No newline at end of file diff --git a/DFFunctions/dfa_lickswrcorr.m b/DFFunctions/dfa_lickswrcorr.m index 14f8d8b0..f6583593 100644 --- a/DFFunctions/dfa_lickswrcorr.m +++ b/DFFunctions/dfa_lickswrcorr.m @@ -34,12 +34,6 @@ end pconf = paramconfig; - -plotfigs = 1; -pausefigs = 0; -savefigs = 1; -savefigas = {'png'}; - bin = .01; tmax = 1; boutNum = 10; @@ -125,6 +119,9 @@ fprintf('lBurst swrs in valid ILI: %d (%.02f pct swrs) \n', sum(swrValidILI), ... sum(swrValidILI)/length(swrStart)*100) swrInBurstStart = swrInBurstStart(swrValidILI); +if isempty(swrInBurstStart) + return +end swrBinILI = swrBinILI(swrValidILI); swrTimeSinceLick = swrTimeSinceLick(swrValidILI); @@ -228,18 +225,18 @@ swrInBurstStart = swrInBurstStart(swrValidILI); swrBinILI = swrBinILI(swrValidILI); swrTimeSinceLick = swrTimeSinceLick(swrValidILI); -% fprintf('swrs with valid ILI: %d (%.02f pct swrs) \n', numel(swrInBurstStart),... -% numel(swrInBurstStart)/length(swrStart)*100) + % fprintf('swrs with valid ILI: %d (%.02f pct swrs) \n', numel(swrInBurstStart),... + % numel(swrInBurstStart)/length(swrStart)*100) % save burst start/end time for each enclosed and valid swr [~,~,swrLickBurstIdx] = histcounts(swrInBurstStart, [burstIntvs(:,1); inf]); swrBurstInterval = burstIntvs(swrLickBurstIdx,:); - + out.swrBinILIShuf{end+1} = swrBinILI; out.swrInBurstStartShuf{end+1} = swrInBurstStart; out.swrTimeSinceLickShuf{end+1} = swrTimeSinceLick; out.swrBurstIntervalShuf{end+1} = swrBurstInterval; - %% Compute shuffle lickDin x swr + %% Compute shuffle lickDin x swr % xcorr norm smooth xc = spikexcorr(sort(swrInBurstStart), lickboutlicks, bin, tmax); normxc = xc.c1vsc2 ./ sqrt(xc.nspikes1 * xc.nspikes2); % normalize xc @@ -254,7 +251,7 @@ p1 = xc.nspikes1/T; p2 = xc.nspikes2/T; % fr in Hz expProb = p1*p2; % per sec. Expected probability - + out.xcShuf{end+1} = xc; out.normxcShuf{end+1} = normxc; out.smthxcShuf{end+1} = smthxc; @@ -286,88 +283,8 @@ out.phasemodShuf{end+1} = phasemod; end end - -%% plot -if plotfigs - [~, fname,~] = fileparts(mfilename('fullpath')); - outdir = sprintf('%s/%s/', pconf.andef{4},fname,animal); - Pp=load_plotting_params({'defaults','dfa_lickswrcorr'}, 'pausefigs', pausefigs, ... - 'savefigs', savefigs); - - %% Xcorr norm smooth, shuff mean/std - sf1 = subaxis(2,2,1,Pp.posparams{:}); - sf1.Tag = 'xcorr'; - - % shuffled xcorr with std ghost trail - xmsh = mean(cell2mat(out.smthxcShuf')); - xstdsh = std(cell2mat(out.smthxcShuf')); %/size(cell2mat(out.smthxcShuf'),1); - plot(out.xc.time, xmsh, 'color', [0 0 1 .2], 'linewidth', 1); - hold on; - fill([out.xc.time'; flipud(out.time')],[xmsh'-xstdsh';flipud(xmsh'+xstdsh')],'b', 'linestyle', ... - 'none', 'facealpha', .2); - % xcorr norm - bar(out.xc.time, out.normxc, 'k', 'facealpha', .2, 'edgealpha', 0) - % xcorr norm smooth - plot(out.time, out.smthxc, 'k') - line([0 0], ylim, 'color', 'k', 'linestyle', '--', 'linewidth', .5) - - ylabel('xcorr'); - xlabel('time from lick s'); - % title('xcorr - hold off; - %% excorr over shuff excorr cdf distr - % relative swr from last lick - sf2 = subaxis(2,2,2); - histogram(cell2mat(out.excorrShuf), 60,'Normalization','probability','edgealpha', 0, 'facecolor', 'k'); - excsort = sort(cell2mat(out.excorrShuf)); - idxsig = round(sigpct*length(out.excorrShuf)); - line([excsort(idxsig) excsort(idxsig)], ylim, 'color', [0 0 0 .8], 'linestyle', '--'); - hold on; - line([out.excorr out.excorr], ylim, 'color', 'r'); - excp = 1-sum(out.excorr>cell2mat(out.excorrShuf))/length(out.excorrShuf); - title(sprintf('excorr %.03f p%.03f', out.excorr, excp)); - ylabel('probability') - xlabel('excess corr') - axis tight - hold off; - %% polar distr phase clustering, swrLickPhase, meanMRVmag - sf3 = subaxis(2,2,3); - % a = polarhistogram(swrLickPhase, 16, 'Normalization', 'pdf', 'edgealpha', 0,... - % 'facealpha', .5); - polarplot([zeros(size(out.swrLickPhase,1),1) out.swrLickPhase]', ... - repmat([0 1],size(out.swrLickPhase,1),1)', 'color', [0 0 0 .4], 'linewidth', 4); - hold on - polarplot([0; out.vecang], [0; out.meanMRVmag], 'color', [1 0 .3], 'linewidth', 4) - grid off - rticks([]) - thetaticks([]) - title('swr ILI-phase') - hold off - axis tight - %% phase mod - sf4 = subaxis(2,2,4); - histogram(cell2mat(out.phasemodShuf), 100, 'Normalization', 'probability', 'edgealpha', 0, 'facecolor', 'k'); - hold on; - mrvsort = sort(cell2mat(out.phasemodShuf)); - idxsig = round(sigpct*length(out.phasemodShuf)); - line([mrvsort(idxsig) mrvsort(idxsig)], ylim, 'color', [0 0 0 .8], 'linestyle', '--'); - hold on - line([out.phasemod out.phasemod], ylim, 'color', 'r'); - modp = 1-sum(out.phasemod>cell2mat(out.phasemodShuf))/length(out.phasemodShuf); - title(sprintf('logMRVmag %.03f p%.03f Rpval%.03f', out.phasemod, modp, pval)); - ylabel('probability') - xlabel('log(Rayleigh Z)') - axis tight - hold off - - %% ---- super axis ----- - sprtit = sprintf('%s %d %d %s', animal, day, epoch, fname(5:end)); - setSuperAxTitle(sprtit) - % ---- pause, save figs ---- - if pausefigs; pause; end - if savefigs; save_figure(outdir, sprtit, 'savefigas', savefigas); end -end end + function out = make_blank() out.index = []; out.time = []; @@ -394,7 +311,8 @@ out.swrPctSinceLick = []; out.vecang = []; out.phasemod = []; -%% shuffle + +%% shuffle out.swrStartShuf = []; % xcTime out.xcShuf = {}; diff --git a/DFFunctions/load_filter_params.m b/DFFunctions/load_filter_params.m index a2c24a2c..321fb147 100644 --- a/DFFunctions/load_filter_params.m +++ b/DFFunctions/load_filter_params.m @@ -159,34 +159,27 @@ % Fp.minstdthresh,'maxvelocity',Fp.maxvelocity,'minvelocity', ... % Fp.minvelocity, 'consensus_numtets',Fp.consensus_numtets,'welldist', ... % Fp.welldist); + %% filter function specific params + case 'dfa_reactivationPSTH' bin = 0.025; % seconds win = [-2 2]; + perPC = 1; - consensus_numtets = 2; % minimum # of tets for consensus event detection - minstdthresh = 3; % STD. how big your ripples are - exclusion_dur = 0; % seconds within which consecutive events are eliminated / ignored - minvelocity = 0; - maxvelocity = 4; 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 + minLickTimeFromSwr = 1; % time duration from closest lick iterator = 'singleepochanal'; datatypes = {'spikes', 'lick', 'ca1rippleskons', 'cellinfo'}; - options = {'cellfilter', cellfilter, 'win', win, 'bin', bin, ... - 'swrTimeFilter', ... - [{{'getconstimes', '($cons == 1)', ... - 'ca1rippleskons', 1,'consensus_numtets',consensus_numtets, ... - 'minstdthresh', minstdthresh,'exclusion_dur',exclusion_dur, ... - 'minvelocity', minvelocity,'maxvelocity',maxvelocity}}, ... - {{'getpriortofirstwell', '($prefirst == 0)'}}], ... + options = {'cellfilter', cellfilter, 'win', win, 'bin', bin, 'perPC', perPC,... + 'templateFilter', ... + {{'get2dstate','($velocity>4)'}}, ... 'burstTimeFilter', ... {{'getLickBout', '($lickBout == 1)', 'maxIntraBurstILI', maxIntraBurstILI, ... 'minBoutLicks', minBoutLicks}}}; - case 'swrlickmod' filtfunction = 'swrlickmod'; diff --git a/DFFunctions/load_plotting_params.m b/DFFunctions/load_plotting_params.m index 8011c32c..88fa493f 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 'Reactivation PerPC' + position = [.1 .1 .5 .5]; + winSE = [-2 2]; + case 'areaTFspect' % area = {{'mec', 'sup'}, {'mec', 'deep'}, {'ca1', 'd'}}; position = [.1 .1 .8 .5]; win = [-.5 .5]; usecolormap = hot; tickFsize = 8; - case 'ReactivationStrength' + case 'Reactivation Full' position = [.1 .1 .5 .5]; winSE = [-2 2]; case 'wetVDryILIphaseSWR' diff --git a/Functions/vals2rgb.m b/Functions/vals2rgb.m new file mode 100644 index 00000000..d32c935a --- /dev/null +++ b/Functions/vals2rgb.m @@ -0,0 +1,63 @@ +function rgb = vals2rgb(vals, colormap, crange) +% Take in a vector of N values and return and return a Nx3 matrix of RGB +% values associated with a given colormap +% +% rgb = AFQ_vals2colormap(vals, [colormap = 'jet'], [crange]) +% +% Inputs: +% vals = A vector of values to map to a colormap or a cell array of +% vectors of values +% colormap = A matlab colormap. Examples: colormap = 'autumn'; +% colormap = 'jet'; colormap = 'hot'; +% crange = The values to map to the minimum and maximum of the colormap. +% Defualts to the full range of values in vals. +% +% Outputs: +% rgb = Nx3 matrix of rgb values mapping each value in vals to the +% corresponding rgb colors. If vals is a cell array then rgb +% will be a cell array of the same length +% +% Example: +% vals = rand(1,100); +% rgb = AFQ_vals2colormap(vals, 'hot'); +% +% Copyright Jason D. Yeatman, June 2012 + +if ~exist('colormap','var') || isempty(colormap) + colormap = 'jet'; +end + +% +if ~iscell(vals) + if ~exist('crange','var') || isempty(crange) + crange = [min(vals) max(vals)]; + end + % Generate the colormap + cmap = eval([colormap '(256)']); + % Normalize the values to be between 1 and 256 + vals(vals < crange(1)) = crange(1); + vals(vals > crange(2)) = crange(2); + valsN = round(((vals - crange(1)) ./ diff(crange)) .* 255)+1; + % Convert any nans to ones + valsN(isnan(valsN)) = 1; + % Convert the normalized values to the RGB values of the colormap + rgb = cmap(valsN, :); +elseif iscell(vals) + if ~exist('crange','var') || isempty(crange) + crange = [min(vertcat(vals{:})) max(vertcat(vals{:}))]; + end + % Generate the colormap + cmap = eval([colormap '(256)']); + for ii = 1:length(vals) + % Normalize the values to be between 1 and 256 for cell ii + valsN = vals{ii}; + valsN(valsN < crange(1)) = crange(1); + valsN(valsN > crange(2)) = crange(2); + valsN = round(((valsN - crange(1)) ./ diff(crange)) .* 255)+1; + % Convert any nans to ones + valsN(isnan(valsN)) = 1; + % Convert the normalized values to the RGB values of the colormap + rgb{ii} = cmap(valsN, :); + end +end +return \ No newline at end of file diff --git a/clusterless_decode/Bon__4_6_fullep.mat b/clusterless_decode/Bon__4_6_fullep.mat new file mode 100644 index 00000000..44e24d3e Binary files /dev/null and b/clusterless_decode/Bon__4_6_fullep.mat differ diff --git a/clusterless_decode/analyze_clusterless_Pos.m b/clusterless_decode/analyze_clusterless_Pos.m new file mode 100644 index 00000000..1a9276b9 --- /dev/null +++ b/clusterless_decode/analyze_clusterless_Pos.m @@ -0,0 +1,1107 @@ + +% Find continuous Theta-associated alternations in Pos decodes +allan = {'Government','Egypt','Chapati','Dave','Higgs','Frank','Bond','Corriander',... + 'Eight','Ten','Conley','Miles'}; +mostan = {'Government','Egypt','Chapati','Dave','Frank','Bond','Corriander',... + 'Eight','Ten','Conley','Miles'}; +core4 = {'Bond','Frank','Dave','Government'}; +anim1st = {'Government','Egypt','Chapati','Dave','Higgs','Frank'}; +anim2nd = {'Bond','Corriander','Eight','Ten','Conley','Miles'}; + +animals_order = {'Government','Egypt','Chapati','Dave','Higgs','Frank','Bond','Corriander',... + 'Iceland','Justice','Kapital','Laplace','Eight','Ten','Conley','Miles'}; + +calculate = 0; +plot_pvalues = 1; +plot_behavior = 1; +report_longest_altpers = 1; + P_thresh_long_altpers = 0.001; +report_prevalence = 0; + + DECODE_RUN = 1; % 0: Uniform, 1: Random walk, 100: (something old?) + EXCURTYPE_toplot = [ 1 ] ; % 1: Pro, 2: Ret, 3: In L, 4: In R + +if plot_pvalues || plot_behavior || report_longest_altpers + animals_toplot = mostan; %mostan; %{'Government'}; %allan; %{'Bond','Frank','Dave','Government'}; %{'Bond','Frank','Dave','Government'}; %'Bond','Frank','Dave', +end + +if calculate + + animals_tocalc = allan; %{'Egypt','Corriander'}; % 'Bond','Frank','Dave','Government','Dave','Corriander','Miles','Eight'}; + manual_dayep = []; + + % Analysis targets % + PRECHOICE_ONLY = 1; % 0: full excur, 1: CP only + EXCURTYPE_tocalc = [1]; % 1: Pro, 2: Ret, 3: In L, 4: In R + PHASECHOICE = 1; % 0: min firing phase + % 1: min LR density ** Choose this + % 2: max LR density + MOVINGPERIODS = 2; % Periods to analyze: 1: vel4, 2: nonimmobile05 (half sec) + LR_prop_thresh = 0.1; % Minimum LR density (vs. C) to calculate LR value + + % Alt period detection params % + ALT_THRESHOLD = 0.1; + NUMSAMP = 10000 + 1; % # of resamples + SAMPMETHOD = 1; % 1: permutation, 2: bootstrap + MAXALTS = 5; + P_THRESH = 0.05; + + % Theta % + nbins = 12; %18; % 36 bins is 10 deg bins, after Jezek 2011 + phasebins = -pi:(2*pi/nbins):pi; + phasebins_c = phasebins(1:(end-1)); +end + +% Saved file directories % +if DECODE_RUN == 1 + decodedir = '/opt/data50/kkay/__Decode/LR_decode_20_4_Randwalk_new'; + thetabindir = '/opt/data50/kkay/__Decode/LR_decode_20_4_Randwalk_new'; + cd(decodedir) +elseif DECODE_RUN == 100 + decodedir = '/opt/data50/kkay/__Decode/Dir_decode_20_4_Uniform'; LR_decode_20_4_Randwalk_new + thetabindir = '/opt/data50/kkay/__Decode/Dir_decode_20_4_Uniform/Thetabins_LR'; + cd(decodedir) +end + +% Internal study figures %%%%%%%%%%%% +plot_Xcorr_LR = 0; +plot_LRprop_hist = 0; +plot_LRprop_phasehist = 0; +plot_altspeed_hist = 0; + +if plot_behavior || plot_pvalues + BEHAVE_P_THRESH = 0.05; + LENGTH_TOPLOT = 3; % [2 3] +end +if plot_pvalues + pbins = 0:.01:.35; +end + + +% Data loading %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +if calculate + + cd(decodedir) + + animals_order = {'Government','Egypt','Chapati','Dave','Higgs','Frank','Bond','Corriander',... + 'Iceland','Justice','Kapital','Laplace','Eight','Ten','Conley','Miles'}; + + for XT = EXCURTYPE_tocalc + + EXCURTYPE = XT; + + out = []; + + for aa = 1:length(animals_tocalc) + + animalname = animals_tocalc{aa}; + animalpref = animalname(1:3); + animalinfo = animaldef(animalname); + daydir = getdaydir(animalname); + an = find(strcmp(animalname,animals_order)); + task = loaddatastruct(animalinfo{2},animalinfo{3},'task'); + + if ~isempty(manual_dayep) + dayeps = manual_dayep; + else + days = 1:length(task); + dayeps = []; + for dd = days + if ~isempty(task{dd}) + eps = wtrackeps(task,dd); + dayeps = [dayeps ; dd * ones(length(eps),1) eps(:)]; + end + end + end + + % Initialize outputs + out.date = date; + out.animalname = animalinfo{3}; + out.dayeps = dayeps; + out.numperms = NUMSAMP; + out.divider1 = '%%%%%%%%%%%%%%%%%%%%%%%%%%%%'; + out.altdist_ep = {}; %nan(NUMSAMP,length(ALTDISTBINS)-1); + out.numaltpers_ep = {}; %nan(NUMSAMP,MAXALTS); + out.altmeans_ep = {}; %nan(1,NUMSAMP); + out.altmedians_ep = {}; %nan(1,NUMSAMP); + out.divider2 = '%%%%%%%%%%%%%%%%%%%%%%%%%%%%'; + out.altperiods = {}; + out.descript = '1: day 2: ep 3: time_a, 4: time_b, 5: tb_a, 6: tb_b, 7: numcyc, 8: xcurtype, 9-11: speed, accel, angspeed, 12: pval'; + + for de = 1:size(dayeps,1) + + d = dayeps(de,1); + ep = dayeps(de,2); + + disp(sprintf('%s: d %d ep %d',animalinfo{3}(1:3),d,ep)) + + % Load Data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + if ~exist('P','var') || isempty(P) || ~strcmp(P.animalname,animalname) || ... + ~all(P.dayep == [d ep]) || EXCURTYPE ~= P.EXTYPE + + % Decode file %%% + filename = sprintf('%s_%d_%d_Pos_%d.mat',animalpref,d,ep,EXCURTYPE); + cd(decodedir); + filedir = dir(filename); + if isempty(filedir) + disp(sprintf('not finding %s, skipping',filename)) + continue + end + load(filedir.name,'P'); + + % Positional % + pos = loaddatastruct(animalinfo{2},animalinfo{3},'pos',P.dayep(1)); + linpos = loaddatastruct(animalinfo{2},animalinfo{3},'linpos',P.dayep(1)); + lindist = linpos{d}{ep}.statematrix.lindist; + CP_buff = choicepointer(linpos{d}{ep}); + CP_bin = ceil(CP_buff); + segind = linpos{d}{ep}.statematrix.segmentIndex; + seg2 = (segind == 2) ; seg3 = (segind == 3) ; + seg4 = (segind == 4) ; seg5 = (segind == 5) ; + Carm = (lindist <= CP_buff) ; % C arm pos times + Larm = ~Carm & (seg2 | seg3) ; % L arm pos times + Rarm = ~Carm & (seg4 | seg5) ; % R arm pos times + postimevec = pos{d}{ep}.data(:,1); + veltrace = pos{d}{ep}.data(:,5); + epstart = postimevec(1); + epend = postimevec(end); + headang = loaddatastruct(animalinfo{2}, animalinfo{3}, 'headang', d); + angveltrace = abs(headang{d}{ep}.angvel) * 180 / pi; % convert to degrees + %[posvec_all] = posvecmaker(lindist,linpos{d}{ep}); + %[posvec_all,outersep] = posvecmaker_compress(posvec_all); + yvec = P.yvec_pos; + outersep = P.CLR_lengths(1) + P.CLR_lengths(2); + + % % Head direction data % + % headdir = linpos{d}{ep}.statematrix.segmentHeadDirection(:,1); % head direction relative to center well -- values > 0 are outbound + % postimevec_nonnan = postimevec(~isnan(headdir)); + % headdir_nonnan = headdir(~isnan(headdir)); + % headdir2 = interp1(postimevec_nonnan,headdir_nonnan,postimevec,'linear'); + % outvec = (headdir2 >= 0); + % invec = (headdir2 < 0); + % disp(num2str(sum(invec & Larm)/30)) + % keyboard + + % By Excursion type, identify Choice Times + C,L,R pos bins %%%% + Choicevec = []; + if EXCURTYPE == 1 % Pro + exflags = [1 3 -11]; + Choicevec = Carm; % Pro Choice Times : Center arm + Cin = yvec <= CP_bin; % Center arm + Lin = (yvec > CP_bin) & (yvec <= outersep); % Left arm + Rin = ~Cin & ~Lin; % Right arm + elseif EXCURTYPE == 2 % Ret + exflags = [2 4 -11]; + Choicevec = Carm; % Ret "Choice" Times : Center arm + Cin = yvec <= CP_bin; % Center arm + Lin = (yvec > CP_bin) & (yvec <= outersep); % Left arm + Rin = ~Cin & ~Lin; % Right arm + elseif EXCURTYPE == 3 % In L + exflags = [-33 4 32]; + Choicevec = Larm; % L inbound "Choice" Times : Left arm + Cin = (yvec > CP_bin) & (yvec <= outersep) ; % Left arm + Lin = yvec > outersep; % Right arm (here the "left side") + Rin = ~Cin & ~Lin; % Right arm + elseif EXCURTYPE == 4 % In R + exflags = [-22 2 23]; + Choicevec = Rarm; % R inbound "Choice" Times: Right arm + Cin = yvec > outersep; % Center arm + Lin = yvec <= CP_bin; % Left arm + Rin = ~Cin & ~Lin; % Right arm + end + + % Identify Excursions (Left, Right, Trackback) %%% + inds = ismember(P.excurlist(:,3),exflags); + excurlist = P.excurlist(inds,1:3); + numexcur = size(excurlist,1); + + % Moving periods % + timefilterscript + if MOVINGPERIODS == 1 + [dummy,~] = evaluatetimefilter(animalinfo{2},animalinfo{3},{ vel4 },[d ep]); + elseif MOVINGPERIODS == 2 + [dummy,~] = evaluatetimefilter(animalinfo{2},animalinfo{3},{ nonimmobile05 },[d ep]); + elseif MOVINGPERIODS == 3 + [dummy,~] = evaluatetimefilter(animalinfo{2},animalinfo{3},{ nonimmobile1 },[d ep]); + end + movingperiods = dummy{d}{ep}; % + movingvec = logical(list2vec(movingperiods,postimevec)); + + % Tets + clust spikes + spikes = loaddatastruct(animalinfo{2},animalinfo{3},'spikes',d); + cellinfo = loaddatastruct(animalinfo{2},animalinfo{3},'cellinfo'); + [adtc_list] = clusteredunits(spikes,cellinfo,an,d,ep,[],1); + maxcell = size(adtc_list,1); + selected_tets = unique(P.selected_tets); + numtets = length(selected_tets); + + % Theta % + thetabins = loaddatastruct(animalinfo{2},animalinfo{3},'thetabins',d); + thetatet = thetabins{d}.thetatet; + if strcmp(animalpref,'Dav') && d == 3 + disp('Dav thetatet = 1, d 3 ep 2'); + thetatet = 1; + %elseif strcmp(animpref,'Bond') && d == 6 + % disp('Bon thetatet = 30, d 6 ep 4'); + % thetatet = 30; + end + theta = loadeegstruct(animalinfo{2},animalinfo{3},'theta',d,ep,thetatet); + thetatrace = theta{d}{ep}{thetatet}.data(:,1); + tph = thetabins{d}.tph{ep}; + tvec_tph = thetabins{d}.timevec_tph{ep}; + %thetaphase = double(theta{d}{ep}{thetatet}.data(:,2))/10000; % hilbert phase + %thetatimevec = geteegtimes(theta{d}{ep}{thetatet}); + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + % I. Construct Tvec and LRdensity, i.e. decode concatenated across excursions + % (Also, histogram LR density to get theta phase at which LR density is minimal (minphase)) %%%%%%%%%%%%%%%%%%%%%% + + Tvec = []; + CLRdec = []; % [row: times // col: L: 1, R: 2 ] + Phasehist_LR = zeros(size(phasebins_c)); + Analysisperiods = [] ; % all analysis pers in epoch + Analysis_totaldur = 0; + + % Get total analysis vector + periods for this epoch + Analysisvec = false(size(postimevec)); + Tvec_collect = []; + for x = 1:numexcur + % Define analysis periods % (Pre-choice point, Excursion times, Movement) + excurvec = list2vec(excurlist(x,[1 2]),postimevec); + + if PRECHOICE_ONLY + intersectvec = Choicevec & excurvec & movingvec; % *** Prior to Choice, This Excursion, Movement **** + else + intersectvec = excurvec & movingvec; % *** This Excursion, Movement **** + end + + Analysisvec(intersectvec) = 1; + % Collect (overlapping) decoding times + Tvec_collect = [Tvec_collect ; P.Tvecs{x}(:)]; + end + Tvec = unique(Tvec_collect); + CLRdec = nan(length(Tvec),3); % [ times x CLR ] + Analysisperiods = vec2list( Analysisvec,postimevec ); + Analysis_totaldur = Analysis_totaldur + round( sum( Analysisperiods(:,2) - Analysisperiods(:,1) )); + + % Iterate excursions to collect + for xx = 1:numexcur + + % Obtain Posterior (C,L,R components) + post_exc = P.Posts{xx}; % Posterior during excursion + tvec_exc = P.Tvecs{xx}; % Time vector of bins during excursion + Cdensity = sum(post_exc(Cin,:),1); + Ldensity = sum(post_exc(Lin,:),1); + Rdensity = sum(post_exc(Rin,:),1); + LRdensity = Ldensity + Rdensity; + LRthreshinds = LRdensity >= LR_prop_thresh; + + tphs = tph(lookup(tvec_exc,tvec_tph - epstart)); % Theta phase of each cut density bin + + % Finally, install onto epoch-master Decode variable + inds_ep = lookup(tvec_exc,Tvec); + CLRdec(inds_ep,:) = [ Cdensity(:) Ldensity(:) Rdensity(:) ]; + + % Also, add LR density to accumulator phase histogram + [~,I] = histc(tphs,phasebins); + for bbb = 1:nbins + %LRdensity_addition = sum(LRdensity( (I' == bbb) & LRthreshinds_cut(:) )); % Correct phase bin & meets threshold + LRdensity_addition = sum(LRdensity( (I' == bbb) )); % Correct phase bin & meets threshold + Phasehist_LR(bbb) = Phasehist_LR(bbb) + LRdensity_addition; % Add to total + end + + % (plot check) Hist of LR densities during this excursion's analysis periods + if plot_LRprop_hist + H = figure('units','normalized','outerposition',[.2 .06 .2 .2]); + hist(LRdensity,20); axis tight + title('LR density','fontweight','bold','fontsize',12) + keyboard + end + + end + + disp(sprintf('%d s of epoch analyzed (EXCURTYPE: %d)',Analysis_totaldur,EXCURTYPE)) + + if Analysis_totaldur == 0 + disp('hmm, skipping') + continue + end + + %%% II. Divide epoch into Theta bins using LR minphase %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % (a) Select phase + divide theta cycles %%%%%%%%%%%%%%%%%%%%%%% + + if PHASECHOICE == 0 + tbins = thetabins{d}.thetabins{ep}; % min firing phase (this is not preferable, since LFP variable) + tbblocks = thetabins{d}.tbblocks{ep}; + minphase = thetabins{d}.minphase; + elseif PHASECHOICE == 1 + [~,I2] = min(Phasehist_LR); % min + minphase = phasebins_c(I2); + [tbins, ~, ~, ~, ~, ~, tbblocks] ... + = thetabinner(tph,tvec_tph,minphase,[epstart epend]); + elseif PHASECHOICE == 2 + [~,I2] = max(Phasehist_LR); % max + minphase = phasebins_c(I2); + [tbins, ~, ~, ~, ~, ~, tbblocks] ... + = thetabinner(tph,tvec_tph,minphase,[epstart epend]); + end + + % (Optional plot check) Theta phase histogram of LR density (ep-wide) %%%%% + if 0 %plot_LRprop_phasehist + H = figure('units','normalized','outerposition',[.1 .3 .2 .2]); + B = bar(phasebins_c,Phasehist_LR,'histc'); hold on + set(B,'facecolor','k') + set(gca,'fontsize',14) + axis tight + plot([minphase minphase],[0 max(Phasehist_LR)],'r-','linewidth',3); hold on + title(sprintf('min phase : %0.3f',minphase)) + keyboard + end + + % (b) Store this in a file (for plotting later) + clear tbins_LR; + tbins_LR.date = date; + tbins_LR.PHASECHOICE = PHASECHOICE; + tbins_LR.minphase = minphase; + tbins_LR.tbins = tbins; + tbins_LR.tbblocks = tbblocks; + cd(thetabindir); + filename_tb = sprintf('Thetabins_LR_%d_%s_%d_%d.mat',EXCURTYPE,animalname(1:3),d,ep); + save(filename_tb,'tbins_LR','-v7.3'); + cd(decodedir) + + % (c) Filter for Theta bins that occur within Analysis periods %%%%%%%%%%%%%%% + filterinds = logical( isExcluded( mean(tbins,2) , Analysisperiods ) ); % Get theta bins that occur analysis periods + Tbins = tbins(filterinds,:); + Tbins_mean = mean(tbins,2); + Numcyc = size(Tbins,1); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + %%% III. Calculate LR Decode + Alt speed + Indicator vectors / Theta bin %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + Tbins; % [ starttime endtime] of theta cycle + Tbins_mean = mean(Tbins,2); % mid-time of theta cycle + Contig_tb = []; % 111122221111: Contig indicator vector for theta bins + LR_tb = nan(Numcyc,1); % R/(R+L) proportion / theta cycle + Alt_tb = nan(Numcyc,1); % LR alternation speed / theta cycle + Excur_tb = []; % Actual excursion-type (L, R, Trackback) the rat was on during the Theta bin + + % Calculate LR index (0: Left, 1: Right) + % Calculate Alt speed + % Construct Indicator vectors + + % Iterate each theta cycle + for tb = 1:Numcyc + + bin_a = lookup(Tbins(tb,1) - epstart,Tvec); % Start index of bin + bin_b = lookup(Tbins(tb,2) - epstart,Tvec); % End index of bin + midtime = mean(Tbins(tb,:)); + %disp(num2str(midtime-epstart)) + + % (Provisional - fix later) Identify the excursion-type of this bin + excur = nan; + for xxxx = size(excurlist,1):-1:1 + if sum(isExcluded(midtime,excurlist(xxxx,[1 2]))) > 0 + excur = [excur ; excurlist(xxxx,3)]; + break + end + end + if isnan(excur) + excur = nan; + disp('nan excur -- should be rare') + %elseif length(unique(excur)) > 1 + % keyboard + %else + % excur = unique(excur); + end + + Excur_tb = [Excur_tb ; excur]; + + % % Check discrepancy due to decoding bin size (e.g. 20 ms) + % decodedur = Tvec(bin_b) - Tvec(bin_a); + % actualdur = Tbins(tb,2) - Tbins(tb,1); + % diffdur = abs(decodedur - actualdur) ; + % disp(num2str(diffdur)); + % if diffdur > 0.008 + % keyboard + % end + + % Post density of C, L, R arms + Cval = sum(CLRdec(bin_a:bin_b,1)); % Center arm density of this theta bin + Lval = sum(CLRdec(bin_a:bin_b,2)); % Left " " + Rval = sum(CLRdec(bin_a:bin_b,3)); % Right " " + + if any(isnan([Cval Lval Rval])) + %keyboard + disp('*********nan CLR val (should only happen occasionally) ****************') + disp(num2str([Cval Lval Rval])) + end + + % Calc LR proportion + LR_prop = (Lval + Rval)/(Cval + Rval + Lval); + + % Calc LR (if exceed min LR density) + if LR_prop > LR_prop_thresh + + LRFLAG = 2; + + if LRFLAG == 1 + % Option 1: Proportion + LRval = Rval / (Lval + Rval); + elseif LRFLAG == 2 + % Option 2: Sign + LRval = Rval / (Lval + Rval); + LRval = sign(LRval - 0.5); % +1: Right, -1: Left + end + + LR_tb(tb) = LRval; + else + LR_tb(tb) = nan; + end + + if tb > 1 + % Check to make sure cycle did not come after discontinuous theta bin gap + if Tbins(tb,1) == Tbins(tb-1,2) % Contiguous + + % Calc Alt speed % + if ~isnan(LR_tb(tb-1)) + Alt_tb(tb) = abs( LR_tb(tb) - LR_tb(tb-1) ); + end + + % Mark Contiguity % + if Contig_tb(end) == 1 + Contig_tb = [Contig_tb ; 1]; % keep + else + Contig_tb = [Contig_tb ; 2]; % keep + end + + else % Discontiguous + + % Mark Discontiguity % + if Contig_tb(end) == 1 + Contig_tb = [Contig_tb ; 2]; % switch + else + Contig_tb = [Contig_tb ; 1]; % switch + end + end + elseif tb == 1 + Contig_tb = [Contig_tb ; 1]; % Initialize w/ 1 + continue + end + + end + + if 0 %(optional) % check alternations + if LRFLAG == 1 + figure; + plot(Tbins_mean-epstart,Alt_tb,'k.'); hold on + plot(Tbins_mean-epstart,Alt_tb,'k-'); hold on + ylim([-0.1 2.1]) + elseif LRFLAG == 2 + plot(Tbins_mean-epstart,Alt_tb,'r.'); hold on + plot(Tbins_mean-epstart,Alt_tb,'r-'); hold on + ylim([-0.1 2.1]) + end + %break + end + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + % IV. Detect Alternation Periods %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % Initialize output + Altpers = []; + Altpers_xc = []; % excursion type of the alt period + Altpers_dur = []; + numaltperiods = nan; % total # of alternation periods this ep + + % Initialize while-loop state variables + cy = 2; % pointer index + altcycdur = 0; % alt per dur + per_a = 1; % alt per start + per_b = nan; % alt per end + + while cy ~= Numcyc + + if Contig_tb(cy) == Contig_tb(cy-1) % Contiguous + + if Alt_tb(cy) > ALT_THRESHOLD % this theta bin qualifies as alternation + altcycdur = altcycdur + 1; % alternation duration (in cycles) + per_b = cy; % update pointer + else % not reach alternation thresh + % Store alternation period + if altcycdur > 2 + Altpers = [Altpers ; per_a per_b ]; + end + % Reset pointers + altcycdur = 0; + per_a = cy; + per_b = nan; + end + + else % Discontiguous + + %Reset pointers + altcycdur = 0; + per_a = cy; + per_b = nan; + + end + + cy = cy + 1; + + end + + numaltperiods = size(Altpers,1); + + if numaltperiods == 0 + disp('no alt periods detected') + continue + end + + % Identify excursion types for alt periods + Altpers_xc = nan(numaltperiods,1); + for vv = 1:numaltperiods + extype = unique( Excur_tb( Altpers(vv,1) : Altpers(vv,2) ) ); + if length(extype) ~= 1 + disp('ALTPERIOD STRADDLES TWO EXCUR: ASSUMING TRACKBACK') + [~,III] = max(abs(extype)); + extype = extype(III); + %keyboard + end + Altpers_xc(vv) = extype; + end + + if isempty(Altpers) + continue + end + + % Tabulate cycle durations + Altpers_dur = [ Altpers(:,2) - Altpers(:,1) ] + 1; % Duration, in theta cycles, of each candidate period + for mm = 1:MAXALTS + out.numaltpers_ep{de}(1,mm) = sum(Altpers_dur >= mm); % Periods with at least mm cycles + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + % V. Permute single alt periods to get P-values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % Initialize output + % [ time_a time_b tb_a tb_b numcyc xcurtype speed accel angspeed pvalue ] + out.altperiods{de} = nan(numaltperiods,10); + + % (i) Preliminary: obtain LRvals for each excursion type (within this epoch) + cycinds_xc = {}; % indices of theta cycles of this excursion type + numcyc_xc = {}; % # of theta cycles of this excursion type + LRvals_xc = {}; % LR values for each of these theta cyccles in this excursion type + for XF_flag = exflags % iterate through each excursion flag value (only positive!) + xf_inds = excurlist(:,3) == XF_flag; + XF = abs(XF_flag); + cycinds_xc{XF} = find( isExcluded(Tbins_mean,excurlist(xf_inds,[1 2])) ); % Right + numcyc_xc{XF} = length(cycinds_xc{XF}); + LRvals_xc{XF} = LR_tb(cycinds_xc{XF}); + end + + % Iterate each alt period + for ll = 1:numaltperiods + + % Basic info % + tb_a = Altpers(ll,1); % theta bin # of first cycle + tb_b = Altpers(ll,2); % theta bin # of last cycle + time_a = Tbins(tb_a,1); % absolute time of start of first cycle + time_b = Tbins(tb_b,2); % absolute time of end of last cycle + numcyc = tb_b - tb_a + 1; % number of theta cycles + xf = abs(Altpers_xc(ll)); % excursion flag (abs makes it indexable) + + if isnan(xf) + disp('xcursion type cant be found') + keyboard + %continue + end + + disp(sprintf('%s: d%d ep%d xc%d altper #%d (%d)',animalinfo{3}(1:3),d,ep,xf,ll,numcyc)) + + % Behavior (calc mean of speed, accel, angular speed) % + pos_a = lookup(time_a,postimevec); + pos_b = lookup(time_b,postimevec); + speed = mean(veltrace(pos_a:pos_b)); % cm/s + accel = mean(diff(veltrace(pos_a:pos_b))); % cm/s^2 + angspeed = mean(angveltrace(pos_a:pos_b)); % deg / s + lowspeedflag = any(veltrace(pos_a:pos_b) < 4); % overlaps with low-speed time + + % Permutation procedure (i, ii) %%%%% + + % (i) Identify indices of the surrounding theta bin contiguous period (or "block period") + % in which this alternation period occurred + % strategy is to use a while loop and step backwards and + % forwards + tbc_a = []; % first bin # of surrounding contig theta period (with decoded LR) + tbc_b = []; % last bin # of surrounding contig theta period (with decoded LR) + + % Initialize % + tbc_a = tb_a; % first bin # of surrounding contig theta period (with decoded LR) + tbc_b = tb_b; % last bin # of surrounding contig theta period (with decoded LR) + blockval = tbblocks(tb_a); % dummy val is 1 or 2: designates a contiguous theta period + + % Find beginning by stepping backwards % + while tbc_a > 0 + if tbblocks(tbc_a) == blockval && ~isnan(LR_tb(tbc_a)) % Same contig theta block && LR is decoded + tbc_a = tbc_a - 1; + else + break % terminate + end + end + + % Find end by stepping forward % + while tbc_b <= Numcyc + if tbblocks(tbc_b) == blockval && ~isnan(LR_tb(tbc_b)) % Same contig theta block && LR is decoded + tbc_b = tbc_b + 1; + else + break % terminate + end + end + + blockdur = tbc_b - tbc_a + 1; % Duration (in cycles) of surrounding block period + + % (ii) Permute + altdetect = nan(1,NUMSAMP); + for rrr = 1:NUMSAMP + + % Resample % + % (here we just go ahead and resample all the LR vals in the epoch) + if SAMPMETHOD == 1 + LRvals_resamp = LRvals_xc{xf}(randperm(numcyc_xc{xf})); % Permutation (order) + elseif SAMPMETHOD == 2 + LRvals_resamp = LRvals_xc{xf}(ceil( numcyc_xc{xf} * rand(1,numcyc_xc{xf}) )) ; % Bootstrap + end + + % Extract out a block of the same size as surrounding contig block period + % (do from the the beginning as this is simplest) + if blockdur < numcyc_xc{xf} % essentially always the case + LR_resamp = LRvals_resamp(1:blockdur); + elseif numcyc_xc{xf} == 0 + disp('Zero numcyc_xc -- should NOT happen but skipping this alternation') + keyboard + continue + elseif blockdur > numcyc_xc{xf} % sometimes (as in trackback) this will happen -- + LR_resamp = LRvals_resamp; + if abs(xf) < 10 + disp('non-trackback excursion type has a lack of permute cycles? unexpected') + keyboard + end + end + + + + % Detect alternations + altspeed_resamp = abs(diff(LR_resamp)); + altvec_resamp = altspeed_resamp > ALT_THRESHOLD; + altpers_resamp = vec2list(altvec_resamp,1:length(altvec_resamp)); + altpers_dur_resamp = altpers_resamp(:,2) - altpers_resamp(:,1) + 1; + + % If produce at least one alternation during this period, + % of at least duration of the actual alternation + if any(altpers_dur_resamp >= numcyc) + altdetect(rrr) = 1; + else + altdetect(rrr) = 0; + end + + end + + % Calculate pvalue + if sum(altdetect) ~= 0 + pval = sum(altdetect) / NUMSAMP; + else + pval = - 1 / NUMSAMP; % negative indicates that the pval is --less than-- the absolute value + end + %%%%%%% end of permute code %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % install outputs + out.altperiods{de}(ll,1) = d; % [ day ep time_a time_b tb_a tb_b numcyc xcurtype speed accel angspeed pvalue ] + out.altperiods{de}(ll,2) = ep; + out.altperiods{de}(ll,3) = time_a - epstart; + out.altperiods{de}(ll,4) = time_b - epstart; + out.altperiods{de}(ll,5) = tb_a; + out.altperiods{de}(ll,6) = tb_b; + out.altperiods{de}(ll,7) = numcyc; + out.altperiods{de}(ll,8) = xf; + out.altperiods{de}(ll,9) = speed; + out.altperiods{de}(ll,10) = accel; + out.altperiods{de}(ll,11) = angspeed; + out.altperiods{de}(ll,12) = lowspeedflag; + out.altperiods{de}(ll,13) = pval; + + end + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % Printed reportage of alt periods >= 3 cycles + NUM_LONGALT_REPORT = 3; + alttimes_long = Tbins(Altpers(find(Altpers_dur >= NUM_LONGALT_REPORT),:))-epstart; + alttimes_long_durs = Altpers_dur(Altpers_dur >= NUM_LONGALT_REPORT); + alttimes_long_sig = out.altperiods{de}(Altpers_dur >= NUM_LONGALT_REPORT,end) < P_THRESH; + alttimes_long_pval = out.altperiods{de}(Altpers_dur >= NUM_LONGALT_REPORT,end); + alttimes_long = [alttimes_long alttimes_long_durs alttimes_long_pval alttimes_long_sig]; + + if ~isempty(alttimes_long) + disp('long alt times:'); + alttimes_long + else + disp('no long alt times') + end + + + end + + cd(decodedir); + filename = sprintf('Analysis_%d_%s',EXCURTYPE,animalpref); + save(filename,'out','-v7.3'); + + end + + end + +end + + + + if plot_pvalues + + % Collect p-values + pvalues_collect = {}; + for L = 2:3 % 2 is exactly 2 alternations, 3 is >= 3 + for M = 1:2 % 1: non-low speed, 2: includes low-speed + pvalues_collect{L}{M} = []; + for aa = 1:length(animals_toplot) + animalname = animals_toplot{aa}; + animalinfo = animaldef(animalname); + cd(decodedir) + filename = dir(sprintf('Analysis_%d_%s.mat',EXCURTYPE_toplot,animalname(1:3))); + load(filename.name,'out') + for de = 1:length(out.altperiods) + if ~isempty(out.altperiods{de}) + if L == 2 + cycinds = out.altperiods{de}(:,7) == 2; + elseif L == 3 + cycinds = out.altperiods{de}(:,7) >= 3; + end + if M == 1 + speedinds = out.altperiods{de}(:,12) == 0; + elseif M == 2 + speedinds = out.altperiods{de}(:,12) == 1; + end + pvalues = abs(out.altperiods{de}(cycinds & speedinds,end)); + pvalues_collect{L}{M} = [pvalues_collect{L}{M} ; pvalues(:)]; + end + end + end + end + end + + %%% p-value scatter +% K = figure('units','normalized','outerposition',[.1 .06 .3 .7]); hold on +% WIDTH = 0.2; +% for L = LENGTH_TOPLOT +% for M = 1:2 % M 1: moving only 2: includes low speed +% if L == 2 +% ptclr = 'k'; +% ptstyle = 'o'; +% ptsize = 50; +% elseif L == 3 +% ptclr = 'k'; +% ptsize = 500; +% ptstyle = '.'; +% end +% if M == 1 +% ptclr = 'k'; +% elseif M == 2 +% ptclr = 'r'; +% end +% for tt = 1:length(pvalues_collect{L}{M}) +% onep = pvalues_collect{L}{M}(tt); +% scatter( WIDTH * rand * [1 1] - WIDTH/2,[onep onep],ptsize,ptclr,ptstyle,'linewidth',2); hold on +% end +% end +% end +% xlim([-1 1]) +% animstring = mat2str(cell2mat(animals_toplot)); +% title(sprintf('%s: p vals, alt pers',animstring),'fontsize',14,'fontweight','bold') + + % p-value histogram (linear version) + if 0 + K = figure('units','normalized','outerposition',[.1 .06 .3 .5]); hold on + for L = LENGTH_TOPLOT + subplot(2,1,L-1); + NNN = []; + for M = 1:2 % M 1: moving only 2: includes low speed + if M == 1 + faceclr = 'k'; + elseif M == 2 + faceclr = 'r'; + end + + NN = histc(pvalues_collect{L}{M},pbins); + NN(end-1) = NN(end-1) + NN(end); + NNN = [NNN ; NN(:)']; + end + + B = bar(bincenterer(pbins),NNN(:,1:(end-1))','stacked'); hold on + %set(B,'facecolor',faceclr); + xlabel('P value','fontsize',14) + ylabel('# of alternation periods','fontsize',14) + + set(gca,'fontsize',14) + + maxvalue = 10*ceil(max(sum(NNN,1))/10); + ylim([0 maxvalue]) + set(gca,'ytick',0:20:maxvalue) + + % p < 0.05 line + plot([0.05 0.05],[0 maxvalue],'k--','linewidth',2) + + % title + animstring = mat2str(cell2mat(animals_toplot)); + title(sprintf('%s: P values, alt pers',animstring),'fontsize',14,'fontweight','bold') + + + + end + + end + + + % p-value histogram ->> log version + K = figure('units','normalized','outerposition',[.1 .06 .3 .5]); hold on + for L = LENGTH_TOPLOT + subplot(2,1,L-1); + NNN = []; + collectvals = []; + for M = 1:2 % M 1: moving only 2: includes low speed + if M == 1 + faceclr = 'k'; + elseif M == 2 + faceclr = 'r'; + end + + logbins = 0:.1:(4+.1); + neglogvals = -log10(abs(pvalues_collect{L}{M})); + + NN = histc(neglogvals,logbins); + NN(end-1) = NN(end-1) + NN(end); + NNN = [NNN ; NN(:)']; + + collectvals = [collectvals ; neglogvals(:)]; + end + + B = bar(bincenterer(logbins),NNN(:,1:(end-1))','stacked'); hold on + %set(B,'facecolor',faceclr); + xlabel('-log10(P value)','fontsize',14) + ylabel('# of alternation periods','fontsize',14) + + set(gca,'fontsize',14) + + xlim([0 max(logbins)]) + maxvalue = 10*ceil(max(sum(NNN,1))/10); + if maxvalue > 20 + set(gca,'ytick',0:10:maxvalue) + else + maxvalue = maxvalue - 5; + set(gca,'ytick',0:5:maxvalue) + end + ylim([0 maxvalue]) + % p < 0.05 line + plot(-log10([0.05 0.05]),[0 maxvalue],'k--','linewidth',2) + + % title + animstring = mat2str(cell2mat(animals_toplot)); + title(sprintf('%s: P values, alt pers',animstring),'fontsize',14,'fontweight','bold') + + + + end + + + + + end + + + if plot_behavior + + % Collect behavior values + behaviorvals = {}; + dayeptimes = {}; + for M = 1:2 + for L = 3 % 2 is exactly 2 alternations, 3 is >= 3 + behaviorvals{L}{M} = []; + dayeptimes{L}{M} = []; + for aa = 1:length(animals_toplot) + animalname = animals_toplot{aa}; + cd(decodedir) + filename = dir(sprintf('Analysis_%d_%s.mat',EXCURTYPE_toplot,animalname(1:3))); + load(filename.name,'out') + for de = 1:length(out.altperiods) + if ~isempty(out.altperiods{de}) + pinds = out.altperiods{de}(:,end) < BEHAVE_P_THRESH; + if L == 2 + cycinds = out.altperiods{de}(:,7) == 2; + elseif L == 3 + cycinds = out.altperiods{de}(:,7) >= 3; + end + if M == 1 + speedinds = out.altperiods{de}(:,12) == 0; + elseif M == 2 + speedinds = out.altperiods{de}(:,12) == 1; + end + % behavior values + behav_vals = out.altperiods{de}(pinds & cycinds & speedinds,[ 9 10 11] ); + behaviorvals{L}{M} = [behaviorvals{L}{M} ; behav_vals]; + % day ep times + dayeptimes{L}{M} = [dayeptimes{L}{M} ; out.altperiods{de}(pinds & cycinds & speedinds,[1 2 3 4 7])]; + end + end + end + end + end + + % behavior scatter plot + K = figure('units','normalized','outerposition',[.5 .06 .3 .5]); hold on + for L = LENGTH_TOPLOT + subplot(2,1,L-1); + for M = 1:2 + + if M == 1 % M = 1: moving only + if L == 2 + ptsize = 500; + ptstyle = '.'; + ptclr = 'k'; %[.85 .85 .85]; + elseif L == 3 + ptsize = 500; + ptstyle = '.'; + ptclr = 'k'; + end + elseif M == 2 % M = 2: low speed too + if L == 2 + ptsize = 50; + ptstyle = 'o'; + ptclr = 'k'; %[.85 .85 .85]; + elseif L == 3 + ptsize = 50; + ptstyle = 'o'; + ptclr = 'k'; + end + end + for tt = 1:size(behaviorvals{L}{M},1) + one_b = behaviorvals{L}{M}(tt,:); + scatter(one_b(1),one_b(3),ptsize,ptclr,ptstyle,'linewidth',2); hold on + %scatter3(one_b(1),one_b(2),one_b(3),ptsize,ptclr,ptstyle,'linewidth',2); hold on + end + end + xlabel('Speed (cm/s)','fontsize',14) + ylabel('Angular speed (deg/s)','fontsize',14) + set(gca,'ytick',0:2:10) + set(gca,'xtick',0:5:40) + %zlabel('accel','fontweight','bold','fontsize',12) + grid on + xlim([0 40]) + ylim([0 10]) + set(gca,'fontsize',14) + end + + animstring = mat2str(cell2mat(animals_toplot)); + subplot(2,1,1); + title(sprintf('%s: Behav vals, alt pers',animstring),'fontsize',14,'fontweight','bold') + + end + + + + + + + if report_longest_altpers + + + % Collect behavior values + altpers_all = []; % 13 columns (same as out.altperiods above) + % 1 2 3 4 5 6 7 8 9 10 11 12 13 14 + % [ a d ep time_a time_b tb_a tb_b durcyc xc speed accel angspeed lowspeedflag pvalue ] + + for M = 1:2 + for aa = 1:length(animals_toplot) + animalname = animals_toplot{aa}; + an = find(strcmp(animalname,animals_order)); + cd(decodedir) + filename = dir(sprintf('Analysis_%d_%s.mat',EXCURTYPE_toplot,animalname(1:3))); + load(filename.name,'out') + for de = 1:length(out.altperiods) + if ~isempty(out.altperiods{de}) && ~isempty(out.altperiods{de}) + pinds = out.altperiods{de}(:,end) < P_thresh_long_altpers; + cycinds = out.altperiods{de}(:,7) >= 3; + if M == 1 + speedinds = out.altperiods{de}(:,12) == 0; + elseif M == 2 + speedinds = out.altperiods{de}(:,12) == 1; + end + % behavior values + allinfo = out.altperiods{de}(pinds & cycinds & speedinds, :); + numalts = size(allinfo,1); + altpers_all = [altpers_all ; an * ones(numalts,1) allinfo]; + end + end + end + end + + altpers_all = sortrows(altpers_all,-8); + + % Report [an day numcyc] + disp('[an day numcyc]') + disp(num2str(altpers_all(:,[1 2 8]))); + + % Report [an day ep numcyc time_a time_b] + disp('[an day ep numcyc time_a time_b]') + disp(num2str(altpers_all(:,[1 2 3 8 4 5]))); + + end + + + + + + + +if report_prevalence + for ee = 1:length(out.altperiods) + d = out.altperiods{ee}(1,1); + ep = out.altperiods{ee}(1,2); + excursions = loaddatastruct(animalinfo{2},animalinfo{3},'excursions',[d ep]); + excurlist = excursions{d}{ep}.excurlist; + inds = ismember(excurlist(:,3),[1 3 -11]); % outbound trajs (1,3) + outbound trackbacks + excurlist = excurlist(inds,1:3); + numoutbound = size(excurlist,1); + + numsigaltpers = sum(out.altperiods{ee}(:,13) < 0.05); + + disp(sprintf('# out excur: %d, # sigaltpers: %d -- %0.2f altper / excur',numoutbound,numsigaltpers,numsigaltpers/numoutbound)) + + end +end + + break + + + diff --git a/clusterless_decode/analyze_clusterless_Pos_old.m b/clusterless_decode/analyze_clusterless_Pos_old.m new file mode 100644 index 00000000..cd2d673b --- /dev/null +++ b/clusterless_decode/analyze_clusterless_Pos_old.m @@ -0,0 +1,856 @@ +datadir = '/opt/data50/kkay/__Decode'; + +calculate = 1; +plot_pvalues = 0; +plot_behavior = 0; +report_prevalence = 0; + + PHASECHOICE = 1; + + if plot_pvalues || plot_behavior + animals_toplot = {'Bond','Frank','Dave','Government'}; %{'Bond','Frank','Dave','Government'}; %'Bond','Frank','Dave', + end + if calculate + animals_tocalc = {'Bond','Frank','Dave','Government'}; %{'Egypt','Corriander'}; % 'Bond','Frank','Dave','Government','Dave','Corriander','Miles','Eight'}; + manual_days = []; + MOVINGPERIODS = 2; % 1: vel4, 2: nonimmobile05 (half sec) + EXCURTYPE = 1; % 1, 2, 3, 4 + % Spatial index + LR_prop_thresh = 0.1; % minimum posterior density that is either in L or R (versus C) + % Alternation detection parameters + ALT_THRESHOLD = 0.5; + NUMSAMP = 10000 + 1; % # of resampling methods + SAMPMETHOD = 1; % 1: permutation, 2: bootstramp + ALTDISTBINS = 0:0.05:1; + MAXALTS = 3; + P_THRESH = 0.001; + end + + % Study figures %%%%%%%%%%%% + plot_Xcorr_LR = 0; + plot_LRprop_hist = 0; + plot_LRprop_phasehist = 0; + plot_altspeed_hist = 0; + + if plot_behavior || plot_pvalues + BEHAVE_P_THRESH = 0.05; + LENGTH_TOPLOT = 2:3; + end + if plot_pvalues + pbins = 0:.01:.35; + end + + + % Data loading %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + if calculate + + cd(datadir) + + % loading %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + animals_order = {'Government','Egypt','Chapati','Dave','Higgs','Frank','Bond','Corriander',... + 'Iceland','Justice','Kapital','Laplace','Eight','Ten','Conley','Miles'}; + + clear out + + for aa = 1:length(animals_tocalc) + + animalname = animals_tocalc{aa}; + animalpref = animalname(1:3); + animalinfo = animaldef(animalname); + daydir = getdaydir(animalname); + an = find(strcmp(animalname,animals_order)); + task = loaddatastruct(animalinfo{2},animalinfo{3},'task'); + if isempty(manual_days) + days = 1:length(task); + else + days = manual_days; + end + dayeps = []; + for dd = days + if ~isempty(task{dd}) + eps = wtrackeps(task,dd); + dayeps = [dayeps ; dd * ones(length(eps),1) eps(:)]; + end + end + + % Initialize outputs + out.date = date; + out.animalname = animalinfo{3}; + out.dayeps = dayeps; + out.numperms = NUMSAMP; + out.ALTDISTBINS = ALTDISTBINS; + out.divider1 = '%%%%%%%%%%%%%%%%%%%%%%%%%%%%'; + out.altdist_ep = {}; %nan(NUMSAMP,length(ALTDISTBINS)-1); + out.numaltpers_ep = {}; %nan(NUMSAMP,MAXALTS); + out.altmeans_ep = {}; %nan(1,NUMSAMP); + out.altmedians_ep = {}; %nan(1,NUMSAMP); + out.divider2 = '%%%%%%%%%%%%%%%%%%%%%%%%%%%%'; + out.altperiods = {}; + out.descript = '1: day 2: ep 3: time_a, 4: time_b, 5: tb_a, 6: tb_b, 7: numcyc, 8: xcurtype, 9-11: speed, accel, angspeed, 12: pval'; + + for de = 1:size(dayeps,1) + + d = dayeps(de,1); + ep = dayeps(de,2); + + disp(sprintf('%s: d %d ep %d',animalinfo{3}(1:3),d,ep)) + + % load Decode file + if ~exist('P','var') || isempty(P) || ~strcmp(P.animalname,animalname) || ... + ~all(P.dayep == [d ep]) + % Decode file %%%%%%%%%%%%%% + filename = sprintf('%s_%d_%d_Pos_1.mat',animalpref,d,ep); + cd(datadir); + filedir = dir(filename); + if isempty(filedir) + disp(sprintf('not finding %s, skipping',filename)) + continue + end + load(filedir.name,'P'); + % Positional %%%%%%%%%%%%%%%%%%% + pos = loaddatastruct(animalinfo{2},animalinfo{3},'pos',P.dayep(1)); + linpos = loaddatastruct(animalinfo{2},animalinfo{3},'linpos',P.dayep(1)); + lindist = linpos{d}{ep}.statematrix.lindist; + CPbuff = choicepointer(linpos{d}{ep}); + CPprevec = lindist < CPbuff; + postimevec = pos{d}{ep}.data(:,1); + veltrace = pos{d}{ep}.data(:,5); + epstart = postimevec(1); + epend = postimevec(end); + headang = loaddatastruct(animalinfo{2}, animalinfo{3}, 'headang', d); + angveltrace = abs(headang{d}{ep}.angvel) * 180 / pi; % convert to degrees + + % Process position data + [posvec_all] = posvecmaker(lindist,linpos{d}{ep}); + [posvec_all,outersep] = posvecmaker_compress(posvec_all); + + % tets, clust spikes + spikes = loaddatastruct(animalinfo{2},animalinfo{3},'spikes',d); + cellinfo = loaddatastruct(animalinfo{2},animalinfo{3},'cellinfo'); + [adtc_list] = clusteredunits(spikes,cellinfo,an,d,ep,[],1); + maxcell = size(adtc_list,1); + selected_tets = unique(P.selected_tets); + numtets = length(selected_tets); + + % LFP %%%%%%%%%%%%%%%%%%%%%%%%%%% + thetabins = loaddatastruct(animalinfo{2},animalinfo{3},'thetabins',d); + thetatet = thetabins{d}.thetatet; + if strcmp(animalpref,'Dav') && d == 3 + disp('Dav thetatet = 1, d 3 ep 2'); + thetatet = 1; + %elseif strcmp(animpref,'Bond') && d == 6 + % disp('Bon thetatet = 30, d 6 ep 4'); + % thetatet = 30; + end + theta = loadeegstruct(animalinfo{2},animalinfo{3},'theta',d,ep,thetatet); + thetatrace = theta{d}{ep}{thetatet}.data(:,1); + tph = thetabins{d}.tph{ep}; + tvec_tph = thetabins{d}.timevec_tph{ep}; + %thetaphase = double(theta{d}{ep}{thetatet}.data(:,2))/10000; % hilbert phase + %thetatimevec = geteegtimes(theta{d}{ep}{thetatet}); + end + + % Identify Moving periods + timefilterscript + if MOVINGPERIODS == 1 + [dummy,~] = evaluatetimefilter(animalinfo{2},animalinfo{3},{ vel4 },[d ep]); + elseif MOVINGPERIODS == 2 + [dummy,~] = evaluatetimefilter(animalinfo{2},animalinfo{3},{ nonimmobile05 },[d ep]); + elseif MOVINGPERIODS == 3 + [dummy,~] = evaluatetimefilter(animalinfo{2},animalinfo{3},{ nonimmobile1 },[d ep]); + end + movingperiods = dummy{d}{ep}; % + movingvec = logical(list2vec(movingperiods,postimevec)); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % Identify excursions %%%%%%%%%%%%%%%%%%%%%%5 + excurlist = []; + if EXCURTYPE == 1 + inds = ismember(P.excurlist(:,3),[1 3 -11]); % Outbound trajs (1,3) + Center well trackbacks + R>>L and L>>R + end + excurlist = P.excurlist(inds,1:3); + numexcur = size(excurlist,1); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % Concatenate across Excurs %%%%%%%%%%%%%%%%% + + SIvec_orig = []; + Altvec_orig = []; + Contigvec = []; % 111122221111: contig vector + + analysis_totaldur = 0; + for xx = 1:numexcur + + % Vector times of Excursion + excurvec = list2vec(excurlist(xx,[1 2]),postimevec); + + % Process posterior into image + Post = P.Posts{xx}; + tvec = P.Tvecs{xx}; + + % Isolate periods to analyze + intersectvec = CPprevec & excurvec & movingvec; % Prior to CP, Excursion, Movement + intersectperiods = vec2list(intersectvec,postimevec); + + % Convert to RGB + yvec = 1:size(Post,1); + Cinds = yvec <= CPbuff; % Center arm + Linds = (yvec > CPbuff) & (yvec <= outersep); % Left arm + Rinds = ~Cinds & ~Linds; % Right arm + + choiceperiodsdur = round(sum(intersectperiods(:,2) - intersectperiods(:,1))); + analysis_totaldur = analysis_totaldur + choiceperiodsdur; + + % Vector times of analysis Periods + choicevec = logical(list2vec(intersectperiods - postimevec(1),tvec)); + Cdensity = sum(Post(Cinds,:),1); + Ldensity = sum(Post(Linds,:),1); + Rdensity = sum(Post(Rinds,:),1); + tvec_cut = tvec(choicevec); % corresponding chopped up times; + + LRdensity = Ldensity + Rdensity; + LRdensity_cut = LRdensity(choicevec); + LRthreshinds = LRdensity >= LR_prop_thresh; + + if plot_LRprop_hist + H = figure('units','normalized','outerposition',[.2 .06 .2 .2]); + hist(LRdensity(LRthreshinds(:) & choicevec(:)),20); axis tight + title('LR density','fontweight','bold','fontsize',12) + keyboard + end + + % Determine min phase for theta (for LR density) + tph; tvec_tph; + tphs = tph(lookup(tvec_cut,tvec_tph - epstart)); % theta phase of each density bin + LRthreshinds2 = LRthreshinds(choicevec); + + nbins = 12; %18; % 36 bins is 10 deg bins, after Jezek 2011 + phasebins = -pi:(2*pi/nbins):pi; + phasebins_c = phasebins(1:(end-1)); + [~,I] = histc(tphs,phasebins); + + phasehist = zeros(size(phasebins_c)); + for bbb = 1:nbins + phasehist(bbb) = sum(LRdensity_cut( (I' == bbb) & LRthreshinds2(:) )); + end + + if PHASECHOICE == 1 + [~,I2] = min(phasehist); % min + minphase = phasebins_c(I2); + elseif PHASECHOICE == 2 + [~,I2] = max(phasehist); %max + minphase = phasebins_c(I2); + end + + % Theta phase histogram of LR density %%%%%%%%%%%%%%%%%%%%%%%% + if plot_LRprop_phasehist + H = figure('units','normalized','outerposition',[.1 .3 .2 .2]); + B = bar(phasebins_c,phasehist,'histc'); hold on + set(B,'facecolor','k') + set(gca,'fontsize',14) + axis tight + plot([minphase minphase],[0 max(phasehist)],'r-','linewidth',3); hold on + title(sprintf('min phase : %0.3f',minphase)) + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%% Identify theta bins (using minphase) %%%%%%%%%%%%%%%%%%%%%%%%%%% + [tbins,~,~,~,~,~,tbblocks] ... + = thetabinner(tph,tvec_tph,minphase,[epstart epend]); + + filterinds = logical(isExcluded(mean(tbins,2),intersectperiods)); % Get bins that occur analysis periods + tbins = tbins(filterinds,:); + tbins_mean = mean(tbins,2); + totalnumcyc = size(tbins,1); + cycletimevec = 1:totalnumcyc; + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % Get Spatial Index (SI) value (0: Left, 1: Right) for each theta bin + SIvec_orig_exc = nan(1,totalnumcyc); + Altvec_orig_exc = nan(1,totalnumcyc); + + for tb = 2:size(tbins,1) + bin_a = lookup(tbins(tb,1)-epstart,tvec); % start index of bin + bin_b = lookup(tbins(tb,2)-epstart,tvec); % end index of bin + Lval = sum(Ldensity(bin_a:bin_b)); + Rval = sum(Rdensity(bin_a:bin_b)); + Cval = sum(Cdensity(bin_a:bin_b)); + LRval = Rval / (Lval + Rval); + % Filter for at least minimum LR posterior density + if (Lval + Rval)/(Cval + Rval + Lval) < LR_prop_thresh + continue + else + % Check to make sure cycle did not come after discontinuous theta bin gap + if tbins(tb,1) == tbins(tb-1,2) + % Store SI value + SIvec_orig_exc(tb) = LRval; + if tb ~= 1 && ~isnan(SIvec_orig_exc(tb-1)) + % Calculate Alt speed + Altvec_orig_exc(tb) = abs( LRval - SIvec_orig_exc(tb-1) ); % alternation speed + end + else + continue + end + end + end + + if isempty(Contigvec) + FLAGVAL = 1; + elseif Contigvec(end) == 1 + FLAGVAL = 2; + else + FLAGVAL = 1; + end + + SIvec_orig = [SIvec_orig ; SIvec_orig_exc(:)]; + Altvec_orig = [Altvec_orig ; Altvec_orig_exc(:)]; + Contigvec = [Contigvec ; FLAGVAL * ones(totalnumcyc,1) ]; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + end + + disp(sprintf('%d sec of excurperiods (%d)',analysis_totaldur,EXCURTYPE)) + + + + % Alternation analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % I. Detecte Alternations throughout epoch %%%%%%%%%%%%%% + + for PM = 1:NUMSAMP + + + + % 2. Calculate Alternation speed %%%%%%%%%%%%%%%%%%%%%%%%%%%%% + for tb = 2:size(tbins,1) + if ~isnan(Altvec_orig(tb)) % if it was in original, it can be calculated here + Altvec_perm(tb) = abs( SIvec_perm(tb) - SIvec_perm(tb-1) ); % alternation speed + end + end + N = histc(Altvec_perm,ALTDISTBINS); + N(end-1) = N(end-1) + N(end); + N(end) = []; + % install in output + out.altdist_ep{de}(PM,:) = N(:); + out.altmeans_ep{de}(PM) = nanmean(Altvec_perm); + out.altmedians_ep{de}(PM) = nanmedian(Altvec_perm); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % 3. Detect Alternation periods %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + altvec = Altvec_perm > ALT_THRESHOLD; + altpers = vec2list(altvec,cycletimevec); + % eliminate alternation periods lasting only 1 cycle + altpers( (altpers(:,2) - altpers(:,1)) == 0 , : ) = []; + + % filter for discontiguous gaps + for p = size(altpers,1):-1:1 + cyclevec = logical(list2vec(altpers(p,:),cycletimevec)); + cyclelist = tbins(cyclevec,:); + % detect any gaps within the alt period + if ~all( cyclelist(1:(end-1),2) == cyclelist(2:end,1)) + %disp('gap in alt period detected, deleting') + altpers(p,:) = []; + continue + end + end + + % Identify / store alternation periods + % (for next analyses + plotting below) + if PM == 1 + Altpers_orig = altpers; % alternation periods in terms of cycle #s + numaltperiods_orig = size(Altpers_orig,1); % total # of alternation periods + Altpers_orig_xc = nan(1,numaltperiods_orig); % excursion type of each alt period + for vv = 1:numaltperiods_orig + % mid time of alternation period + altperiods_midtime = tbins_mean(round(mean(altpers(vv,:),2))); + % excursion # + excurnum = find( excurlist(:,1) < altperiods_midtime & ... + excurlist(:,2) > altperiods_midtime ); + % excursion type + if isempty(excurnum) + Altpers_orig_xc(vv) = nan; + disp('***************** excurtype not found, setting to nan') + else + Altpers_orig_xc(vv) = excurlist(excurnum,3); + end + end + end + + % Tabulate cycle durations + altpers_dur = [altpers(:,2) - altpers(:,1)] + 1; % duration, in theta cycles, of each candidate period + for mm = 1:MAXALTS + out.numaltpers_ep{de}(PM,mm) = sum(altpers_dur >= mm); % periods with at least 4 cycles + end + if PM == 1 + Altpers_dur_orig = altpers_dur; + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + end + + % II. Detect/Permute @ single altperiod level %%%%%%%%%%%%%% + + % Preliminarily, obtain SIvals for each xcursion type for this epoch + clear cycinds_xc numcyc_xc sivals_xc + cycinds_xc{1} = find( isExcluded(tbins_mean,excurlist(excurlist(:,3) == 1,[1 2])) ); % Right + cycinds_xc{3} = find( isExcluded(tbins_mean,excurlist(excurlist(:,3) == 3,[1 2])) ); % Left + cycinds_xc{11} = find( isExcluded(tbins_mean,excurlist(excurlist(:,3) == -11,[1 2])) ); % Trackback + numcyc_xc{1} = length(cycinds_xc{1}); + numcyc_xc{3} = length(cycinds_xc{3}); + numcyc_xc{11} = length(cycinds_xc{11}); + sivals_xc{1} = SIvec_orig(cycinds_xc{1}); + sivals_xc{3} = SIvec_orig(cycinds_xc{3}); + sivals_xc{11} = SIvec_orig(cycinds_xc{11}); + + % initialize output + out.altperiods{de} = nan(numaltperiods_orig,10); % [ time_a time_b tb_a tb_b numcyc xcurtype speed accel angspeed pvalue ] + + % Iterate through each alternation period + for ll = 1:numaltperiods_orig + + % Identify basic information for this alternation period + tb_a = Altpers_orig(ll,1); % theta bin # of first cycle + tb_b = Altpers_orig(ll,2); % theta bin # of last cycle + time_a = tbins(tb_a,1); % absolute time of start of first cycle + time_b = tbins(tb_b,2); % absolute time of end of last cycle + numcyc = tb_b-tb_a+1; % number of theta cycles + xc = Altpers_orig_xc(ll); % excursion type + + if xc == -11 + xc = 11; % trackbacks >> can't index with negative value + end + if isnan(xc) + disp('xcursion type cant be found, skipping') + continue + end + + % Identify behavior (speed, accel, angular speed) + % take mean + pos_a = lookup(time_a,postimevec); + pos_b = lookup(time_b,postimevec); + speed = mean(veltrace(pos_a:pos_b)); % cm/s + accel = mean(diff(veltrace(pos_a:pos_b))); % cm/s^2 + angspeed = mean(angveltrace(pos_a:pos_b)); % deg / s + lowspeedflag = any(veltrace(pos_a:pos_b) < 4); + + % Permute each event %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % First need to identify indices of the surrounding continuous thetabin period ("block period") in which this occurred + tbc_a = nan; % "theta bin contig" first bin # of contig theta period (with decoded SI) + tbc_b = nan; % "theta bin contig" last bin # of contig theta period (with decoded SI) + blockval = tbblocks(tb_a); % dummy val is 1 or 2 -- designates a contiguous theta period + % initialize + tbc_a = tb_a; + tbc_b = tb_b; + % now find the beginning end by stepping + while tbc_a > 0 + if tbblocks(tbc_a) == blockval && ~isnan(SIvec_orig(tbc_a)) % same contiguous theta block && SIvec is decoded + tbc_a = tbc_a - 1; + else + break + end + end + while tbc_b <= totalnumcyc + if tbblocks(tbc_b) == blockval && ~isnan(SIvec_orig(tbc_b)) % same contiguous theta block && SIvec is decoded + tbc_b = tbc_b + 1; + else + break + end + end + blockdur = tbc_b - tbc_a; % duration of surrounding block period + + altdetect = nan(1,NUMSAMP); % maximum alternation is 19 + for rrr = 1:NUMSAMP + + % Resampling method (here just resample all the sivals) + if SAMPMETHOD == 1 + sivals_resamp = sivals_xc{xc}(randperm(numcyc_xc{xc})); % Permute order + elseif SAMPMETHOD == 2 + sivals_resamp = sivals_xc{xc}(ceil( numcyc_xc{xc} * rand(1,numcyc_xc{xc}) )) ; % Bootstrap + end + + % Take out the block period (from the beginning, since simplest) + si_resamp = sivals_resamp(1:blockdur); + + % Detect alternations + altspeed_resamp = abs(diff(si_resamp)); + altvec_resamp = altspeed_resamp > ALT_THRESHOLD; + altpers_resamp = vec2list(altvec_resamp,1:length(altvec_resamp)); + altpers_dur_resamp = altpers_resamp(:,2) - altpers_resamp(:,1) + 1; + + % If produce at least one alternation during this period, + % of at least duration of the actual alternation + if any(altpers_dur_resamp >= numcyc) + altdetect(rrr) = 1; + else + altdetect(rrr) = 0; + end + end + % Calculate pvalue + if sum(altdetect) ~= 0 + pval = sum(altdetect) / NUMSAMP; + else + pval = - 1 / NUMSAMP; % negative indicates that the pval is --less than-- the absolute value + end + %%%%%%%%%%%%% end of permute code %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % install outputs + out.altperiods{de}(ll,1) = d; % [ day ep time_a time_b tb_a tb_b numcyc xcurtype speed accel angspeed pvalue ] + out.altperiods{de}(ll,2) = ep; + out.altperiods{de}(ll,3) = time_a - epstart; + out.altperiods{de}(ll,4) = time_b - epstart; + out.altperiods{de}(ll,5) = tb_a; + out.altperiods{de}(ll,6) = tb_b; + out.altperiods{de}(ll,7) = numcyc; + out.altperiods{de}(ll,8) = xc; + out.altperiods{de}(ll,9) = speed; + out.altperiods{de}(ll,10) = accel; + out.altperiods{de}(ll,11) = angspeed; + out.altperiods{de}(ll,12) = lowspeedflag; + out.altperiods{de}(ll,13) = pval; + + end + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % 2. Detect Alternation periods %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + altvec = Altvec_perm > ALT_THRESHOLD; + altpers = vec2list(altvec,cycletimevec); + % eliminate alternation periods lasting only 1 cycle + altpers( (altpers(:,2) - altpers(:,1)) == 0 , : ) = []; + + % filter for discontiguous gaps + for p = size(altpers,1):-1:1 + cyclevec = logical(list2vec(altpers(p,:),cycletimevec)); + cyclelist = tbins(cyclevec,:); + % detect any gaps within the alt period + if ~all( cyclelist(1:(end-1),2) == cyclelist(2:end,1)) + %disp('gap in alt period detected, deleting') + altpers(p,:) = []; + continue + end + end + + % Store original Alternation periods (for plotting below) + if PM == 1 + Altpers_orig = altpers; + end + + % Tabulate cycle durations + altpers_dur = [altpers(:,2) - altpers(:,1)] + 1; % duration, in theta cycles, of each candidate period + for mm = 1:MAXALTS + out.numaltpers_ep{de}(PM,mm) = sum(altpers_dur >= mm); % periods with at least 4 cycles + end + if PM == 1 + Altpers_dur_orig = altpers_dur; + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % Printed reportage of alt periods >= 3 cycles + NUM_LONGALT_REPORT = 3; + alttimes_long = tbins(Altpers_orig(find(Altpers_dur_orig >= NUM_LONGALT_REPORT),:))-epstart; + alttimes_long_durs = Altpers_dur_orig(Altpers_dur_orig >= NUM_LONGALT_REPORT); + alttimes_long_sig = out.altperiods{de}(Altpers_dur_orig >= NUM_LONGALT_REPORT,end) < P_THRESH; + alttimes_long_pval = out.altperiods{de}(Altpers_dur_orig >= NUM_LONGALT_REPORT,end); + alttimes_long = [alttimes_long alttimes_long_durs alttimes_long_pval alttimes_long_sig]; + + if ~isempty(alttimes_long) + disp('long alt times:'); + alttimes_long + else + disp('no long alt times') + end + + + + end + + cd(datadir); + filename = sprintf('Analysis_%s',animalpref); + save(filename,'out','-v7.3'); + + end + + + end + + + + if plot_pvalues + + % Collect p-values + pvalues_collect = {}; + for L = 2:3 % 2 is exactly 2 alternations, 3 is >= 3 + for M = 1:2 % 1: non-low speed, 2: includes low-speed + pvalues_collect{L}{M} = []; + cd(datadir) + for aa = 1:length(animals_toplot) + animalname = animals_toplot{aa}; + animalinfo = animaldef(animalname); + filename = dir(sprintf('Analysis_%s.mat',animalname(1:3))); + load(filename.name,'out') + for de = 1:length(out.altperiods) + if ~isempty(out.altperiods{de}) + if L == 2 + cycinds = out.altperiods{de}(:,7) == 2; + elseif L == 3 + cycinds = out.altperiods{de}(:,7) >= 3; + end + if M == 1 + speedinds = out.altperiods{de}(:,12) == 0; + elseif M == 2 + speedinds = out.altperiods{de}(:,12) == 1; + end + pvalues = abs(out.altperiods{de}(cycinds & speedinds,end)); + pvalues_collect{L}{M} = [pvalues_collect{L}{M} ; pvalues(:)]; + end + end + end + end + end + + %%% p-value scatter +% K = figure('units','normalized','outerposition',[.1 .06 .3 .7]); hold on +% WIDTH = 0.2; +% for L = LENGTH_TOPLOT +% for M = 1:2 % M 1: moving only 2: includes low speed +% if L == 2 +% ptclr = 'k'; +% ptstyle = 'o'; +% ptsize = 50; +% elseif L == 3 +% ptclr = 'k'; +% ptsize = 500; +% ptstyle = '.'; +% end +% if M == 1 +% ptclr = 'k'; +% elseif M == 2 +% ptclr = 'r'; +% end +% for tt = 1:length(pvalues_collect{L}{M}) +% onep = pvalues_collect{L}{M}(tt); +% scatter( WIDTH * rand * [1 1] - WIDTH/2,[onep onep],ptsize,ptclr,ptstyle,'linewidth',2); hold on +% end +% end +% end +% xlim([-1 1]) +% animstring = mat2str(cell2mat(animals_toplot)); +% title(sprintf('%s: p vals, alt pers',animstring),'fontsize',14,'fontweight','bold') + + % p-value histogram + if 0 + K = figure('units','normalized','outerposition',[.1 .06 .3 .5]); hold on + for L = LENGTH_TOPLOT + subplot(2,1,L-1); + NNN = []; + for M = 1:2 % M 1: moving only 2: includes low speed + if M == 1 + faceclr = 'k'; + elseif M == 2 + faceclr = 'r'; + end + + NN = histc(pvalues_collect{L}{M},pbins); + NN(end-1) = NN(end-1) + NN(end); + NNN = [NNN ; NN(:)']; + end + + B = bar(bincenterer(pbins),NNN(:,1:(end-1))','stacked'); hold on + %set(B,'facecolor',faceclr); + xlabel('P value','fontsize',14) + ylabel('# of alternation periods','fontsize',14) + + set(gca,'fontsize',14) + + maxvalue = 10*ceil(max(sum(NNN,1))/10); + ylim([0 maxvalue]) + set(gca,'ytick',0:20:maxvalue) + + % p < 0.05 line + plot([0.05 0.05],[0 maxvalue],'k--','linewidth',2) + + % title + animstring = mat2str(cell2mat(animals_toplot)); + title(sprintf('%s: P values, alt pers',animstring),'fontsize',14,'fontweight','bold') + + + + end + + end + % p-value histogram ->> log version + K = figure('units','normalized','outerposition',[.1 .06 .3 .5]); hold on + for L = LENGTH_TOPLOT + subplot(2,1,L-1); + NNN = []; + collectvals = []; + for M = 1:2 % M 1: moving only 2: includes low speed + if M == 1 + faceclr = 'k'; + elseif M == 2 + faceclr = 'r'; + end + + logbins = 0:.1:(4+.1); + neglogvals = -log10(abs(pvalues_collect{L}{M})); + + NN = histc(neglogvals,logbins); + NN(end-1) = NN(end-1) + NN(end); + NNN = [NNN ; NN(:)']; + + collectvals = [collectvals ; neglogvals(:)]; + end + + B = bar(bincenterer(logbins),NNN(:,1:(end-1))','stacked'); hold on + %set(B,'facecolor',faceclr); + xlabel('-log10(P value)','fontsize',14) + ylabel('# of alternation periods','fontsize',14) + + set(gca,'fontsize',14) + + xlim([0 max(logbins)]) + maxvalue = 10*ceil(max(sum(NNN,1))/10); + if maxvalue > 20 + set(gca,'ytick',0:10:maxvalue) + else + maxvalue = maxvalue - 5; + set(gca,'ytick',0:5:maxvalue) + end + ylim([0 maxvalue]) + % p < 0.05 line + plot(-log10([0.05 0.05]),[0 maxvalue],'k--','linewidth',2) + + % title + animstring = mat2str(cell2mat(animals_toplot)); + title(sprintf('%s: P values, alt pers',animstring),'fontsize',14,'fontweight','bold') + + + + end + + + + + end + + + if plot_behavior + + % Collect behavior values + behaviorvals = {}; + dayeptimes = {}; + for M = 1:2 + for L = 2:3 % 2 is exactly 2 alternations, 3 is >= 3 + behaviorvals{L}{M} = []; + dayeptimes{L}{M} = []; + for aa = 1:length(animals_toplot) + animalname = animals_toplot{aa}; + filename = dir(sprintf('Analysis_%s.mat',animalname(1:3))); + load(filename.name,'out') + for de = 1:length(out.altperiods) + if ~isempty(out.altperiods{de}) + pinds = out.altperiods{de}(:,end) < BEHAVE_P_THRESH; + if L == 2 + cycinds = out.altperiods{de}(:,7) == 2; + elseif L == 3 + cycinds = out.altperiods{de}(:,7) >= 3; + end + if M == 1 + speedinds = out.altperiods{de}(:,12) == 0; + elseif M == 2 + speedinds = out.altperiods{de}(:,12) == 1; + end + % behavior values + behav_vals = out.altperiods{de}(pinds & cycinds & speedinds,[ 9 10 11] ); + behaviorvals{L}{M} = [behaviorvals{L}{M} ; behav_vals]; + % day ep times + dayeptimes{L}{M} = [dayeptimes{L}{M} ; out.altperiods{de}(pinds & cycinds & speedinds,[1 2 3 4 7])]; + end + end + end + end + end + + % behavior scatter plot + K = figure('units','normalized','outerposition',[.5 .06 .3 .5]); hold on + for L = LENGTH_TOPLOT + subplot(2,1,L-1); + for M = 1:2 + + if M == 1 % M = 1: moving only + if L == 2 + ptsize = 500; + ptstyle = '.'; + ptclr = 'k'; %[.85 .85 .85]; + elseif L == 3 + ptsize = 500; + ptstyle = '.'; + ptclr = 'k'; + end + elseif M == 2 % M = 2: low speed too + if L == 2 + ptsize = 50; + ptstyle = 'o'; + ptclr = 'k'; %[.85 .85 .85]; + elseif L == 3 + ptsize = 50; + ptstyle = 'o'; + ptclr = 'k'; + end + end + for tt = 1:size(behaviorvals{L}{M},1) + one_b = behaviorvals{L}{M}(tt,:); + scatter(one_b(1),one_b(3),ptsize,ptclr,ptstyle,'linewidth',2); hold on + %scatter3(one_b(1),one_b(2),one_b(3),ptsize,ptclr,ptstyle,'linewidth',2); hold on + end + end + xlabel('Speed (cm/s)','fontsize',14) + ylabel('Angular speed (deg/s)','fontsize',14) + set(gca,'ytick',0:2:10) + %zlabel('accel','fontweight','bold','fontsize',12) + grid on + xlim([0 60]) + ylim([0 10]) + set(gca,'fontsize',14) + end + + animstring = mat2str(cell2mat(animals_toplot)); + subplot(2,1,1); + title(sprintf('%s: Behav vals, alt pers',animstring),'fontsize',14,'fontweight','bold') + + end + + dayeptimes{3}{1} + + +if report_prevalence + for ee = 1:length(out.altperiods) + d = out.altperiods{ee}(1,1); + ep = out.altperiods{ee}(1,2); + excursions = loaddatastruct(animalinfo{2},animalinfo{3},'excursions',[d ep]); + excurlist = excursions{d}{ep}.excurlist; + inds = ismember(excurlist(:,3),[1 3 -11]); % outbound trajs (1,3) + outbound trackbacks + excurlist = excurlist(inds,1:3); + numoutbound = size(excurlist,1); + + numsigaltpers = sum(out.altperiods{ee}(:,13) < 0.05); + + disp(sprintf('# out excur: %d, # sigaltpers: %d -- %0.2f altper / excur',numoutbound,numsigaltpers,numsigaltpers/numoutbound)) + + end +end + + break + + + diff --git a/clusterless_decode/analyze_clusterless_sfn.m b/clusterless_decode/analyze_clusterless_sfn.m new file mode 100644 index 00000000..9c7323fe --- /dev/null +++ b/clusterless_decode/analyze_clusterless_sfn.m @@ -0,0 +1,874 @@ +datadir = '/opt/data50/kkay/__Decode'; + +calculate = 0; +plot_pvalues = 1; +plot_behavior = 1; +report_prevalence = 1; + + PHASECHOICE = 1; + + if plot_pvalues || plot_behavior + animals_toplot = {'Bond','Frank','Dave','Government'}; %{'Bond','Frank','Dave','Government'}; %'Bond','Frank','Dave', + end + if calculate + animals_tocalc = {'Bond','Frank','Dave','Government'}; %{'Egypt','Corriander'}; % 'Bond','Frank','Dave','Government','Dave','Corriander','Miles','Eight'}; + manual_days = []; + MOVINGPERIODS = 2; % 1: vel4, 2: nonimmobile05 (half sec) + EXCURPERIODS = 1; % 1: Excursions + CPbuff) & (yvec <= outersep); % Left arm + Rinds = ~Cinds & ~Linds; % Right arm + + excurperiods = []; + + % Identify excursions + if EXCURPERIODS == 1 + % Excursion + before CP + excursions = loaddatastruct(animalinfo{2},animalinfo{3},'excursions',P.dayep(1)); + excurlist = excursions{d}{ep}.excurlist; + inds = ismember(excurlist(:,3),[1 3 -11 23 32]); % Outbound trajs (1,3) + Center well trackbacks + R>>L and L>>R + excurlist = excurlist(inds,1:3); + excurvec = list2vec(excurlist(:,[1 2]),postimevec); + intersectvec = CPprevec & excurvec & movingvec; + excurperiods = vec2list(intersectvec,postimevec); + else + keyboard + end + + choiceperiodsdur = round(sum(excurperiods(:,2) - excurperiods(:,1))); + disp(sprintf('%d sec of excurperiods (%d)',choiceperiodsdur,EXCURPERIODS)) + + choicevec = logical(list2vec(excurperiods - postimevec(1),tvec)); + + Cdensity = sum(postfull(Cinds,:),1); + Ldensity = sum(postfull(Linds,:),1); + Rdensity = sum(postfull(Rinds,:),1); + tvec_cut = tvec(choicevec); % corresponding chopped up times; + + % XCorr Left vs. Right %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + if plot_Xcorr_LR + [Xcorr_LR,tlag] = xcorr(Rdensity(choicevec),Ldensity(choicevec)); + H = figure('units','normalized','outerposition',[.1 .6 .2 .2]); + plot(tlag/1000,Xcorr_LR,'k-','linewidth',2); + xlim([-0.5 0.5]) + set(gca,'fontsize',14) + end + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % Theta phase histogram of L vs. R density >> minphase %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + Ldensity; Rdensity; + LRdensity = Ldensity + Rdensity; + LRthreshinds = LRdensity >= LR_prop_thresh; + + if plot_LRprop_hist + H = figure('units','normalized','outerposition',[.2 .06 .2 .2]); + hist(LRdensity(LRthreshinds(:) & choicevec(:)),20); axis tight + title('LR density','fontweight','bold','fontsize',12) + end + + % Determine min phase for theta (for LR density) + thetaphase; + thetatimevec; + tphs = thetaphase(lookup(tvec_cut,thetatimevec - epstart)); % theta phase of each density bin + LRthreshinds2 = LRthreshinds(lookup(tvec_cut,tvec)); + nbins = 12; %18; % 36 bins is 10 deg bins, after Jezek 2011 + phasebins = -pi:(2*pi/nbins):pi; + phasebins_c = phasebins(1:(end-1)); + [~,I] = histc(tphs,phasebins); + + phasehist = zeros(size(phasebins_c)); + for bbb = 1:nbins + phasehist(bbb) = sum(LRdensity( (I == bbb) & LRthreshinds2(:) )); + end + + if PHASECHOICE == 1 + [~,I2] = min(phasehist); % min + minphase = phasebins_c(I2); + elseif PHASECHOICE == 2 + [~,I2] = max(phasehist); %max + minphase = phasebins_c(I2); + end + + % Theta phase histogram of LR density %%%%%%%%%%%%%%%%%%%%%%%% + if plot_LRprop_phasehist + H = figure('units','normalized','outerposition',[.1 .3 .2 .2]); + B = bar(phasebins_c,phasehist,'histc'); hold on + set(B,'facecolor','k') + set(gca,'fontsize',14) + axis tight + plot([minphase minphase],[0 max(phasehist)],'r-','linewidth',3); hold on + title(sprintf('min phase : %0.3f',minphase)) + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%% Identify theta bins (using minphase) %%%%%%%%%%%%%%%%%%%%%%%%%%% + [tbins,~,~,~,~,~,tbblocks] ... + = thetabinner(thetaphase,thetatimevec,minphase,[epstart epend]); + % Filter out bins that don't occur during excurperiods (outbound excursion, and prior to CP) + filterinds = logical(isExcluded(mean(tbins,2),excurperiods)); + tbins = tbins(filterinds,:); + tbins_mean = mean(tbins,2); + totalnumcyc = size(tbins,1); + cycletimevec = 1:totalnumcyc; + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % Get Spatial Index (SI) value (0: Left, 1: Right) for each theta bin + SIvec_orig = nan(1,totalnumcyc); + Altvec_orig = nan(1,totalnumcyc); + + for tb = 2:size(tbins,1) + bin_a = lookup(tbins(tb,1)-epstart,tvec); % start index of bin + bin_b = lookup(tbins(tb,2)-epstart,tvec); % end index of bin + Lval = sum(Ldensity(bin_a:bin_b)); + Rval = sum(Rdensity(bin_a:bin_b)); + Cval = sum(Cdensity(bin_a:bin_b)); + LRval = Rval / (Lval + Rval); + % Filter for at least minimum LR posterior density + if (Lval + Rval)/(Cval + Rval + Lval) < LR_prop_thresh + continue + else + % Check to make sure cycle did not come after discontinuous theta bin gap + if tbins(tb,1) == tbins(tb-1,2) + % Store SI value + SIvec_orig(tb) = LRval; + if tb ~= 1 && ~isnan(SIvec_orig(tb-1)) + % Calculate Alt speed + Altvec_orig(tb) = abs( LRval - SIvec_orig(tb-1) ); % alternation speed + end + else + continue + end + end + end + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % Alternation analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % I. Detect/Permute @ epoch level %%%%%%%%%%%%%% + + for PM = 1:NUMSAMP + + if mod(PM,100) == 0 + disp(['Permute #: ' num2str(PM)]) + end + + % 1. Permute %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % Initialize permuted output + SIvec_perm = SIvec_orig; % copy SI values + Altvec_perm = Altvec_orig; % + + % Iterate through excursion type + if PM > 1 + % must identify non-NaN cycles prior to permutation + nonnancyc = ~isnan(SIvec_orig) ; + for XTYPE = [1 3 -11] % 1: L, 2: R, -11: Trackback + exinds = excurlist(:,3) == XTYPE; + % identify all thetabins that fall in this excursion type + excurtype_pers = excurlist(exinds,[1 2]); + if ~isempty(excurtype_pers) + excurinds = logical(isExcluded(tbins_mean,excurtype_pers)); + cycinds = find(excurinds & nonnancyc(:)); +% cycinds = find( isExcluded(tbins_mean,excurtype_pers) ); + numcyc = length(cycinds); + + % Resampling method + if SAMPMETHOD == 1 + % Permute order + sampord = randperm(numcyc); + elseif SAMPMETHOD == 2 + % Bootstrap + sampord = ceil( rand(1,numcyc) * numcyc) ; + end + + % Apply sampling method to SI values + SIvec_perm(cycinds) = SIvec_orig(cycinds(sampord)); + end + end + end + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % 2. Calculate Alternation speed %%%%%%%%%%%%%%%%%%%%%%%%%%%%% + for tb = 2:size(tbins,1) + if ~isnan(Altvec_orig(tb)) % if it was in original, it can be calculated here + Altvec_perm(tb) = abs( SIvec_perm(tb) - SIvec_perm(tb-1) ); % alternation speed + end + end + N = histc(Altvec_perm,ALTDISTBINS); + N(end-1) = N(end-1) + N(end); + N(end) = []; + % install in output + out.altdist_ep{de}(PM,:) = N(:); + out.altmeans_ep{de}(PM) = nanmean(Altvec_perm); + out.altmedians_ep{de}(PM) = nanmedian(Altvec_perm); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % 3. Detect Alternation periods %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + altvec = Altvec_perm > ALT_THRESHOLD; + altpers = vec2list(altvec,cycletimevec); + % eliminate alternation periods lasting only 1 cycle + altpers( (altpers(:,2) - altpers(:,1)) == 0 , : ) = []; + + % filter for discontiguous gaps + for p = size(altpers,1):-1:1 + cyclevec = logical(list2vec(altpers(p,:),cycletimevec)); + cyclelist = tbins(cyclevec,:); + % detect any gaps within the alt period + if ~all( cyclelist(1:(end-1),2) == cyclelist(2:end,1)) + %disp('gap in alt period detected, deleting') + altpers(p,:) = []; + continue + end + end + + % Identify / store alternation periods + % (for next analyses + plotting below) + if PM == 1 + Altpers_orig = altpers; % alternation periods in terms of cycle #s + numaltperiods_orig = size(Altpers_orig,1); % total # of alternation periods + Altpers_orig_xc = nan(1,numaltperiods_orig); % excursion type of each alt period + for vv = 1:numaltperiods_orig + % mid time of alternation period + altperiods_midtime = tbins_mean(round(mean(altpers(vv,:),2))); + % excursion # + excurnum = find( excurlist(:,1) < altperiods_midtime & ... + excurlist(:,2) > altperiods_midtime ); + % excursion type + if isempty(excurnum) + Altpers_orig_xc(vv) = nan; + disp('***************** excurtype not found, setting to nan') + else + Altpers_orig_xc(vv) = excurlist(excurnum,3); + end + end + end + + % Tabulate cycle durations + altpers_dur = [altpers(:,2) - altpers(:,1)] + 1; % duration, in theta cycles, of each candidate period + for mm = 1:MAXALTS + out.numaltpers_ep{de}(PM,mm) = sum(altpers_dur >= mm); % periods with at least 4 cycles + end + if PM == 1 + Altpers_dur_orig = altpers_dur; + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + end + + % II. Detect/Permute @ single altperiod level %%%%%%%%%%%%%% + + % Preliminarily, obtain SIvals for each xcursion type for this epoch + clear cycinds_xc numcyc_xc sivals_xc + cycinds_xc{1} = find( isExcluded(tbins_mean,excurlist(excurlist(:,3) == 1,[1 2])) ); % Right + cycinds_xc{3} = find( isExcluded(tbins_mean,excurlist(excurlist(:,3) == 3,[1 2])) ); % Left + cycinds_xc{11} = find( isExcluded(tbins_mean,excurlist(excurlist(:,3) == -11,[1 2])) ); % Trackback + numcyc_xc{1} = length(cycinds_xc{1}); + numcyc_xc{3} = length(cycinds_xc{3}); + numcyc_xc{11} = length(cycinds_xc{11}); + sivals_xc{1} = SIvec_orig(cycinds_xc{1}); + sivals_xc{3} = SIvec_orig(cycinds_xc{3}); + sivals_xc{11} = SIvec_orig(cycinds_xc{11}); + + % initialize output + out.altperiods{de} = nan(numaltperiods_orig,10); % [ time_a time_b tb_a tb_b numcyc xcurtype speed accel angspeed pvalue ] + + % Iterate through each alternation period + for ll = 1:numaltperiods_orig + + % Identify basic information for this alternation period + tb_a = Altpers_orig(ll,1); % theta bin # of first cycle + tb_b = Altpers_orig(ll,2); % theta bin # of last cycle + time_a = tbins(tb_a,1); % absolute time of start of first cycle + time_b = tbins(tb_b,2); % absolute time of end of last cycle + numcyc = tb_b-tb_a+1; % number of theta cycles + xc = Altpers_orig_xc(ll); % excursion type + + if xc == -11 + xc = 11; % trackbacks >> can't index with negative value + end + if isnan(xc) + disp('xcursion type cant be found, skipping') + continue + end + + % Identify behavior (speed, accel, angular speed) + % take mean + pos_a = lookup(time_a,postimevec); + pos_b = lookup(time_b,postimevec); + speed = mean(veltrace(pos_a:pos_b)); % cm/s + accel = mean(diff(veltrace(pos_a:pos_b))); % cm/s^2 + angspeed = mean(angveltrace(pos_a:pos_b)); % deg / s + lowspeedflag = any(veltrace(pos_a:pos_b) < 4); + + % Permute each event %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % First need to identify indices of the surrounding continuous thetabin period ("block period") in which this occurred + tbc_a = nan; % "theta bin contig" first bin # of contig theta period (with decoded SI) + tbc_b = nan; % "theta bin contig" last bin # of contig theta period (with decoded SI) + blockval = tbblocks(tb_a); % dummy val is 1 or 2 -- designates a contiguous theta period + % initialize + tbc_a = tb_a; + tbc_b = tb_b; + % now find the beginning end by stepping + while tbc_a > 0 + if tbblocks(tbc_a) == blockval && ~isnan(SIvec_orig(tbc_a)) % same contiguous theta block && SIvec is decoded + tbc_a = tbc_a - 1; + else + break + end + end + while tbc_b <= totalnumcyc + if tbblocks(tbc_b) == blockval && ~isnan(SIvec_orig(tbc_b)) % same contiguous theta block && SIvec is decoded + tbc_b = tbc_b + 1; + else + break + end + end + blockdur = tbc_b - tbc_a; % duration of surrounding block period + + altdetect = nan(1,NUMSAMP); % maximum alternation is 19 + for rrr = 1:NUMSAMP + + % Resampling method (here just resample all the sivals) + if SAMPMETHOD == 1 + sivals_resamp = sivals_xc{xc}(randperm(numcyc_xc{xc})); % Permute order + elseif SAMPMETHOD == 2 + sivals_resamp = sivals_xc{xc}(ceil( numcyc_xc{xc} * rand(1,numcyc_xc{xc}) )) ; % Bootstrap + end + + % Take out the block period (from the beginning, since simplest) + si_resamp = sivals_resamp(1:blockdur); + + % Detect alternations + altspeed_resamp = abs(diff(si_resamp)); + altvec_resamp = altspeed_resamp > ALT_THRESHOLD; + altpers_resamp = vec2list(altvec_resamp,1:length(altvec_resamp)); + altpers_dur_resamp = altpers_resamp(:,2) - altpers_resamp(:,1) + 1; + + % If produce at least one alternation during this period, + % of at least duration of the actual alternation + if any(altpers_dur_resamp >= numcyc) + altdetect(rrr) = 1; + else + altdetect(rrr) = 0; + end + end + % Calculate pvalue + if sum(altdetect) ~= 0 + pval = sum(altdetect) / NUMSAMP; + else + pval = - 1 / NUMSAMP; % negative indicates that the pval is --less than-- the absolute value + end + %%%%%%%%%%%%% end of permute code %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % install outputs + out.altperiods{de}(ll,1) = d; % [ day ep time_a time_b tb_a tb_b numcyc xcurtype speed accel angspeed pvalue ] + out.altperiods{de}(ll,2) = ep; + out.altperiods{de}(ll,3) = time_a - epstart; + out.altperiods{de}(ll,4) = time_b - epstart; + out.altperiods{de}(ll,5) = tb_a; + out.altperiods{de}(ll,6) = tb_b; + out.altperiods{de}(ll,7) = numcyc; + out.altperiods{de}(ll,8) = xc; + out.altperiods{de}(ll,9) = speed; + out.altperiods{de}(ll,10) = accel; + out.altperiods{de}(ll,11) = angspeed; + out.altperiods{de}(ll,12) = lowspeedflag; + out.altperiods{de}(ll,13) = pval; + + end + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % 2. Detect Alternation periods %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + altvec = Altvec_perm > ALT_THRESHOLD; + altpers = vec2list(altvec,cycletimevec); + % eliminate alternation periods lasting only 1 cycle + altpers( (altpers(:,2) - altpers(:,1)) == 0 , : ) = []; + + % filter for discontiguous gaps + for p = size(altpers,1):-1:1 + cyclevec = logical(list2vec(altpers(p,:),cycletimevec)); + cyclelist = tbins(cyclevec,:); + % detect any gaps within the alt period + if ~all( cyclelist(1:(end-1),2) == cyclelist(2:end,1)) + %disp('gap in alt period detected, deleting') + altpers(p,:) = []; + continue + end + end + + % Store original Alternation periods (for plotting below) + if PM == 1 + Altpers_orig = altpers; + end + + % Tabulate cycle durations + altpers_dur = [altpers(:,2) - altpers(:,1)] + 1; % duration, in theta cycles, of each candidate period + for mm = 1:MAXALTS + out.numaltpers_ep{de}(PM,mm) = sum(altpers_dur >= mm); % periods with at least 4 cycles + end + if PM == 1 + Altpers_dur_orig = altpers_dur; + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % Printed reportage of alt periods >= 3 cycles + NUM_LONGALT_REPORT = 3; + alttimes_long = tbins(Altpers_orig(find(Altpers_dur_orig >= NUM_LONGALT_REPORT),:))-epstart; + alttimes_long_durs = Altpers_dur_orig(Altpers_dur_orig >= NUM_LONGALT_REPORT); + alttimes_long_sig = out.altperiods{de}(Altpers_dur_orig >= NUM_LONGALT_REPORT,end) < P_THRESH; + alttimes_long_pval = out.altperiods{de}(Altpers_dur_orig >= NUM_LONGALT_REPORT,end); + alttimes_long = [alttimes_long alttimes_long_durs alttimes_long_pval alttimes_long_sig]; + + if ~isempty(alttimes_long) + disp('long alt times:'); + alttimes_long + else + disp('no long alt times') + end + + + + end + + cd(datadir); + filename = sprintf('Analysis_%s',animalpref); + save(filename,'out','-v7.3'); + + end + + + end + + + + if plot_pvalues + + % Collect p-values + pvalues_collect = {}; + for L = 2:3 % 2 is exactly 2 alternations, 3 is >= 3 + for M = 1:2 % 1: non-low speed, 2: includes low-speed + pvalues_collect{L}{M} = []; + cd(datadir) + for aa = 1:length(animals_toplot) + animalname = animals_toplot{aa}; + animalinfo = animaldef(animalname); + filename = dir(sprintf('Analysis_%s.mat',animalname(1:3))); + load(filename.name,'out') + for de = 1:length(out.altperiods) + if ~isempty(out.altperiods{de}) + if L == 2 + cycinds = out.altperiods{de}(:,7) == 2; + elseif L == 3 + cycinds = out.altperiods{de}(:,7) >= 3; + end + if M == 1 + speedinds = out.altperiods{de}(:,12) == 0; + elseif M == 2 + speedinds = out.altperiods{de}(:,12) == 1; + end + pvalues = abs(out.altperiods{de}(cycinds & speedinds,end)); + pvalues_collect{L}{M} = [pvalues_collect{L}{M} ; pvalues(:)]; + end + end + end + end + end + + %%% p-value scatter +% K = figure('units','normalized','outerposition',[.1 .06 .3 .7]); hold on +% WIDTH = 0.2; +% for L = LENGTH_TOPLOT +% for M = 1:2 % M 1: moving only 2: includes low speed +% if L == 2 +% ptclr = 'k'; +% ptstyle = 'o'; +% ptsize = 50; +% elseif L == 3 +% ptclr = 'k'; +% ptsize = 500; +% ptstyle = '.'; +% end +% if M == 1 +% ptclr = 'k'; +% elseif M == 2 +% ptclr = 'r'; +% end +% for tt = 1:length(pvalues_collect{L}{M}) +% onep = pvalues_collect{L}{M}(tt); +% scatter( WIDTH * rand * [1 1] - WIDTH/2,[onep onep],ptsize,ptclr,ptstyle,'linewidth',2); hold on +% end +% end +% end +% xlim([-1 1]) +% animstring = mat2str(cell2mat(animals_toplot)); +% title(sprintf('%s: p vals, alt pers',animstring),'fontsize',14,'fontweight','bold') + + % p-value histogram + if 0 + K = figure('units','normalized','outerposition',[.1 .06 .3 .5]); hold on + for L = LENGTH_TOPLOT + subplot(2,1,L-1); + NNN = []; + for M = 1:2 % M 1: moving only 2: includes low speed + if M == 1 + faceclr = 'k'; + elseif M == 2 + faceclr = 'r'; + end + + NN = histc(pvalues_collect{L}{M},pbins); + NN(end-1) = NN(end-1) + NN(end); + NNN = [NNN ; NN(:)']; + end + + B = bar(bincenterer(pbins),NNN(:,1:(end-1))','stacked'); hold on + %set(B,'facecolor',faceclr); + xlabel('P value','fontsize',14) + ylabel('# of alternation periods','fontsize',14) + + set(gca,'fontsize',14) + + maxvalue = 10*ceil(max(sum(NNN,1))/10); + ylim([0 maxvalue]) + set(gca,'ytick',0:20:maxvalue) + + % p < 0.05 line + plot([0.05 0.05],[0 maxvalue],'k--','linewidth',2) + + % title + animstring = mat2str(cell2mat(animals_toplot)); + title(sprintf('%s: P values, alt pers',animstring),'fontsize',14,'fontweight','bold') + + + + end + + end + % p-value histogram ->> log version + K = figure('units','normalized','outerposition',[.1 .06 .3 .5]); hold on + for L = LENGTH_TOPLOT + subplot(2,1,L-1); + NNN = []; + collectvals = []; + for M = 1:2 % M 1: moving only 2: includes low speed + if M == 1 + faceclr = 'k'; + elseif M == 2 + faceclr = 'r'; + end + + logbins = 0:.1:(4+.1); + neglogvals = -log10(abs(pvalues_collect{L}{M})); + + NN = histc(neglogvals,logbins); + NN(end-1) = NN(end-1) + NN(end); + NNN = [NNN ; NN(:)']; + + collectvals = [collectvals ; neglogvals(:)]; + end + + B = bar(bincenterer(logbins),NNN(:,1:(end-1))','stacked'); hold on + %set(B,'facecolor',faceclr); + xlabel('-log10(P value)','fontsize',14) + ylabel('# of alternation periods','fontsize',14) + + set(gca,'fontsize',14) + + xlim([0 max(logbins)]) + maxvalue = 10*ceil(max(sum(NNN,1))/10); + if maxvalue > 20 + set(gca,'ytick',0:10:maxvalue) + else + maxvalue = maxvalue - 5; + set(gca,'ytick',0:5:maxvalue) + end + ylim([0 maxvalue]) + % p < 0.05 line + plot(-log10([0.05 0.05]),[0 maxvalue],'k--','linewidth',2) + + % title + animstring = mat2str(cell2mat(animals_toplot)); + title(sprintf('%s: P values, alt pers',animstring),'fontsize',14,'fontweight','bold') + + + + end + + + + + end + + + if plot_behavior + + % Collect behavior values + behaviorvals = {}; + dayeptimes = {}; + for M = 1:2 + for L = 2:3 % 2 is exactly 2 alternations, 3 is >= 3 + behaviorvals{L}{M} = []; + dayeptimes{L}{M} = []; + for aa = 1:length(animals_toplot) + animalname = animals_toplot{aa}; + filename = dir(sprintf('Analysis_%s.mat',animalname(1:3))); + load(filename.name,'out') + for de = 1:length(out.altperiods) + if ~isempty(out.altperiods{de}) + pinds = out.altperiods{de}(:,end) < BEHAVE_P_THRESH; + if L == 2 + cycinds = out.altperiods{de}(:,7) == 2; + elseif L == 3 + cycinds = out.altperiods{de}(:,7) >= 3; + end + if M == 1 + speedinds = out.altperiods{de}(:,12) == 0; + elseif M == 2 + speedinds = out.altperiods{de}(:,12) == 1; + end + % behavior values + behav_vals = out.altperiods{de}(pinds & cycinds & speedinds,[ 9 10 11] ); + behaviorvals{L}{M} = [behaviorvals{L}{M} ; behav_vals]; + % day ep times + dayeptimes{L}{M} = [dayeptimes{L}{M} ; out.altperiods{de}(pinds & cycinds & speedinds,[1 2 3 4 7])]; + end + end + end + end + end + + % behavior scatter plot + K = figure('units','normalized','outerposition',[.5 .06 .3 .5]); hold on + for L = LENGTH_TOPLOT + subplot(2,1,L-1); + for M = 1:2 + + if M == 1 % M = 1: moving only + if L == 2 + ptsize = 500; + ptstyle = '.'; + ptclr = 'k'; %[.85 .85 .85]; + elseif L == 3 + ptsize = 500; + ptstyle = '.'; + ptclr = 'k'; + end + elseif M == 2 % M = 2: low speed too + if L == 2 + ptsize = 50; + ptstyle = 'o'; + ptclr = 'k'; %[.85 .85 .85]; + elseif L == 3 + ptsize = 50; + ptstyle = 'o'; + ptclr = 'k'; + end + end + for tt = 1:size(behaviorvals{L}{M},1) + one_b = behaviorvals{L}{M}(tt,:); + scatter(one_b(1),one_b(3),ptsize,ptclr,ptstyle,'linewidth',2); hold on + %scatter3(one_b(1),one_b(2),one_b(3),ptsize,ptclr,ptstyle,'linewidth',2); hold on + end + end + xlabel('Speed (cm/s)','fontsize',14) + ylabel('Angular speed (deg/s)','fontsize',14) + set(gca,'ytick',0:2:10) + %zlabel('accel','fontweight','bold','fontsize',12) + grid on + xlim([0 60]) + ylim([0 10]) + set(gca,'fontsize',14) + end + + animstring = mat2str(cell2mat(animals_toplot)); + subplot(2,1,1); + title(sprintf('%s: Behav vals, alt pers',animstring),'fontsize',14,'fontweight','bold') + + end + + dayeptimes{3}{1} + + +if report_prevalence + for ee = 1:length(out.altperiods) + d = out.altperiods{ee}(1,1); + ep = out.altperiods{ee}(1,2); + excursions = loaddatastruct(animalinfo{2},animalinfo{3},'excursions',[d ep]); + excurlist = excursions{d}{ep}.excurlist; + inds = ismember(excurlist(:,3),[1 3 -11]); % outbound trajs (1,3) + outbound trackbacks + excurlist = excurlist(inds,1:3); + numoutbound = size(excurlist,1); + + numsigaltpers = sum(out.altperiods{ee}(:,13) < 0.05); + + disp(sprintf('# out excur: %d, # sigaltpers: %d -- %0.2f altper / excur',numoutbound,numsigaltpers,numsigaltpers/numoutbound)) + + end +end + + break + + + diff --git a/clusterless_decode/bond04_clusterless_alltet.m b/clusterless_decode/bond04_clusterless_alltet.m new file mode 100755 index 00000000..df807e64 --- /dev/null +++ b/clusterless_decode/bond04_clusterless_alltet.m @@ -0,0 +1,885 @@ +%load('bond04_vecLF_fold.mat'); +%load('bond04_vecLF.mat'); + +load('/opt/data13/kkay/Bon/bonlinpos04.mat'); + +%%% basic data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +postimevec=linpos{4}{2}.statematrix.time; +poslin=vecLF(:,2); + poslin2 = poslin'; + +%%% basic values: xdel, xbins, dt %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +xdel=0.1; +xbins=min(poslin):xdel:max(poslin); % -3 : 0.1 : 3 + numlinbins = length(xbins); +dt=postimevec(2)-postimevec(1); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% make State M %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +stateM=[]; +[~,binind]=histc(poslin,xbins); % histogram of linearized positions +xbinstate=[binind(1:end-1) binind(2:end)]; % [ ] +%by column = departuring +for bb = 1:numlinbins + % nextstates is the next states; + nextstates = xbinstate( xbinstate(:,1)==bb , 2 ); + if ~isempty(nextstates) + % histograms the nextstates and normalizes the distribution (why?) + stateM(:,bb) = histc( nextstates , linspace(1,numlinbins,numlinbins))./size(nextstates,1 ); + elseif isempty(nextstates) + stateM(:,bb) = zeros(1,numlinbins); + end +end +K = inline('exp(-(x.^2+y.^2)/2/sig^2)'); %gaussian +[dx,dy] = meshgrid([-1:1]); +sig = 0.5; +weight = K(sig,dx,dy)/sum(sum(K(sig,dx,dy))); %normalizing weights +stateM_gaus = conv2(stateM,weight,'same'); %gaussian smoothed +stateM_gausnorm = stateM_gaus*diag(1./sum(stateM_gaus,1)); %normalized to confine probability to 1 + + + + +%%%%%% encode per tetrode +load('/mnt/vortex_data/mkarlsso/backup/bond/bond04/01-149/bond04-01_params.mat'); +%ind_t1=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)&filedata.params(:,8)>10&(filedata.params(:,2)>100|filedata.params(:,3)>100|filedata.params(:,4)>100|filedata.params(:,5)>100)); + % identify spikes that occur in transcribed epoch times +ind_t1= find( filedata.params(:,1)/10000 >= postimevec(1) & filedata.params(:,1)/10000 <= postimevec(end) ); + % spike timestamps +time_t1=filedata.params(ind_t1,1); + % mark vector (4 channel amplitudes) +mark0_t1=[filedata.params(ind_t1,2) filedata.params(ind_t1,3) filedata.params(ind_t1,4) filedata.params(ind_t1,5)]; + % spike timestamps in sec +time2_t1=time_t1/10000; + +spikeT0_t1=time2_t1; +[procInd0_t1,procInd1_t1]=histc(spikeT0_t1,postimevec); +procInd_t1=find(procInd0_t1); +spikeT_t1=postimevec(procInd_t1); +spike_t1=procInd0_t1'; +[~,rawInd0_t1]=histc(spikeT0_t1,time2_t1); +markAll_t1(:,1)=procInd1_t1;markAll_t1(:,2:5)=mark0_t1(rawInd0_t1(rawInd0_t1~=0),:); +mdel=20; ms=min(min(markAll_t1(:,2:5))):mdel:max(max(markAll_t1(:,2:5))); +sxker=2*xdel; smker = mdel; T=size(postimevec,1); +occ=normpdf(xbins'*ones(1,T),ones(length(xbins),1)*poslin2,sxker)*ones(T,length(ms)); +%occ: columns are identical; occupancy based on position; denominator +Xnum_t1=normpdf(xbins'*ones(1,length(spikeT0_t1)),ones(length(xbins),1)*poslin2(procInd1_t1),sxker); +%Xnum: Gaussian kernel estimators for position +Lint_t1=sum(Xnum_t1,2)./occ(:,1)./dt; %integral +Lint_t1=Lint_t1./sum(Lint_t1); + +load('/mnt/vortex_data/mkarlsso/backup/bond/bond04/02-151/bond04-02_params.mat'); +%ind_t2=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)&filedata.params(:,8)>10&(filedata.params(:,2)>100|filedata.params(:,3)>100|filedata.params(:,4)>100|filedata.params(:,5)>100)); +ind_t2=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)); +time_t2=filedata.params(ind_t2,1); +mark0_t2=[filedata.params(ind_t2,2) filedata.params(ind_t2,3) filedata.params(ind_t2,4) filedata.params(ind_t2,5)]; +time2_t2=time_t2/10000; +spikeT0_t2=time2_t2; +[procInd0_t2,procInd1_t2]=histc(spikeT0_t2,postimevec); +procInd_t2=find(procInd0_t2); +spikeT_t2=postimevec(procInd_t2); +spike_t2=procInd0_t2'; +[~,rawInd0_t2]=histc(spikeT0_t2,time2_t2); +markAll_t2(:,1)=procInd1_t2;markAll_t2(:,2:5)=mark0_t2(rawInd0_t2(rawInd0_t2~=0),:); +mdel=20; ms=min(min(markAll_t2(:,2:5))):mdel:max(max(markAll_t2(:,2:5))); +sxker=2*xdel; smker = mdel; T=size(postimevec,1); +occ=normpdf(xbins'*ones(1,T),ones(length(xbins),1)*poslin2,sxker)*ones(T,length(ms)); +%occ: columns are identical; occupancy based on position; denominator +Xnum_t2=normpdf(xbins'*ones(1,length(spikeT0_t2)),ones(length(xbins),1)*poslin2(procInd1_t2),sxker); +%Xnum: Gaussian kernel estimators for position +Lint_t2=sum(Xnum_t2,2)./occ(:,1)./dt; %integral +Lint_t2=Lint_t2./sum(Lint_t2); + +load('/mnt/vortex_data/mkarlsso/backup/bond/bond04/04-107/bond04-04_params.mat'); +%ind_t4=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)&filedata.params(:,8)>10&(filedata.params(:,2)>100|filedata.params(:,3)>100|filedata.params(:,4)>100|filedata.params(:,5)>100)); +ind_t4=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)); +time_t4=filedata.params(ind_t4,1); +mark0_t4=[filedata.params(ind_t4,2) filedata.params(ind_t4,3) filedata.params(ind_t4,4) filedata.params(ind_t4,5)]; +time2_t4=time_t4/10000; +spikeT0_t4=time2_t4; +[procInd0_t4,procInd1_t4]=histc(spikeT0_t4,postimevec); +procInd_t4=find(procInd0_t4); +spikeT_t4=postimevec(procInd_t4); +spike_t4=procInd0_t4'; +[~,rawInd0_t4]=histc(spikeT0_t4,time2_t4); +markAll_t4(:,1)=procInd1_t4;markAll_t4(:,2:5)=mark0_t4(rawInd0_t4(rawInd0_t4~=0),:); +mdel=20; ms=min(min(markAll_t4(:,2:5))):mdel:max(max(markAll_t4(:,2:5))); +sxker=2*xdel; smker = mdel; T=size(postimevec,1); +occ=normpdf(xbins'*ones(1,T),ones(length(xbins),1)*poslin2,sxker)*ones(T,length(ms)); +%occ: columns are identical; occupancy based on position; denominator +Xnum_t4=normpdf(xbins'*ones(1,length(spikeT0_t4)),ones(length(xbins),1)*poslin2(procInd1_t4),sxker); +%Xnum: Gaussian kernel estimators for position +Lint_t4=sum(Xnum_t4,2)./occ(:,1)./dt; %integral +Lint_t4=Lint_t4./sum(Lint_t4); + +load('/mnt/vortex_data/mkarlsso/backup/bond/bond04/05-105/bond04-05_params.mat'); +%ind_t5=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)&filedata.params(:,8)>10&(filedata.params(:,2)>100|filedata.params(:,3)>100|filedata.params(:,4)>100|filedata.params(:,5)>100)); +ind_t5=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)); +time_t5=filedata.params(ind_t5,1); +mark0_t5=[filedata.params(ind_t5,2) filedata.params(ind_t5,3) filedata.params(ind_t5,4) filedata.params(ind_t5,5)]; +time2_t5=time_t5/10000; +spikeT0_t5=time2_t5; +[procInd0_t5,procInd1_t5]=histc(spikeT0_t5,postimevec); +procInd_t5=find(procInd0_t5); +spikeT_t5=postimevec(procInd_t5); +spike_t5=procInd0_t5'; +[~,rawInd0_t5]=histc(spikeT0_t5,time2_t5); +markAll_t5(:,1)=procInd1_t5;markAll_t5(:,2:5)=mark0_t5(rawInd0_t5(rawInd0_t5~=0),:); +mdel=20; ms=min(min(markAll_t5(:,2:5))):mdel:max(max(markAll_t5(:,2:5))); +sxker=2*xdel; smker = mdel; T=size(postimevec,1); +occ=normpdf(xbins'*ones(1,T),ones(length(xbins),1)*poslin2,sxker)*ones(T,length(ms)); +%occ: columns are identical; occupancy based on position; denominator +Xnum_t5=normpdf(xbins'*ones(1,length(spikeT0_t5)),ones(length(xbins),1)*poslin2(procInd1_t5),sxker); +%Xnum: Gaussian kernel estimators for position +Lint_t5=sum(Xnum_t5,2)./occ(:,1)./dt; %integral +Lint_t5=Lint_t5./sum(Lint_t5); + +load('/mnt/vortex_data/mkarlsso/backup/bond/bond04/07-162/bond04-07_params.mat'); +%ind_t7=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)&filedata.params(:,8)>10&(filedata.params(:,2)>100|filedata.params(:,3)>100|filedata.params(:,4)>100|filedata.params(:,5)>100)); +ind_t7=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)); +time_t7=filedata.params(ind_t7,1); +mark0_t7=[filedata.params(ind_t7,2) filedata.params(ind_t7,3) filedata.params(ind_t7,4) filedata.params(ind_t7,5)]; +time2_t7=time_t7/10000; +spikeT0_t7=time2_t7; +[procInd0_t7,procInd1_t7]=histc(spikeT0_t7,postimevec); +procInd_t7=find(procInd0_t7); +spikeT_t7=postimevec(procInd_t7); +spike_t7=procInd0_t7'; +[~,rawInd0_t7]=histc(spikeT0_t7,time2_t7); +markAll_t7(:,1)=procInd1_t7;markAll_t7(:,2:5)=mark0_t7(rawInd0_t7(rawInd0_t7~=0),:); +mdel=20; ms=min(min(markAll_t7(:,2:5))):mdel:max(max(markAll_t7(:,2:5))); +sxker=2*xdel; smker = mdel; T=size(postimevec,1); +occ=normpdf(xbins'*ones(1,T),ones(length(xbins),1)*poslin2,sxker)*ones(T,length(ms)); +%occ: columns are identical; occupancy based on position; denominator +Xnum_t7=normpdf(xbins'*ones(1,length(spikeT0_t7)),ones(length(xbins),1)*poslin2(procInd1_t7),sxker); +%Xnum: Gaussian kernel estimators for position +Lint_t7=sum(Xnum_t7,2)./occ(:,1)./dt; %integral +Lint_t7=Lint_t7./sum(Lint_t7); + +load('/mnt/vortex_data/mkarlsso/backup/bond/bond04/10-146/bond04-10_params.mat'); +%ind_t10=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)&filedata.params(:,8)>10&(filedata.params(:,2)>100|filedata.params(:,3)>100|filedata.params(:,4)>100|filedata.params(:,5)>100)); +ind_t10=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)); +time_t10=filedata.params(ind_t10,1); +mark0_t10=[filedata.params(ind_t10,2) filedata.params(ind_t10,3) filedata.params(ind_t10,4) filedata.params(ind_t10,5)]; +time2_t10=time_t10/10000; +spikeT0_t10=time2_t10; +[procInd0_t10,procInd1_t10]=histc(spikeT0_t10,postimevec); +procInd_t10=find(procInd0_t10); +spikeT_t10=postimevec(procInd_t10); +spike_t10=procInd0_t10'; +[~,rawInd0_t10]=histc(spikeT0_t10,time2_t10); +markAll_t10(:,1)=procInd1_t10;markAll_t10(:,2:5)=mark0_t10(rawInd0_t10(rawInd0_t10~=0),:); +mdel=20; ms=min(min(markAll_t10(:,2:5))):mdel:max(max(markAll_t10(:,2:5))); +sxker=2*xdel; smker = mdel; T=size(postimevec,1); +occ=normpdf(xbins'*ones(1,T),ones(length(xbins),1)*poslin2,sxker)*ones(T,length(ms)); +%occ: columns are identical; occupancy based on position; denominator +Xnum_t10=normpdf(xbins'*ones(1,length(spikeT0_t10)),ones(length(xbins),1)*poslin2(procInd1_t10),sxker); +%Xnum: Gaussian kernel estimators for position +Lint_t10=sum(Xnum_t10,2)./occ(:,1)./dt; %integral +Lint_t10=Lint_t10./sum(Lint_t10); + +load('/mnt/vortex_data/mkarlsso/backup/bond/bond04/11-096/bond04-11_params.mat'); +%ind_t11=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)&filedata.params(:,8)>10&(filedata.params(:,2)>100|filedata.params(:,3)>100|filedata.params(:,4)>100|filedata.params(:,5)>100)); +ind_t11=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)); +time_t11=filedata.params(ind_t11,1); +mark0_t11=[filedata.params(ind_t11,2) filedata.params(ind_t11,3) filedata.params(ind_t11,4) filedata.params(ind_t11,5)]; +time2_t11=time_t11/10000; +spikeT0_t11=time2_t11; +[procInd0_t11,procInd1_t11]=histc(spikeT0_t11,postimevec); +procInd_t11=find(procInd0_t11); +spikeT_t11=postimevec(procInd_t11); +spike_t11=procInd0_t11'; +[~,rawInd0_t11]=histc(spikeT0_t11,time2_t11); +markAll_t11(:,1)=procInd1_t11;markAll_t11(:,2:5)=mark0_t11(rawInd0_t11(rawInd0_t11~=0),:); +mdel=20; ms=min(min(markAll_t11(:,2:5))):mdel:max(max(markAll_t11(:,2:5))); +sxker=2*xdel; smker = mdel; T=size(postimevec,1); +occ=normpdf(xbins'*ones(1,T),ones(length(xbins),1)*poslin2,sxker)*ones(T,length(ms)); +%occ: columns are identical; occupancy based on position; denominator +Xnum_t11=normpdf(xbins'*ones(1,length(spikeT0_t11)),ones(length(xbins),1)*poslin2(procInd1_t11),sxker); +%Xnum: Gaussian kernel estimators for position +Lint_t11=sum(Xnum_t11,2)./occ(:,1)./dt; %integral +Lint_t11=Lint_t11./sum(Lint_t11); + +load('/mnt/vortex_data/mkarlsso/backup/bond/bond04/12-136/bond04-12_params.mat'); +%ind_t12=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)&filedata.params(:,8)>10&(filedata.params(:,2)>100|filedata.params(:,3)>100|filedata.params(:,4)>100|filedata.params(:,5)>100)); +ind_t12=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)); +time_t12=filedata.params(ind_t12,1); +mark0_t12=[filedata.params(ind_t12,2) filedata.params(ind_t12,3) filedata.params(ind_t12,4) filedata.params(ind_t12,5)]; +time2_t12=time_t12/10000; +spikeT0_t12=time2_t12; +[procInd0_t12,procInd1_t12]=histc(spikeT0_t12,postimevec); +procInd_t12=find(procInd0_t12); +spikeT_t12=postimevec(procInd_t12); +spike_t12=procInd0_t12'; +[~,rawInd0_t12]=histc(spikeT0_t12,time2_t12); +markAll_t12(:,1)=procInd1_t12;markAll_t12(:,2:5)=mark0_t12(rawInd0_t12(rawInd0_t12~=0),:); +mdel=20; ms=min(min(markAll_t12(:,2:5))):mdel:max(max(markAll_t12(:,2:5))); +sxker=2*xdel; smker = mdel; T=size(postimevec,1); +occ=normpdf(xbins'*ones(1,T),ones(length(xbins),1)*poslin2,sxker)*ones(T,length(ms)); +%occ: columns are identical; occupancy based on position; denominator +Xnum_t12=normpdf(xbins'*ones(1,length(spikeT0_t12)),ones(length(xbins),1)*poslin2(procInd1_t12),sxker); +%Xnum: Gaussian kernel estimators for position +Lint_t12=sum(Xnum_t12,2)./occ(:,1)./dt; %integral +Lint_t12=Lint_t12./sum(Lint_t12); + +load('/mnt/vortex_data/mkarlsso/backup/bond/bond04/13-094/bond04-13_params.mat'); +%ind_t13=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)&filedata.params(:,8)>10&(filedata.params(:,2)>100|filedata.params(:,3)>100|filedata.params(:,4)>100|filedata.params(:,5)>100)); +ind_t13=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)); +time_t13=filedata.params(ind_t13,1); +mark0_t13=[filedata.params(ind_t13,2) filedata.params(ind_t13,3) filedata.params(ind_t13,4) filedata.params(ind_t13,5)]; +time2_t13=time_t13/10000; +spikeT0_t13=time2_t13; +[procInd0_t13,procInd1_t13]=histc(spikeT0_t13,postimevec); +procInd_t13=find(procInd0_t13); +spikeT_t13=postimevec(procInd_t13); +spike_t13=procInd0_t13'; +[~,rawInd0_t13]=histc(spikeT0_t13,time2_t13); +markAll_t13(:,1)=procInd1_t13;markAll_t13(:,2:5)=mark0_t13(rawInd0_t13(rawInd0_t13~=0),:); +mdel=20; ms=min(min(markAll_t13(:,2:5))):mdel:max(max(markAll_t13(:,2:5))); +sxker=2*xdel; smker = mdel; T=size(postimevec,1); +occ=normpdf(xbins'*ones(1,T),ones(length(xbins),1)*poslin2,sxker)*ones(T,length(ms)); +%occ: columns are identical; occupancy based on position; denominator +Xnum_t13=normpdf(xbins'*ones(1,length(spikeT0_t13)),ones(length(xbins),1)*poslin2(procInd1_t13),sxker); +%Xnum: Gaussian kernel estimators for position +Lint_t13=sum(Xnum_t13,2)./occ(:,1)./dt; %integral +Lint_t13=Lint_t13./sum(Lint_t13); + +load('/mnt/vortex_data/mkarlsso/backup/bond/bond04/14-096/bond04-14_params.mat'); +%ind_t14=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)&filedata.params(:,8)>10&(filedata.params(:,2)>100|filedata.params(:,3)>100|filedata.params(:,4)>100|filedata.params(:,5)>100)); +ind_t14=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)); +time_t14=filedata.params(ind_t14,1); +mark0_t14=[filedata.params(ind_t14,2) filedata.params(ind_t14,3) filedata.params(ind_t14,4) filedata.params(ind_t14,5)]; +time2_t14=time_t14/10000; +spikeT0_t14=time2_t14; +[procInd0_t14,procInd1_t14]=histc(spikeT0_t14,postimevec); +procInd_t14=find(procInd0_t14); +spikeT_t14=postimevec(procInd_t14); +spike_t14=procInd0_t14'; +[~,rawInd0_t14]=histc(spikeT0_t14,time2_t14); +markAll_t14(:,1)=procInd1_t14;markAll_t14(:,2:5)=mark0_t14(rawInd0_t14(rawInd0_t14~=0),:); +mdel=20; ms=min(min(markAll_t14(:,2:5))):mdel:max(max(markAll_t14(:,2:5))); +sxker=2*xdel; smker = mdel; T=size(postimevec,1); +occ=normpdf(xbins'*ones(1,T),ones(length(xbins),1)*poslin2,sxker)*ones(T,length(ms)); +%occ: columns are identical; occupancy based on position; denominator +Xnum_t14=normpdf(xbins'*ones(1,length(spikeT0_t14)),ones(length(xbins),1)*poslin2(procInd1_t14),sxker); +%Xnum: Gaussian kernel estimators for position +Lint_t14=sum(Xnum_t14,2)./occ(:,1)./dt; %integral +Lint_t14=Lint_t14./sum(Lint_t14); + +load('/mnt/vortex_data/mkarlsso/backup/bond/bond04/17-108/bond04-17_params.mat'); +%ind_t17=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)&filedata.params(:,8)>10&(filedata.params(:,2)>100|filedata.params(:,3)>100|filedata.params(:,4)>100|filedata.params(:,5)>100)); +ind_t17=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)); +time_t17=filedata.params(ind_t17,1); +mark0_t17=[filedata.params(ind_t17,2) filedata.params(ind_t17,3) filedata.params(ind_t17,4) filedata.params(ind_t17,5)]; +time2_t17=time_t17/10000; +spikeT0_t17=time2_t17; +[procInd0_t17,procInd1_t17]=histc(spikeT0_t17,postimevec); +procInd_t17=find(procInd0_t17); +spikeT_t17=postimevec(procInd_t17); +spike_t17=procInd0_t17'; +[~,rawInd0_t17]=histc(spikeT0_t17,time2_t17); +markAll_t17(:,1)=procInd1_t17;markAll_t17(:,2:5)=mark0_t17(rawInd0_t17(rawInd0_t17~=0),:); +mdel=20; ms=min(min(markAll_t17(:,2:5))):mdel:max(max(markAll_t17(:,2:5))); +sxker=2*xdel; smker = mdel; T=size(postimevec,1); +occ=normpdf(xbins'*ones(1,T),ones(length(xbins),1)*poslin2,sxker)*ones(T,length(ms)); +%occ: columns are identical; occupancy based on position; denominator +Xnum_t17=normpdf(xbins'*ones(1,length(spikeT0_t17)),ones(length(xbins),1)*poslin2(procInd1_t17),sxker); +%Xnum: Gaussian kernel estimators for position +Lint_t17=sum(Xnum_t17,2)./occ(:,1)./dt; %integral +Lint_t17=Lint_t17./sum(Lint_t17); + +load('/mnt/vortex_data/mkarlsso/backup/bond/bond04/18-138/bond04-18_params.mat'); +%ind_t18=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)&filedata.params(:,8)>10&(filedata.params(:,2)>100|filedata.params(:,3)>100|filedata.params(:,4)>100|filedata.params(:,5)>100)); +ind_t18=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)); +time_t18=filedata.params(ind_t18,1); +mark0_t18=[filedata.params(ind_t18,2) filedata.params(ind_t18,3) filedata.params(ind_t18,4) filedata.params(ind_t18,5)]; +time2_t18=time_t18/10000; +spikeT0_t18=time2_t18; +[procInd0_t18,procInd1_t18]=histc(spikeT0_t18,postimevec); +procInd_t18=find(procInd0_t18); +spikeT_t18=postimevec(procInd_t18); +spike_t18=procInd0_t18'; +[~,rawInd0_t18]=histc(spikeT0_t18,time2_t18); +markAll_t18(:,1)=procInd1_t18;markAll_t18(:,2:5)=mark0_t18(rawInd0_t18(rawInd0_t18~=0),:); +mdel=20; ms=min(min(markAll_t18(:,2:5))):mdel:max(max(markAll_t18(:,2:5))); +sxker=2*xdel; smker = mdel; T=size(postimevec,1); +occ=normpdf(xbins'*ones(1,T),ones(length(xbins),1)*poslin2,sxker)*ones(T,length(ms)); +%occ: columns are identical; occupancy based on position; denominator +Xnum_t18=normpdf(xbins'*ones(1,length(spikeT0_t18)),ones(length(xbins),1)*poslin2(procInd1_t18),sxker); +%Xnum: Gaussian kernel estimators for position +Lint_t18=sum(Xnum_t18,2)./occ(:,1)./dt; %integral +Lint_t18=Lint_t18./sum(Lint_t18); + +load('/mnt/vortex_data/mkarlsso/backup/bond/bond04/19-130/bond04-19_params.mat'); +%ind_t19=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)&filedata.params(:,8)>10&(filedata.params(:,2)>100|filedata.params(:,3)>100|filedata.params(:,4)>100|filedata.params(:,5)>100)); +ind_t19=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)); +time_t19=filedata.params(ind_t19,1); +mark0_t19=[filedata.params(ind_t19,2) filedata.params(ind_t19,3) filedata.params(ind_t19,4) filedata.params(ind_t19,5)]; +time2_t19=time_t19/10000; +spikeT0_t19=time2_t19; +[procInd0_t19,procInd1_t19]=histc(spikeT0_t19,postimevec); +procInd_t19=find(procInd0_t19); +spikeT_t19=postimevec(procInd_t19); +spike_t19=procInd0_t19'; +[~,rawInd0_t19]=histc(spikeT0_t19,time2_t19); +markAll_t19(:,1)=procInd1_t19;markAll_t19(:,2:5)=mark0_t19(rawInd0_t19(rawInd0_t19~=0),:); +mdel=20; ms=min(min(markAll_t19(:,2:5))):mdel:max(max(markAll_t19(:,2:5))); +sxker=2*xdel; smker = mdel; T=size(postimevec,1); +occ=normpdf(xbins'*ones(1,T),ones(length(xbins),1)*poslin2,sxker)*ones(T,length(ms)); +%occ: columns are identical; occupancy based on position; denominator +Xnum_t19=normpdf(xbins'*ones(1,length(spikeT0_t19)),ones(length(xbins),1)*poslin2(procInd1_t19),sxker); +%Xnum: Gaussian kernel estimators for position +Lint_t19=sum(Xnum_t19,2)./occ(:,1)./dt; %integral +Lint_t19=Lint_t19./sum(Lint_t19); + +load('/mnt/vortex_data/mkarlsso/backup/bond/bond04/20-132/bond04-20_params.mat'); +%ind_t20=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)&filedata.params(:,8)>10&(filedata.params(:,2)>100|filedata.params(:,3)>100|filedata.params(:,4)>100|filedata.params(:,5)>100)); +ind_t20=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)); +time_t20=filedata.params(ind_t20,1); +mark0_t20=[filedata.params(ind_t20,2) filedata.params(ind_t20,3) filedata.params(ind_t20,4) filedata.params(ind_t20,5)]; +time2_t20=time_t20/10000; +spikeT0_t20=time2_t20; +[procInd0_t20,procInd1_t20]=histc(spikeT0_t20,postimevec); +procInd_t20=find(procInd0_t20); +spikeT_t20=postimevec(procInd_t20); +spike_t20=procInd0_t20'; +[~,rawInd0_t20]=histc(spikeT0_t20,time2_t20); +markAll_t20(:,1)=procInd1_t20;markAll_t20(:,2:5)=mark0_t20(rawInd0_t20(rawInd0_t20~=0),:); +mdel=20; ms=min(min(markAll_t20(:,2:5))):mdel:max(max(markAll_t20(:,2:5))); +sxker=2*xdel; smker = mdel; T=size(postimevec,1); +occ=normpdf(xbins'*ones(1,T),ones(length(xbins),1)*poslin2,sxker)*ones(T,length(ms)); +%occ: columns are identical; occupancy based on position; denominator +Xnum_t20=normpdf(xbins'*ones(1,length(spikeT0_t20)),ones(length(xbins),1)*poslin2(procInd1_t20),sxker); +%Xnum: Gaussian kernel estimators for position +Lint_t20=sum(Xnum_t20,2)./occ(:,1)./dt; %integral +Lint_t20=Lint_t20./sum(Lint_t20); + +load('/mnt/vortex_data/mkarlsso/backup/bond/bond04/22-171/bond04-22_params.mat'); +%ind_t22=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)&filedata.params(:,8)>10&(filedata.params(:,2)>100|filedata.params(:,3)>100|filedata.params(:,4)>100|filedata.params(:,5)>100)); +ind_t22=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)); +time_t22=filedata.params(ind_t22,1); +mark0_t22=[filedata.params(ind_t22,2) filedata.params(ind_t22,3) filedata.params(ind_t22,4) filedata.params(ind_t22,5)]; +time2_t22=time_t22/10000; +spikeT0_t22=time2_t22; +[procInd0_t22,procInd1_t22]=histc(spikeT0_t22,postimevec); +procInd_t22=find(procInd0_t22); +spikeT_t22=postimevec(procInd_t22); +spike_t22=procInd0_t22'; +[~,rawInd0_t22]=histc(spikeT0_t22,time2_t22); +markAll_t22(:,1)=procInd1_t22;markAll_t22(:,2:5)=mark0_t22(rawInd0_t22(rawInd0_t22~=0),:); +mdel=20; ms=min(min(markAll_t22(:,2:5))):mdel:max(max(markAll_t22(:,2:5))); +sxker=2*xdel; smker = mdel; T=size(postimevec,1); +occ=normpdf(xbins'*ones(1,T),ones(length(xbins),1)*poslin2,sxker)*ones(T,length(ms)); +%occ: columns are identical; occupancy based on position; denominator +Xnum_t22=normpdf(xbins'*ones(1,length(spikeT0_t22)),ones(length(xbins),1)*poslin2(procInd1_t22),sxker); +%Xnum: Gaussian kernel estimators for position +Lint_t22=sum(Xnum_t22,2)./occ(:,1)./dt; %integral +Lint_t22=Lint_t22./sum(Lint_t22); + +load('/mnt/vortex_data/mkarlsso/backup/bond/bond04/23-152/bond04-23_params.mat'); +%ind_t23=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)&filedata.params(:,8)>10&(filedata.params(:,2)>100|filedata.params(:,3)>100|filedata.params(:,4)>100|filedata.params(:,5)>100)); +ind_t23=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)); +time_t23=filedata.params(ind_t23,1); +mark0_t23=[filedata.params(ind_t23,2) filedata.params(ind_t23,3) filedata.params(ind_t23,4) filedata.params(ind_t23,5)]; +time2_t23=time_t23/10000; +spikeT0_t23=time2_t23; +[procInd0_t23,procInd1_t23]=histc(spikeT0_t23,postimevec); +procInd_t23=find(procInd0_t23); +spikeT_t23=postimevec(procInd_t23); +spike_t23=procInd0_t23'; +[~,rawInd0_t23]=histc(spikeT0_t23,time2_t23); +markAll_t23(:,1)=procInd1_t23;markAll_t23(:,2:5)=mark0_t23(rawInd0_t23(rawInd0_t23~=0),:); +mdel=20; ms=min(min(markAll_t23(:,2:5))):mdel:max(max(markAll_t23(:,2:5))); +sxker=2*xdel; smker = mdel; T=size(postimevec,1); +occ=normpdf(xbins'*ones(1,T),ones(length(xbins),1)*poslin2,sxker)*ones(T,length(ms)); +%occ: columns are identical; occupancy based on position; denominator +Xnum_t23=normpdf(xbins'*ones(1,length(spikeT0_t23)),ones(length(xbins),1)*poslin2(procInd1_t23),sxker); +%Xnum: Gaussian kernel estimators for position +Lint_t23=sum(Xnum_t23,2)./occ(:,1)./dt; %integral +Lint_t23=Lint_t23./sum(Lint_t23); + +load('/mnt/vortex_data/mkarlsso/backup/bond/bond04/27-158/bond04-27_params.mat'); +%ind_t27=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)&filedata.params(:,8)>10&(filedata.params(:,2)>100|filedata.params(:,3)>100|filedata.params(:,4)>100|filedata.params(:,5)>100)); +ind_t27=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)); +time_t27=filedata.params(ind_t27,1); +mark0_t27=[filedata.params(ind_t27,2) filedata.params(ind_t27,3) filedata.params(ind_t27,4) filedata.params(ind_t27,5)]; +time2_t27=time_t27/10000; +spikeT0_t27=time2_t27; +[procInd0_t27,procInd1_t27]=histc(spikeT0_t27,postimevec); +procInd_t27=find(procInd0_t27); +spikeT_t27=postimevec(procInd_t27); +spike_t27=procInd0_t27'; +[~,rawInd0_t27]=histc(spikeT0_t27,time2_t27); +markAll_t27(:,1)=procInd1_t27;markAll_t27(:,2:5)=mark0_t27(rawInd0_t27(rawInd0_t27~=0),:); +mdel=20; ms=min(min(markAll_t27(:,2:5))):mdel:max(max(markAll_t27(:,2:5))); +sxker=2*xdel; smker = mdel; T=size(postimevec,1); +occ=normpdf(xbins'*ones(1,T),ones(length(xbins),1)*poslin2,sxker)*ones(T,length(ms)); +%occ: columns are identical; occupancy based on position; denominator +Xnum_t27=normpdf(xbins'*ones(1,length(spikeT0_t27)),ones(length(xbins),1)*poslin2(procInd1_t27),sxker); +%Xnum: Gaussian kernel estimators for position +Lint_t27=sum(Xnum_t27,2)./occ(:,1)./dt; %integral +Lint_t27=Lint_t27./sum(Lint_t27); + +load('/mnt/vortex_data/mkarlsso/backup/bond/bond04/29-118/bond04-29_params.mat'); +%ind_t29=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)&filedata.params(:,8)>10&(filedata.params(:,2)>100|filedata.params(:,3)>100|filedata.params(:,4)>100|filedata.params(:,5)>100)); +ind_t29=find(filedata.params(:,1)/10000>=postimevec(1)&filedata.params(:,1)/10000<=postimevec(end)); +time_t29=filedata.params(ind_t29,1); +mark0_t29=[filedata.params(ind_t29,2) filedata.params(ind_t29,3) filedata.params(ind_t29,4) filedata.params(ind_t29,5)]; +time2_t29=time_t29/10000; +spikeT0_t29=time2_t29; +[procInd0_t29,procInd1_t29]=histc(spikeT0_t29,postimevec); +procInd_t29=find(procInd0_t29); +spikeT_t29=postimevec(procInd_t29); +spike_t29=procInd0_t29'; +[~,rawInd0_t29]=histc(spikeT0_t29,time2_t29); +markAll_t29(:,1)=procInd1_t29;markAll_t29(:,2:5)=mark0_t29(rawInd0_t29(rawInd0_t29~=0),:); +mdel=20; ms=min(min(markAll_t29(:,2:5))):mdel:max(max(markAll_t29(:,2:5))); +sxker=2*xdel; smker = mdel; T=size(postimevec,1); +occ=normpdf(xbins'*ones(1,T),ones(length(xbins),1)*poslin2,sxker)*ones(T,length(ms)); +%occ: columns are identical; occupancy based on position; denominator +Xnum_t29=normpdf(xbins'*ones(1,length(spikeT0_t29)),ones(length(xbins),1)*poslin2(procInd1_t29),sxker); +%Xnum: Gaussian kernel estimators for position +Lint_t29=sum(Xnum_t29,2)./occ(:,1)./dt; %integral +Lint_t29=Lint_t29./sum(Lint_t29); + + + +%% +time0=[time_t1;time_t2;time_t4;time_t5;time_t7;time_t10;time_t11;time_t12;time_t13;time_t14;time_t17;time_t18;time_t19;time_t20;time_t22;time_t23;time_t27;time_t29]; +[time,timeInd]=sort(time0); +mark0=[mark0_t1;mark0_t2;mark0_t4;mark0_t5;mark0_t7;mark0_t10;mark0_t11;mark0_t12;mark0_t13;mark0_t14;mark0_t17;mark0_t18;mark0_t19;mark0_t20;mark0_t22;mark0_t23;mark0_t27;mark0_t29]; +mark0=mark0(timeInd,:); +procInd1=[procInd1_t1;procInd1_t2;procInd1_t4;procInd1_t5;procInd1_t7;procInd1_t10;procInd1_t11;procInd1_t12;procInd1_t13;procInd1_t14;procInd1_t17;procInd1_t18;procInd1_t19;procInd1_t20;procInd1_t22;procInd1_t23;procInd1_t27;procInd1_t29]; +procInd1=procInd1(timeInd,:); + +len_t1=length(time_t1);len_t2=length(time_t2);len_t4=length(time_t4);len_t5=length(time_t5); +len_t7=length(time_t7);len_t11=length(time_t11);len_t10=length(time_t10);len_t12=length(time_t12); +len_t13=length(time_t13);len_t14=length(time_t14);len_t17=length(time_t17);len_t18=length(time_t18); +len_t19=length(time_t19);len_t20=length(time_t20);len_t22=length(time_t22);len_t23=length(time_t23); +len_t27=length(time_t27);len_t29=length(time_t29); + +tet_ind=zeros(length(time),5); +%indicator matrix +%row: time point; column: which tetrode spikes +for i=1:length(time0) + if timeInd(i)>=1 && timeInd(i)<=len_t1 + tet_ind(i,1)=1; %tet1 + elseif timeInd(i)>=len_t1+1 && timeInd(i)<=len_t1+len_t2 + tet_ind(i,2)=1; %tet2 + elseif timeInd(i)>=len_t1+len_t2+1 && timeInd(i)<=len_t1+len_t2+len_t4 + tet_ind(i,3)=1; %tet4 + elseif timeInd(i)>=len_t1+len_t2+len_t4+1 && timeInd(i)<=len_t1+len_t2+len_t4+len_t5 + tet_ind(i,4)=1; %tet5 + elseif timeInd(i)>=len_t1+len_t2+len_t4+len_t5+1 && timeInd(i)<=len_t1+len_t2+len_t4+len_t5+len_t7 + tet_ind(i,5)=1; %tet7 + elseif timeInd(i)>=len_t1+len_t2+len_t4+len_t5+len_t7+1 && timeInd(i)<=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10 + tet_ind(i,6)=1; %tet10 + elseif timeInd(i)>=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+1 && timeInd(i)<=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11 + tet_ind(i,7)=1; %tet11 + elseif timeInd(i)>=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+1 && timeInd(i)<=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12 + tet_ind(i,8)=1; %tet12 + elseif timeInd(i)>=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12+1 && timeInd(i)<=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12+len_t13 + tet_ind(i,9)=1; %tet13 + elseif timeInd(i)>=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12+len_t13+1 && timeInd(i)<=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12+len_t13+len_t14 + tet_ind(i,10)=1; %tet14 + elseif timeInd(i)>=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12+len_t13+len_t14+1 && timeInd(i)<=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12+len_t13+len_t14+len_t17 + tet_ind(i,11)=1; %tet17 + elseif timeInd(i)>=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12+len_t13+len_t14+len_t17+1 && timeInd(i)<=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12+len_t13+len_t14+len_t17+len_t18 + tet_ind(i,12)=1; %tet18 + elseif timeInd(i)>=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12+len_t13+len_t14+len_t17+len_t18+1 && timeInd(i)<=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12+len_t13+len_t14+len_t17+len_t18+len_t19 + tet_ind(i,13)=1; %tet19 + elseif timeInd(i)>=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12+len_t13+len_t14+len_t17+len_t18+len_t19+1 && timeInd(i)<=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12+len_t13+len_t14+len_t17+len_t18+len_t19+len_t20 + tet_ind(i,14)=1; %tet20 + elseif timeInd(i)>=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12+len_t13+len_t14+len_t17+len_t18+len_t19+len_t20+1 && timeInd(i)<=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12+len_t13+len_t14+len_t17+len_t18+len_t19+len_t20+len_t22 + tet_ind(i,15)=1; %tet22 + elseif timeInd(i)>=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12+len_t13+len_t14+len_t17+len_t18+len_t19+len_t20+len_t22+1 && timeInd(i)<=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12+len_t13+len_t14+len_t17+len_t18+len_t19+len_t20+len_t22+len_t23 + tet_ind(i,16)=1; %tet23 + elseif timeInd(i)>=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12+len_t13+len_t14+len_t17+len_t18+len_t19+len_t20+len_t22+len_t23+1 && timeInd(i)<=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12+len_t13+len_t14+len_t17+len_t18+len_t19+len_t20+len_t22+len_t23+len_t27 + tet_ind(i,17)=1; %tet27 + elseif timeInd(i)>=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12+len_t13+len_t14+len_t17+len_t18+len_t19+len_t20+len_t22+len_t23+len_t27+1 && timeInd(i)<=len_t1+len_t2+len_t4+len_t5+len_t7+len_t10+len_t11+len_t12+len_t13+len_t14+len_t17+len_t18+len_t19+len_t20+len_t22+len_t23+len_t27+len_t29 + tet_ind(i,18)=1; %tet29 + end +end + +tet_sum=tet_ind.*cumsum(tet_ind,1); %row: time point; column: index of spike per tetrode + + +%% +mdel=20; ms=100:mdel:max(mark0(:)); +sxker=2*xdel; smker = mdel; T=size(postimevec,1); +occ=normpdf(xbins'*ones(1,T),ones(length(xbins),1)*poslin2,sxker)*ones(T,length(ms)); +%occ: columns are identical; occupancy based on position; denominator +Xnum=normpdf(xbins'*ones(1,length(time0)),ones(length(xbins),1)*poslin2(procInd1),sxker); +%Xnum: Gaussian kernel estimators for position +Lint=sum(Xnum,2)./occ(:,1)./dt; %integral +Lint=Lint./sum(Lint); +%Lint: conditional intensity function for the unmarked case + + + +%%%% Ripples %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +ripoption = 2; % == 1 for old ripples, == 2 for ripplescons + + +%% +if ripoption == 1 + + load('/opt/data13/kkay/Bon/bonripples04.mat'); + %ind_t=[1 2 4 5 7 10 11 12 13 14 17 18 19 20 22 23 27 29]; + ex=4;ep=2; + + ti=round(postimevec(1)*1000):1:round(postimevec(end)*1000); % 1-ms time vector + %position time; from 2461011 to 3404979; length(ti)~~max(spike_tim) + + % Make ind_ripple + ind_ripple=zeros(size(ti,2),size(ind_t,2)); %