% This script analyzes preprocessed data generated by the CellExplorer
% framework (Petersen et al., 2021, Neuron). It can be run on the dataset
% containing recordings of dorsal CA1 activity during a continuous
% alternation task. The dataset can be found in the Buzsaki Lab Databank:
% https://buzsakilab.com/wp/projects/entry/67125/. All access paths are
% local (specific to the computer on which the analysis was run) thus,
% pathnames have to be changed before running the script.
% Dependencies: timespikeshuffler.m; placefieldidentifierr.m;
% PFidentifierSimple.m; placeTMdecoder.m; placeTMdecoderR.m; 

old = 'no';
% sessnum = input('Give a session number: '); % don't give "1"
% sessnum = sessnums(pm);
sessnum = 24;
disp(sessnum)
% lfp_chan = input('LFP channel for spike-phase coupling: ');
sessionname = {'ham11_34-36',... %1
    'ham12_82-84','ham12_85_88-89','ham12_92-95','ham12_96-98',... %5
    'ham12_102-104','ham12_110-113','ham12_119-121','ham12_125-127','ham12_139-140_142',... %10
    'ham12_153-155',... %11
    'ham21_18-22',... %12
    'ham21_27-29','ham21_33-35','ham21_92-94',... %15
    'ham21_95-97','ham21_102-105','ham21_106-108',... %18
    'ham21_109-111','ham21_113-115','ham21_117-119','ham21_123-125',... %22
    'ham25_35-37','ham25_44-46','ham25_47-49','ham25_51-53','ham25_58-60',... %27
    'ham31_24-27',... %28
    'ham31_28-30',... %29
    'ham31_31-34',... %30
    'ham31_42-45',... %31
    'ham34_32-34','ham34_37-39','ham26_33-35','ham21_98-100','ham25_38-40','ham25_41-43','ham34_48-52','ham34_40-43',... %39
    'ham34_120-122','ham34_100-102','ham34_123-125','ham34_126-128','ham34_132-134','ham34_135-137','ham34_138-140',... %46
    'ham34_141-143','ham34_144-146','ham34_147-149','ham34_150-152','ham34_153-155',... %51
    'ham34_74-76','ham20_16-18','ham34_53-56','ham34_68-70','ham34_81-83','ham34_87-90','ham34_91-93','ham34_97-99',... %59
    'ham34_156-158','ham34_162-165','ham34_166-170','ham34_171-174','ham34_129-131','ham34_159-161'}; %65

%     'ham34_77-80','ham34_81-83','ham34_84-86','ham34_87-90','ham34_91-93','ham34_97-99','ham34_171-174',...

% sns = input('Save (yes) or not to save (no)?','s');
sns = 'yes';
% shgen = input('Generate place cell surrogate data? /yes OR no/?','s');
shgen = 'no';
disp(['Loading data...',sessionname{sessnum}])
% datadir = ['E:\HAMalysis\Peter_preprocessed\' sessionname{sessnum}(1:5) '\' sessionname{sessnum} '_amp'];
datadir = ['D:\vargav\Hamalysis\splitters\datasets_allnostim\' sessionname{sessnum} '_amp'];
datadirBehav = 'D:\vargav\Hamalysis\splitters\full_maze_behavioral_files\';
resdir = 'D:\vargav\Hamalysis\splitters\PFA_results';
mkdir(resdir,sessionname{sessnum})
resdirf = [resdir '\' sessionname{sessnum}];
respath = fullfile(resdir,sessionname{sessnum});
%SUPER IMPORTANT: ONLY THE CONTROL SESSIONS WERE CHECKED AND VERIFIED!!!!
pc_ID_full = {[1 5 8 9 10 12 13 14 15 16 18 21 22 23 24 25],... %1
    [1 2 3 4 5 7 9 10],[1 2 4 5 7 8 9 10 11 12 14 15 18 19 20],[2 3 6 7 8 9 10 11 12 14 15],[1 2 3 4 5 6 8 9 10 11 14 15 16 17 18 19 20],... %5
    [2 3 4 5 9 10 11 13 14 15],[2 3 5 7 8 9 11 16 17 18],[1 3 5 6 7 8 9 10 11 12 13 14 15],[3 4 5 6 7 8 9 10 12 13 14 15],[1 3 4 5 6 7 10 11 12 13 14 15],... %10
    [2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 20 25 27 28 30 31 32 33 34 35 36 37 38 39 40],... %11
    [2 3 4 5 8 9 10 11 12 13 14 15 16 17 18 19 20],... %12
    [2 3 4 5 6 7 8 9 10 13 14 15 26 28 29 30 31 32 33 34],[2 3 4 5 14 15 16 17 18 19 20 21 22],[1 5 6 8 10 11 13 14 15 17 20 22 23 24 25],... %15
    [1 2 3 4 5 8 9 12 13 16 17 18 19 20],[2 3 5 6 7 8 9 10 11 12 13 14 16 17 18 19 22 23 24 25 26 27 28 30],[3 5 6 7 8 11 14 15],... %18
    [1 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 26 27],[1 2 3 4 5 7 8 9 10 12 14 17 18 19 20],[2 3 4 5 6 7 8 10 11 12 13 14 15 17 19 20],[1 2 3 4 5 6 8 11 12 13 14 15],... %22
    [3 4 5 6 7 8 9 10 11 15 16 17 18],[2 3 4 5 6 7 8 9 11 12 13 14 15],[4 5 6 7 8 9 10 11 12 13 14 15 17 18],[2 3 4 5 6 7 8 11 12 13 14 15 16],[1 2 3 4 5 6 10 11 12 13 14 15 16],... %27
    [3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56],... %28
    [4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57],... %29
    [2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52],... %30
    [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 25 26 27 28 29 30 31 32 33 34 35 36],... %31
    [1 4 5 6 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27],[1 4 5 6 7 8 10 11 12 13],[1 2 3 4 5 6 7 8 9 10 11 12 13 14 16 17],[1 4 5 6 7 8 9 10],[1 3 4 5 6 7 12 13],[3 4 5 6 7 8 11 12 13],... %37
    [1 2 3 4 7 8 9 10],[1 2 3 4 6 7 8 9 10 11 12],[2 3 5 7 8 9 10 11 12 13 14 15 16 19 20],... %40
    [2 3 4 5 6 7 8 9 10 11 12 13 16 17 18],[2 3 4 5 6 7 9 10 12 13 14 15 16 18 20 21 22 23 24 25 26],... %42
    [1 2 3 4 9 10 11 12 13 14 15 16 17 18 19 20],[1 2 3 4 5 6 8 9 10 11 12 13 15 16 17],... %44
    [1 2 3 4 5 6 7 9 13 15 16 17],... %45
    [1 2 3 4 5 6 8 12 14 17 18 19 20 21 22 23],... %46
    [3 4 5 6 7 8 10 11 12 13 14 16 17 21 22 23 24],... %47
    [1 2 5 6 7 8 9 10 12 13 14 16 17 18 19 20 24 25 26 27 28 29 30 31],... %48
    [4 5 6 7 8 9 10 11 13 14 15 16 19 20 21 22 23],[2 3 5 6 7 8 9 10 11 15 16 18],... %50
    [3 4 5 6 7 8 9 10 11 12 13 14 16],[5 6 7 8 9 10 11 12 13 14 16],... %52
    [1 2 3 4 6 7 8 9 10],[2 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28],[1 2 4 5 6 7],... %55
    [3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21],[3 4 5 6 7 8 9 10 11 12 13 14 15 16],[3 4 5 6 7 8 9 10 11 12 13 14 15 16 17],... %58
    [4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 24 25 26],... %59
    [1 2 3 4 5 6 7 11 12 13 14 16 17 18 19 20],[2 3 4 5 6 7 9 11 12 13 14 17 18 19 20 21 22 23 24],... %61
    [5 6 7 8 9 11 12 13 14 15 16 17 18 19 20],[4 5 6 7 8 9 11 12 13 15 16 17 18 19 20 21 22 23],[2 3 4 5 6 7 10 13 14 15 16],... %64
    [2 3 4 5 6 7 9 10 11 12 13 14 15 16 17 21 22 23 24 25 26]}; %65
IN_ID_full = {[],[],[],[],[],[],[],[],[],[],[],... %11
    [1 6 7 21],...%C-12
    [],[],[],[],[],[],...%18
    [2 24 25],...%C-19
    [],[],[],...%22
    [1 2 12 13 14],[1 10],[1 2 3 16 19],[1 9 10],[7 8 9],...%C-27
    [1 2 38],[1 2 3 35],[1 33 34],[21 22 23 24],...%C-31
    [2 3 7 8 9 10],[],15,[2 3],[2 8 9 10 11],[1 2 9 10],[5 6],[],...%C-39
    [1 4 6 17 18],[1 14 15],[1 8 11 17 19],[5 6 7 8 21],[7 14],[5 10 11 12],[7 9 10 11 13 15 16],...%C-46
    [1 2 11 15 18 19 20 25],[3 4 11 15 21 22 23],[1 2 3 12 17 18],[1 4 12 13 14 17],[1 2 15],...%C-51
    [],5,[],3,[],[],[],[],...%C-59
    [8 9 10 15],[1 8 10 15 16 25],[1 2 3 4 10],[1 2 3 10 14],[1 8 9 11 12],[1 8 18 19 20]}; %65
pc_ID = pc_ID_full{sessnum}; %principal cells
IN_ID = IN_ID_full{sessnum}; %interneurons

%%
%Section 2: Extraction and definition of necessary variables
disp('Preaparations...')

load([datadir '\' sessionname{sessnum} '_amp.spikes.cellinfo.mat']);
load([datadirBehav '\' sessionname{sessnum} '_amp.trials.behavior.mat']);
load([datadir '\' sessionname{sessnum} '_amp.InstantaneousTheta.lfp.mat']);
load([datadirBehav sessionname{sessnum} '_amp.animal.behavior.mat']);
load([datadir '\' sessionname{sessnum} '_amp.spikes.linpos.recalculated.mat']);
%     load([datadir '\' sessionname{sessnum} '_amp.optogenetics.manipulation.mat']);

%definition of zone boundaries
stem = [20,62]; %changed from 16 to 20
transitionZ = [0,20]; %changed from 16 to 20
leftrim = [62,134];
rightrim = [135,205];
rewardpos = [119 193];
rewardrange = [118 192;120 194];

sro = 20000;
fps = 100; %but check for each session (can be 120)!
alltime = animal.time; %animal times on the maze

linpos = spike_linpos.PC_linpos_recalculated; %linearized position of PC spikes (not full coverage, discarded)
fullpos = spikes.pos(pc_ID); %circular position of PC spikes
INpos = spike_linpos.IN_linpos_recalculated; %linear position of IN spikes
trialnum = spike_linpos.PC_trials_recalculated; %trial numbers assigned to PC spikes
trialnumIN = spike_linpos.IN_trials_recalculated; %trial numbers assigned to IN spikes


otimes = spikes.times(pc_ID); %spike times of PCs
otimeIN = spikes.times(IN_ID); %spike times of INs
times = cell(1,length(otimes));
timeIN = cell(1,length(otimeIN));
numspk = zeros(1,length(otimes));
for ti = 1 : length(otimes)
    times{ti} = otimes{ti}(otimes{ti} >= alltime(1) & otimes{ti} <= alltime(end));
    linpos{ti} = linpos{ti}(otimes{ti} >= alltime(1) & otimes{ti} <= alltime(end));
    numspk(ti) = size(times{ti},1);
end
times = times(numspk >= 100);
pc_ID(numspk < 100) = [];
linpos = linpos(numspk >= 100);
fullpos = fullpos(numspk >= 100);
otimes = otimes(numspk >= 100);
trialnum = trialnum(numspk >= 100);

for fp1 = 1 : length(fullpos)
    fullpos{fp1} = fullpos{fp1}(:,otimes{fp1} >= alltime(1) & otimes{fp1} <= alltime(end));
end

numINspk = zeros(1,length(otimeIN));
for tj = 1 : length(otimeIN)
    timeIN{tj} = otimeIN{tj}(otimeIN{tj} >= alltime(1) & otimeIN{tj} <= alltime(end));
    INpos{tj} = INpos{tj}(otimeIN{tj} >= alltime(1) & otimeIN{tj} <= alltime(end));
    numINspk(tj) = size(timeIN{tj},1);
end
timeIN = timeIN(numINspk >= 1000);
IN_ID(numINspk < 1000) = [];
INpos = INpos(numINspk >= 1000);
otimeIN = otimeIN(numINspk >= 1000);
trialnumIN = trialnumIN(numINspk >= 1000);

spikephase = spikes.theta_phase(pc_ID); %phase of PC spikes
spikephaseIN = spikes.theta_phase(IN_ID); %phase of PC spikes

for tm = 1 : length(times)
    trialnum{tm} = trialnum{tm}(otimes{tm} >= alltime(1) & otimes{tm} <= alltime(end));
    ttn = trialnum{tm};
    nanfinder = zeros(1,size(ttn,1));
    nanfinder(isnan(ttn)) = 1;
    dnf = diff(nanfinder);
    bgn = find(dnf == 1) + 1;
    fin = find(dnf == -1);
    if nanfinder(1) == 1
        bgn = [1,bgn];
    end
    if nanfinder(end) == 1
        fin = [fin,length(nanfinder)];
    end
    for tmm = 1 : length(bgn)
        if bgn(tmm) == 1
            ttn(bgn(tmm) : fin(tmm)) = 1;
        end
        if fin(tmm) < size(ttn,1) && tmm > 1
            ttn(bgn(tmm) : fin(tmm)) = ttn(bgn(tmm) - 1);
        end
        if fin(tmm) == size(ttn,1) && bgn(tmm) ~= 1
            ttn(bgn(tmm) : fin(tmm)) = ttn(bgn(tmm) - 1);
        end
    end
    trialnum{tm} = ttn;
end

for tn = 1 : length(timeIN)
    trialnumIN{tn} = trialnumIN{tn}(otimeIN{tn} >= alltime(1) & otimeIN{tn} <= alltime(end));
    ttIN = trialnumIN{tn};
    naninf = zeros(1,size(ttIN,1));
    naninf(isnan(ttIN)) = 1;
    dinf = diff(naninf);
    bgin = find(dinf == 1) + 1;
    fiin = find(dinf == -1);
    if naninf(1) == 1
        bgin = [1,bgin];
    end
    if naninf(end) == 1
        fiin = [fiin,length(naninf)];
    end
    for tnn = 1 : length(bgin)
        if bgin(tnn) == 1
            ttIN(bgin(tnn) : fiin(tnn)) = 1;
        end
        if fiin(tnn) < size(ttIN,1) && tnn > 1
            ttIN(bgin(tnn) : fiin(tnn)) = ttIN(bgin(tnn) - 1);
        end
        if fiin(tnn) == size(ttIN,1)
            ttIN(bgin(tnn) : fiin(tnn)) = ttIN(bgin(tnn) - 1);
        end
    end
    trialnumIN{tn} = ttIN;
end

unitarm = cell(1,length(linpos)); %cell array for stem spike linear positions
unitarm_full = cell(1,length(linpos)); %cell array for stem spike circular positions
unitPTP_arm = cell(1,length(linpos)); %cell array for spike times in stem
unitDP_arm = cell(1,length(linpos)); %cell array for ISIs in stem
unitarm_phase = cell(1,length(linpos)); %cell array for spike phase in stem

unitranz = cell(1,length(linpos)); %cell array for TZ spike linear positions
unitranz_full = cell(1,length(linpos)); %cell array for TZ spike circular positions
unitPTP_tranz = cell(1,length(linpos)); %cell array for spike times in stem
unitDP_tranz = cell(1,length(linpos)); %cell array for ISIs in stem
unitranz_phase = cell(1,length(linpos)); %cell array for spike phase in stem

unitPTP_full = cell(1,length(linpos)); %cell array for spike times on the maze (full)
unitDP_full = cell(1,length(linpos)); %cell array for ISIs on the maze (full)
unit_fullphase = cell(1,length(linpos)); %cell array for spike phase on the maze (full)

unitranz_fullin = cell(1,length(linpos)); %cell array for IN spike circular positions on the maze
INtfull = cell(1,length(INpos)); %IN spike times on the maze
INfullpos = cell(1,length(INpos)); %IN spike positions on the maze

unitarm_fullin = cell(1,length(linpos)); %cell array for IN spike stem circular positions
INtarm = cell(1,length(INpos)); %IN spike times in stem
INarmpos = cell(1,length(INpos)); %IN spike positions in the stem
INarmphase = cell(1,length(INpos)); %IN spike phases in the stem

INtranz = cell(1,length(INpos)); %IN spike times in TZ
INtranzpos = cell(1,length(INpos)); %IN spike positions in TZ
INtranzphase = cell(1,length(INpos)); %IN spike phases in TZ

X_lim = [min(animal.pos(1,:)) max(animal.pos(1,:))];
Y_lim = [min(animal.pos(2,:)) max(animal.pos(2,:))];
allpos = animal.pos_linearized; %animal's position
speed = animal.speed; %animal's speed

trnm = trials.trials2;
nanfind = zeros(1,size(trnm,2));
nanfind(isnan(trnm)) = 1;
dnt = diff(nanfind);
bgin = find(dnt == 1) + 1;
finn = find(dnt == -1);
if nanfind(1) == 1
    bgin = [1,bgin];
end
if nanfind(end) == 1
    finn = [finn,length(nanfind)];
end
for tmm = 1 : length(bgin)
    if bgin(tmm) == 1
        trnm(bgin(tmm) : finn(tmm)) = 1;
    end
    if finn(tmm) < size(trnm,2) && tmm > 1
        trnm(bgin(tmm) : finn(tmm)) = trnm(bgin(tmm) - 1);
    end
    if finn(tmm) == size(trnm,2)
        trnm(bgin(tmm) : finn(tmm)) = trnm(bgin(tmm) - 1);
    end
end
trialnums = trnm; %number of trials
trials_act = trialnums; %number of trials after removing NaNs
speedact = speed;
allpos_act = allpos; %animal's position w/o NaNs
alltime_act = alltime; %animal time w/o NaNs
alltimepos = [alltime_act;allpos_act]; %time and position merged

trials_stem = trials_act(allpos_act >= 20 & allpos_act < 62);
trials_TZS = trials_act(allpos_act >= 0 & allpos_act < 62);
fullanimpos = [animal.pos;animal.pos_linearized];
% fullanimpos(:,isnan(trialnums)) = [];
fullanimposTZS = fullanimpos(:,fullanimpos(4,:) >= 0 & fullanimpos(4,:) < 62);
fullanimposStem = fullanimpos(:,fullanimpos(4,:) >= 20 & fullanimpos(4,:) < 62);
startime = alltime(1); %time of session start
alltimepos_stem = alltimepos(:,alltimepos(2,:) >= 20 & alltimepos(2,:) < 62); %animal time / position in stem
alltimepos_stem(1,:) = (alltimepos_stem(1,:) - startime) * sro; %start on stem converted to SRo
% stemstimpeaks = (stemstimpeaks - startime) * sro; %stimulation peaks on stem (redundant?)

trialsides = trials.stat; %side codes for trials
diffsides = diff(trialsides);
% trials2drop = 2; %dropping first N trials but just for the population vector analysis!
% trialsides = trialsides(trials2drop + 1 : end);
lefts = find(trialsides == 1); %left trials
rights = find(trialsides == 2); %right trials
LR_prop = length(lefts) / length(rights); %L/R ratio

errortrials = find(diffsides == 0) + 1; %error trials
errorlefts =  lefts; %left errors
errorights = rights; %right errors
for et1 = 1 : length(lefts)
    terl = find(lefts(et1) == errortrials);
    if isempty(terl) == 1
        errorlefts(et1) = 0;
    end
end
for et2 = 1 : length(rights)
    terr = find(rights(et2) == errortrials);
    if isempty(terr) == 1
        errorights(et2) = 0;
    end
end
errorlefts(errorlefts == 0) = []; %cleaned error lefts
errorights(errorights == 0) = []; %cleaned error rights
purelefts = lefts;
purerights = rights;
for et3 = 1 : length(purelefts)
    tet1 = find(purelefts(et3) == errorlefts);
    if isempty(tet1) == 0
        purelefts(et3) = 0;
    end
end
for et4 = 1 : length(purerights)
    tet2 = find(purerights(et4) == errorights);
    if isempty(tet2) == 0
        purerights(et4) = 0;
    end
end
purelefts(purelefts == 0) = []; %cleaned correct lefts
purerights(purerights == 0) = [];%cleaned correct rights
if length(purelefts) > length(purerights)
    purelefts(end) = [];
elseif length(purelefts) < length(purerights)
    purerights(end) = [];
end

checktrialnum = zeros(1,length(linpos));
for cj = 1 : length(linpos)
    checktrialnum(cj) = max(trialnum{cj});
end
constrialnum = max(checktrialnum); % number of trials

animTrialPosCode = zeros(1,length(trials_act));
for ap = 1 : length(trials_act)
    if ~isempty(find(purelefts == trials_act(ap)))
        animTrialPosCode(ap) = 1;
    elseif ~isempty(find(purerights == trials_act(ap)))
        animTrialPosCode(ap) = 2;
    elseif ~isempty(find(errorlefts == trials_act(ap)))
        animTrialPosCode(ap) = 3;
    elseif ~isempty(find(errorights == trials_act(ap)))
        animTrialPosCode(ap) = 4;
    end
end

cleftAnimPos = fullanimpos(1:2,animTrialPosCode == 1);
crightAnimPos = fullanimpos(1:2,animTrialPosCode == 2);
eleftAnimPos = fullanimpos(1:2,animTrialPosCode == 3);
erightAnimPos = fullanimpos(1:2,animTrialPosCode == 4);

% trials_act_trunc = trials_act(trials_act > trials2drop);
% alltimepos_trunc = alltimepos(:,trials_act > trials2drop);
% trialfull = zeros(1,length(trialnums));
% trialfull(isnan(trialnums)) = 1;
% dtf = diff(trialfull);
% tbg = find(dtf == 1);
% tfns = find(dtf == -1);
% trialbegins = tbg + 1;
% trialends = tfns;
% if trialfull(1) == 1
%     trialbegins = [1 trialbegins];
% end
% if trialfull(end) == 1
%     trialends = [tfns length(trialnums)];
% end
% trialfull(trialbegins(1) : trialends(1)) = 1;
% for tf1 = 1 : constrialnum
%     trialfull(trialbegins(tf1 + 1) : trialends(tf1 + 1)) = tf1;
% end
% t2add = trialnums;
% t2add(isnan(t2add)) = 0;
% trialfull = trialfull + t2add;
% trialextstem = trialfull();

PFA_trial_parameters = struct('correct_lefts',purelefts,'correct_rights',purerights,'error_lefts',errorlefts,'error_rights',errorights);
filenameTP = [sessionname{sessnum} '_trial_parameters.mat'];
save(fullfile(resdir,sessionname{sessnum},filenameTP),'PFA_trial_parameters')

%%
%Section 3: Trial separation based on sides

%stem = [20,62];
%transitionZ = [0,20];
%leftrim = [62,134];
%rightrim = [135,205];

alltrials = cell(1,max(trials_act));
alltrialtimes = cell(1,max(trials_act));
alltrialtimestem = cell(1,max(trials_act));
alltrialtimetranz = cell(1,max(trials_act));
trialspeeds = cell(1,max(trials_act));
for b0 = 1 : max(trials_act)
    t0 = find(trials_act == b0);
    alltrials{b0} = allpos_act(t0); %animal position is broken down into trials
    alltrialtimes{b0} = alltime_act(t0); %time on maze is divided based on trials
    alltrialtimestem{b0} = alltrialtimes{b0}(alltrials{b0} >= 20 & alltrials{b0} < 62); %time in stem
    alltrialtimetranz{b0} = alltrialtimes{b0}(alltrials{b0} >= 0 & alltrials{b0} < 20); %time in transition zone
    trialspeeds{b0} = speed(t0); %instantaneous speed in trials
end

fulltrialpos = cell(1,constrialnum);
for bx = 1 : constrialnum
    tx0 = find(trials_TZS == bx);
    fulltrialpos{bx} = fullanimposTZS(:,tx0);
end

%mean speed in trials on the stem
armspeeds = cell(1,max(trials_act));
tranzspeed = cell(1,max(trials_act));
avearmspeeds = zeros(1,max(trials_act));
avetranzspeeds = zeros(1,max(trials_act));
avespeeds = zeros(1,max(trials_act));
for ba1 = 1 : max(trials_act)
    armspeeds{ba1} = trialspeeds{ba1}(alltrials{ba1} >= 20 & alltrials{ba1} < 62);
    tranzspeed{ba1} = trialspeeds{ba1}(alltrials{ba1} >= 0 & alltrials{ba1} < 20);
    avearmspeeds(ba1) = mean(armspeeds{ba1});
    avetranzspeeds(ba1) = mean(tranzspeed{ba1});
    avespeeds(ba1) = mean(trialspeeds{ba1});
end

%animal position in left trials
leftrials = cell(1,length(lefts));
leftfullpos = cell(1,length(lefts));
cfulleftspeed = cell(1,length(lefts));
for b1 = 1 : length(lefts)
    tl = find(trials_act == lefts(b1));
    leftrials{b1} = allpos_act(tl);
    leftfullpos{b1} = fullanimpos(:,tl);
    cfulleftspeed{b1} = speedact(tl);
end

%animal position in right trials
rightrials = cell(1,length(rights));
rightfullpos = cell(1,length(rights));
cfullrightspeed = cell(1,length(rights));
for b2 = 1 : length(rights)
    tr = find(trials_act == rights(b2));
    rightrials{b2} = allpos_act(tr);
    rightfullpos{b2} = fullanimpos(:,tr);
    cfullrightspeed{b2} = speedact(tr);
end

%animal position and speed in correct left trials
cleftrials = cell(1,length(purelefts));
cleftspeed = cell(1,length(purelefts));
cleftarmspeed = cell(1,length(purelefts));
for b1 = 1 : length(purelefts)
    ctpl = find(trials_act == purelefts(b1));
    cleftrials{b1} = fullanimpos(:,ctpl);
    cleftspeed{b1} = speedact(ctpl);
    cleftarmspeed{b1} = cleftspeed{b1}(cleftrials{b1}(4,:) < 62);
end

%animal position and speed in correct right trials
crightrials = cell(1,length(purerights));
crightspeed = cell(1,length(purerights));
crightarmspeed = cell(1,length(purerights));
for b2 = 1 : length(purerights)
    ctpr = find(trials_act == purerights(b2));
    crightrials{b2} = fullanimpos(:,ctpr);
    crightspeed{b2} = speedact(ctpr);
    crightarmspeed{b2} = crightspeed{b2}(crightrials{b2}(4,:) < 62);
end

%stem positions in left trials
leftcentrals = cell(1,length(lefts));
leftarms = cell(1,length(lefts));
for b3 = 1 : length(lefts)
    leftcentrals{b3} = leftrials{b3}(leftrials{b3} < 62);
    leftarms{b3} = leftrials{b3}(leftrials{b3} >= 62 & leftrials{b3} <= 134);
end

%stem positions in right trials
rightcentrals = cell(1,length(rights));
rightarms = cell(1,length(rights));
for b4 = 1 : length(rights)
    rightcentrals{b4} = rightrials{b4}(rightrials{b4} < 62);
    rightarms{b4} = rightrials{b4}(rightrials{b4} > 134 & rightrials{b4} <= 205);
end

%%
%Section 4: distribution of occupancy time in spatial bins
disp('Calculation of occupancy...')

spatbins = 20 : 2 : 62; %stem bins: changed from 16, see placefieldanalyzer_batch
tranzbins = 0 : 2 : 20; %transition zone bins
TZSbins = 0 : 2 : 62;
armbins = 62 : 2 : 190;
leftbins = 62 : 2 : 134;
rightbins = 135 : 2 : 205;
fullbins = 0 : 2 : 205;
fulldistr = zeros(size(alltrials,2),length(fullbins)-1);
alldistr = zeros(size(alltrials,2),length(spatbins)-1);
armdistr = zeros(size(alltrials,2),length(armbins)-1);
tranzdistr = zeros(size(alltrials,2),length(tranzbins)-1);
fullarmposbini = cell(1,constrialnum);
armtzposbini = cell(1,constrialnum);
for dj0 = 1 : size(alltrials,2)
    fulldistr(dj0,:) = histcounts(alltrials{dj0},fullbins);
    alldistr(dj0,:) = histcounts(alltrials{dj0},spatbins);
    armdistr(dj0,:) = histcounts(alltrials{dj0},armbins);
    tranzdistr(dj0,:) = histcounts(alltrials{dj0},tranzbins);
    tobedistr = fulltrialpos{dj0}(4,:);
    [ps1,ps2,ps3] = histcounts(tobedistr,spatbins);
    fullarmposbini{dj0} = ps3;
    [pstz1,pstz2,pstz3] = histcounts(tobedistr,TZSbins);
    armtzposbini{dj0} = pstz3;
end

leftdistr = zeros(length(lefts),length(spatbins)-1);
leftarmdistr = zeros(length(lefts),length(leftbins)-1);
for dj1 = 1 : length(lefts)
    leftdistr(dj1,:) = histcounts(leftcentrals{dj1},spatbins);
    leftarmdistr(dj1,:) = histcounts(leftarms{dj1},leftbins);
end

rightdistr = zeros(length(rights),length(spatbins)-1);
rightarmdistr = zeros(length(rights),length(rightbins)-1);
for dj2 = 1 : length(rights)
    rightdistr(dj2,:) = histcounts(rightcentrals{dj2},spatbins);
    rightarmdistr(dj2,:) = histcounts(rightarms{dj2},rightbins);
end

%%
%Section 5: Left-Right trajectory distance in the stem and TZ
disp('Left-right trajectory separation...')

leftbini = armtzposbini(lefts);
rightbini = armtzposbini(rights);
leftpos = fulltrialpos(lefts);
leftpos = cat(2,leftpos{:});
leftX = leftpos(1,:);
rightpos = fulltrialpos(rights);
rightpos = cat(2,rightpos{:});
rightX = rightpos(1,:);
leftbinsd = armtzposbini(lefts);
leftbinsd = cat(2,leftbinsd{:});
rightbinsd = armtzposbini(rights);
rightbinsd = cat(2,rightbinsd{:});
LRdistonstem = zeros(1,length(TZSbins));
LRdsd = zeros(1,length(TZSbins));
LRdistcell = cell(1,length(TZSbins));
for sx1 = 1 : length(TZSbins) - 1
    Lfb = leftX(leftbinsd == sx1);
    Rfb = rightX(rightbinsd == sx1);
    allsubt = zeros(length(Lfb),length(Rfb));
    for sx2 = 1 : length(Lfb)
        for sx3 = 1 : length(Rfb)
            allsubt(sx2,sx3) = abs(Lfb(sx2) - Rfb(sx3));
        end
    end
    LRdistcell{sx1} = reshape(allsubt,1,size(allsubt,1) * size(allsubt,2));
    LRdistonstem(sx1) = mean(mean(allsubt));
    LRdsd(sx1) = std(LRdistcell{sx1});
end

pureleftpos = fulltrialpos(purelefts);
pureleftpos = cat(2,pureleftpos{:});
pureleftX = pureleftpos(1,:);
purerightpos = fulltrialpos(purerights);
purerightpos = cat(2,purerightpos{:});
purerightX = purerightpos(1,:);
pureleftbinsd = armtzposbini(purelefts);
pureleftbinsd = cat(2,pureleftbinsd{:});
purerightbinsd = armtzposbini(purerights);
purerightbinsd = cat(2,purerightbinsd{:});
cLRdistonstem = zeros(1,length(TZSbins));
cLRdsd = zeros(1,length(TZSbins));
cLRdistcell = cell(1,length(TZSbins));
for sx1 = 1 : length(TZSbins) - 1
    pLfb = pureleftX(pureleftbinsd == sx1);
    pRfb = purerightX(purerightbinsd == sx1);
    allpsubt = zeros(length(pLfb),length(pRfb));
    for sx2 = 1 : length(pLfb)
        for sx3 = 1 : length(pRfb)
            allpsubt(sx2,sx3) = abs(pLfb(sx2) - pRfb(sx3));
        end
    end
    cLRdistcell{sx1} = reshape(allpsubt,1,size(allpsubt,1) * size(allpsubt,2));
    cLRdistonstem(sx1) = mean(mean(allpsubt));
    cLRdsd(sx1) = std(cLRdistcell{sx1});
end

% if strcmp(sns,'yes') == 1
%     LRdistruct = struct('allLRdistances',LRdistcell,'allLRdistave',LRdistonstem,'correctLRdistances',cLRdistcell,'correctLRdistave',cLRdistonstem);
%     save([resdir '\' sessionname{sessnum} '_LeftRightDistance.mat'],'LRdistruct')
% end

PFA_trajdistr_raw = struct('LR_distances_all',LRdistcell,'LR_distances_mean',LRdistonstem,'LR_distances_correct_all',cLRdistcell,'LR_distances_correct_mean',cLRdistonstem);
PFA_trajdistr = struct('LR_distances_mean',LRdistonstem,'LR_distances_std',LRdsd,'LR_distances_correct_mean',cLRdistonstem,'LR_distances_correct_std',cLRdsd);
filenameTRD = [sessionname{sessnum} '_trajectory_divergence.mat'];
save(fullfile(resdir,sessionname{sessnum},filenameTRD),'PFA_trajdistr')

%%
%Section 6: shuffling IN spiketimes
disp('Shuffling interneuron spiketimes...')

timeshIN = cell(1,length(timeIN));
for sq1 = 1 : length(timeIN)
%     disp(sq1)
    selectedspikes = timeIN{sq1};
    [shiftedtrain] = timespikeshuffler(selectedspikes);
    timeshIN{sq1} = shiftedtrain;
end
% shuffIN = struct('shuffled IN times',timeshIN);

%%
%Section 7: loading spike positions and spiketimes into cell array
disp('Calculation of smoothed tuning curves: position per trial vectors...')

mazestart = alltime(1);
mazend = alltime(end);
%full maze (PCs)
for j1 = 1 : length(pc_ID)
    tnum = trialnum{j1}(times{j1} >= mazestart & times{j1} <= mazend);
    farmspikes = cell(1,constrialnum);
    linmazespikes = cell(1,constrialnum);
    mspkt = cell(1,constrialnum);
    mspkdp = cell(1,constrialnum);
    mspkp = cell(1,constrialnum);
    for j2 = 1 : constrialnum
        qvect = fullpos{j1}([1,2],tnum == j2);
        linvect = linpos{j1}(tnum == j2);
        tivect = times{j1}(tnum == j2);
        pivect = spikephase{j1}(tnum == j2);
        if isempty(qvect) == 0
            farmspikes{j2} = qvect;
            linmazespikes{j2} = linvect;
            mspkt{j2} = tivect;
            mspkp{j2} = pivect;
            mspkdp{j2} = diff(tivect);
        else
            farmspikes{j2} = NaN;
            linmazespikes{j2} = NaN;
            mspkt{j2} = NaN;
            mspkdp{j2} = NaN;
            mspkp{j2} = NaN;
        end
    end
%     farmspikes(isempty(farmspikes) == 1) = [];
%     linmazespikes(isempty(linmazespikes) == 1) = [];
%     mspkt(isempty(mspkt) == 1) = [];
%     mspkp(isempty(mspkp) == 1) = [];
    unitarm_full{j1} = farmspikes;
    unitarm_fullin{j1} = linmazespikes;
    unitPTP_full{j1} = mspkt;
    unitDP_full{j1} = mspkdp;
    unit_fullphase{j1} = mspkp;
end

%only stem (PCs)
for i1 = 1 : length(pc_ID)
    tnums = trialnum{i1}(times{i1} >= mazestart & times{i1} <= mazend);
    armspikes = cell(1,constrialnum);
    armspkt = cell(1,constrialnum);
    armdp = cell(1,constrialnum);
    armphase = cell(1,constrialnum);
    for i2 = 1 : constrialnum
        ttvect = linpos{i1}(tnums == i2);
        tvect = ttvect(ttvect >= 20 & ttvect < 62); % only central arm, chenged from 16 to 20
        phasevect = spikephase{i1}(tnums == i2);
        phasevect = phasevect(ttvect >= 20 & ttvect < 62);
        ativect = times{i1}(tnums == i2);
        ativect = ativect(ttvect >= 20 & ttvect < 62);
        if isempty(tvect) == 0
            armspikes{i2} = tvect;
            armspkt{i2} = ativect;
            armdp{i2} = diff(tvect(isnan(tvect) == 0));
            armphase{i2} = phasevect;
        else
            armspikes{i2} = NaN;
            armspkt{i2} = NaN;
            armdp{i2} = NaN;
            armphase{i2} = NaN;
        end
    end
%     armspikes(isempty(armspikes) == 1) = [];
%     armspkt(isempty(armspkt) == 1) = [];
    unitarm{i1} = armspikes;
    unitPTP_arm{i1} = armspkt;
    unitDP_arm{i1} = armdp;
    unitarm_phase{i1} = armphase;
end

%only TZ (PCs)
for i1 = 1 : length(pc_ID)
    tnumz = trialnum{i1}(times{i1} >= mazestart & times{i1} <= mazend);
    TZspikes = cell(1,constrialnum);
    TZspkt = cell(1,constrialnum);
    TZdp = cell(1,constrialnum);
    TZphase = cell(1,constrialnum);
    for i2 = 1 : constrialnum
        ttzvect = linpos{i1}(tnumz == i2);
        tzvect = ttzvect(ttzvect >= 0 & ttzvect < 20); % only central arm, chenged from 16 to 20
        phaseTZvect = spikephase{i1}(tnumz == i2);
        phaseTZvect = phaseTZvect(ttzvect >= 0 & ttzvect < 20);
        atzivect = times{i1}(tnumz == i2);
        atzivect = atzivect(ttzvect >= 0 & ttzvect < 20);
        if isempty(tzvect) == 0
            TZspikes{i2} = tzvect;
            TZspkt{i2} = atzivect;
            TZdp{i2} = diff(tzvect(isnan(tzvect) == 0));
            TZphase{i2} = phasevect;
        else
            TZspikes{i2} = NaN;
            TZspkt{i2} = NaN;
            TZdp{i2} = NaN;
            TZphase{i2} = NaN;
        end
    end
%     armspikes(isempty(armspikes) == 1) = [];
%     armspkt(isempty(armspkt) == 1) = [];
    unitranz{i1} = TZspikes;
    unitPTP_tranz{i1} = TZspkt;
    unitDP_tranz{i1} = TZdp;
    unitranz_phase{i1} = TZphase;
end

%full maze (INs)
for nj1 = 1 : length(IN_ID)
    tnumIN = trialnumIN{nj1}(timeIN{nj1} >= mazestart & timeIN{nj1} <= mazend);
    inspikepos = cell(1,constrialnum);
    inspiketimes = cell(1,constrialnum);
    for nj2 = 1 : constrialnum
        intivect = timeIN{nj1}(tnumIN == nj2);
        inivect = INpos{nj1}(tnumIN == nj2);
        if isempty(inivect) == 0
            inspiketimes{nj2} = intivect;
            inspikepos{nj2} = inivect;
        else
            inspiketimes{nj2} = NaN;
            inspikepos{nj2} = NaN;
        end
    end
    INtfull{nj1} = inspiketimes;
    INfullpos{nj1} = inspikepos;
end

%only stem (INs)
for nj3 = 1 : length(IN_ID)
    tnumsIN = trialnumIN{nj3}(timeIN{nj3} >= mazestart & timeIN{nj3} <= mazend);
    inarmspikepos = cell(1,constrialnum);
    inarmspiketimes = cell(1,constrialnum);
    inarmspikephase = cell(1,constrialnum);
    for nj4 = 1 : constrialnum
        inativect = timeIN{nj3}(tnumsIN == nj4);
        tiniavect = INpos{nj3}(tnumsIN == nj4);
        inphasevect = spikephaseIN{nj3}(tnumsIN == nj4);
        iniavect = tiniavect(tiniavect >= 20 & tiniavect < 62);
        inativect = inativect(tiniavect >= 20 & tiniavect < 62);
        inphasevect = inphasevect(tiniavect >= 20 & tiniavect < 62);
        if isempty(inivect) == 0
            inarmspiketimes{nj4} = inativect;
            inarmspikepos{nj4} = iniavect;
            inarmspikephase{nj4} = inphasevect;
        else
            inarmspiketimes{nj4} = NaN;
            inarmspikepos{nj4} = NaN;
            inarmspikephase{nj4} = NaN;
        end
    end
    INtarm{nj3} = inarmspiketimes;
    INarmpos{nj3} = inarmspikepos;
    INarmphase{nj3} = inarmspikephase;
end

%only TZ (INs)
for nj3 = 1 : length(IN_ID)
    tnumzIN = trialnumIN{nj3}(timeIN{nj3} >= mazestart & timeIN{nj3} <= mazend);
    inTZspikepos = cell(1,constrialnum);
    inTZspiketimes = cell(1,constrialnum);
    inTZspikephase = cell(1,constrialnum);
    for nj4 = 1 : constrialnum
        inatzivect = timeIN{nj3}(tnumzIN == nj4);
        tiniavectz = INpos{nj3}(tnumzIN == nj4);
        inphasevectz = spikephaseIN{nj3}(tnumzIN == nj4);
        iniavectz = tiniavectz(tiniavectz >= 0 & tiniavectz < 20);
        inatzivect = inatzivect(tiniavectz >= 0 & tiniavectz < 20);
        inphasevectz = inphasevectz(tiniavectz >= 0 & tiniavectz < 20);
        if isempty(iniavectz) == 0
            inTZspiketimes{nj4} = inatzivect;
            inTZspikepos{nj4} = iniavectz;
            inTZspikephase{nj4} = inphasevectz;
        else
            inTZspiketimes{nj4} = NaN;
            inTZspikepos{nj4} = NaN;
            inTZspikephase{nj4} = NaN;
        end
    end
    INtranz{nj3} = inTZspiketimes;
    INtranzpos{nj3} = inTZspikepos;
    INtranzphase{nj3} = inTZspikephase;
end

if exist('timeshIN','var')
%shuffled IN spiketimes in stem
    shINtarm = cell(1,length(timeshIN));
    for nq1 = 1 : length(timeshIN)
        trialshc = cell(1,constrialnum);
        for nq2 = 1 : constrialnum
%             disp([num2str(nq1) '-' num2str(nq2)])
            tshc = cell(1,size(timeshIN{nq1},1));
            for nq3 = 1 : size(timeshIN{nq1},1)
                if ~isempty(alltrialtimestem{nq2})
                    tshc{nq3} = timeshIN{nq1}(nq3,timeshIN{nq1}(nq3,:) >= alltrialtimestem{nq2}(1) & timeshIN{nq1}(nq3,:) <= alltrialtimestem{nq2}(end));
                else
                    tshc{nq3} = NaN;
                end
            end
            trialshc{nq2} = tshc;
        end
        shINtarm{nq1} = trialshc;
    end
%shuffled IN spiketimes in TZ
    shINtz = cell(1,length(timeshIN));
    for np1 = 1 : length(timeshIN)
        trialshcz = cell(1,constrialnum);
        for np2 = 1 : constrialnum
%             disp([num2str(np1) '-' num2str(np2)])
            tshcz = cell(1,size(timeshIN{np1},1));
            for np3 = 1 : size(timeshIN{np1},1)
                if ~isempty(alltrialtimetranz{np2})
                    tshcz{np3} = timeshIN{np1}(np3,timeshIN{np1}(np3,:) >= alltrialtimetranz{np2}(1) & timeshIN{np1}(np3,:) <= alltrialtimetranz{np2}(end));
                else
                    tshcz{np3} = NaN;
                end
            end
            trialshcz{np2} = tshcz;
        end
        shINtz{np1} = trialshc;
    end
end

PFA_PC_trialspikepos = struct('full_maze_spikepos',unitarm_full,'full_maze_spikelinpos',unitarm_fullin,'full_maze_spiketimes',unitPTP_full,'full_maze_spikephase',unit_fullphase,...
    'stem_spikelinpos',unitarm,'stem_spiketimes',unitPTP_arm,'stem_spikephase',unitarm_phase,...
    'TZ_spikelinpos',unitranz,'TZ_spiketimes',unitPTP_tranz,'TZ_spikephase',unitranz_phase);
PFA_IN_trialspikepos = struct('IN_full_spikelinpos',INfullpos,'IN_full_spiketimes',INtfull,...
    'IN_stem_spikelinpos',INarmpos,'IN_stem_spiketimes',INtarm,'IN_stem_spikephase',INarmphase,...
    'IN_TZ_spikelinpos',INtranzpos,'IN_TZ_spiketimes',INtranz,'IN_TZ_spikephase',INtranzphase,...
    'IN_shuffled_stem_spiketimes',shINtarm,'IN_shuffled_TZ_spiketimes',shINtz);
filenamePC = [sessionname{sessnum} '_spike_positions_PC_trials.mat'];
filenameIN = [sessionname{sessnum} '_spike_positions_IN_trials.mat'];
disp('Saving spike position data...(just for checking if saving is really happening)')
save(fullfile(resdir,sessionname{sessnum},filenamePC),'PFA_PC_trialspikepos','-v7.3')
save(fullfile(resdir,sessionname{sessnum},filenameIN),'PFA_IN_trialspikepos','-v7.3')

%%
%Section 8: calculation of tuning curves: counting spike numbers in spatial bins (all trials counted) in the stem
disp('Calculation of smoothed tuning curves: firing rate vectors...')

% 2 cm bins
spatbins2 = 20:2:62; %chenged from 16 to 20
spatfullbins = 0 : 2 : 205;
spatvects2 = cell(1,length(pc_ID));
spatFRvects2 = cell(1,length(pc_ID));
tranzFRvects2 = cell(1,length(pc_ID));
fullFRvects2 = cell(1,length(pc_ID)); % full maze raw tuning curves, 2 cm bins
fullFRS = zeros(1,length(pc_ID));
avenumspk2 = zeros(length(pc_ID),constrialnum);
for i3 = 1 : length(pc_ID)
    spatmtx2 = zeros(length(spatbins2)-1,constrialnum);
    spatmtx2FR = zeros(length(spatbins2)-1,constrialnum);
    tranzmtx2FR = zeros(length(tranzbins)-1,constrialnum);
    spatmtx2FRfull = zeros(length(spatfullbins)-1,constrialnum);
    for i5 = 1 : length(unitarm_fullin{i3})
        tfn2 = histcounts(unitarm_fullin{i3}{i5},spatfullbins);
        tfnz = histcounts(unitranz{i3}{i5},tranzbins);
        if isempty(tfn2)
            spatmtx2FRfull(:,i5) = spatmtx2FRfull(:,i5);
        else
            tdfr = fulldistr(i5,:);
            if ~isequal(size(tfn2,1),size(tdfr,1))
                tdfr = tdfr';
            end
            spatmtx2FRfull(:,i5) = (tfn2 ./ tdfr) * fps;
        end
        if isempty(tfnz)
            tranzmtx2FR(:,i5) = tranzmtx2FR(:,i5);
        else
            tzfr = tranzdistr(i5,:);
            if ~isequal(size(tfnz,1),size(tzfr,1))
                tzfr = tzfr';
            end
            tranzmtx2FR(:,i5) = (tfnz ./ tzfr) * fps;
        end
    end
    for i4 = 1 : length(unitarm{i3})
        tn2 = histcounts(unitarm{i3}{i4},spatbins2);
        if isempty(tn2)
            spatmtx2(:,i4) = spatmtx2(:,i4);
            spatmtx2FR(:,i4) = spatmtx2FR(:,i4);
        else
            tdr = alldistr(i4,:);
            if ~isequal(size(tn2,1),size(tdr,1))
                tdr = tdr';
            end
            spatmtx2FR(:,i4) = (tn2 ./ tdr) * fps;
        end
        avenumspk2(i3,i4) = mean(spatmtx2(:,i4));
    end
    spatvects2{i3} = spatmtx2; %tuning curves w 2 cm bins (spike count)
    spatFRvects2{i3} = spatmtx2FR; %tuning curves w 2 cm bins (FR)
    tranzFRvects2{i3} = tranzmtx2FR;
    fullFRvects2{i3} = spatmtx2FRfull;
    fullFRS(i3) = mean(mean(spatmtx2FR)); %whole session firing rate
end

%IN tuning curves
inSPKvects2 = cell(1,length(INarmpos));
inFRvects2 = cell(1,length(INarmpos));
intzFRvects2 = cell(1,length(INarmpos));
infullFRvect = cell(1,length(INarmpos));
for in1 = 1 : length(INarmpos)
    %stem
    INmtx2 = zeros(constrialnum,length(spatbins2)-1);
    INfrmtx2 = zeros(constrialnum,length(spatbins2)-1);
    for in2 = 1 : constrialnum
        tin2 = histcounts(INarmpos{in1}{in2},spatbins2);
        if size(tin2,1) > size(tin2,2)
            tin2 = tin2';
        end
        if isempty(tin2) == 1
            INmtx2(in2,:) = INmtx2(in2,:);
            INfrmtx2(in2,:) = INfrmtx2(in2,:);
        else
            INmtx2(in2,:) = tin2;
            tdr = alldistr(in2,:);
%             if isequal(size(tin2,2),size(tdr,2)) == 0
%                 tdr = tdr';
%             end
            INfrmtx2(in2,:) = (tin2 ./ tdr) * fps;
        end
    end
    inSPKvects2{in1} = INmtx2;
    inFRvects2{in1} = INfrmtx2;
    %TZ
    INtzfrmtx2 = zeros(constrialnum,length(tranzbins)-1);
    for in4 = 1 : constrialnum
        tinz = histcounts(INtranzpos{in1}{in4},tranzbins);
        if size(tinz,1) > size(tinz,2)
            tinz = tinz';
        end
        if isempty(tinz)
            INtzfrmtx2(in4,:) = INtzfrmtx2(in4,:);
        else
            tzdr = tranzdistr(in4,:);
%             if isequal(size(tin2,2),size(tdr,2)) == 0
%                 tdr = tdr';
%             end
            INtzfrmtx2(in4,:) = (tinz ./ tzdr) * fps;
        end
    end
    intzFRvects2{in1} = INtzfrmtx2;
    %full
    INfullmtx = zeros(constrialnum,length(spatfullbins)-1);
    for in3 = 1 : constrialnum
        tinf = histcounts(INfullpos{in1}{in3},spatfullbins);
        if size(tinf,1) > size(tinf,2)
            tinf = tinf';
        end
        if isempty(tinf) == 1
            INfullmtx(in3,:) = INfullmtx(in3,:);
        else
            INfullmtx(in3,:) = tinf;
            tindfr = fulldistr(in3,:);
%             if isequal(size(tin2,2),size(tdr,2)) == 0
%                 tdr = tdr';
%             end
            INfullmtx(in3,:) = (tinf ./ tindfr) * fps;
        end
    end
    infullFRvect{in1} = INfullmtx;
end

%counting spike numbers in spatial bins (all trials counted) on the sides
%(only 2 cm bins)
spatvects2_leftarms = cell(1,length(pc_ID));
spatvects2_rightarms = cell(1,length(pc_ID));
spatFRvects2_leftarms = cell(1,length(pc_ID));
spatFRvects2_rightarms = cell(1,length(pc_ID));
for i5 = 1 : length(unitarm_fullin)
    spatmtx2_left = zeros(length(leftbins)-1,length(lefts));
    spatmtx2_right = zeros(length(rightbins)-1,length(rights));
    spatmtx2_FRleft = zeros(length(leftbins)-1,length(lefts));
    spatmtx2_FRright = zeros(length(rightbins)-1,length(rights));
    leftspikes = unitarm_fullin{i5}(lefts);
    rightspikes = unitarm_fullin{i5}(rights);
    for i6 = 1 : length(leftspikes)
        spkleft = leftspikes{i6};
        if isempty(spkleft) == 1
            spatmtx2_left(:,i6) = spatmtx2_left(:,i6);
            spatmtx2_FRleft(:,i6) = spatmtx2_FRleft(:,i6);
        else
            spatmtx2_left(:,i6) = histcounts(spkleft,leftbins);
            spatmtx2_FRleft(:,i6) = spatmtx2_left(:,i6) ./ leftarmdistr(i6,:)';
        end
    end
    spatvects2_leftarms{i5} = spatmtx2_left;
    spatFRvects2_leftarms{i5} = spatmtx2_FRleft;
    
    for i7 = 1 : length(rightspikes)
        spkright = rightspikes{i7};
        if isempty(spkright) == 1
            spatmtx2_right(:,i7) = spatmtx2_right(:,i7);
            spatmtx2_FRright(:,i7) = spatmtx2_FRright(:,i7);
        else
            spatmtx2_right(:,i7) = histcounts(spkright,rightbins);
            spatmtx2_FRright(:,i7) = spatmtx2_right(:,i7) ./ rightarmdistr(i7,:)';
        end
    end
    spatvects2_rightarms{i5} = spatmtx2_right; %spike count tuning curve (2 cm bins)
    spatFRvects2_rightarms{i5} = spatmtx2_FRright; %FR tuning curves (2 cm bins)
end

PFA_PC_FRvectors = struct('session_FR',fullFRS,'stem_FR_vectors',spatFRvects2,'TZ_FR_vectors',tranzFRvects2,'full_maze_FR_vectors',fullFRvects2);
PFA_IN_FRvectors = struct('IN_stem_FR_vectors',inFRvects2,'IN_TZ_FR_vectors',intzFRvects2,'IN_full_maze_FR_vectors',infullFRvect);
filenameFRvPC = [sessionname{sessnum} '_PC_FR_vectors.mat'];
filenameFRvIN = [sessionname{sessnum} '_IN_FR_vectors.mat'];
save(fullfile(resdir,sessionname{sessnum},filenameFRvPC),'PFA_PC_FRvectors')
save(fullfile(resdir,sessionname{sessnum},filenameFRvIN),'PFA_IN_FRvectors')

%%
% Section 9: Gaussian smoothing of (putative) place cell tuning curves and elimination
% of units with too few spikes
disp('Calculation of smoothed tuning curves: smoothing firing rate vectors...')

spkcountS = zeros(1,length(spatvects2));
spkcountL = zeros(1,length(spatvects2));
spkcountR = zeros(1,length(spatvects2));
upspatcells = cell(5,length(spatvects2));
spatwin = gausswin(500);
for ui0 = 1 : length(spatFRvects2)
    spkcountS(ui0) = size(cat(1,unitarm{ui0}{:}),1);
    tscL = unitarm_fullin{ui0}(lefts);
    spkcountL(ui0) = size(cat(1,tscL{:}),1);
    tscR = unitarm_fullin{ui0}(rights);
    spkcountR(ui0) = size(cat(1,tscR{:}),1);
    upspatmtxFull = zeros(size(fullFRvects2{ui0},2),size(fullFRvects2{ui0},1) - 1);
    upspatmtxStem = zeros(size(spatFRvects2{ui0},2),size(spatFRvects2{ui0},1) - 1);
    upspatmtxTranz = zeros(size(tranzFRvects2{ui0},2),size(tranzFRvects2{ui0},1) - 1);
    upspatmtxLeft = zeros(length(lefts),size(spatFRvects2_leftarms{ui0},1));
    upspatmtxRight = zeros(length(rights),size(spatFRvects2_rightarms{ui0},1));
    %full maze
    tmtxFull = fullFRvects2{ui0}(1 : end-1,:)'; %can be spike count- or FR-based
    for uS = 1 : constrialnum
        uptvFull = tmtxFull(uS,:);
        uptvFull(isnan(uptvFull)) = 0;
        uptvFull(isinf(uptvFull)) = 0;
        orixFull = 1 : length(uptvFull);
        newxFull = 1 : 1 / 100 : length(uptvFull);
        upvectFull = interp1(orixFull,uptvFull,newxFull);
        convectFull = conv(upvectFull,spatwin);
        convectFull = convectFull(length(spatwin) / 2 + 1 : length(upvectFull) + length(spatwin) / 2);
        upspatmtxFull(uS,:) = convectFull(1 : 100 : end);
    end
    %stem
    tmtxStem = spatFRvects2{ui0}(1 : end-1,:)'; %can be spike count- or FR-based
    for uS = 1 : constrialnum
        uptvStem = tmtxStem(uS,:);
        uptvStem(isnan(uptvStem)) = 0;
        uptvStem(isinf(uptvStem)) = 0;
        orixStem = 1 : length(uptvStem);
        newxStem = 1 : 1 / 100 : length(uptvStem);
        upvectStem = interp1(orixStem,uptvStem,newxStem);
        convectStem = conv(upvectStem,spatwin);
        convectStem = convectStem(length(spatwin) / 2 + 1 : length(upvectStem) + length(spatwin) / 2);
        upspatmtxStem(uS,:) = convectStem(1 : 100 : end);
    end
    %TZ
    tmtranz = tranzFRvects2{ui0}(1 : end-1,:)'; %can be spike count- or FR-based
    for uT = 1 : constrialnum
        uptranz = tmtranz(uT,:);
        uptranz(isnan(uptranz)) = 0;
        uptranz(isinf(uptranz)) = 0;
        orixTranz = 1 : length(uptranz);
        newxTranz = 1 : 1 / 100 : length(uptranz);
        upvecTranz = interp1(orixTranz,uptranz,newxTranz);
        convecTranz = conv(upvecTranz,spatwin);
        convecTranz = convecTranz(length(spatwin) / 2 + 1 : length(upvecTranz) + length(spatwin) / 2);
        upspatmtxTranz(uT,:) = convecTranz(1 : 100 : end);
    end
    %left arm
    tmtxLeft = spatFRvects2_leftarms{ui0}'; %can be spike count- or FR-based
    for uL = 1 : length(lefts)
        uptvLeft = tmtxLeft(uL,:);
        uptvLeft(isnan(uptvLeft)) = 0;
        uptvLeft(isinf(uptvLeft)) = 0;
        orixLeft = 1 : length(uptvLeft);
        newxLeft = 1 : 1 / 100 : length(uptvLeft);
        upvectLeft = interp1(orixLeft,uptvLeft,newxLeft);
        convectLeft = conv(upvectLeft,spatwin);
        convectLeft = convectLeft(length(spatwin) / 2 + 1 : length(upvectLeft) + length(spatwin) / 2);
        tutL = convectLeft(1 : 100 : end);
        if length(tutL) > size(upspatmtxLeft,2)
            upspatmtxLeft(uL,:) = tutL(1 : end - 1);
        end
        if isequal(length(tutL),size(upspatmtxLeft,2)) == 1
            upspatmtxLeft(uL,:) = tutL;
        end
    end
    %right arm
    tmtxRight = spatFRvects2_rightarms{ui0}'; %can be spike count- or FR-based
    for uR = 1 : length(rights)
        uptvRight = tmtxRight(uR,:);
        uptvRight(isnan(uptvRight)) = 0;
        uptvRight(isinf(uptvRight)) = 0;
        orixRight = 1 : length(uptvRight);
        newxRight = 1 : 1 / 100 : length(uptvRight);
        upvectRight = interp1(orixRight,uptvRight,newxRight);
        convectRight = conv(upvectRight,spatwin);
        convectRight = convectRight(length(spatwin) / 2 + 1 : length(upvectRight) + length(spatwin) / 2);
        tutR = convectRight(1 : 100 : end);
        if length(tutR) > size(upspatmtxRight,2)
            upspatmtxRight(uR,:) = tutR(1 : end - 1);
        elseif isequal(length(tutR),size(upspatmtxRight,2)) == 1
            upspatmtxRight(uR,:) = tutR;
        end
    end
    upspatcells{1,ui0} = upspatmtxStem; %stem smoothed tuning curves
    upspatcells{2,ui0} = upspatmtxLeft; %left smoothed tuning curves
    upspatcells{3,ui0} = upspatmtxRight; %right smoothed tuning curves
    upspatcells{4,ui0} = upspatmtxFull; %full maze smoothed tuning curves
    upspatcells{5,ui0} = upspatmtxTranz; %TZ smoothed tuning curves
end

%Gaussian smoothing of (putative) interneuron tuning curves
INcountS = zeros(1,length(inFRvects2));
INcountL = zeros(1,length(inFRvects2));
INcountR = zeros(1,length(inFRvects2));
upINspatcells = cell(2,length(inFRvects2));
for ui0 = 1 : length(inFRvects2)
    INcountS(ui0) = size(cat(1,INarmpos{ui0}{:}),1);
    tinscL = INfullpos{ui0}(lefts);
    INcountL(ui0) = size(cat(1,tinscL{:}),1);
    tinscR = INfullpos{ui0}(rights);
    INcountR(ui0) = size(cat(1,tinscR{:}),1);
    upINspatmtxFull = zeros(size(infullFRvect{ui0},1),size(infullFRvect{ui0},2) - 1);
    upINspatmtxStem = zeros(size(inFRvects2{ui0},1),size(inFRvects2{ui0},2) - 1);
    upINspatmtxTZ = zeros(size(intzFRvects2{ui0},1),size(intzFRvects2{ui0},2) - 1);
%     upINspatmtxLeft = zeros(length(lefts),size(spatFRvects2_leftarms{ui0},1));
%     upINspatmtxRight = zeros(length(rights),size(spatFRvects2_rightarms{ui0},1));
    %full maze
    tmtxFullin = infullFRvect{ui0}(:,1:end-1); %can be spike count- or FR-based
    for uS = 1 : constrialnum
        uptvFullin = tmtxFullin(uS,:);
        uptvFullin(isnan(uptvFullin)) = 0;
        uptvFullin(isinf(uptvFullin)) = 0;
        orixFullin = 1 : length(uptvFullin);
        newxFullin = 1 : 1 / 100 : length(uptvFullin);
        upvectFullin = interp1(orixFullin,uptvFullin,newxFullin);
        convectFullin = conv(upvectFullin,spatwin);
        convectFullin = convectFullin(length(spatwin) / 2 + 1 : length(upvectFullin) + length(spatwin) / 2);
        upINspatmtxFull(uS,:) = convectFullin(1 : 100 : end);
    end
    %stem
    tmtxStemin = inFRvects2{ui0}(:,1:end-1); %can be spike count- or FR-based
    for uS = 1 : constrialnum
        uptvStemin = tmtxStemin(uS,:);
        uptvStemin(isnan(uptvStemin)) = 0;
        uptvStemin(isinf(uptvStemin)) = 0;
        orixStemin = 1 : length(uptvStemin);
        newxStemin = 1 : 1 / 100 : length(uptvStemin);
        upvectStemin = interp1(orixStemin,uptvStemin,newxStemin);
        convectStemin = conv(upvectStemin,spatwin);
        convectStemin = convectStemin(length(spatwin) / 2 + 1 : length(upvectStemin) + length(spatwin) / 2);
        upINspatmtxStem(uS,:) = convectStemin(1 : 100 : end);
    end
    %TC
    tmTZin = intzFRvects2{ui0}(:,1:end-1); %can be spike count- or FR-based
    for uT = 1 : constrialnum
        upTZin = tmTZin(uT,:);
        upTZin(isnan(upTZin)) = 0;
        upTZin(isinf(upTZin)) = 0;
        orixTZin = 1 : length(upTZin);
        newxTZin = 1 : 1 / 100 : length(upTZin);
        upvecTZin = interp1(orixTZin,upTZin,newxTZin);
        convecTZin = conv(upvecTZin,spatwin);
        convecTZin = convecTZin(length(spatwin) / 2 + 1 : length(upvecTZin) + length(spatwin) / 2);
        upINspatmtxTZ(uT,:) = convecTZin(1 : 100 : end);
    end
    %left arm
%     tmtxLeft = spatFRvects2_leftarms{ui0}'; %can be spike count- or FR-based
%     for uL = 1 : length(lefts)
%         uptvLeft = tmtxLeft(uL,:);
%         uptvLeft(isnan(uptvLeft)) = 0;
%         uptvLeft(isinf(uptvLeft)) = 0;
%         orixLeft = 1 : length(uptvLeft);
%         newxLeft = 1 : 1 / 100 : length(uptvLeft);
%         upvectLeft = interp1(orixLeft,uptvLeft,newxLeft);
%         convectLeft = conv(upvectLeft,spatwin);
%         convectLeft = convectLeft(length(spatwin) / 2 + 1 : length(upvectLeft) + length(spatwin) / 2);
%         tutL = convectLeft(1 : 100 : end);
%         if length(tutL) > size(upspatmtxLeft,2)
%             upspatmtxLeft(uL,:) = tutL(1 : end - 1);
%         end
%         if isequal(length(tutL),size(upspatmtxLeft,2)) == 1
%             upspatmtxLeft(uL,:) = tutL;
%         end
%     end
    %right arm
%     tmtxRight = spatFRvects2_rightarms{ui0}'; %can be spike count- or FR-based
%     for uR = 1 : length(rights)
%         uptvRight = tmtxRight(uR,:);
%         uptvRight(isnan(uptvRight)) = 0;
%         uptvRight(isinf(uptvRight)) = 0;
%         orixRight = 1 : length(uptvRight);
%         newxRight = 1 : 1 / 100 : length(uptvRight);
%         upvectRight = interp1(orixRight,uptvRight,newxRight);
%         convectRight = conv(upvectRight,spatwin);
%         convectRight = convectRight(length(spatwin) / 2 + 1 : length(upvectRight) + length(spatwin) / 2);
%         tutR = convectRight(1 : 100 : end);
%         if length(tutR) > size(upspatmtxRight,2)
%             upspatmtxRight(uR,:) = tutR(1 : end - 1);
%         elseif isequal(length(tutR),size(upspatmtxRight,2)) == 1
%             upspatmtxRight(uR,:) = tutR;
%         end
%     end
    upINspatcells{1,ui0} = upINspatmtxStem; %stem smoothed tuning curves
%     upINspatcells{2,ui0} = upINspatmtxLeft; %left smoothed tuning curves
%     upINspatcells{3,ui0} = upINspatmtxRight; %right smoothed tuning curves
    upINspatcells{4,ui0} = upINspatmtxFull; %full maze smoothed tuning curves
    upINspatcells{5,ui0} = upINspatmtxTZ; %TZ maze smoothed tuning curves
end

%sorting out pyramidal cells with not enough spikes (limit = 300): NOT
%USED!
clusteredPCS = find(spkcountS > 300);
clusteredPCL = find(spkcountL > 300);
clusteredPCR = find(spkcountR > 300);

mtx2Scorr = upspatcells(1,:);
mtx2Full = upspatcells(4,:);
mtx2TZ = upspatcells(5,:);
selectedFRvects = spatFRvects2;
selectedFRfull = fullFRvects2;
selFRS = fullFRS;

%NOT USED:-----------------------------------------------------------------
abstrmmRatio = zeros(1,length(mtx2Scorr));
for aj1 = 1 : length(mtx2Scorr)
    if min(mean(mtx2Scorr{aj1})) == 0
        mime = 1;
    else
        mime = min(mean(mtx2Scorr{aj1}));
    end
    abstrmmRatio(aj1) = max(mean(mtx2Scorr{aj1})) / mime;
end
ihavePF = find(abstrmmRatio >= 3);
%--------------------------------------------------------------------------

ihmtx2Scorr = mtx2Scorr;
ihmtx2TZcorr = mtx2TZ;
ihmtx2full = mtx2Full;

PFA_smoothed_PC_FR_vectors = struct('PC_vectors',upspatcells);
PFA_smoothed_IN_FR_vectors = struct('IN_vectors',upINspatcells);
filenamePCsFRv = [sessionname{sessnum} '_smoothed_PC_FR_vectors.mat'];
filenameINsFRv = [sessionname{sessnum} '_smoothed_IN_FR_vectors.mat'];
save(fullfile(resdir,sessionname{sessnum},filenamePCsFRv),'PFA_smoothed_PC_FR_vectors')
save(fullfile(resdir,sessionname{sessnum},filenameINsFRv),'PFA_smoothed_IN_FR_vectors')

%%
%Place field identification
disp('Searching for place fields...')
mazepart = 'full';
inX = mtx2Full;
idX = 1 : length(mtx2Full);
idXX = 1 : length(mtx2Full);
numoT = constrialnum;
spaBin = spatfullbins;
spktimes = unitPTP_full;
spkpos = unitarm_fullin;
rawTC = fullFRvects2;

[placesupercell,pFR,allSfields,~,~] = placefieldidentifierr(pc_ID,inX,...
    idX,idXX,numoT,spktimes,spkpos,rawTC,spaBin,mazepart);

[PFregionS,cellswithfields] = PFidentifierSimple(mtx2Full,fullFRvects2);

Pplaces = placesupercell;
PpFR = pFR;
PallSfields = allSfields;
PPFregionS = PFregionS(cellswithfields ~= 0);

%Selection of place cells with stem fields & their fields in the stem
stnst1 = zeros(1,length(PPFregionS));
stemregions = cell(1,length(PPFregionS));
mtx2Field = mtx2Full(cellswithfields > 0);
for sf1 = 1 : length(PPFregionS)
    stnst2 = zeros(1,size(PPFregionS{sf1},2));
    for sf2 = 1 : size(PPFregionS{sf1},2)
        tmn = nanmean(mtx2Field{sf1}(PPFregionS{sf1}(3,sf2) : PPFregionS{sf1}(4,sf2),PPFregionS{sf1}(1,sf2) : PPFregionS{sf1}(2,sf2)));
        PFp = PPFregionS{sf1}(1,sf2) + find(tmn == max(tmn));
        if ((PPFregionS{sf1}(1,sf2) >= stem(1)/2 - 1 && PPFregionS{sf1}(1,sf2) <= stem(2)/2 + 1) && (PFp > stem(1)/2 - 1 && PFp < stem(2)/2 + 1)) || ...
                ((PPFregionS{sf1}(2,sf2) >= stem(1)/2 - 1 && PPFregionS{sf1}(2,sf2) <= stem(2)/2 + 1) && (PFp > stem(1)/2 - 1 && PFp < stem(2)/2 + 1))
            stnst2(sf2) = 1;
        end
    end
    stemregions{sf1} = PPFregionS{sf1}(:,stnst2 == 1);
    if sum(stnst2) > 0
        stnst1(sf1) = 1;
    end
end

stemPC_ID = find(stnst1 == 1);
stemregions = stemregions(~cellfun('isempty',stemregions));

PFA_place_fields = struct('Place_fields_in_trials',placesupercell,'FR_threshold',pFR);
PFA_place_field_zones = struct('Plce_field_zones',PFregionS,'PC_with_PF',cellswithfields);
PFA_stem_place_field_zones = struct('Place_cells_with_stem_PF',stemPC_ID,'Stem_PF_regions',stemregions);
filenamePF = [sessionname{sessnum} '_place_fields.mat'];
filenamePFZ = [sessionname{sessnum} '_place_field_zones.mat'];
filenamePFSZ = [sessionname{sessnum} '_place_field_stem_zones.mat'];
save(fullfile(resdir,sessionname{sessnum},filenamePF),'PFA_place_fields')
save(fullfile(resdir,sessionname{sessnum},filenamePFZ),'PFA_place_field_zones')
save(fullfile(resdir,sessionname{sessnum},filenamePFSZ),'PFA_stem_place_field_zones')

%%
%Section 11. Crosscorrelation of tuning curves in the stem arranged by time
disp('Calculating tuning curve correlations...')

%stem segment of smoothed tuning curves with stem PFs

ccStem = mtx2Field(stemPC_ID);
cccell = cell(1,length(ccStem));
for cj1 = 1 : length(ccStem)
    ccSmtx = ccStem{cj1}(:,stem(1)/2 : stem(2)/2);
    [Xc,Xp] = corr(ccSmtx','Type','spearman');
    cccell{cj1} = cat(3,Xc,abs(log10(Xp)));
end

%TZ segment of all smoothed tuning curves
ccTZ = mtx2Field;
ccTZcell = cell(1,length(ccTZ));
for cj1 = 1 : length(ccTZ)
%     disp(cj1)
    ccTZmtx = ccTZ{cj1}(:,1 : stem(1)/2);
    [XcTZ,XpTZ] = corr(ccTZmtx','Type','spearman');
    ccTZcell{cj1} = cat(3,XcTZ,abs(log10(XpTZ)));
end

%off-diagonal correlations: stem
offdimtx = zeros(constrialnum - 1,constrialnum,length(cccell));
offdimeans = zeros(length(cccell),constrialnum);
for oi = 1 : length(cccell)
    odc = NaN(constrialnum - 1,constrialnum - 1);
    tecm = squeeze(cccell{oi}(:,:,1));
    trcell = cell(1,constrialnum);
    for oi1 = 1 : constrialnum - 1
        for oi2 = 1 : constrialnum - oi1
            odc(oi1,oi2 + oi1) = tecm(oi2,oi2 + oi1);
        end
    end
    offdimtx(:,:,oi) = odc;
    offdimeans(oi,:) = nanmean(odc);
end

%off-diagonal correlations: TZ
offdiTZmtx = zeros(constrialnum - 1,constrialnum,length(ccTZcell));
offdiTZmeans = zeros(length(ccTZcell),constrialnum);
for oi = 1 : length(ccTZcell)
    odTZ = NaN(constrialnum - 1,constrialnum - 1);
    teTZm = squeeze(ccTZcell{oi}(:,:,1));
    tzcell = cell(1,constrialnum);
    for oi1 = 1 : constrialnum - 1
        for oi2 = 1 : constrialnum - oi1
            odTZ(oi1,oi2 + oi1) = teTZm(oi2,oi2 + oi1);
        end
    end
    offdiTZmtx(:,:,oi) = odTZ;
    offdiTZmeans(oi,:) = nanmean(odTZ);
end

%Crosscorrelation of tuning curves in the stem arranged by context (left or
%right): only for stem
ccLRcell = cell(1,length(ccStem));
contcatlens = [1,length(purelefts),length(purerights),length(errorlefts),length(errorights)];
catcompmtx = zeros(4,4,length(ccStem));
for cq1 = 1 : length(ccStem)
    LRmtx = [ccStem{cq1}(purelefts,:);ccStem{cq1}(purerights,:);...
        ccStem{cq1}(errorlefts,:);ccStem{cq1}(errorights,:)];
    ccLRcell{cq1} = corr(LRmtx','Type','spearman');
    catcomp = zeros(4,4);
    startc = 0;
    stopc = 0;
    for cq2 = 1 : 4
        startc = startc + contcatlens(cq2);
        stopc = stopc + contcatlens(cq2 + 1);
        startr = 0;
        stopr = 0;
        for cq3 = 1 : 4
            startr = startr + contcatlens(cq3);
            stopr = stopr + contcatlens(cq3 + 1);
            catcomp(cq2,cq3) = mean(nanmean(ccLRcell{cq1}(startc : stopc,startr : stopr)));
        end
    end
    catcompmtx(:,:,cq1) = catcomp;
end

%off-diagonal correlations: only for stem
offdiLRcell = cell(1,length(ccLRcell));
offdiLRmeans = cell(1,length(ccLRcell));
for oi = 1 : length(ccLRcell)
    odclr = NaN(size(ccLRcell{oi},1) - 1,size(ccLRcell{oi},1) - 1);
    teclr = ccLRcell{oi};
    tlrcell = cell(1,size(ccLRcell{oi},1));
    for oi1 = 1 : size(ccLRcell{oi},1) - 1
        for oi2 = 1 : size(ccLRcell{oi},2) - oi1
            odclr(oi1,oi2 + oi1) = teclr(oi2,oi2 + oi1);
        end
    end
    offdiLRcell{oi} = odclr;
    offdiLRmeans{oi} = nanmean(odclr);
end

%within and across context shuffled correlations without the error trials
% purecont005 = zeros(6,length(shuffcorrs));
% for fc1 = 1 : length(shuffcorrs)
%     shpLm = zeros(1,numrep);
%     shpRm = zeros(1,numrep);
%     shpLR = zeros(1,numrep);
%     for fc2 = 1 : size(shuffcorrs{fc1},3)
%         shuffplefts = shuffcorrs{fc1}(purelefts,purelefts,fc2);
%         shuffprights = shuffcorrs{fc1}(purerights,purerights,fc2);
%         shuffpLR = shuffcorrs{fc1}(purelefts,purerights,fc2);
%         shpLm(fc2) = nanmean(nanmean(shuffplefts));
%         shpRm(fc2) = nanmean(nanmean(shuffprights));
%         shpLR(fc2) = nanmean(nanmean(shuffpLR));
%     end
%     purecont005(:,fc1) = [quantile(shpLm,[0.025 0.975])';...
%         quantile(shpRm,[0.025 0.975])';...
%         quantile(shpLR,[0.025 0.975])'];
% end

%within and across context shuffled correlations with the error trials
%(prospective code)
% prospcont005 = zeros(6,length(shuffcorrs));
% for fc3 = 1 : length(shuffcorrs)
%     shLm = zeros(1,numrep);
%     shRm = zeros(1,numrep);
%     shLR = zeros(1,numrep);
%     for fc4 = 1 : size(shuffcorrs{fc3},3)
%         shufflefts = shuffcorrs{fc3}(lefts,lefts,fc4);
%         shuffrights = shuffcorrs{fc3}(rights,rights,fc4);
%         shuffLR = shuffcorrs{fc3}(lefts,rights,fc4);
%         shLm(fc4) = nanmean(nanmean(shufflefts));
%         shRm(fc4) = nanmean(nanmean(shuffrights));
%         shLR(fc4) = nanmean(nanmean(shuffLR));
%     end
%     prospcont005(:,fc3) = [quantile(shLm,[0.025 0.975])';...
%         quantile(shRm,[0.025 0.975])';...
%         quantile(shLR,[0.025 0.975])'];
% end

%within and across context correlations: retrospective and prospective coding
%stem
LRmemcross = zeros(8,length(cccell));
for co3 = 1 : length(cccell)
    withincL = cccell{co3}(purelefts,purelefts);
    withincR = cccell{co3}(purerights,purerights);
    withineL = cccell{co3}(errorlefts,errorlefts);
    withineR = cccell{co3}(errorights,errorights);
    leftmem = cccell{co3}(purerights,errorlefts);
    rightmem = cccell{co3}(purelefts,errorights);
    leftprosp = cccell{co3}(purelefts,errorlefts);
    rightprosp = cccell{co3}(purerights,errorights);
    LRmemcross(1,co3) = nanmean(nanmean(withincL)); %correct left-left correlation
    LRmemcross(2,co3) = nanmean(nanmean(withincR)); %correct right-right correlation
    LRmemcross(3,co3) = nanmean(nanmean(withineL)); %error left-left correlation
    LRmemcross(4,co3) = nanmean(nanmean(withineR)); %error right-right correlation
    LRmemcross(5,co3) = nanmean(nanmean(leftmem)); %retrospective left
    LRmemcross(6,co3) = nanmean(nanmean(rightmem)); %retrospective right
    LRmemcross(7,co3) = nanmean(nanmean(leftprosp)); %left prospective
    LRmemcross(8,co3) = nanmean(nanmean(rightprosp)); %right prospective
end

%TZ
LRTZmemcross = zeros(8,length(ccTZcell));
for co3 = 1 : length(ccTZcell)
    withincLTZ = ccTZcell{co3}(purelefts,purelefts);
    withincRTZ = ccTZcell{co3}(purerights,purerights);
    withineLTZ = ccTZcell{co3}(errorlefts,errorlefts);
    withineRTZ = ccTZcell{co3}(errorights,errorights);
    leftmemTZ = ccTZcell{co3}(purerights,errorlefts);
    rightmemTZ = ccTZcell{co3}(purelefts,errorights);
    leftprospTZ = ccTZcell{co3}(purelefts,errorlefts);
    rightprospTZ = ccTZcell{co3}(purerights,errorights);
    LRTZmemcross(1,co3) = nanmean(nanmean(withincLTZ)); %correct left-left correlation
    LRTZmemcross(2,co3) = nanmean(nanmean(withincRTZ)); %correct right-right correlation
    LRTZmemcross(3,co3) = nanmean(nanmean(withineLTZ)); %error left-left correlation
    LRTZmemcross(4,co3) = nanmean(nanmean(withineRTZ)); %error right-right correlation
    LRTZmemcross(5,co3) = nanmean(nanmean(leftmemTZ)); %retrospective left
    LRTZmemcross(6,co3) = nanmean(nanmean(rightmemTZ)); %retrospective right
    LRTZmemcross(7,co3) = nanmean(nanmean(leftprospTZ)); %left prospective
    LRTZmemcross(8,co3) = nanmean(nanmean(rightprospTZ)); %right prospective
end

PFA_tuningCurve_correlations = struct('TC_stem_correlations',cccell,'Off_diagonal_stem_correlations',offdimtx,'Off_diagonal_mean_stem_correlations',offdimeans,...
    'Off_diagonal_TZ_correlations',offdiTZmtx,'Off_diagonal_TZ_mean_correlations',offdiTZmeans,'Trail_by_trial_mean_correlations',catcompmtx,...
    'Off_diagonal_LeftRight_correlations',offdiLRcell,'Off_diagonal_LeftRight_mean_correlations',offdiLRmeans,...
    'RetroProsp_stem_correlations',LRmemcross,'RetroProsp_TZ_correlations',LRTZmemcross);
PFA_tuningCurve_TZ_correlations = struct('TC_TZ_correlations',ccTZcell);
filenameTCx = [sessionname{sessnum} '_tuningCurve_correlations.mat'];
filenameTZTCx = [sessionname{sessnum} '_TZ_tuningCurve_correlations.mat'];
save(fullfile(resdir,sessionname{sessnum},filenameTCx),'PFA_tuningCurve_correlations')
save(fullfile(resdir,sessionname{sessnum},filenameTZTCx),'PFA_tuningCurve_TZ_correlations')

%%
%PF centers
disp('PF center position...')

bigPFpeakcell = cell(1,length(stemregions));
PFposcell = cell(1,length(stemregions));
for pf1 = 1 : length(stemregions)
    sourcePFmtx = mtx2Field{stemPC_ID(pf1)}(:,1 : 33); %33 = stem + 1 bin
    controlPFmtx = mtx2Field{stemPC_ID(pf1)}(:,1 : 33); %33 = stem + 1 bin
    for fp1 = 1 : size(stemregions{pf1},2)
        for fp2 = 1 : size(controlPFmtx,1)
            controlPFmtx(fp2,stemregions{pf1}(1,fp1) : stemregions{pf1}(2,fp1)) = NaN;
        end
    end
    PFpcell = cell(1,size(stemregions{pf1},2));
    PFLRpos = zeros(2,size(stemregions{pf1},2));
    for pf2 = 1 : size(stemregions{pf1},2)
        PFLRpeak = zeros(2,length(purelefts));
        for pf3 = 1 : length(purelefts)
            controlL = nanmean(controlPFmtx(purelefts(pf3),:));
            if controlL == 0
                controlL = 0.1;
            end
            controlR = nanmean(controlPFmtx(purerights(pf3),:));
            if controlR == 0
                controlR = 0.1;
            end
            startune = stemregions{pf1}(1,pf2);
            stoptune = stemregions{pf1}(2,pf2);
            if  stoptune > 33
               stoptune = 33;
            end
            ttunesegmL = sourcePFmtx(purelefts(pf3),startune : stoptune);
            ttunesegmR = sourcePFmtx(purerights(pf3),startune : stoptune);
            tpfL = find(ttunesegmL == max(ttunesegmL)) + startune;
            tpfR = find(ttunesegmR == max(ttunesegmR)) + startune;
            PFLRpeak(1,pf3) = tpfL(1);
            PFLRpeak(2,pf3) = tpfR(1);
            checkLpeak = ttunesegmL(ttunesegmL == max(ttunesegmL));
            checkRpeak = ttunesegmR(ttunesegmR == max(ttunesegmR));
            if checkLpeak / controlL < 2
                PFLRpeak(1,pf3) = NaN;
            end
            if checkRpeak / controlR < 2
                PFLRpeak(2,pf3) = NaN;
            end
        end
        PFpcell{pf2} = PFLRpeak;
        PFLRpos(1,pf2) = nanmean(PFLRpeak(1,:));
        PFLRpos(2,pf2) = nanmean(PFLRpeak(2,:));
    end
    bigPFpeakcell{pf1} = PFpcell;
    PFposcell{pf1} = PFLRpos;
end

PFA_PF_center_position = struct('All_PF_peaks',bigPFpeakcell,'PF_peak_mean_positions',PFposcell);
filename_PFCP = [sessionname{sessnum} '_PC_PF_center_position.mat'];
save(fullfile(resdir,sessionname{sessnum},filename_PFCP),'PFA_PF_center_position')

%%
%MUA Firing rate variability
alltrialisi = cell(1,constrialnum);
alltrialpos = cell(1,constrialnum);
isinum = zeros(1,constrialnum);
posnum = zeros(1,constrialnum);
for vj1 = 1 : constrialnum
    pertrialisi = cell(1,length(unitDP_full));
    pertrialpos = cell(1,length(unitDP_full));
    for vj2 = 1 : length(unitDP_full)
        pertrialisi{vj2} = unitDP_full{vj2}{vj1};
        pertrialpos{vj2} = unitarm_fullin{vj2}{vj1};
        pertrialpos{vj2} = pertrialpos{vj2}(2:end);
    end
    pertrialisi = pertrialisi(~cellfun('isempty',pertrialisi));
    pertrialpos = pertrialpos(~cellfun('isempty',pertrialpos));
    alltrialisi{vj1} = cat(1,pertrialisi{:});
    alltrialpos{vj1} = cat(1,pertrialpos{:});
    isinum(vj1) = length(alltrialisi{vj1});
    posnum(vj1) = length(alltrialpos{vj1});
    if isinum(vj1) - posnum(vj1) > 0
        alltrialisi{vj1}(randperm(isinum(vj1),isinum(vj1) - posnum(vj1))) = [];
    end
    if isinum(vj1) - posnum(vj1) < 0
        alltrialpos{vj1}(randperm(isinum(vj1),abs(isinum(vj1) - posnum(vj1)))) = [];
    end
end
allallpos = cat(1,alltrialpos{:});
allallisi = cat(1,alltrialisi{:});
[av1,av2,av3] = histcounts(allallpos',spatfullbins);
spatspkvar = cell(1,length(spatfullbins) - 1);
ssv = zeros(1,length(spatfullbins) - 1);
for vj3 = 1 : length(spatfullbins)
    spatspkvar{vj3} = allallisi(av3 == vj3);
    ssv(vj3) = std(spatspkvar{vj3},'omitnan').^2;
end

PFA_MUA_FR_var = struct('MUA_FR_variability',ssv);
filename_MUA_var = [sessionname{sessnum} '_MUA_FR_var.mat'];
save(fullfile(resdir,sessionname{sessnum},filename_MUA_var),'PFA_MUA_FR_var')

%%
%Correlation of pyramidal cell and interneuron tuning curves
pleft = purelefts;
pright = purerights;
if length(pleft) == length(pright) + 1
    pleft(1) = [];
elseif length(pright) == length(pleft) + 1
    pright(1) = [];
end

%full
%trial-by-trial
INxPtcMcell = cell(1,size(upINspatcells,2));
INxPtcMmtx = zeros(4,length(mtx2Field),size(upINspatcells,2));
for ip1 = 1 : size(upINspatcells,2)
    INxPtc = cell(1,length(mtx2Field));
    INxPtcv = zeros(4,length(mtx2Field));
    for ip2 = 1 : length(mtx2Field)
        txtc = corr(upINspatcells{4,ip1}',mtx2Field{ip2}','Type','spearman');
        txtc1 = txtc(pleft,pleft);
        txtc2 = txtc(pright,pright);
        txtc3 = txtc(pleft,pright);
        txtc4 = txtc(pright,pleft);
        INxPtc{ip2} = cat(3,txtc1,txtc2,txtc3,txtc4);
        INxPtcv(:,ip2) = [nanmean(nanmean(txtc1));nanmean(nanmean(txtc2));...
            nanmean(nanmean(txtc3));nanmean(nanmean(txtc4))];
    end
    INxPtcMcell{ip1} = INxPtc;
    INxPtcMmtx(:,:,ip1) = INxPtcv;
end

%bin-by-bin
INxPloc = cell(1,size(upINspatcells,2));
for ip3 = 1 : size(upINspatcells,2)
    ixpRP = zeros(size(mtx2Field{ip3},2),2,size(mtx2Field,2));
    for ip4 = 1 : size(mtx2Field,2)
        [ixpR,ixpP] = corr(upINspatcells{4,ip3},mtx2Field{ip4},'Type','spearman');
        ixpRP(:,:,ip4) = [diag(ixpR),diag(ixpP)];
    end
    INxPloc{ip3} = ixpRP;
end

PFA_PC_IN_TC_correlations = struct('Trial_by_trial_PC_IN_correlations',INxPtcMcell,'Trial_by_trial_PC_IN_mean_correlations',INxPtcMmtx,...
    'Bin_by_bin_PC_IN_correlations',INxPloc);
filenamePCxIN = [sessionname{sessnum} '_PC_IN_TC_correlations.mat'];
save(fullfile(resdir,sessionname{sessnum},filenamePCxIN),'PFA_PC_IN_TC_correlations')

%%
%Difference of max-normalized IN and PC tuning curves
disp('Calculating PC-IN TC difference...')

INpPCmc = cell(1,size(upINspatcells,2));
for sp1 = 1 : size(upINspatcells,2)
    trialaveIN = max(upINspatcells{4,sp1},[],2);
    mnINmtx = upINspatcells{4,sp1} ./ repmat(trialaveIN,1,size(upINspatcells{4,sp1},2));
    mnINmtxL = mnINmtx(purelefts,:);
    mnINmtxR = mnINmtx(purerights,:);
    INpPCmtx = zeros(2,size(mnINmtx,2),length(mtx2Field));
    for sp2 = 1 : length(mtx2Field)
        trialavePC = max(mtx2Field{sp2},[],2);
        mnPCmtx = mtx2Field{sp2} ./ repmat(trialavePC,1,size(mtx2Field{sp2},2));
        mnPCmtxL = mnPCmtx(purelefts,:);
        mnPCmtxR = mnPCmtx(purerights,:);
        INpPCL = zeros(1,size(mnINmtx,2));
        INpPCR = zeros(1,size(mnINmtx,2));
        for sp3 = 1 : size(mnINmtx,2)
            compINvL = mnINmtxL(:,sp3);
            compPCvL = mnPCmtxL(:,sp3);
            compINvR = mnINmtxR(:,sp3);
            compPCvR = mnPCmtxR(:,sp3);
            compINvL(isnan(compINvL)) = [];
            compPCvL(isnan(compPCvL)) = [];
            compINvR(isnan(compINvR)) = [];
            compPCvR(isnan(compPCvR)) = [];
            if ~isempty(compINvL) && ~isempty(compPCvL)
                INpPCL(sp3) = abs(log10(ranksum(compINvL,compPCvL)));
            end
            if ~isempty(compINvR) && ~isempty(compPCvR)
                INpPCR(sp3) = abs(log10(ranksum(compINvR,compPCvR)));
            end
        end
        INpPCmtx(:,:,sp2) = [INpPCL;INpPCR];
    end
    INpPCmc{sp1} = INpPCmtx;
end

PFA_PC_IN_TC_difference = struct('PC_IN_TC_difference',INpPCmc);
filenamedPCIN = [sessionname{sessnum} '_PC_IN_TC_difference.mat'];
save(fullfile(resdir,sessionname{sessnum},filenamedPCIN),'PFA_PC_IN_TC_difference')

%%
%template matching-based decoding using cosine similarity for correlation
disp('Decoding position by template matching...')

L1O = 'yes';
SoFM = 'f';
how2DCm = {'P','L','R'};
propm = [0.5 1 1];
titp = InstantaneousTheta.signal_phase;

%Choosing the trial set for building the spatial PV template resulting the
%smallest decoding error
if length(purelefts) >= 10 && length(purerights) >= 10
    decrepeats = 100;
elseif length(purelefts) >= 5 && length(purerights) >= 5
    decrepeats = 20;
else
    decrepeats = 1;
end

mainDE = zeros(length(how2DCm),decrepeats);
fChoT = cell(1,length(how2DCm));
for dq1 = 1 : length(how2DCm)
    how2DC = how2DCm{dq1};
    seltimes = times;
    selphase = spikephase;
    selspktimes = unitPTP_full;
    selvects = fullFRvects2;
    preselvects = fullFRvects2;
    ChoT = cell(1,decrepeats);
    for dq2 = 1 : decrepeats
        disp([how2DC,': ',num2str(dq2)])
        SoFM = 'f';
        celltype = 'PC';
        allChoTrials = NaN;
        LChoTrials = NaN;
        RChoTrials = NaN;
        [PCT,PCP,rearrthetatimes,firstspiketimes,ofst,thetatimeorders,percyclespknum,...
        ttarm,alltempos,trialcycles,posincycles,poscycleinbins,thetaspeed,thetavectors,...
        spatempmtx,spatemp,ethP_CS,de_CS,nnde_CS,qde_CS,estPv_CS,realposvect,timeposdistr,...
        randAlltrials,randLtrials,randRtrials] = placeTMdecoder(SoFM,how2DC,alltime,titp,seltimes,selphase,selspktimes,selvects,preselvects,...
        constrialnum,purelefts,purerights,alltrialtimes,alltrials,trialspeeds,spatbins2,spatfullbins,allChoTrials,LChoTrials,RChoTrials);
        mainDE(dq1,dq2) = median(abs(cat(2,de_CS{:})),'omitnan');
        if strcmp(how2DC,'P')
            ChoT{dq2} = randAlltrials;
        end
        if strcmp(how2DC,'L')
            ChoT{dq2} = randLtrials;
        end
        if strcmp(how2DC,'R')
            ChoT{dq2} = randRtrials;
        end
    end
    fChoT{dq1} = ChoT(mainDE(dq1,:) == min(mainDE(dq1,:)));
end

%running the decoding using the spatial template composed of the
%aforechosen trials

celltype = {'PC','IN','All'};
for mj = 1 : length(how2DCm)
PCTc = cell(1,3);
PCPc = cell(1,3);
Ratt = cell(1,3);
Fspktt = cell(1,3);
ofstc = cell(1,3);
TTO = cell(1,3);
pCspkn = cell(1,3);
tCyc = cell(1,3);
PinC = cell(1,3);
PinCb = cell(1,3);
Thspeed = cell(1,3);
Thvect = cell(1,3);
esThpos = cell(1,3);
DecErr = cell(1,3);
nnDecErr = cell(1,3);
qDecErr = cell(1,3);
ePv = cell(1,3);
rPv = cell(1,3);
tPd = cell(1,3);
if strcmp(how2DCm{mj},'P')
    Cht = fChoT{mj};
    allChoTrials = Cht{1};
end
if strcmp(how2DCm{mj},'L')
    Cht = fChoT{mj};
    LChoTrials = Cht{1};
end
if strcmp(how2DCm{mj},'R')
    Cht = fChoT{mj};
    RChoTrials = Cht{1};
end
for ck = 1 : length(celltype)
    disp(celltype{ck})
    cT = celltype{ck};
    SoFM = 'f';
    how2DC = how2DCm{mj};
    prop = propm(mj);
    if strcmp(cT,'PC')
        seltimes = times;
        selphase = spikephase;
        selspktimes = unitPTP_full;
        selvects = fullFRvects2;
        preselvects = fullFRvects2;
    elseif strcmp(cT,'IN')
        seltimes = timeIN;
        selphase = spikephaseIN;
        selspktimes = INtfull;
        preselvects = infullFRvect;
        for oi = 1 : length(preselvects)
            preselvects{oi} = preselvects{oi}';
        end
        selvects = preselvects;
    elseif strcmp(cT,'All')
        seltimes = cat(2,times,timeIN);
        selphase = cat(2,spikephase,spikephaseIN);
        selspktimes = cat(2,unitPTP_full,INtfull);
        prepre = infullFRvect;
        for oj = 1 : length(prepre)
            prepre{oj} = prepre{oj}';
        end
        selvects = cat(2,fullFRvects2,prepre);
        preselvects = selvects;
    end
    [PCT,PCP,rearrthetatimes,firstspiketimes,ofst,thetatimeorders,percyclespknum,...
        ttarm,alltempos,trialcycles,posincycles,poscycleinbins,thetaspeed,thetavectors,...
        spatempmtx,spatemp,ethP_CS,de_CS,nnde_CS,qde_CS,estPv_CS,realposvect,timeposdistr,...
        randAlltrials,randLtrials,randRtrials] = placeTMdecoder(SoFM,how2DC,alltime,titp,seltimes,selphase,selspktimes,selvects,preselvects,...
        constrialnum,purelefts,purerights,alltrialtimes,alltrials,trialspeeds,spatbins2,spatfullbins,allChoTrials,LChoTrials,RChoTrials);
    PCTc{ck} = PCT;
    PCPc{ck} = PCP;
    Ratt{ck} = rearrthetatimes;
    Fspktt{ck} = firstspiketimes;
    ofstc{ck} = ofst;
    TTO{ck} = thetatimeorders;
    pCspkn{ck} = percyclespknum;
    tCyc{ck} = trialcycles;
    PinC{ck} = posincycles;
    PinCb{ck} = poscycleinbins;
    Thspeed{ck} = thetaspeed;
    Thvect{ck} = thetavectors;
%     esThpos{ck} = esthetapos;
    DecErr{ck} = de_CS;
    nnDecErr{ck} = nnde_CS;
    qDecErr{ck} = qde_CS;
    ePv{ck} = estPv_CS;
    rPv{ck} = realposvect;
    tPd{ck} = timeposdistr;
    PFA_decoding = struct('PerCycleTime',PCTc,'PerCyclePhase',PCPc,'RearrangedThetaTimes',Ratt,'FirstSpikeInCycleTimes',Fspktt,...
        'SortedFirstSpikeTimes',ofstc,'OrderOfSpikesInCycles',TTO,'PerCycleSpikeNumber',pCspkn,'CycleTrialIdentity',tCyc,...
        'PositionInCycles',PinC,'PositionBinInCycles',PinCb,'SpeedInCycles',Thspeed,'ThetaCyclePopVectors',Thvect,...
        'DecodingError',DecErr,'NaNfreeDecodingError',nnDecErr,'QuantileOfDecodingError',qDecErr,...
        'DecodedPosition',ePv,'RealPosition',rPv,'PositionDistribution',tPd);
end
save([resdirf '\' 'PFA_' how2DCm{mj} '_decodingresults.mat'],'PFA_decoding')
end

%%
%contribution of each neuron to spatial coding (leave one out): here, the
%spatial template is built from the trials determined in the previous
%section (by minimizing decoding error)
disp('Contribution of each neuron to the population place code...')

SoFM = 'f';
seltimes = cat(2,times,timeIN);
selphase = cat(2,spikephase,spikephaseIN);
selspktimes = cat(2,unitPTP_full,INtarm);
prepre = infullFRvect;
for oj = 1 : length(prepre)
    prepre{oj} = prepre{oj}';
end
selvects = cat(2,fullFRvects2,prepre);
DCerr = cell(1,length(seltimes));
NNDCerr = cell(1,length(seltimes));
qDCerr = cell(1,length(seltimes));
EPv = cell(1,length(seltimes));
RPv = cell(1,length(seltimes));
contrib = zeros(1,length(seltimes));
for di = 1 : length(how2DCm)
    how2DC = how2DCm{di};
    load([resdirf '\' 'PFA_' how2DCm{di} '_decodingresults.mat']);
    DErr = PFA_decoding(1).DecodingError;
    if strcmp(how2DCm{di},'P')
        Cht = fChoT{di};
        allChoTrials = Cht{1};
    end
    if strcmp(how2DCm{di},'L')
        Cht = fChoT{di};
        LChoTrials = Cht{1};
    end
    if strcmp(how2DCm{di},'R')
        Cht = fChoT{di};
        RChoTrials = Cht{1};
    end
    for l1 = 1 : length(seltimes)
        disp([how2DCm{di},'_contribution_of_cell_',num2str(l1)])
        rseltimes = seltimes;
        rseltimes(l1) = [];
        rselphase = selphase;
        rselphase(l1) = [];
        rselspktimes = selspktimes;
        rselspktimes(l1) = [];
        rselvects = selvects;
        rselvects(l1) = [];
        rpreselvects = rselvects;
        [PCT,PCP,rearrthetatimes,firstspiketimes,ofst,thetatimeorders,percyclespknum,...
    ttarm,alltempos,trialcycles,posincycles,poscycleinbins,thetaspeed,thetavectors,...
    spatempmtx,spatemp,ethP_CS,de_CS,nnde_CS,qde_CS,estPv_CS,realposvect,timeposdistr] =...
    placeTMdecoderR(SoFM,how2DC,alltime,titp,rseltimes,rselphase,rselspktimes,rpreselvects,...
    constrialnum,alltrialtimes,alltrials,trialspeeds,spatbins2,spatfullbins,purelefts,purerights,allChoTrials,LChoTrials,RChoTrials);
        DCerr{l1} = de_CS;
        NNDCerr{l1} = nnde_CS;
        qDCerr{l1} = qde_CS;
        EPv{l1} = ethP_CS;
        RPv{l1} = realposvect;
        contrib(l1) = (mean(abs(cat(2,DCerr{l1}{:})),'omitnan') - mean(abs(cat(2,DErr{:})),'omitnan')) / (mean(abs(cat(2,DCerr{l1}{:})),'omitnan') + mean(abs(cat(2,DErr{:})),'omitnan'));
    end
    PFA_contribution = struct('DecodingError',DCerr,'NonNanDecodingError',NNDCerr,'QuantDecodingError',qDCerr,'EstimatedPosition',EPv,'RealPosition',RPv,'Contribution',contrib);
    save([resdirf '\' 'PFA_' how2DCm{di} '_contrib2coding.mat'],'PFA_contribution')
end

%%
%Splitter analysis: con't from before the population analysis

%Decide if a cell is correlated more with the future or with the past based
%on TC difference
disp('Analyzing context-dependent firing...')

choicecorrmtx = zeros(4,constrialnum,length(mtx2Field));
choiceBasedTemplates = zeros(4,length(spatbins2),length(mtx2Field));
choices = zeros(length(mtx2Field),constrialnum,2); %1: correct left, 2: correct right, 3: error left, 4: error right
choiceprop = zeros(length(mtx2Field),4);
covarcorr = zeros(1,length(mtx2Field));
for sm1 = 1 : length(mtx2Field)
    stemTC = mtx2Field{sm1}(:,stem(1)/2 : stem(2)/2);
    if length(errorlefts) >= 5 && length(errorights) >= 5
        choiceBasedTemplates(:,:,sm1) = [mean(stemTC(purelefts,:));...
        mean(stemTC(purerights,:));mean(stemTC(errorlefts,:));...
        mean(stemTC(errorights,:))];
    elseif length(errorlefts) >= 5 && length(errorights) < 5
        choiceBasedTemplates(:,:,sm1) = [mean(stemTC(purelefts,:));...
        mean(stemTC(purerights,:));mean(stemTC(errorlefts,:));...
        mean(stemTC([],:))];
    elseif length(errorlefts) < 5 && length(errorights) >= 5
        choiceBasedTemplates(:,:,sm1) = [mean(stemTC(purelefts,:));...
        mean(stemTC(purerights,:));mean(stemTC([],:));...
        mean(stemTC(errorights,:))];
    elseif length(errorlefts) < 5 && length(errorights) < 5
        choiceBasedTemplates(:,:,sm1) = [mean(stemTC(purelefts,:));...
        mean(stemTC(purerights,:));mean(stemTC([],:));...
        mean(stemTC([],:))];
    end
    tpltcorr = zeros(4,constrialnum);
    for sm2 = 1 : constrialnum
        pttplt = stemTC(sm2,:);
        tpltcorr(1,sm2) = corr(pttplt',choiceBasedTemplates(1,:,sm1)','type','spearman');
        tpltcorr(2,sm2) = corr(pttplt',choiceBasedTemplates(2,:,sm1)','type','spearman');
        tpltcorr(3,sm2) = corr(pttplt',choiceBasedTemplates(3,:,sm1)','type','spearman');
        tpltcorr(4,sm2) = corr(pttplt',choiceBasedTemplates(4,:,sm1)','type','spearman');
        tcm = find(tpltcorr(:,sm2) == max(tpltcorr(:,sm2)));
        if ~isempty(tcm)
            tcm = tcm(1);
            choices(sm1,sm2,1) = tpltcorr(tcm,sm2);
            choices(sm1,sm2,2) = tcm;
        else
            choices(sm1,sm2,1) = NaN;
            choices(sm1,sm2,2) = NaN;
        end
    end
    choicecorrmtx(:,:,sm1) = tpltcorr;
    choiceprop(sm1,1) = length(choices(sm1,choices(sm1,:,2) == 1,2)) / constrialnum;
    choiceprop(sm1,2) = length(choices(sm1,choices(sm1,:,2) == 2,2)) / constrialnum;
    choiceprop(sm1,3) = length(choices(sm1,choices(sm1,:,2) == 3,2)) / constrialnum;
    choiceprop(sm1,4) = length(choices(sm1,choices(sm1,:,2) == 4,2)) / constrialnum;
    covarcorr(sm1) = nanstd(choices(sm1,:,1)) / nanmean(choices(sm1,:,1));
end

%bin-by-bin tuning curve difference: choice-based without error trials in
%stem
stemTCc = mtx2Field(stemPC_ID);
eoLRd = zeros(length(stemTCc),length(spatbins2));
for bj1 = 1 : length(stemTCc)
    stemTCmtx = stemTCc{bj1}(:,stem(1)/2 : stem(2)/2);
    onlyleft = stemTCmtx(purelefts,:);
    onlyright = stemTCmtx(purerights,:);
    extofdiff = zeros(1,size(stemTCmtx,2));
    for bj2 = 1 : length(spatbins2)
        extofdiff(bj2) = abs(log10(ranksum(onlyleft(:,bj2),onlyright(:,bj2))));
    end
    eoLRd(bj1,:) = extofdiff;
end

%in TZ
tzTCc = mtx2Field;
tzLRd = zeros(length(tzTCc),10);
for bj1 = 1 : length(tzTCc)
    tzTCmtx = tzTCc{bj1}(:,1 : stem(1)/2);
    tzoleft = tzTCmtx(purelefts,:);
    tzoright = tzTCmtx(purerights,:);
    tzofdiff = zeros(1,size(tzTCmtx,2));
    for bj2 = 1 : size(tzTCmtx,2)
        tzofdiff(bj2) = abs(log10(ranksum(tzoleft(:,bj2),tzoright(:,bj2))));
    end
    tzLRd(bj1,:) = tzofdiff;
end

%bin-by-bin tuning curve difference of INs: choice-based without error
%trials in stem
inLRd = zeros(size(upINspatcells,2),length(spatbins2));
for ib1 = 1 : size(upINspatcells,2)
    instem = upINspatcells{4,ib1}(:,stem(1)/2 : stem(2)/2);
    inleft = instem(purelefts,:);
    inright = instem(purerights,:);
    inLRdi = zeros(1,size(instem,2));
    for ib2 = 1 : size(instem,2)
        inLRdi(ib2) = abs(log10(ranksum(inleft(:,ib2),inright(:,ib2))));
    end
    inLRd(ib1,:) = inLRdi;
end

%in TZ
intzLRd = zeros(size(upINspatcells,2),10);
for ib1 = 1 : size(upINspatcells,2)
    inTZ = upINspatcells{4,ib1}(:,1 : stem(1)/2);
    intzleft = inTZ(purelefts,:);
    intzright = inTZ(purerights,:);
    intzLRdi = zeros(1,size(inTZ,2));
    for ib2 = 1 : size(inTZ,2)
        intzLRdi(ib2) = abs(log10(ranksum(intzleft(:,ib2),intzright(:,ib2))));
    end
    intzLRd(ib1,:) = intzLRdi;
end

%bin-by-bin stem TC difference of PCs and INs: only L-R error trials
eoerrLRd = zeros(length(stemTCc),length(spatbins2));
eoerrINLRd = zeros(size(upINspatcells,2),length(spatbins2));
if length(errorlefts) >= 5 && length(errorights) >= 5
for bj1 = 1 : length(stemTCc)
    stemTCmtx = stemTCc{bj1}(:,stem(1)/2 : stem(2)/2);
    errleft = stemTCmtx(errorlefts,:);
    erright = stemTCmtx(errorights,:);
    errtofdiff = zeros(1,size(stemTCmtx,2));
    for bj2 = 1 : size(stemTCmtx,2)
        errtofdiff(bj2) = abs(log10(ranksum(errleft(:,bj2),erright(:,bj2))));
    end
    eoerrLRd(bj1,:) = errtofdiff;
end
for bj1i = 1 : size(upINspatcells,2)
    instem = upINspatcells{4,bj1i}(:,stem(1)/2 : stem(2)/2);
    errINleft = instem(errorlefts,:);
    errINright = instem(errorights,:);
    errINdiff = zeros(1,size(instem,2));
    for bj2i = 1 : size(instem,2)
        errINdiff(bj2i) = abs(log10(ranksum(errINleft(:,bj2i),errINright(:,bj2i))));
    end
    eoerrINLRd(bj1i,:) = errINdiff;
end
else
    eoerrLRd = NaN(size(stemTCc,2),length(spatbins2));
    eoerrINLRd = NaN(size(upINspatcells,2),length(spatbins2));
end

%bin-by-bin TZ TC difference of PCs and INs: only L-R error trials
eoerrtzLRd = zeros(length(stemTCc),10);
eoerrINtzLRd = zeros(size(upINspatcells,2),10);
if length(errorlefts) >= 5 && length(errorights) >= 5
for bj1 = 1 : length(stemTCc)
    tzTCmtx = stemTCc{bj1}(:,1 : stem(1)/2);
    errtzleft = tzTCmtx(errorlefts,:);
    errtzright = tzTCmtx(errorights,:);
    errtztofdiff = zeros(1,size(tzTCmtx,2));
    for bj2 = 1 : size(tzTCmtx,2)
        errtztofdiff(bj2) = abs(log10(ranksum(errtzleft(:,bj2),errtzright(:,bj2))));
    end
    eoerrtzLRd(bj1,:) = errtztofdiff;
end
for bj1i = 1 : size(upINspatcells,2)
    inTZ = upINspatcells{4,bj1i}(:,1 : stem(1)/2);
    errtzINleft = inTZ(errorlefts,:);
    errtzINright = inTZ(errorights,:);
    errtzINdiff = zeros(1,size(inTZ,2));
    for bj2i = 1 : size(inTZ,2)
        errtzINdiff(bj2i) = abs(log10(ranksum(errtzINleft(:,bj2i),errtzINright(:,bj2i))));
    end
    eoerrINtzLRd(bj1i,:) = errtzINdiff;
end
else
    eoerrtzLRd = NaN(size(stemTCc,2),size(stemTCc{1},2));
    eoerrINtzLRd = NaN(size(upINspatcells,2),size(upINspatcells{4,1},2));    
end

%bin-by-bin stem TC difference: L-R error vs. correct trials: PCs
correrrLd = zeros(length(stemTCc),size(spatbins2,2));
correrrRd = zeros(length(stemTCc),size(spatbins2,2));
correrrLRd = zeros(length(stemTCc),size(spatbins2,2));
correrrRLd = zeros(length(stemTCc),size(spatbins2,2));
for bj1 = 1 : length(stemTCc)
    stemTCmtx = stemTCc{bj1}(:,stem(1)/2 : stem(2)/2);
    onlyleft = stemTCmtx(purelefts,:);
    onlyright = stemTCmtx(purerights,:);
    errleft = stemTCmtx(errorlefts,:);
    erright = stemTCmtx(errorights,:);
    correrrLdiff = zeros(1,size(stemTCmtx,2));
    correrrRdiff = zeros(1,size(stemTCmtx,2));
    correrrLRdiff = zeros(1,size(stemTCmtx,2));
    correrrRLdiff = zeros(1,size(stemTCmtx,2));
    for bj2 = 1 : size(stemTCmtx,2)
        if length(errorlefts) >= 5
            correrrLdiff(bj2) = abs(log10(ranksum(onlyleft(:,bj2),errleft(:,bj2))));
        else
            correrrLdiff(bj2) = NaN;
        end
        if length(errorights) >= 5
            correrrRdiff(bj2) = abs(log10(ranksum(onlyright(:,bj2),erright(:,bj2))));
        else
            correrrRdiff(bj2) = NaN;
        end
        if length(errorights) >= 5
            correrrLRdiff(bj2) = abs(log10(ranksum(onlyleft(:,bj2),erright(:,bj2))));
        else
            correrrLRdiff(bj2) = NaN;
        end
        if length(errorlefts) >= 5
            correrrRLdiff(bj2) = abs(log10(ranksum(onlyright(:,bj2),errleft(:,bj2))));
        else
            correrrRLdiff(bj2) = NaN;
        end
    end
    correrrLd(bj1,:) = correrrLdiff;
    correrrRd(bj1,:) = correrrRdiff;
    correrrLRd(bj1,:) = correrrLRdiff;
    correrrRLd(bj1,:) = correrrRLdiff;
end

%INs
correrrINLd = zeros(size(upINspatcells,2),length(spatbins2));
correrrINRd = zeros(size(upINspatcells,2),length(spatbins2));
correrrINLRd = zeros(size(upINspatcells,2),length(spatbins2));
correrrINRLd = zeros(size(upINspatcells,2),length(spatbins2));
for bj1 = 1 : size(upINspatcells,2)
    instem = upINspatcells{4,bj1}(:,stem(1)/2 : stem(2)/2);
    onlyleftIN = instem(purelefts,:);
    onlyrightIN = instem(purerights,:);
    errleftIN = instem(errorlefts,:);
    errightIN = instem(errorights,:);
    correrrLdiffIN = zeros(1,size(instem,2));
    correrrRdiffIN = zeros(1,size(instem,2));
    correrrLRdiffIN = zeros(1,size(instem,2));
    correrrRLdiffIN = zeros(1,size(instem,2));
    for bj2 = 1 : size(instem,2)
        if length(errorlefts) >= 5
            correrrLdiffIN(bj2) = abs(log10(ranksum(onlyleftIN(:,bj2),errleftIN(:,bj2))));
        else
            correrrLdiffIN(bj2) = NaN;
        end
        if length(errorights) >= 5
            correrrRdiffIN(bj2) = abs(log10(ranksum(onlyrightIN(:,bj2),errightIN(:,bj2))));
        else
            correrrRdiffIN(bj2) = NaN;
        end
        if length(errorights) >= 5
            correrrLRdiffIN(bj2) = abs(log10(ranksum(onlyleftIN(:,bj2),errightIN(:,bj2))));
        else
            correrrLRdiffIN(bj2) = NaN;
        end
        if length(errorlefts) >= 5
            correrrRLdiffIN(bj2) = abs(log10(ranksum(onlyrightIN(:,bj2),errleftIN(:,bj2))));
        else
            correrrRLdiffIN(bj2) = NaN;
        end
    end
    correrrINLd(bj1,:) = correrrLdiffIN;
    correrrINRd(bj1,:) = correrrRdiffIN;
    correrrINLRd(bj1,:) = correrrLRdiffIN;
    correrrINRLd(bj1,:) = correrrRLdiffIN;
end

%bin-by-bin TZ TC difference: L-R error vs. correct trials: PCs
correrrtzLd = zeros(length(stemTCc),10);
correrrtzRd = zeros(length(stemTCc),10);
correrrtzLRd = zeros(length(stemTCc),10);
correrrtzRLd = zeros(length(stemTCc),10);
for bj1 = 1 : length(stemTCc)
    tzTCmtx = stemTCc{bj1}(:,1 : stem(1)/2);
    onlytzleft = tzTCmtx(purelefts,:);
    onlytzright = tzTCmtx(purerights,:);
    errtzleft = tzTCmtx(errorlefts,:);
    errtzright = tzTCmtx(errorights,:);
    correrrtzLdiff = zeros(1,size(tzTCmtx,2));
    correrrtzRdiff = zeros(1,size(tzTCmtx,2));
    correrrtzLRdiff = zeros(1,size(tzTCmtx,2));
    correrrtzRLdiff = zeros(1,size(tzTCmtx,2));
    for bj2 = 1 : size(tzTCmtx,2)
        if length(errorlefts) >= 5
            correrrtzLdiff(bj2) = abs(log10(ranksum(onlytzleft(:,bj2),errtzleft(:,bj2))));
        else
            correrrtzLdiff(bj2) = NaN;
        end
        if length(errorights) >= 5
            correrrtzRdiff(bj2) = abs(log10(ranksum(onlytzright(:,bj2),errtzright(:,bj2))));
        else
            correrrtzRdiff(bj2) = NaN;
        end
        if length(errorights) >= 5
            correrrtzLRdiff(bj2) = abs(log10(ranksum(onlytzleft(:,bj2),errtzright(:,bj2))));
        else
            correrrtzLRdiff(bj2) = NaN;
        end
        if length(errorlefts) >= 5
            correrrtzRLdiff(bj2) = abs(log10(ranksum(onlytzright(:,bj2),errtzleft(:,bj2))));
        else
            correrrtzRLdiff(bj2) = NaN;
        end
    end
    correrrtzLd(bj1,:) = correrrtzLdiff;
    correrrtzRd(bj1,:) = correrrtzRdiff;
    correrrtzLRd(bj1,:) = correrrtzLRdiff;
    correrrtzRLd(bj1,:) = correrrtzRLdiff;
end

%INs
correrrtzINLd = zeros(size(upINspatcells,2),10);
correrrtzINRd = zeros(size(upINspatcells,2),10);
correrrtzINLRd = zeros(size(upINspatcells,2),10);
correrrtzINRLd = zeros(size(upINspatcells,2),10);
for bj1 = 1 : size(upINspatcells,2)
    inTZ = upINspatcells{4,bj1}(:,1 : stem(1)/2);
    onlyleftzIN = inTZ(purelefts,:);
    onlyrightzIN = inTZ(purerights,:);
    errleftzIN = inTZ(errorlefts,:);
    errightzIN = inTZ(errorights,:);
    correrrtzLdiffIN = zeros(1,size(inTZ,2));
    correrrtzRdiffIN = zeros(1,size(inTZ,2));
    correrrtzLRdiffIN = zeros(1,size(inTZ,2));
    correrrtzRLdiffIN = zeros(1,size(inTZ,2));
    for bj2 = 1 : size(inTZ,2)
        if length(errorlefts) >= 5
            correrrtzLdiffIN(bj2) = abs(log10(ranksum(onlyleftzIN(:,bj2),errleftzIN(:,bj2))));
        else
            correrrtzLdiffIN(bj2) = NaN;
        end
        if length(errorights) >= 5
            correrrtzRdiffIN(bj2) = abs(log10(ranksum(onlyrightzIN(:,bj2),errightzIN(:,bj2))));
        else
            correrrtzRdiffIN(bj2) = NaN;
        end
        if length(errorights) >= 5
            correrrtzLRdiffIN(bj2) = abs(log10(ranksum(onlyleftzIN(:,bj2),errightzIN(:,bj2))));
        else
            correrrtzLRdiffIN(bj2) = NaN;
        end
        if length(errorlefts) >= 5
            correrrtzRLdiffIN(bj2) = abs(log10(ranksum(onlyrightzIN(:,bj2),errleftzIN(:,bj2))));
        else
            correrrtzRLdiffIN(bj2) = NaN;
        end
    end
    correrrtzINLd(bj1,:) = correrrtzLdiffIN;
    correrrtzINRd(bj1,:) = correrrtzRdiffIN;
    correrrtzINLRd(bj1,:) = correrrtzLRdiffIN;
    correrrtzINRLd(bj1,:) = correrrtzRLdiffIN;
end

PFA_TC_differences = struct('Choice_context_selectivity',choices,'Choice_context_correlation',choicecorrmtx,'Choice_proportions',choiceprop,'Choice_ratios',covarcorr,...
    'Correct_LeftRight_stem_PC_TC_difference',eoLRd,'Correct_LeftRight_TZ_PC_TC_difference',tzLRd,...
    'Correct_LeftRight_stem_IN_TC_difference',inLRd,'Correct_LeftRight_TZ_IN_TC_difference',intzLRd,...
    'Error_LeftRight_stem_PC_TC_difference',eoerrLRd,'Error_LeftRight_stem_IN_TC_difference',eoerrINLRd,...
    'Error_LeftRight_TZ_PC_TC_difference',eoerrtzLRd,'Error_LeftRight_TZ_IN_TC_difference',eoerrINtzLRd,...
    'Correct_Error_stem_PC_TC_Ldiff',correrrLd,'Correct_Error_stem_PC_TC_Rdiff',correrrRd,...
    'Correct_Error_stem_PC_TC_LRdiff',correrrLRd,'Correct_Error_stem_PC_TC_RLdiff',correrrRLd,...
    'Correct_Error_stem_IN_TC_Ldiff',correrrINLd,'Correct_Error_stem_IN_TC_Rdiff',correrrINRd,...
    'Correct_Error_stem_IN_TC_LRdiff',correrrINLRd,'Correct_Error_stem_IN_TC_RLdiff',correrrINRLd,...
    'Correct_Error_TZ_PC_TC_Ldiff',correrrtzLd,'Correct_Error_TZ_PC_TC_Rdiff',correrrtzRd,...
    'Correct_Error_TZ_PC_TC_LRdiff',correrrtzLRd,'Correct_Error_TZ_PC_TC_RLdiff',correrrtzRLd,...
    'Correct_Error_TZ_IN_TC_Ldiff',correrrtzINLd,'Correct_Error_TZ_IN_TC_Rdiff',correrrtzINRd,...
    'Correct_Error_TZ_IN_TC_LRdiff',correrrtzINLRd,'Correct_Error_TZ_IN_TC_RLdiff',correrrtzINRLd);
filename_TCd = [sessionname{sessnum} '_TC_differences.mat'];
save(fullfile(resdir,sessionname{sessnum},filename_TCd),'PFA_TC_differences')

%%
%%Population correlation analysis (both place fielded and place fieldless
%%cells are used)
disp('Population correlation analysis...')

%trials (both sides) is >= 5.
%Templates
pthres = abs(log10(0.05));
bigFRmtx = cat(3,mtx2Full{:});
% bigFRmtx = permute(bigFRmtx,[2,1,3]);
leftFRmtx = bigFRmtx(purelefts,:,:);
rightFRmtx = bigFRmtx(purerights,:,:);
popX = NaN(3,size(bigFRmtx,1));

locXmtx = zeros(size(bigFRmtx,2),size(bigFRmtx,2),size(bigFRmtx,1));
locPmtx = zeros(size(bigFRmtx,2),size(bigFRmtx,2),size(bigFRmtx,1));
popCC = zeros(2,size(bigFRmtx,1));
for cq1 = 1 : size(bigFRmtx,1)
    locmtx = squeeze(bigFRmtx(cq1,:,:))';
    [locXmtx(:,:,cq1),locPmtx(:,:,cq1)] = corr(locmtx,'type','spearman');
    popCC(1,cq1) = nanmean(nanmean(locXmtx(:,:,cq1)));
    popCC(2,cq1) = abs(log10(nanmean(nanmean(locPmtx(:,:,cq1)))));
end
popXScores = [mean(popCC(1,:)),mean(popCC(1,purelefts)),mean(popCC(1,purerights));...
    mean(popCC(2,:)),mean(popCC(2,purelefts)),mean(popCC(2,purerights))];

PFA_population_correlation = struct('Pop_correlation_matrix',locXmtx,'Pop_corrrelation_signif_matrix',locPmtx,...
    'Mean_correlation_values',popCC,'Session_correlation_scores',popXScores);
filename_PVx = [sessionname{sessnum} '_population_correlation.mat'];
save(fullfile(resdir,sessionname{sessnum},filename_PVx),'PFA_population_correlation')

%%
%trial-by-trial firing rate (tuning curve and raw)
disp('Calculation of selectivity indices...')

% sFRv = selectedFRvects(stemPC_ID);

sFRt = fullFRvects2(cellswithfields > 0);
sFRt = sFRt(stemPC_ID);
sFRv = cell(1,length(sFRt));
for ik = 1 : length(sFRt)
    sFRv{ik} = sFRt{ik}(1:31,:);
end
trialvects = zeros(length(stemTCc),constrialnum);
trialFRSvect = zeros(length(stemTCc),constrialnum);
trialFRTZvect = zeros(length(stemTCc),constrialnum);
for i1 = 1 : length(stemTCc)
trialvects(i1,:) = sum(stemTCc{i1}(:,1:end-1),2); %Note: it's full maze!!!!
trialFRSvect(i1,:) = nanmean(sFRv{i1}(9:end,:));
trialFRTZvect(i1,:) = nanmean(sFRv{i1}(1:9,:));
end

%FR comparison in left/right trials: stem
lrFRcmp = zeros(1,length(stemTCc));
lfFR = zeros(2,length(stemTCc));
lfFRratio = zeros(1,length(stemTCc));
for lr1 = 1 : length(stemTCc)
    leftFR = trialFRSvect(lr1,purelefts);
    leftFR(isinf(leftFR)) = NaN;
    rightFR = trialFRSvect(lr1,purerights);
    rightFR(isinf(rightFR)) = NaN;
    lfFR(:,lr1) = [nanmean(leftFR);nanmean(rightFR)];
    lrFRcmp(lr1) = ranksum(leftFR(~isnan(leftFR)),rightFR(~isnan(rightFR)));
    lfFRratio(lr1) = (lfFR(1,lr1) - lfFR(2,lr1)) / (lfFR(1,lr1) + lfFR(2,lr1));
end

%FR comparison in left/right trials: TZ
lrFRTZcmp = zeros(1,length(stemTCc));
lfFRTZ = zeros(2,length(stemTCc));
lfFRTZratio = zeros(1,length(stemTCc));
for lr1 = 1 : length(stemTCc)
    leftzFR = trialFRTZvect(lr1,purelefts);
    leftzFR(isinf(leftzFR)) = NaN;
    rightzFR = trialFRTZvect(lr1,purerights);
    rightzFR(isinf(rightzFR)) = NaN;
    lfFRTZ(:,lr1) = [nanmean(leftzFR);nanmean(rightzFR)];
    lrFRTZcmp(lr1) = ranksum(leftzFR(~isnan(leftzFR)),rightzFR(~isnan(rightzFR)));
    lfFRTZratio(lr1) = (lfFRTZ(1,lr1) - lfFRTZ(2,lr1)) / (lfFRTZ(1,lr1) + lfFRTZ(2,lr1));
end

%trial by trial splitterness map defined as sum(FR(i)) - sum(FR(i-1)) / sum(FR(i)) +
%sum(FR(i-1))
zonepopside = cell(1,length(sFRv));
fullpopside = zeros(length([purelefts,purerights]),size(sFRv{1},1),length(sFRv));
popside = zeros(fix(length([purelefts,purerights])/2),size(sFRv{1},1),length(sFRv));
splitmag = zeros(fix(length([purelefts,purerights])/2),size(sFRv{1},1),length(sFRv));
for j1 = 1 : length(sFRv)
    %only correct trials counted
    sspFR = sFRv{j1}';
    sspFR([errorlefts,errorights],:) = [];
    evenodds = [2 : 2 : size(sspFR,1);1 : 2 : size(sspFR,1) - 1];
    for j2 = 1 : length(evenodds)
        tstore = sspFR(evenodds(:,j2),:);
        for j3 = 1 : size(sspFR,2)
            popside(j2,j3,j1) = (tstore(1,j3) - tstore(2,j3)) / (tstore(1,j3) + tstore(2,j3));
            splitmag(j2,j3,j1) = abs(popside(j2,j3,j1));
            if isnan(popside(j2,j3,j1))
                popside(j2,j3,j1) = -1.2;
                splitmag(j2,j3,j1) = -1.2;
            end
        end
    end
    %all trials counted
    ccpFR = sFRv{j1}';
    for j4 = 1 : size(sspFR,1) - 1
        for j5 = 1 : size(sspFR,2)
            fullpopside(j4,j5,j1) = (ccpFR(j4 + 1,j5) - ccpFR(j4,j5)) / (ccpFR(j4 + 1,j5) + ccpFR(j4,j5));
            if isnan(fullpopside(j4,j5,j1))
                fullpopside(j4,j5,j1) = -1.1;
            end
        end
    end
    tps = squeeze(popside(:,:,j1));
    tzps = tps(:,1 : 9);
    stemps = tps(:,10 : end);
    zps = zeros(size(popside,1),size(stemregions{j1},2));
    for j6 = 1 : size(stemregions{j1},2)
        actzone = [stemregions{j1}(1,j6) stemregions{j1}(2,j6)];
        actzone = actzone - 10;
        if actzone(1) <= 0
            actzone(1) = 1;
        end
        if actzone(2) > 22
            actzone(2) = 22;
        end
        if ~isnan(actzone(1))
            for j7 = 1 : size(popside,1)
                tps(tps == -1.2) = NaN;
                zps(j7,j6) = nanmean(stemps(j7,actzone(1) : actzone(2)));
            end
        else
            zps(j7,j6) = NaN;
        end
    end
    zonepopside{j1} = [nanmean(tzps,2),zps]; %SI place field zones of correct trials
end

%trail-by-trial place field zone splitterness
PFZsplitterness = cell(2,length(stemregions));
for u1 = 1 : length(stemregions)
    triPFsplt  = cell(1,size(stemregions{u1},2));
    zoneSplt = cell(1,size(stemregions{u1},2));
    for u2 = 1 : size(stemregions{u1},2)
        tFRmtx = sFRv{u1}';
        tFRmtx = tFRmtx(:,10 : end);
        if ~isempty(stemregions{u1})
            trialrange = (stemregions{u1}(4,u2) - stemregions{u1}(3,u2)) + 1;
            startrial = stemregions{u1}(3,u2);
            spatzones = [stemregions{u1}(1,u2) stemregions{u1}(2,u2)];
            spatzones = spatzones - 10;
            if spatzones(1) <= 0
                spatzones(1) = 1;
            end
            if spatzones(2) > 22
                spatzones(2) = 22;
            end
            pfzone = tFRmtx(stemregions{u1}(3,u2) : stemregions{u1}(4,u2),spatzones(1) : spatzones(2));
            zonesplitness = NaN(2,size(pfzone,1));
            for u3 = 1 : size(pfzone,1) - 1
                 zonesplitness(1,u3) = (nanmean(pfzone(u3 + 1,:)) - nanmean(pfzone(u3,:))) / ...
                     (nanmean(pfzone(u3 + 1,:)) + nanmean(pfzone(u3,:)));
                zonesplitness(2,u3) = abs(zonesplitness(1,u3));
            end
%             zonesplitness(:,isnan(zonesplitness(1,:))) = [];
            zoneSplt{u2} = zonesplitness;
        else
            zoneSplt{u2} = NaN;
        end
    end
    tzt = sFRv{u1}';
    tzzone = tzt(:,1 : 9);
    tzsplitness = NaN(2,size(tzzone,1));
    for u4 = 1 : size(tzzone,1) - 1
        tzsplitness(1,u4) = (nanmean(tzzone(u4 + 1,:)) - nanmean(tzzone(u4,:))) / ...
             (nanmean(tzzone(u4 + 1,:)) + nanmean(tzzone(u4,:)));
        tzsplitness(2,u4) = abs(tzsplitness(1,u4));
    end
    PFZsplitterness{1,u1} = zoneSplt;
    PFZsplitterness{2,u1} = tzsplitness;
end

%assign splitterness to spikes / trial in stem
sspos = unitarm(cellswithfields > 0);
sspos = sspos(stemPC_ID);
sstimes = unitPTP_arm(cellswithfields > 0);
sstimes = sstimes(stemPC_ID);
tzpos = unitranz(cellswithfields > 0);
tzpos = tzpos(stemPC_ID);
tztimes = unitPTP_tranz(cellswithfields > 0);
tztimes = tztimes(stemPC_ID);
fullspksides = cell(1,length(sstimes));
sidespkratio = zeros(1,length(sstimes));
fulltzspksides = cell(1,length(sstimes));
sidespkratzio = zeros(1,length(tztimes));
for uz1 = 1 : length(PFZsplitterness)
    tspk = sspos{uz1}';
    tspkt = sstimes{uz1}';
    tzspk = tzpos{uz1}';
    tzspkt = tztimes{uz1}';
    trialsidecell = cell(1,constrialnum);
    tzrialsidecell = cell(1,constrialnum);
    for uz2 = 1 : constrialnum
        trialspkp = tspk{uz2}';
        trialspkt = tspkt{uz2}';
        tzrialspkp = tzspk{uz2}';
        tzrialspkt = tzspkt{uz2}';
        trialspkside = zeros(1,length(trialspkp));
        for uz3 = 1 : size(stemregions{uz1},2)
            serange = [stemregions{uz1}(1,uz3) stemregions{uz1}(2,uz3)];
            terange = [stemregions{uz1}(3,uz3) stemregions{uz1}(4,uz3)];
            if uz2 >= terange(1) && uz2 < terange(2)
                trialspkside(trialspkp >= serange(1)*2 & trialspkp <= serange(2)*2) = PFZsplitterness{1,uz1}{uz3}(1,uz2 - terange(1) + 1);
            elseif uz2 == terange(2)
                trialspkside(trialspkp >= serange(1)*2 & trialspkp <= serange(2)*2) = PFZsplitterness{1,uz1}{uz3}(1,uz2 - terange(1));
            else
                trialspkside(trialspkp >= serange(1)*2 & trialspkp <= serange(2)*2) = NaN;
            end
        end
        trialsidecell{uz2} = trialspkside;
        %TZ
        tzrialsidecell{uz2} = ones(1,length(tzrialspkp)) * PFZsplitterness{2,uz1}(1,uz2);
    end
    fullspksides{uz1} = trialsidecell;
    sidespkratio(uz1) = length(find(cat(2,fullspksides{uz1}{:}) ~= 0)) / length(cat(2,fullspksides{uz1}{:}));
    fulltzspksides{uz1} = tzrialsidecell;
    sidespkratzio(uz1) = length(find(cat(2,fulltzspksides{uz1}{:}) ~= 0)) / length(cat(2,fulltzspksides{uz1}{:}));
end

%SI for interneurons (only correct trials)
intzstemFR = cell(1,length(infullFRvect));
for ip = 1 : length(infullFRvect)
    intzstemFR{ip} = infullFRvect{ip}(:,1 : 31);
end
popINside = zeros(fix(length([purelefts,purerights])/2),size(intzstemFR{1},2),length(inFRvects2));
splitINmag = zeros(fix(length([purelefts,purerights])/2),size(intzstemFR{1},2),length(inFRvects2));
for js1 = 1 : length(infullFRvect)
    INfr = intzstemFR{js1};
    INfr([errorlefts,errorights],:) = [];
    if mod(size(INfr,1),2) > 0
        INfr(end,:) = [];
    end
    eoin = [2 : 2 : size(INfr,1);1 : 2 : size(INfr,1)];
    for js2 = 1 : length(eoin)
        tinstore = INfr(eoin(:,js2),:);
        for js3 = 1 : size(intzstemFR{js1},2)
            popINside(js2,js3,js1) = (tinstore(1,js3) - tinstore(2,js3)) / (tinstore(1,js3) + tinstore(2,js3));
            splitINmag(js2,js3,js1) = abs(popINside(js2,js3,js1));
            if isnan(popINside(js2,js3,js1))
                popINside(js2,js3,js1) = -1.2;
                splitINmag(js2,js3,js1) = -1.2;
            end
        end
    end
end

PFA_SI = struct('LeftRight_stem_FR_comp',lrFRcmp,'LeftRight_stem_FR_ratio',lfFRratio,'LeftRight_TZ_FR_comp',lrFRTZcmp,'LeftRight_TZ_FR_ratio',lfFRTZratio,...
    'Loc_correct_trial_SI',popside,'Loc_correct_trial_SI_magnitude',splitmag,'Loc_all_trial_SI',fullpopside,'PFzone_correct_trial_SI',zonepopside);
PFA_PFZ_SI1 = struct('PFzone_all_trial_SI',PFZsplitterness);
PFA_PFZ_SI2 = struct('Spike_stem_SI',fullspksides,'Spike_TZ_SI',fulltzspksides,'Spike_stem_side_ratio',sidespkratio,'Spike_TZ_side_ratio',sidespkratzio,...
    'Loc_correct_trial_IN_SI',popINside,'Loc_correct_trial_IN_SI_magnitude',splitINmag);
filename_SI = [sessionname{sessnum} '_selectivity_indices.mat'];
filename_PFZSI1 = [sessionname{sessnum} '_PFZ_selectivity_indices1.mat'];
filename_PFZSI2 = [sessionname{sessnum} '_PFZ_selectivity_indices2.mat'];
save(fullfile(resdir,sessionname{sessnum},filename_SI),'PFA_SI')
save(fullfile(resdir,sessionname{sessnum},filename_PFZSI1),'PFA_PFZ_SI1')
save(fullfile(resdir,sessionname{sessnum},filename_PFZSI2),'PFA_PFZ_SI2')

%%
%splitterness in theta cycles in stem
disp('Theta-scale selectivity index...')

selspktimes = unitPTP_arm(cellswithfields > 0);
selspktimes = selspktimes(stemPC_ID);
thetacyclesplit = cell(1,constrialnum);
trialcyclesplit = zeros(1,constrialnum);
thetahalfcyclesplit = cell(1,constrialnum);
for jt1 = 1 : constrialnum
    %SI per theta cycle
    tavesplit = zeros(1,length(trialcycles{jt1}) - 1);
    for jt2 = 1 : length(trialcycles{jt1}) - 1
        percyclesplit = cell(1,length(fullspksides));
        for jt3 = 1 : length(fullspksides)
            percyclesplit{jt3} = fullspksides{jt3}{jt1}(selspktimes{jt3}{jt1} >= trialcycles{jt1}(jt2) & selspktimes{jt3}{jt1} < trialcycles{jt1}(jt2 + 1));
        end
        tavesplit(jt2) = nanmean(cat(2,percyclesplit{:}));
    end
    thetacyclesplit{jt1} = tavesplit;
    trialcyclesplit(jt1) = nanmean(thetacyclesplit{jt1});
    %SI per half theta cycle
    trialhalfcycles = zeros(1,length(trialcycles{jt1}) * 2 - 1);
    intercyc = diff(trialcycles{jt1}) / 2;
    trialhalfcycles(1 : 2 : end) = trialcycles{jt1};
    trialhalfcycles(2 : 2 : end - 1) = trialcycles{jt1}(1 : end - 1) + intercyc;
    tavehalfsplit = zeros(1,length(trialhalfcycles) - 1);
    for jt2 = 1 : length(trialhalfcycles) - 1
        perhalfcyclesplit = cell(1,length(fullspksides));
        for jt3 = 1 : length(fullspksides)
            perhalfcyclesplit{jt3} = fullspksides{jt3}{jt1}(selspktimes{jt3}{jt1} >= trialhalfcycles(jt2) & selspktimes{jt3}{jt1} < trialhalfcycles(jt2 + 1));
        end
        tavehalfsplit(jt2) = nanmean(cat(2,perhalfcyclesplit{:}));
    end
    thetahalfcyclesplit{jt1} = tavehalfsplit;
end

allthsplt = cat(2,thetacyclesplit{:});
allthsplt = allthsplt(~isnan(allthsplt));
sessionsplit = length(allthsplt(allthsplt ~= 0)) / length(allthsplt);

%splitterness in theta cycles in TZ
seltzspktimes = unitPTP_tranz(cellswithfields > 0);
seltzspktimes = seltzspktimes(stemPC_ID);
thetacyclesplitz = cell(1,constrialnum);
trialcyclesplitz = zeros(1,constrialnum);
thetahalfcyclesplitz = cell(1,constrialnum);
for jt1 = 1 : constrialnum
    %SI per theta cycle
    tavesplitz = zeros(1,length(trialcycles{jt1}) - 1);
    for jt2 = 1 : length(trialcycles{jt1}) - 1
        percyclesplitz = cell(1,length(fulltzspksides));
        for jt3 = 1 : length(fulltzspksides)
            percyclesplitz{jt3} = fulltzspksides{jt3}{jt1}(seltzspktimes{jt3}{jt1} >= trialcycles{jt1}(jt2) & seltzspktimes{jt3}{jt1} < trialcycles{jt1}(jt2 + 1));
        end
        tavesplitz(jt2) = nanmean(cat(2,percyclesplitz{:}));
    end
    thetacyclesplitz{jt1} = tavesplitz;
    trialcyclesplitz(jt1) = nanmean(thetacyclesplitz{jt1});
    %SI per half theta cycle
    trialhalfcyclesz = zeros(1,length(trialcycles{jt1}) * 2 - 1);
    intercycz = diff(trialcycles{jt1}) / 2;
    trialhalfcyclesz(1 : 2 : end) = trialcycles{jt1};
    trialhalfcyclesz(2 : 2 : end - 1) = trialcycles{jt1}(1 : end - 1) + intercycz;
    tavehalfsplitz = zeros(1,length(trialhalfcyclesz) - 1);
    for jt2 = 1 : length(trialhalfcyclesz) - 1
        perhalfcyclesplitz = cell(1,length(fulltzspksides));
        for jt3 = 1 : length(fulltzspksides)
            perhalfcyclesplitz{jt3} = fulltzspksides{jt3}{jt1}(seltzspktimes{jt3}{jt1} >= trialhalfcyclesz(jt2) & seltzspktimes{jt3}{jt1} < trialhalfcyclesz(jt2 + 1));
        end
        tavehalfsplitz(jt2) = nanmean(cat(2,perhalfcyclesplitz{:}));
    end
    thetahalfcyclesplitz{jt1} = tavehalfsplitz;
end

allthspltz = cat(2,thetacyclesplitz{:});
allthspltz = allthspltz(~isnan(allthspltz));
sessionsplitz = length(allthspltz(allthspltz ~= 0)) / length(allthspltz);

PFA_theta_SI = struct('Theta_half_cycle_SI_stem',thetahalfcyclesplit,'Full_theta_cycle_SI_stem',allthsplt,'Session_stem_SI',sessionsplit,...
    'Theta_half_cycle_SI_TZ',thetahalfcyclesplitz,'Full_theta_cycle_SI_TZ',allthspltz,'Session_TZ_SI',sessionsplitz);
filename_thSI = [sessionname{sessnum} '_selectivity_in_theta_cycles.mat'];
save(fullfile(resdir,sessionname{sessnum},filename_thSI),'PFA_theta_SI')

%%
%correlation of splitterness and choice trail-by-trial: consistency of the
%SI code: stem
disp('Consistency of selectivity...')

splithresh = 0 : 0.1 : 0.9;
trialsplitterness = zeros(length(sFRv),constrialnum);
trialsplittzerness = zeros(length(sFRv),constrialnum);
choiceprobs = cell(1,length(sFRv));
choiceprobtz = cell(1,length(sFRv));
maxchoiceprob = zeros(length(sFRv),8);
maxchoiceprobtz = zeros(length(sFRv),8);
for z1 = 1 : length(sFRv)
    trialsplt = zeros(1,constrialnum);
    trialspltz = zeros(1,constrialnum);
    for z2 = 1 : constrialnum - 1
        trialsplt(z2) = (trialFRSvect(z1,z2 + 1) - trialFRSvect(z1,z2)) / (trialFRSvect(z1,z2 + 1) + trialFRSvect(z1,z2));
        trialspltz(z2) = (trialFRTZvect(z1,z2 + 1) - trialFRTZvect(z1,z2)) / (trialFRTZvect(z1,z2 + 1) + trialFRTZvect(z1,z2));
    end
    trialsplitterness(z1,:) = trialsplt;
    trialsplittzerness(z1,:) = trialspltz;
    checkchoice2 = zeros(1,constrialnum);
    checkchoice2(purelefts) = 1;
    checkchoice2(purerights) = -2;
    checkchoice2(errorlefts) = 4;
    checkchoice2(errorights) = -5;
    probvect = zeros(length(splithresh),8);
    probvectz = zeros(length(splithresh),8);
    for z3 = 1 : length(splithresh)
        checkchoice1 = zeros(1,constrialnum);
        checkchoice1tz = zeros(1,constrialnum);
        negs = find(trialsplt < -splithresh(z3));
        poss = find(trialsplt > splithresh(z3));
        negtz = find(trialspltz < -splithresh(z3));
        postz = find(trialspltz > splithresh(z3));
        checkchoice1(negs) = 1;
        checkchoice1(poss) = -1;
        checkchoice1tz(negtz) = 1;
        checkchoice1tz(postz) = -1;
        sumchoice = checkchoice1 + checkchoice2;
        sumchoicetz = checkchoice1tz + checkchoice2;
        matNchoiceL = find(sumchoice == 2);
        matNchoiceR = find(sumchoice == -1);
        matNchoiceeL = find(sumchoice == 5);
        matNchoiceeR = find(sumchoice == -4);
        matPchoiceL = find(sumchoice == 0);
        matPchoiceR = find(sumchoice == -3);
        matPchoiceeL = find(sumchoice == 3);
        matPchoiceeR = find(sumchoice == -6);
        matNchoiceLtz = find(sumchoicetz == 2);
        matNchoiceRtz = find(sumchoicetz == -1);
        matNchoiceeLtz = find(sumchoicetz == 5);
        matNchoiceeRtz = find(sumchoicetz == -4);
        matPchoiceLtz = find(sumchoicetz == 0);
        matPchoiceRtz = find(sumchoicetz == -3);
        matPchoiceeLtz = find(sumchoicetz == 3);
        matPchoiceeRtz = find(sumchoicetz == -6);
        probNL = length(matNchoiceL) / length(purelefts);
        probNR = length(matNchoiceR) / length(purerights);
        probNeL = length(matNchoiceeL) / length(errorlefts);
        probNeR = length(matNchoiceeR) / length(errorights);
        probPL = length(matPchoiceL) / length(purelefts);
        probPR = length(matPchoiceR) / length(purerights);
        probPeL = length(matPchoiceeL) / length(errorlefts);
        probPeR = length(matPchoiceeR) / length(errorights);
        probtzNL = length(matNchoiceLtz) / length(purelefts);
        probtzNR = length(matNchoiceRtz) / length(purerights);
        probtzNeL = length(matNchoiceeLtz) / length(errorlefts);
        probtzNeR = length(matNchoiceeRtz) / length(errorights);
        probtzPL = length(matPchoiceLtz) / length(purelefts);
        probtzPR = length(matPchoiceRtz) / length(purerights);
        probtzPeL = length(matPchoiceeLtz) / length(errorlefts);
        probtzPeR = length(matPchoiceeRtz) / length(errorights);
        probvect(z3,:) = [probNL probNR probNeL probNeR probPL probPR probPeL probPeR];
        probvectz(z3,:) = [probtzNL probtzNR probtzNeL probtzNeR probtzPL probtzPR probtzPeL probtzPeR];
    end
    choiceprobs{z1} = probvect;
    choiceprobtz{z1} = probvectz;
    for z4 = 1 : size(probvect,2)
        tmcp = choiceprobs{z1}(choiceprobs{z1}(:,z4) == max(choiceprobs{z1}(:,z4)),z4);
        if isempty(tmcp) == 0
            maxchoiceprob(z1,z4) = tmcp(1);
        else
            maxchoiceprob(z1,z4) = NaN;
        end
    end
    for z5 = 1 : size(probvectz,2)
        tzmcp = choiceprobtz{z1}(choiceprobtz{z1}(:,z5) == max(choiceprobtz{z1}(:,z5)),z5);
        if isempty(tzmcp) == 0
            maxchoiceprobtz(z1,z5) = tzmcp(1);
        else
            maxchoiceprobtz(z1,z5) = NaN;
        end
    end
    if length(errorlefts) < 5
        maxchoiceprob(:,[3 7]) = NaN;
        maxchoiceprobtz(:,[3 7]) = NaN;
    elseif length(errorights) < 5
        maxchoiceprob(:,[4 8]) = NaN;
        maxchoiceprobtz(:,[4 8]) = NaN;
    end
end

%correlation of splitterness and choice trail-by-trial: consistency of the
%SI code: stem

%correlation of splitterness and choice in stem PF regions
choiceAllPFprobs = cell(1,length(sFRv));
maxchoiceAllPFprob = cell(1,length(sFRv));
choiceprobPerTrial = cell(1,length(sFRv));
for zf1 = 1 : length(sFRv)
    choicePFprobs = cell(1,size(stemregions{zf1},2));
    maxchoicePFprob = zeros(size(stemregions{zf1},2),8);
    allspfc = cell(1,size(stemregions{zf1},2));
    for zf2 = 1 : size(stemregions{zf1},2)
        PFrange = (stemregions{zf1}(4,zf2) - stemregions{zf1}(3,zf2)) + 1;
        pureLpf = purelefts(purelefts >= stemregions{zf1}(3,zf2) & purelefts <= stemregions{zf1}(4,zf2)) - (stemregions{zf1}(3,zf2) - 1);
        pureRpf = purerights(purerights >= stemregions{zf1}(3,zf2) & purerights <= stemregions{zf1}(4,zf2)) - (stemregions{zf1}(3,zf2) - 1);
        errorLpf = errorlefts(errorlefts >= stemregions{zf1}(3,zf2) & errorlefts <= stemregions{zf1}(4,zf2)) - (stemregions{zf1}(3,zf2) - 1);
        errorRpf = errorights(errorights >= stemregions{zf1}(3,zf2) & errorights <= stemregions{zf1}(4,zf2)) - (stemregions{zf1}(3,zf2) - 1);
        checkPFchoice2 = zeros(1,size(PFZsplitterness{1,zf1}{zf2},2));
        checkPFchoice2(pureLpf) = 1;
        checkPFchoice2(pureRpf) = -2;
        checkPFchoice2(errorLpf) = 4;
        checkPFchoice2(errorRpf) = -5;
        probPFvect = zeros(length(splithresh),8);
        sPFchoice = zeros(length(splithresh),size(PFZsplitterness{1,zf1}{zf2},2));
        for zf3 = 1 : length(splithresh)
            checkPFchoice1 = zeros(1,size(PFZsplitterness{1,zf1}{zf2},2));
            pfnegs = find(PFZsplitterness{1,zf1}{zf2}(1,:) < -splithresh(zf3));
            pfposs = find(PFZsplitterness{1,zf1}{zf2}(1,:) > splithresh(zf3));
            checkPFchoice1(pfnegs) = 1;
            checkPFchoice1(pfposs) = -1;
            sumPFchoice = checkPFchoice1 + checkPFchoice2;
            sPFchoice(zf3,:) = sumPFchoice;
            matNPFchoiceL = find(sumPFchoice == 2);
            matNPFchoiceR = find(sumPFchoice == -1);
            matNPFchoiceeL = find(sumPFchoice == 5);
            matNPFchoiceeR = find(sumPFchoice == -4);
            matPPFchoiceL = find(sumPFchoice == 0);
            matPPFchoiceR = find(sumPFchoice == -3);
            matPPFchoiceeL = find(sumPFchoice == 3);
            matPPFchoiceeR = find(sumPFchoice == -6);
            probNLpf = length(matNPFchoiceL) / length(pureLpf);
            probNRpf = length(matNPFchoiceR) / length(pureRpf);
            probNeLpf = length(matNPFchoiceeL) / length(errorLpf);
            probNeRpf = length(matNPFchoiceeR) / length(errorRpf);
            probPLpf = length(matPPFchoiceL) / length(pureLpf);
            probPRpf = length(matPPFchoiceR) / length(pureRpf);
            probPeLpf = length(matPPFchoiceeL) / length(errorLpf);
            probPeRpf = length(matPPFchoiceeR) / length(errorRpf);
            probPFvect(zf3,:) = [probNLpf probNRpf probNeLpf probNeRpf probPLpf probPRpf probPeLpf probPeRpf];
        end
        choicePFprobs{zf2} = probPFvect;
        allspfc{zf2} = sPFchoice(1,:);
        for zf4 = 1 : size(probPFvect,2)
            tmcpf = choicePFprobs{zf2}(choicePFprobs{zf2}(:,zf4) == max(choicePFprobs{zf2}(:,zf4)),zf4);
            if isempty(tmcpf) == 0
                maxchoicePFprob(zf2,zf4) = tmcpf(1);
            else
                maxchoicePFprob(zf2,zf4) = NaN;
            end
        end
        if length(errorlefts) < 5
            maxchoicePFprob(:,[3 7]) = NaN;
        elseif length(errorights) < 5
            maxchoicePFprob(:,[4 8]) = NaN;
        end
    end
    choiceAllPFprobs{zf1} = choicePFprobs;
    maxchoiceAllPFprob{zf1} = maxchoicePFprob;
    choiceprobPerTrial{zf1} = allspfc;
end

PFA_choices = struct('Choice_stem_SI_consistency',maxchoiceprob,'Choice_TZ_SI_consistency',maxchoiceprobtz,...
    'PFzone_choice_probabilities',choiceAllPFprobs,'PFzone_choice_SI_consistency',maxchoiceAllPFprob,'Trialwise_choice_SI_consistency',choiceprobPerTrial);
filename_choices = [sessionname{sessnum} '_choice_selectivity_correlation.mat'];
save(fullfile(resdir,sessionname{sessnum},filename_choices),'PFA_choices')

%%
%FR comparison in PF zones
disp('Firing rate difference in place fields between left and right trials...')

FR_PF_comp = cell(2,length(stemregions));
for fq1 = 1 : length(stemregions)
    actmtx = sFRv{fq1}';
    FRintracmp = zeros(6,size(stemregions{fq1},2));
    frcL = cell(1,size(stemregions{fq1},2));
    frcR = cell(1,size(stemregions{fq1},2));
    freL = cell(1,size(stemregions{fq1},2));
    freR = cell(1,size(stemregions{fq1},2));
    for fq2 = 1 : size(stemregions{fq1},2)
        trLrange = purelefts(purelefts >= stemregions{fq1}(3,fq2) & purelefts <= stemregions{fq1}(4,fq2));
        trRrange = purerights(purerights >= stemregions{fq1}(3,fq2) & purerights <= stemregions{fq1}(4,fq2));
        treLrange = errorlefts(errorlefts >= stemregions{fq1}(3,fq2) & errorlefts <= stemregions{fq1}(4,fq2));
        treRrange = errorights(errorights >= stemregions{fq1}(3,fq2) & errorights <= stemregions{fq1}(4,fq2));
        sprange = stemregions{fq1}(1,fq2) - 10 : stemregions{fq1}(2,fq2) - 10;
        sprange(sprange <= 0) = 1;
        sprange(sprange > 22) = 22;
        frcL{fq2} = nanmean(actmtx(trLrange,sprange),2);
        nncL = length(frcL{fq2}(~isnan(frcL{fq2})));
        frcR{fq2} = nanmean(actmtx(trRrange,sprange),2);
        nncR = length(frcR{fq2}(~isnan(frcR{fq2})));
        freL{fq2} = nanmean(actmtx(treLrange,sprange),2);
        nneL = length(freL{fq2}(~isnan(freL{fq2})));
        freR{fq2} = nanmean(actmtx(treRrange,sprange),2);
        nneR = length(freR{fq2}(~isnan(freR{fq2})));
        if (length(treLrange) >= 5 && length(treRrange) >= 5) && (nncL >= 5 && nncR >= 5 && nneL >= 5 && nneR >= 5)
            FRintracmp(:,fq2) = abs(log10([ranksum(frcL{fq2},frcR{fq2}),ranksum(freL{fq2},freR{fq2}),ranksum(frcL{fq2},freL{fq2}),ranksum(frcL{fq2},freR{fq2}),...
                ranksum(frcR{fq2},freR{fq2}),ranksum(frcR{fq2},freL{fq2})]));
        elseif length(treLrange) >= 5 && length(treRrange) < 5 && (nncL >= 5 && nncR >= 5 && nneL >= 5)
            FRintracmp(:,fq2) = abs(log10([ranksum(frcL{fq2},frcR{fq2}),NaN,ranksum(frcL{fq2},freL{fq2}),NaN,...
                NaN,ranksum(frcR{fq2},freL{fq2})]));
        elseif length(treLrange) < 5 && length(treRrange) >= 5 && (nncL >= 5 && nncR >= 5 && nneR >= 5)
            FRintracmp(:,fq2) = abs(log10([ranksum(frcL{fq2},frcR{fq2}),NaN,NaN,ranksum(frcL{fq2},freR{fq2}),...
                ranksum(frcR{fq2},freR{fq2}),NaN]));
        elseif length(treLrange) < 5 && length(treRrange) < 5 && (nncL >= 5 && nncR >= 5)
            FRintracmp(:,fq2) = abs(log10([ranksum(frcL{fq2},frcR{fq2}),NaN,NaN,NaN,...
                NaN,NaN]));
        end
    end

    if size(stemregions{fq1},2) > 1
        numcomp = factorial(size(stemregions{fq1},2)) / (2 * factorial(size(stemregions{fq1},2) - 2));
        FRintercmp = zeros(4,numcomp);
        counter = 0;
        for fq3 = 1 : length(frcL) - 1
            for fq4 = fq3 + 1 : length(frcL)
                counter = counter + 1;
                if (length(treLrange) >= 5 && length(treRrange) >= 5) && (~isempty(freL{fq3}) && ~isempty(freL{fq4})) && (~isempty(freR{fq3}) && ~isempty(freR{fq4})) &&...
                        size(find(frcL{fq3} == ~isnan(frcL{fq3})),1) >= 5 && size(find(frcR{fq3} == ~isnan(frcR{fq3})),1) >= 5 &&...
                        size(find(frcL{fq4} == ~isnan(frcL{fq4})),1) >= 5 && size(find(frcR{fq4} == ~isnan(frcR{fq4})),1) >= 5 &&...
                        size(find(freL{fq3} == ~isnan(freL{fq3})),1) >= 5 && size(find(freR{fq3} == ~isnan(freR{fq3})),1) >= 5 &&...
                        size(find(freL{fq4} == ~isnan(freL{fq4})),1) >= 5 && size(find(freR{fq4} == ~isnan(freR{fq4})),1) >= 5
                    FRintercmp(:,counter) = abs(log10([ranksum(frcL{fq3},frcL{fq4}),ranksum(frcR{fq3},frcR{fq4}),ranksum(freL{fq3},freL{fq4}),ranksum(freR{fq3},freR{fq4})]));
                elseif (length(treLrange) >= 5 && length(treRrange) < 5) && ~isempty(freL{fq3}) && ~isempty(freL{fq4}) &&...
                        size(find(frcL{fq3} == ~isnan(frcL{fq3})),1) >= 5 && size(find(frcR{fq3} == ~isnan(frcR{fq3})),1) >= 5 &&...
                        size(find(frcL{fq4} == ~isnan(frcL{fq4})),1) >= 5 && size(find(frcR{fq4} == ~isnan(frcR{fq4})),1) >= 5 &&...
                        size(find(freL{fq3} == ~isnan(freL{fq3})),1) >= 5 && size(find(freL{fq4} == ~isnan(freL{fq4})),1) >= 5
                    FRintercmp(:,counter) = abs(log10([ranksum(frcL{fq3},frcL{fq4}),ranksum(frcR{fq3},frcR{fq4}),ranksum(freL{fq3},freL{fq4}),NaN]));
                elseif (length(treLrange) < 5 && length(treRrange) >= 5) && ~isempty(freR{fq3}) && ~isempty(freR{fq4}) &&...
                        size(find(frcL{fq3} == ~isnan(frcL{fq3})),1) >= 5 && size(find(frcR{fq3} == ~isnan(frcR{fq3})),1) >= 5 &&...
                        size(find(frcL{fq4} == ~isnan(frcL{fq4})),1) >= 5 && size(find(frcR{fq4} == ~isnan(frcR{fq4})),1) >= 5 &&...
                        size(find(freR{fq3} == ~isnan(freR{fq3})),1) >= 5 && size(find(freR{fq4} == ~isnan(freR{fq4})),1) >= 5
                    FRintercmp(:,counter) = abs(log10([ranksum(frcL{fq3},frcL{fq4}),ranksum(frcR{fq3},frcR{fq4}),NaN,ranksum(freR{fq3},freR{fq4})]));
                elseif length(treLrange) < 5 && length(treRrange) < 5 && size(find(frcL{fq3} == ~isnan(frcL{fq3})),1) >= 5 &&...
                        size(find(frcR{fq3} == ~isnan(frcR{fq3})),1) >= 5 && size(find(frcL{fq4} == ~isnan(frcL{fq4})),1) >= 5 &&...
                        size(find(frcR{fq4} == ~isnan(frcR{fq4})),1) >= 5
                    FRintercmp(:,counter) = abs(log10([ranksum(frcL{fq3},frcL{fq4}),ranksum(frcR{fq3},frcR{fq4}),NaN,NaN]));
                end
            end
        end
    else
        FRintercmp = NaN;
    end
    FR_PF_comp{1,fq1} = FRintracmp;
    FR_PF_comp{2,fq1} = FRintercmp;
end

%Splitterness of interneurons calculated in 2cm spatial bins trial by trial using the
%FR-based tuning curves
% selFRvects = spatFRvects2(clusteredPCS);

% popinside = zeros(length(inFRvects2),constrialnum - 1);
% splitinmag = zeros(length(inFRvects2),constrialnum - 1);
% for j1 = 1 : length(inFRvects2)
%     tinstore = nanmean(inFRvects2{j1}(:,1:end-1),2);
%     for j2 = 2 : constrialnum
%         popinside(j1,j2) = (tinstore(j2 - 1) - tinstore(j2)) / (tinstore(j2 - 1) + tinstore(j2));
%         splitinmag(j1,j2) = abs(popinside(j1,j2));
%     end
% end

% selFRvects = inFRvects2;
% splmtx = zeros(constrialnum,length(spatbins2)-1,length(selFRvects));
% for l1 = 1 : length(selFRvects)
%     temtx = selFRvects{l1};
%     for l2 = 1 : constrialnum - 1
%         tevect_a = temtx(l2 + 1,1 : end - 1);
%         tevect_a(isnan(tevect_a)) = 0;
%         tevect_b = temtx(l2,1 : end - 1);
%         tevect_b(isnan(tevect_b)) = 0;
%         splmtx(l2,:,l1) = (tevect_a - tevect_b) ./ (tevect_a + tevect_b);
%     end
% end

PFA_FiringRate_comparison = struct('PF_zone_FR_comp',FR_PF_comp);
filename_FRc = [sessionname{sessnum} '_FiringRate_comparison.mat'];
save(fullfile(resdir,sessionname{sessnum},filename_FRc),'PFA_FiringRate_comparison')

%%
%choice decoding from theta cycle PVs
disp('Decoding choice from theta cycle population vectors...')

% poscycleinbins
% thetacyclesplit
% thetahalfcyclesplit
posinTC = PinCb{1};
spattempSI = NaN(constrialnum,length(spatbins2) - 1);
actbins = find(spatfullbins >= 20 & spatfullbins < 62);
for tr = 1 : constrialnum
    selPinTC = posinTC{tr}(posinTC{tr} > 10 & posinTC{tr} <= 31);
%     thhalfcsplt = reshape(thetahalfcyclesplit{tr},2,size(thetahalfcyclesplit{tr},2) / 2);
    for ts = 2 : length(actbins)
        spattempSI(tr,ts - 1) = nanmean(thetacyclesplit{tr}(selPinTC >= actbins(ts - 1) & selPinTC < actbins(ts)));
        ssSI = spattempSI;
        if isnan(ssSI(tr,ts))
            ssSI(tr,ts) = -10;
        end
    end
end
avessSI = abs(nanmean(spattempSI)); %SI from theta cycle averages along the spatial axis

%distribution of SI along the spatial axis in the stem...
SIbins = -1 : 0.1 : 1;
ntrts = zeros(length(SIbins) - 1,length(actbins));
ntrtshalves = zeros(length(SIbins) - 1,length(actbins),2);
for st = 1 : length(actbins)
    SIspattrial = cell(1,constrialnum);
    SI2spattrial = cell(1,constrialnum);
    for rt = 1 : constrialnum
        tSI2 = reshape(thetahalfcyclesplit{rt},2,size(thetahalfcyclesplit{rt},2) / 2);
        SIspattrial{rt} = thetacyclesplit{rt}(posinTC{rt} == actbins(st));
        SI2spattrial{rt} = tSI2(:,posinTC{rt} == actbins(st));
    end
    tsptrial = cat(2,SIspattrial{:});
    t2sptrial = cat(2,SI2spattrial{:});
    if ~isempty(tsptrial)
        ntrts(:,st) = histcounts(tsptrial,SIbins,'normalization','probability')'; %full theta cycle
    end
    if ~isempty(t2sptrial)
        ntrtshalves(:,st,1) = histcounts(t2sptrial(1,:),SIbins,'normalization','probability')'; %1st half of theta cycles
        ntrtshalves(:,st,2) = histcounts(t2sptrial(2,:),SIbins,'normalization','probability')'; %2nd half of theta cycles
    end
end

%...and in the TZ
ntzrts = zeros(length(SIbins) - 1,10);
ntzrtshalves = zeros(length(SIbins) - 1,10,2);
for st = 1 : 8
%     TZPinTC = posinTC{tr}(posinTC{tr} >= 1 & posinTC{tr} < 9);
    SIspattrialz = cell(1,constrialnum);
    SI2spattrialz = cell(1,constrialnum);
    for rt = 1 : constrialnum
        tzSI2 = reshape(thetahalfcyclesplitz{rt},2,size(thetahalfcyclesplitz{rt},2) / 2);
        SIspattrialz{rt} = thetacyclesplitz{rt}(posinTC{rt} == st);
        SI2spattrialz{rt} = tzSI2(:,posinTC{rt} == st);
    end
    tsptrialz = cat(2,SIspattrialz{:});
    t2sptrialz = cat(2,SI2spattrialz{:});
    if ~isempty(tsptrialz)
        ntzrts(:,st) = histcounts(tsptrialz,SIbins,'normalization','probability')'; %full theta cycle
    end
    if ~isempty(t2sptrialz)
        ntzrtshalves(:,st,1) = histcounts(t2sptrialz(1,:),SIbins,'normalization','probability')'; %1st half of theta cycles
        ntzrtshalves(:,st,2) = histcounts(t2sptrialz(2,:),SIbins,'normalization','probability')'; %2nd half of theta cycles
    end
end

%choice decoding based on correlation with templates in stem
%1. generation of theta cycle templates bin-by-bin
actbinsTZS = find(spatfullbins > 0 & spatfullbins < 62);
thetavect = Thvect{1};
spattempPV = NaN(size(thetavect{1},1),length(actbinsTZS),constrialnum);
for tr = 1 : constrialnum
    selPinTC = posinTC{tr}(posinTC{tr} > 0 & posinTC{tr} <= 31);
    for ts = 1 : length(actbinsTZS)
        spattempPV(:,ts,tr) = nanmean(thetavect{tr}(:,selPinTC == actbinsTZS(ts)),2);
    end
end
cLtplt = nanmean(spattempPV(:,:,purelefts),3);
cRtplt = nanmean(spattempPV(:,:,purerights),3);
eLtplt = nanmean(spattempPV(:,:,errorlefts),3);
eRtplt = nanmean(spattempPV(:,:,errorights),3);

allthetaPV = cat(2,thetavect{:});
%bin-by-bin and trial-by-trial correlation of theta cycle PVs and context-dep. PV
%templates
thetacyclecorr = cell(1,constrialnum);
tcycsides = cell(1,constrialnum);
for tr = 1 : constrialnum
    tspatcyc = cell(1,length(actbinsTZS));
    for ts = 1 : length(actbinsTZS)
%         tcyclecorr(tr,ts,:) = [corr(spattempPV(:,ts,tr),cLtplt(:,ts),'type','spearman'),...
%             corr(spattempPV(:,ts,tr),cRtplt(:,ts),'type','spearman')];
        if ~isempty(poscycleinbins{tr}(poscycleinbins{tr} == ts))
            tcycc = zeros(2,length(poscycleinbins{tr}(poscycleinbins{tr} == ts)));
            for tss = 1 : size(poscycleinbins{tr}(poscycleinbins{tr} == ts),2)
                ttPV = thetavect{tr}(:,poscycleinbins{tr} == ts);
                tcycc(:,tss) = [corr(ttPV(:,tss),cLtplt(:,ts));corr(ttPV(:,tss),cRtplt(:,ts))];
            end
            tspatcyc{ts} = tcycc;
        end
    end
    tspatvect = cat(2,tspatcyc{:});
    tcycsides{tr} = tspatvect(1,:) > tspatvect(2,:);
    thetacyclecorr{tr} = tspatvect;
end

tccv = cat(2,thetacyclecorr{:});
tcsv = cat(2,tcycsides{:});
tccLR = [tccv(1,tcsv == 1),tccv(2,tcsv == 0)];
sideidx = [find(tcsv == 1),find(tcsv == 0)];
[as,bs] = sort(sideidx);
tccS = tccLR(bs);

corrbins = -0.5 : 0.05 : 1;
LTCxd = histcounts(tccv(1,tcsv == 1),corrbins,'normalization','probability');
RTCxd = histcounts(tccv(2,tcsv == 0),corrbins,'normalization','probability');

% figure,
% plot(find(tcsv == 1),tccS(tcsv == 1),'b.','markersize',12),
% hold on
% plot(find(tcsv == 0),tccS(tcsv == 0),'r.','markersize',12),
% plot(tccS,'k')

PFA_choicecode = struct('Theta_cycle_stem_SI_distr',spattempSI,'Theta_half_cycle_stem_SI',ntrtshalves,...
    'Theta_half_cycle_TZ_SI',ntzrtshalves,'Spatial_template_source_mtx',spattempPV,'Theta_cycle_spatial_tplt_corr',thetacyclecorr,...
    'Theta_cycle_sidedness',tcycsides,'Theta_cycle_correlation_with_sidedness',tccS,...
    'Left_correlation_values',LTCxd,'Right_correlation_values',RTCxd);
filename_chc = [sessionname{sessnum} '_choice_decoding.mat'];
save(fullfile(resdir,sessionname{sessnum},filename_chc),'PFA_choicecode')

%%
%sparsity and spatial info of PCs
disp('Calculation of spatial metrics...')

stemPCFR = fullFRvects2(cellswithfields > 0);
stemPCFR = stemPCFR(stemPC_ID);
sparsityPC = zeros(length(stemPCFR),constrialnum);
spkinfoPC = zeros(length(stemPCFR),constrialnum);
sparsall = zeros(1,length(stemPCFR));
spkinfall = zeros(1,length(stemPCFR));
for to1 = 1 : length(stemPCFR)
    for to2 = 1 : constrialnum
        trialoccup = fulldistr(to2,:) / sum(fulldistr(to2,:)); %p(x)
        tspv = stemPCFR{to1}';
        trialspatfr = (tspv(to2,:) ./ fulldistr(to2,:)) * fps; %spatial FR
        trialrate = nanmean(trialspatfr); % mean FR of trial
        sparsityPC(to1,to2) = (nansum(trialoccup .* trialspatfr)).^2 / nansum(trialoccup .* trialspatfr.^2);
        spkinfoPC(to1,to2) = nansum((trialspatfr / trialrate) .* log2(trialspatfr / trialrate) .* trialoccup);
    end
    sparsall(to1) = nanmean(sparsityPC(to1,:));
    spkinfall(to1) = nanmean(spkinfoPC(to1,:));
end
sparsityPCL = sparsityPC(:,purelefts);
spkinfoPCL = spkinfoPC(:,purelefts);
sparsityPCR = sparsityPC(:,purerights);
spkinfoPCR = spkinfoPC(:,purerights);

%LR comparison (PC)
sparsLR = zeros(1,length(stemPCFR));
spkinfLR = zeros(1,length(stemPCFR));
for at1 = 1 : length(stemPCFR)
    sparscL = sparsityPCL(at1,:);
    sparscR = sparsityPCR(at1,:);
    sparscL(isnan(sparscL)) = [];
    sparscR(isnan(sparscR)) = [];
    spkinfL = spkinfoPCL(at1,:);
    spkinfR = spkinfoPCR(at1,:);
    spkinfL(isnan(spkinfL)) = [];
    spkinfR(isnan(spkinfR)) = [];
    if ~isempty(sparscL) && ~isempty(sparscR)
        sparsLR(at1) = abs(log10(ranksum(sparscL,sparscR)));
    end
    if ~isempty(spkinfL) && ~isempty(spkinfR)
        spkinfLR(at1) = abs(log10(ranksum(spkinfL,spkinfR)));
    end
end

%sparsity and spatial info of INs
sparsIN = zeros(length(inFRvects2),constrialnum);
spkinfIN = zeros(length(inFRvects2),constrialnum);
sparsallIN = zeros(1,length(inFRvects2));
spkinfallIN = zeros(1,length(inFRvects2));
for to1 = 1 : length(inFRvects2)
    for to2 = 1 : constrialnum
        trialoccup = alldistr(to2,:) / sum(alldistr(to2,:)); %p(x)
        tspvin = inFRvects2{to1};
        trialspatfrin = (tspvin(to2,:) ./ alldistr(to2,:)) * fps; %spatial FR
        trialINrate = nanmean(trialspatfrin); % mean FR of trial
        sparsIN(to1,to2) = (nansum(trialoccup .* trialspatfrin)).^2 / nansum(trialoccup .* trialspatfrin.^2);
        spkinfIN(to1,to2) = nansum((trialspatfrin / trialINrate) .* log2(trialspatfrin / trialINrate) .* trialoccup);
    end
    sparsallIN(to1) = nanmean(sparsIN(to1,:));
    spkinfallIN(to1) = nanmean(spkinfIN(to1,:));
end
sparsityINL = sparsIN(:,purelefts);
spkinfoINL = spkinfIN(:,purelefts);
sparsityINR = sparsIN(:,purerights);
spkinfoINR = spkinfIN(:,purerights);

%LR comparison (IN)
sparsinLR = zeros(1,length(inFRvects2));
spkinfinLR = zeros(1,length(inFRvects2));
for at1 = 1 : length(inFRvects2)
    sparscINL = sparsityINL(at1,:);
    sparscINR = sparsityINR(at1,:);
    sparscINL(isnan(sparscINL)) = [];
    sparscINR(isnan(sparscINR)) = [];
    spkinfINL = spkinfoINL(at1,:);
    spkinfINR = spkinfoINR(at1,:);
    spkinfINL(isnan(spkinfINL)) = [];
    spkinfINR(isnan(spkinfINR)) = [];
    if ~isempty(sparscINL) && ~isempty(sparscINR)
        sparsinLR(at1) = abs(log10(ranksum(sparscINL,sparscINR)));
    end
    if ~isempty(spkinfINL) && ~isempty(spkinfINR)
        spkinfinLR(at1) = abs(log10(ranksum(spkinfINL,spkinfINR)));
    end
end

%spatial coherence: pyramidal cells
%selection of FR tuning curves based on spike number in stem
%spatial coherence
spatcohFR = zeros(length(stemPCFR),constrialnum);
for sc1 = 1 : length(stemPCFR)
    scR = zeros(1,constrialnum);
    for sc2 = 1 : constrialnum
        zpspv = [0;stemPCFR{sc1}(:,sc2);0];
        zpspv(isnan(zpspv)) = 0;
        movaveTC = zeros(1,size(stemPCFR{sc1},1));
        for sc3 = 2 : size(stemPCFR{sc1},1) + 1
            movaveTC(sc3 - 1) = mean(zpspv(sc3 - 1 : sc3 + 1));
        end
        rawFRv = stemPCFR{sc1}(:,sc2);
        rawFRv(isnan(rawFRv)) = 0;
        RscR = corr(movaveTC',rawFRv,'Type','spearman');
        scR(sc2) = 0.5 * log((1 + RscR) / (1 - RscR));
    end
    spatcohFR(sc1,:) = scR;
end
spatcohFRL = spatcohFR(:,purelefts);
spatcohFRR = spatcohFR(:,purerights);

%LR comparison (PC)
spacohdLR = zeros(1,length(stemPCFR));
for at1 = 1 : length(stemPCFR)
    cohL = spatcohFRL(at1,:);
    cohR = spatcohFRR(at1,:);
    cohL(isnan(cohL)) = [];
    cohR(isnan(cohR)) = [];
    if ~isempty(cohL) && ~isempty(cohR)
        spacohdLR(at1) = abs(log10(ranksum(cohL,cohR)));
    end
end

%spatial coherence: interneurons
%selection of FR tuning curves based on spike number in stem
INspatFRvects2 = infullFRvect;
%spatial coherence
spatcohINFR = zeros(length(INspatFRvects2),constrialnum);
for sc1 = 1 : length(INspatFRvects2)
    scinR = zeros(1,constrialnum);
    INspfrv = INspatFRvects2{sc1}';
    for sc2 = 1 : constrialnum
        inspv = [0;INspfrv(:,sc2);0];
        inspv(isnan(inspv)) = 0;
        movaveIN = zeros(1,size(INspfrv,1));
        for sc3 = 2 : size(INspfrv,1) + 1
            movaveIN(sc3 - 1) = mean(inspv(sc3 - 1 : sc3 + 1));
        end
        rawinFRv = INspfrv(:,sc2);
        rawinFRv(isnan(rawinFRv)) = 0;
        RscinR = corr(movaveIN',rawinFRv,'Type','spearman');
        scinR(sc2) = 0.5 * log((1 + RscinR) / (1 - RscinR));
    end
    spatcohINFR(sc1,:) = scinR;
end
spatcohINFRL = spatcohINFR(:,purelefts);
spatcohINFRR = spatcohINFR(:,purerights);

%LR comparison (IN)
spacohdINLR = zeros(1,length(INspatFRvects2));
for at1 = 1 : length(INspatFRvects2)
    cohinL = spatcohINFR(at1,:);
    cohinR = spatcohINFR(at1,:);
    cohinL(isnan(cohinL)) = [];
    cohinR(isnan(cohinR)) = [];
    if ~isempty(cohinL) && ~isempty(cohinR)
        spacohdINLR(at1) = abs(log10(ranksum(cohinL,cohinR)));
    end
end

PFA_spatial_metrics = struct('PC_sparsity',sparsityPC,'PC_spike_info',spkinfoPC,'Mean_PC_sparsity',sparsall,...
    'Mean_spike_info',spkinfall,'PC_left_sparsity',sparsityPCL,'PC_right_sparsity',sparsityPCR,...
    'PC_left_spike_info',spkinfoPCL,'PC_right_spike_info',spkinfoPCR,'Sparsity_LeftRight_comparison',sparsLR,...
    'PC_spike_info_LeftRight_comparison',spkinfLR,...
    'IN_sparsity',sparsIN,'IN_spike_info',spkinfIN,'Mean_IN_sparsity',sparsallIN,'Mean_IN_spike_info',spkinfallIN,...
    'IN_left_sparsity',sparsityINL,'IN_right_sparsity',sparsityINR,...
    'IN_left_spike_info',spkinfoINL,'IN_right_spike_info',spkinfoINR,'IN_sparsity_LeftRight_comparison',sparsinLR,...
    'IN_spike_info_LeftRight_comparison',spkinfinLR,...
    'PC_spatial_coherence',spatcohFR,'PC_left_spatial_coherence',spatcohFRL,'PC_right_spatial_coherence',spatcohFRR,...
    'PC_spatial_coherence_LeftRight_comp',spacohdLR,...
    'IN_spatial_coherence',spatcohINFR,'IN_left_spatial_coherence',spatcohINFRL,...
    'IN_right_spatial_coherence',spatcohINFRR,'IN_spatial_coherence_LeftRight_comp',spacohdINLR);
filename_SM = [sessionname{sessnum} '_spatial_metrics.mat'];
save(fullfile(resdir,sessionname{sessnum},filename_SM),'PFA_spatial_metrics')

%%
% Analysis of theta sequences

disp('Analysis of theta sequences:...')
thetavectors = PFA_decoding.ThetaCyclePopVectors;
cyclespeed = PFA_decoding.SpeedInCycles;
leftspeed = cyclespeed(purelefts);
rightspeed = cyclespeed(purerights);
lesv = cat(2,leftspeed{:});
risv = cat(2,rightspeed{:});
leftvectors = thetavectors(purelefts);
rightvectors = thetavectors(purerights);
leftv = cat(2,leftvectors{:});
rightv = cat(2,rightvectors{:});
checkLv = sum(leftv);
checkRv = sum(rightv);
Lev = leftv(:,checkLv > 0);
Riv = rightv(:,checkRv > 0);
LevInd = find(checkLv > 0);
RivInd = find(checkRv > 0);
lesvnz = lesv(checkLv > 0);
risvnz = risv(checkRv > 0);
allrv = [Lev';Riv'];
allsv = [lesvnz,risvnz];
numclu = 2 : 100; %iteration from 2 to 100 clusters
avesil = zeros(1,length(numclu));
for kj = 1 : length(numclu)
    disp(['Testing cluster number: ',num2str(numclu(kj))])
    idx = kmeans(allrv,numclu(kj),'Distance','correlation','MaxIter',100,'Replicates',20);
    sil = silhouette(allrv,idx,'correlation');
    avesil(kj) = mean(sil);
end
optimalCluNum = find(avesil == max(avesil));
%final clustering step
fidx = kmeans(allrv,optimalCluNum,'Distance','correlation','MaxIter',100,'Replicates',20);
leftFidx = fidx(1 : size(Lev,2));
rightFidx = fidx(size(Lev,2) + 1 : end);
%separation of theta cycle PV clusters based on position
ThetaCyclePositions = PFA_decoding.PositionBinInCycles;
leftTCP = ThetaCyclePositions(purelefts);
rightTCP = ThetaCyclePositions(purerights);
leftTCP = cat(2,leftTCP{:});
rightTCP = cat(2,rightTCP{:});
leftTCP = leftTCP(checkLv > 0);
rightTCP = rightTCP(checkRv > 0);
allTCPnz = [leftTCP,rightTCP];
%calculation of cluster sidedness, location and speed
spatfullbins = 0 : 2 : 205;
cluspeed = cell(1,optimalCluNum);
cluLR = zeros(3,optimalCluNum);
clulocbins = 1 : length(spatfullbins);
speedbins = 0 : 2 : 100;
tcld = zeros(optimalCluNum,length(clulocbins) - 1);
aveTCPVs = zeros(optimalCluNum,size(allrv,2));
sdTCPVs = zeros(optimalCluNum,size(allrv,2));
TC_PV_clu = cell(1,optimalCluNum);
clsd = zeros(optimalCluNum,length(speedbins) - 1);
medcluspeed = zeros(1,optimalCluNum);
for qc = 1 : optimalCluNum
    cluLR(1,qc) = length(leftFidx(leftFidx == qc));
    cluLR(2,qc) = length(rightFidx(rightFidx == qc));
    cluLR(3,qc) = cluLR(1,qc) / cluLR(2,qc);
    tcluloc = allTCPnz(fidx == qc);
    tcld(qc,:) = histcounts(tcluloc,clulocbins,'Normalization','probability');
    TC_PV_clu{qc} = allrv(fidx == qc,:);
    aveTCPVs(qc,:) = mean(allrv(fidx == qc,:),'omitnan')';
    sdTCPVs(qc,:) = std(allrv(fidx == qc,:),'omitnan')';
    cluspeed{qc} = allsv(fidx == qc);
    clsd(qc,:) = histcounts(cluspeed{qc},speedbins,'Normalization','probability');
    medcluspeed(qc) = median(cluspeed{qc},'omitnan');
end
%Theta cycle participation of PCs
PC_TC_particip = zeros(size(allrv,2),optimalCluNum,2);
PC_TC_domclu = cell(2,size(allrv,2));
for qd = 1 : size(allrv,2)
    for qdd = 1 : length(TC_PV_clu)
        PC_TC_particip(qd,qdd,1) = length(find(TC_PV_clu{qdd}(:,qd) > 1));
        PC_TC_particip(qd,qdd,2) = length(find(TC_PV_clu{qdd}(:,qd) > 1)) / size(TC_PV_clu{qdd},1);
    end
    PC_TC_domclu{2,qd} = find(PC_TC_particip(qd,:,2) == max(PC_TC_particip(qd,:,2)));
    PC_TC_domclu{1,qd} = PC_TC_particip(qd,PC_TC_domclu{2,qd},2);
end

thetaCycleCoFiring = zeros(size(allrv,2),size(allrv,2));
for cf1 = 1 : size(allrv,2) - 1
    for cf2 = cf1 + 1 : size(allrv,2)
        coFfinder = zeros(1,size(allrv,1));
        for cf3 = 1 : size(allrv,1)
            if allrv(cf3,cf1) > 1 && allrv(cf3,cf2) > 1
                coFfinder(cf3) = 1;
            end
        end
        thetaCycleCoFiring(cf1,cf2) = (length(coFfinder(coFfinder == 1)) / length(coFfinder)) * 100;
    end
end

PFA_TC_PV_matrix = struct('TC_population_vectors',allrv,'Left_vectors',Lev,'Right_vectors',Riv);
PFA_TC_PV_clusters = struct('TC_population_vector_clusters',TC_PV_clu);
PFA_TC_PV_speeds = struct('speed_in_PV_clusters',cluspeed);
PFA_TC_PV_cluster_charact = struct('silhouette_scores',avesil,'optimal_cluster_number',optimalCluNum,'cluster_identifiers',fidx,...
    'left_identifiers',leftFidx,'right_identifiers',rightFidx,'PV_cluster_context_bias',cluLR,'PV_cluster_spatial_distribution',tcld,...
    'the_clusters',aveTCPVs,'speed_distr_in_clusters',clsd,'median_speed_in_clusters',medcluspeed,'participation_in_PV_clusters',PC_TC_particip,...
    'ThCycleCoFiring',thetaCycleCoFiring);
PFA_TC_PV_particip_cluident = struct('Dominant_TC_PV_cluster_perCell',PC_TC_domclu);
filenameTCPVm = [sessionname{sessnum} '_TC_PVs.mat'];
filenameTCPVclu = [sessionname{sessnum} '_TC_PV_clusters.mat'];
filenameTCPVcluSpeed = [sessionname{sessnum} '_TC_PV_clusterspeed.mat'];
filenameTCPVcluCharact = [sessionname{sessnum} '_TC_PV_cluster_characteristics.mat'];
filenameTCPVcluIdent = [sessionname{sessnum} '_TC_PV_domcluster_identity.mat'];
save(fullfile(resdir,sessionname{sessnum},filenameTCPVm),'PFA_TC_PV_matrix')
save(fullfile(resdir,sessionname{sessnum},filenameTCPVclu),'PFA_TC_PV_clusters')
save(fullfile(resdir,sessionname{sessnum},filenameTCPVcluSpeed),'PFA_TC_PV_speeds')
save(fullfile(resdir,sessionname{sessnum},filenameTCPVcluCharact),'PFA_TC_PV_cluster_charact')
save(fullfile(resdir,sessionname{sessnum},filenameTCPVcluIdent),'PFA_TC_PV_particip_cluident')

%%
%autocorrelation of splitter trials: place cells
disp('Autocorrelation-based measures...')

maxlag = 1; %second
lag = 0.005; %seconds
nfft = 512;
nsfft = 2^15;
fss = 200;
freqres = 200 / nfft;
freqsres = sro / nsfft;
AClen = (2 * maxlag) / lag;
spectwin = blackmanharris(1000); %50 ms kernel for SDF
ACwin = blackmanharris(1000); %50 ms kernel for AC
afc = fir1(128,2/100,'high');
splitAC = cell(1,length(sstimes));
splitfAC = cell(1,length(sstimes));
splitsAC = cell(1,length(sstimes));
splitSpectr = cell(1,length(sstimes));
spkrasters = cell(1,length(sstimes));
for sa3 = 1 : length(sstimes)
    PFR_AC = cell(1,size(stemregions{sa3},2));
    PFR_fAC = cell(1,size(stemregions{sa3},2));
    PFR_spect = cell(1,size(stemregions{sa3},2));
    domfreqL = zeros(3,size(stemregions{sa3},2));
    domfreqR = zeros(3,size(stemregions{sa3},2));
    domfreqLR = zeros(3,size(stemregions{sa3},2),2);
    spikerastercell = cell(1,size(stemregions{sa3},2));
    for sa4 = 1 : size(stemregions{sa3},2)
        %left trials
        pureLpf = purelefts(purelefts >= stemregions{sa3}(3,sa4) & purelefts <= stemregions{sa3}(4,sa4)) - (stemregions{sa3}(3,sa4) - 1);
        pureRpf = purerights(purerights >= stemregions{sa3}(3,sa4) & purerights <= stemregions{sa3}(4,sa4)) - (stemregions{sa3}(3,sa4) - 1);
        LPFspikes = sstimes{sa3}(pureLpf);
        RPFspikes = sstimes{sa3}(pureRpf);
        LPFspikespos = sspos{sa3}(pureLpf);
        RPFspikespos = sspos{sa3}(pureRpf);
        sprialacL = zeros(length(pureLpf),AClen);
        trialLPFspect = zeros(length(pureLpf),100);
        trialRPFspect = zeros(length(pureRpf),100);
        LRspikecell = cell(1,2);
        Lspikeraster = zeros(length(pureLpf),2 * sro);
        for sa5 = 1 : length(pureLpf)
            splitarmspk1 = LPFspikes{sa5}(LPFspikespos{sa5} >= stemregions{sa3}(1,sa4)*2 & LPFspikespos{sa5} <= stemregions{sa3}(2,sa4)*2);
            if ~isempty(splitarmspk1)
                tLspk = round((splitarmspk1 - splitarmspk1(1)) * sro);
                tLspk(1) = 1;
                Lspikeraster(sa5,tLspk) = 1;
            end
            
%             spectLspikes = fix(splitarmspk1 * sro);
%             if length(spectLspikes) >= 3
%                 ssLvect = zeros(1,spectLspikes(end) - spectLspikes(1));
%                 ssLcv = zeros(1,spectLspikes(end) - spectLspikes(1));
%                 ssLcv(fix(length(ssLcv) / 2)) = 1;
%                 spectLspikes = spectLspikes - spectLspikes(1);
%                 spectLspikes(1) = 1;
%                 ssLvect(spectLspikes) = 1;
%                 sdfL = conv(ssLvect,spectwin);
%                 ppxL = pwelch(sdfL,[],[],nsfft);
%                 trialLPFspect(sa5,:) = ppxL(1 : 100);
%             else
%                 trialLPFspect(sa5,:) = NaN;
%             end
            if length(splitarmspk1) >= 3
                saccountL = zeros(length(splitarmspk1),AClen);
                for sa6 = 1 : length(splitarmspk1)
                    for sa7 = 1 : AClen
                        saccountL(sa6,sa7) = length(find(splitarmspk1 >= splitarmspk1(sa6) - maxlag + (sa7-1) * lag & splitarmspk1 < splitarmspk1(sa6) - maxlag + sa7 * lag));
                    end
                end
                if isnan(splitarmspk1) == 0 & size(splitarmspk1,1) > 5
                    sprialacL(sa5,:) = sum(saccountL) / size(LPFspikes{sa5},1);
                else
                    sprialacL(sa5,:) = zeros(1,AClen);
                end
            else
                sprialacL(sa5,:) = zeros(1,AClen);
            end
        end
        
        orixACL = 1 : size(sprialacL,2);
        newxACL = 1 : 1 / 100 : size(sprialacL,2);
        upvectACL = interp1(orixACL,sprialacL',newxACL);
        upvectACL = upvectACL';
        smoothACL = zeros(size(upvectACL,1),size(upvectACL,2));
        for sb = 1 : size(upvectACL,1)
            tACL = conv(upvectACL(sb,:),ACwin);
            smoothACL(sb,:) = tACL(length(ACwin) / 2 + 1 : size(upvectACL,2) + length(ACwin) / 2);
        end
        smoothACL = smoothACL(:,1 : 100 : end);
        fsmoothACL = filtfilt(afc,1,smoothACL');
        fsmoothACL = fsmoothACL';
        
        if size(sprialacL,1) > 1
            ttacL = mean(sprialacL);
        else
            ttacL = sprialacL;
        end
        ttacL(201) = 0;
        fttacL = filtfilt(afc,1,ttacL);
        pxL = pwelch(fttacL,[],[],nfft);
        psxL = pwelch(fsmoothACL',[],[],nfft);
        psxL = psxL';
        ACLtrialspect = zeros(size(psxL,1),2);
        for sc = 1 : size(psxL,1)
            tacloc = find(psxL(sc,:) == max(psxL(sc,:)));
            if ~isempty(tacloc)
                tacloc = tacloc(1);
                ACLtrialspect(sc,:) = [tacloc * fss / nfft,...
                   psxL(sc,tacloc) / sum(psxL(sc,:))];
            else
                ACLtrialspect(sc,:) = [NaN,NaN];
            end
        end
        
        maxfLloc = find(pxL == max(pxL));
        maxfLloc = maxfLloc(1);
        domfreqL(:,sa4) = [maxfLloc * freqres;pxL(maxfLloc);pxL(maxfLloc) / sum(pxL)];

        %right trials
        sprialacR = zeros(length(pureRpf),AClen);
        Rspikeraster = zeros(length(pureRpf),2 * sro);
        for sa8 = 1 : length(pureRpf)
            splitarmspk2 = RPFspikes{sa8}(RPFspikespos{sa8} >= stemregions{sa3}(1,sa4)*2 & RPFspikespos{sa8} <= stemregions{sa3}(2,sa4)*2);
            if ~isempty(splitarmspk2)
                tRspk = round((splitarmspk2 - splitarmspk2(1)) * sro);
                tRspk(1) = 1;
                Rspikeraster(sa8,tRspk') = 1;
            end
            
%             spectRspikes = fix(splitarmspk2 * sro);
%             if length(spectRspikes) >= 3
%                 ssRvect = zeros(1,spectRspikes(end) - spectRspikes(1));
%                 ssRcv = zeros(1,spectRspikes(end) - spectRspikes(1));
%                 ssRcv(fix(length(ssRcv) / 2)) = 1;
%                 spectRspikes = spectRspikes - spectRspikes(1);
%                 spectRspikes(1) = 1;
%                 ssRvect(spectRspikes) = 1;
%                 sdfR = conv(ssRvect,spectwin);
%                 ppxR = pwelch(sdfR,[],[],nsfft);
%                 trialRPFspect(sa8,:) = ppxR(1 : 100);
%             else
%                 trialRPFspect(sa8,:) = NaN;
%             end
            if length(splitarmspk2) >= 3
                saccountR = zeros(length(splitarmspk2),AClen);
                for sa9 = 1 : length(splitarmspk2)
                    for sa10 = 1 : AClen
                        saccountR(sa9,sa10) = length(find(splitarmspk2 >= splitarmspk2(sa9) - maxlag + (sa10-1) * lag & splitarmspk2 < splitarmspk2(sa9) - maxlag + sa10 * lag));
                    end
                end
                if ~isnan(splitarmspk2) & size(splitarmspk2,1) > 5
                    sprialacR(sa8,:) = sum(saccountR) / size(RPFspikes{sa8},1);
                else
                    sprialacR(sa8,:) = zeros(1,AClen);
                end
            else
                sprialacR(sa8,:) = zeros(1,AClen);
            end
        end
        
        orixACR = 1 : size(sprialacR,2);
        newxACR = 1 : 1 / 100 : size(sprialacR,2);
        upvectACR = interp1(orixACL,sprialacR',newxACR);
        upvectACR = upvectACR';
        smoothACR = zeros(size(upvectACR,1),size(upvectACR,2));
        for sb = 1 : size(upvectACR,1)
            tACR = conv(upvectACR(sb,:),ACwin);
            smoothACR(sb,:) = tACR(length(ACwin) / 2 + 1 : size(upvectACR,2) + length(ACwin) / 2);
        end
        smoothACR = smoothACR(:,1 : 100 : end);
        fsmoothACR = filtfilt(afc,1,smoothACR');
        fsmoothACR = fsmoothACR';
        
        if size(sprialacR,1)
            ttacR = mean(sprialacR);
        else
            ttacR = sprialacR;
        end
        ttacR(201) = 0;
        fttacR = filtfilt(afc,1,ttacR);
        pxR = pwelch(fttacR,[],[],nfft);
        psxR = pwelch(fsmoothACR',[],[],nfft);
        psxR = psxR';
        ACRtrialspect = zeros(size(psxR,1),2);
        for sc = 1 : size(psxR,1)
            tacloc = find(psxR(sc,:) == max(psxR(sc,:)));
            if ~isempty(tacloc)
                tacloc = tacloc(1);
                ACRtrialspect(sc,:) = [tacloc * fss / nfft,...
                   psxR(sc,tacloc) / sum(psxR(sc,:))];
            else
                ACRtrialspect(sc,:) = [NaN,NaN];
            end
        end
        
        maxfRloc = find(pxR == max(pxR));
        maxfRloc = maxfRloc(1);
        domfreqR(:,sa4) = [maxfRloc * freqres;pxR(maxfRloc);pxR(maxfRloc) / sum(pxR)];
        domfreqLR(:,:,1) = domfreqL;
        domfreqLR(:,:,2) = domfreqR;
        LRspikecell{1} = Lspikeraster;
        LRspikecell{2} = Rspikeraster;
        spikerastercell{sa4} = LRspikecell;
        
        PFR_AC{sa4} = [ttacL;ttacR];
        PFR_fAC{sa4} = [fttacL;fttacR];
        PFR_spect{sa4} = [pxL,pxR];
    end
    splitAC{sa3} = PFR_AC;
    splitfAC{sa3} = PFR_fAC;
    splitsAC{sa3} = PFR_spect;
    splitSpectr{sa3} = domfreqLR;
    spkrasters{sa3} = spikerastercell;
end

%autocorrelation analysis of interneurons

IN_AC_mtx = cell(2,length(INtarm));
IN_AC_vect = zeros(2,AClen,length(INtarm));
IN_AC_spect = zeros(2,nfft / 2 + 1,length(INtarm));
for ia3 = 1 : length(INtarm)
    %left trials
    leftINt = INtarm{ia3}(purelefts);
    INLac = zeros(length(purelefts),AClen);
    for ia4 = 1 : length(purelefts)
        leftINspk = leftINt{ia4};
        if length(leftINspk) >= 3
                INaccountL = zeros(length(leftINspk),AClen);
                for ia5 = 1 : length(leftINspk)
                    for ia6 = 1 : AClen
                        INaccountL(ia5,ia6) = length(find(leftINspk >= leftINspk(ia5) - maxlag + (ia6-1) * lag & leftINspk < leftINspk(ia5) - maxlag + ia6 * lag));
                    end
                end
                if ~isnan(leftINspk) & size(leftINspk,1) > 5
                    INLac(ia4,:) = sum(INaccountL) / size(leftINspk,1);
                else
                    INLac(ia4,:) = zeros(1,AClen);
                end
            else
                INLac(ia4,:) = zeros(1,AClen);
        end
    end
    INLac(:,201) = 0;
    IN_AC_mtx{1,ia3} = INLac;
    
    mINLac = nanmean(INLac);
    fmINLac = filtfilt(afc,1,mINLac);
    pINL = pwelch(fmINLac,[],[],nfft);
    IN_AC_vect(1,:,ia3) = fmINLac;
    IN_AC_spect(1,:,ia3) = pINL;
    
    %right trials
    rightINt = INtarm{ia3}(purerights);
    INRac = zeros(length(purerights),AClen);
    for ia4 = 1 : length(purerights)
        rightINspk = rightINt{ia4};
        if length(rightINspk) >= 3
                INaccountR = zeros(length(rightINspk),AClen);
                for ia5 = 1 : length(rightINspk)
                    for ia6 = 1 : AClen
                        INaccountR(ia5,ia6) = length(find(rightINspk >= rightINspk(ia5) - maxlag + (ia6-1) * lag & rightINspk < rightINspk(ia5) - maxlag + ia6 * lag));
                    end
                end
                if ~isnan(rightINspk) & size(rightINspk,1) > 5
                    INRac(ia4,:) = sum(INaccountR) / size(rightINspk,1);
                else
                    INRac(ia4,:) = zeros(1,AClen);
                end
            else
                INRac(ia4,:) = zeros(1,AClen);
        end
    end
    INRac(:,201) = 0;
    IN_AC_mtx{2,ia3} = INRac;
    
    mINRac = nanmean(INRac);
    fmINRac = filtfilt(afc,1,mINRac);
    pINR = pwelch(fmINRac,[],[],nfft);
    IN_AC_vect(2,:,ia3) = fmINRac;
    IN_AC_spect(2,:,ia3) = pINR;
end

PFA_autocorrelation = struct('PC_PFZone_AC',splitAC,'PC_PFZone_filtered_AC',splitfAC,...
    'PC_PFZone_AC_spectra',splitsAC,'PC_dominant_mod_frequency',splitSpectr,...
    'PC_AC_raster',spkrasters,'IN_AC',IN_AC_vect,'IN_AC_spectra',IN_AC_spect);
filename_AC = [sessionname{sessnum} '_autocorrelation.mat'];
save(fullfile(resdir,sessionname{sessnum},filename_AC),'PFA_autocorrelation','-v7.3')

%%
%crosscorrelation analysis of IN - Pyr coupling in PF zones
%...global CC
disp('Crosscorrelation-based connectivity...')

maxlag = 1; %second
maxxlag = 0.1; %second
lag = 0.001; %seconds
fss = 1000;
ACxlen = (2 * maxxlag) / lag;
afcx = fir1(256,2/500,'high');
INxPmaster = cell(1,length(INtarm));
signifmaster = cell(1,length(INtarm));

cellsWithZones = cellswithfields;
cellsWithZones(cellsWithZones == 0) = [];
cellsWithStemZones = cellsWithZones(stemPC_ID);
sftimes = times(cellsWithStemZones);
masterCC = zeros(length(sftimes),ACxlen,length(INtarm));
mup99 = zeros(length(sftimes),ACxlen,length(INtarm));
for mc1 = 1 : length(INtarm)
    INspikes = timeIN{mc1};
    PIm = zeros(length(sftimes),ACxlen);
    upper99 = zeros(length(sftimes),ACxlen);
    for mc2 = 1 : length(sftimes)
        disp(['Correlating: ','PC',num2str(mc2),' with IN',num2str(mc1)])
        PCspikes = sftimes{mc2};
        PCxIN = zeros(length(sftimes{mc2}),ACxlen);
        for mc3 = 1 : length(sftimes{mc2})
            CCxbins = PCspikes(mc3)-maxxlag:lag:PCspikes(mc3)+maxxlag;
            PCxIN(mc3,:) = histcounts(INspikes,CCxbins);
%             for mc4 = 1 : ACxlen
%                 PCxIN(mc3,mc4) = length(find(INspikes >= PCspikes(mc3) - maxxlag + (mc4-1) * lag & INspikes < PCspikes(mc3) - maxxlag + mc4 * lag));
%             end
        end
        shCC = zeros(size(timeshIN{mc1},1),ACxlen);
        for me1 = 1 : size(timeshIN{mc1},1)
            shINspikes = timeshIN{mc1}(me1,:);
            tempCC = zeros(length(sftimes{mc2}),ACxlen);
            for me2 = 1 : length(sftimes{mc2})
                CCbins = PCspikes(me2)-maxxlag:lag:PCspikes(me2)+maxxlag;
                tempCC(me2,:) = histcounts(shINspikes,CCbins);
%                 for me3 = 1 : ACxlen
%                     tempCC(me2,me3) = length(find(shINspikes >= PCspikes(me2) - maxxlag + (me3-1) * lag & shINspikes < PCspikes(me2) - maxxlag + me3 * lag));
%                 end
            end
            shCC(me1,:) = mean(tempCC,'omitnan');
        end
        for me4 = 1 : ACxlen
            upper99(mc2,me4) = quantile(shCC(:,me4),0.99);
        end
        PIm(mc2,:) = mean(PCxIN,'omitnan');
    end
    masterCC(:,:,mc1) = PIm;
    mup99(:,:,mc1) = upper99;
end

PFA_PC_IN_fullCC = struct('IN_PC_full_crosscorrelations',masterCC,'shuffling_based_significance',mup99);
filename = [sessionname{sessnum} '_PC_IN_full_CC.mat'];
save(fullfile(resdir,sessionname{sessnum},filename),'PFA_PC_IN_fullCC')

%%
%...PF CC
for is1 = 1 : length(INtarm)
%     disp(['IN ',num2str(is1)])
    INxPCppc = cell(1,length(sstimes));
    PFsgnl = cell(1,length(sstimes));
    for is2 = 1 : length(sstimes)
%         disp(['PC ',num2str(is2)])
        INxPLpf = zeros(size(stemregions{is2},2),ACxlen);
        INxPRpf = zeros(size(stemregions{is2},2),ACxlen);
        PFsignlevels = zeros(4,ACxlen,size(stemregions{is2},2));
        if ~isempty(stemregions{is2})
        for is3 = size(stemregions{is2},2)
%             disp(['PF ',num2str(is3)])
            pureLpf = purelefts(purelefts >= stemregions{is2}(3,is3) & purelefts <= stemregions{is2}(4,is3)) - (stemregions{is2}(3,is3) - 1);
            pureRpf = purerights(purerights >= stemregions{is2}(3,is3) & purerights <= stemregions{is2}(4,is3)) - (stemregions{is2}(3,is3) - 1);
            Leftpos = alltrials(pureLpf);
            Rightpos = alltrials(pureRpf);
            Leftimes = alltrialtimes(pureLpf);
            Rightimes = alltrialtimes(pureRpf);
            LPFspikes = sstimes{is2}(pureLpf);
            RPFspikes = sstimes{is2}(pureRpf);
            LPFspikespos = sspos{is2}(pureLpf);
            RPFspikespos = sspos{is2}(pureRpf);
            INLspikes = INtarm{is1}(pureLpf);
            INRspikes = INtarm{is1}(pureRpf);
            INLspikepos = INarmpos{is1}(pureLpf);
            INRspikepos = INarmpos{is1}(pureRpf);
            shINLspikes = shINtarm{is1}(pureLpf);
            shINRspikes = shINtarm{is1}(pureRpf);
            INxPLm = zeros(length(pureLpf),ACxlen);
            signrangesL = zeros(2,ACxlen,length(pureLpf));
            signrangesR = zeros(2,ACxlen,length(pureRpf));
            %left trials
            for is4 = 1 : length(pureLpf) %trial cycle
                %defining time boundaries of place field zones
                PFLtimebounds = Leftimes{is4}(Leftpos{is4} >= stemregions{is2}(1,is3)*2 & Leftpos{is4} <= stemregions{is2}(2,is3)*2);
                splitarmspkL = LPFspikes{is4}(LPFspikespos{is4} >= stemregions{is2}(1,is3)*2 & LPFspikespos{is4} <= stemregions{is2}(2,is3)*2);
                INLpfspk = INLspikes{is4}(INLspikepos{is4} >= stemregions{is2}(1,is3)*2 & INLspikepos{is4} <= stemregions{is2}(2,is3)*2);
                if length(splitarmspkL) >= 3 & length(INLpfspk) >= 30
                    %real
                    INxPL = zeros(length(splitarmspkL),ACxlen);
                    %CC calculation
                    for is5 = 1 : length(splitarmspkL)
                        for is6 = 1 : ACxlen
                            INxPL(is5,is6) = length(find(INLpfspk >= splitarmspkL(is5) - maxxlag + (is6-1) * lag & INLpfspk < splitarmspkL(is5) - maxxlag + is6 * lag));
                        end
                    end
                    %shuffled
                    shINxPL = zeros(length(splitarmspkL),ACxlen,size(shINLspikes{is4},2));
                    for sis1 = 1 : size(shINLspikes{is4},2)
                        shINLtrain = shINLspikes{is4}{sis1}(shINLspikes{is4}{sis1} >= PFLtimebounds(1) & shINLspikes{is4}{sis1} <= PFLtimebounds(end));
                        %CC calculation
                        for sis2 = 1 : length(splitarmspkL)
                            for sis3 = 1 : ACxlen
                                shINxPL(sis2,sis3,sis1) = length(find(shINLtrain >= splitarmspkL(sis2) - maxxlag + (sis3-1) * lag & shINLtrain < splitarmspkL(sis2) - maxxlag + sis3 * lag));
                            end
                        end
                    end
                    tsrL = nanmean(shINxPL,1);
                    tsrL = squeeze(tsrL)';
                    signrangesL(:,:,is4) = quantile(tsrL,[0.01 0.99],1);
                    if isempty(splitarmspkL(isnan(splitarmspkL) == 1)) == 1 & size(splitarmspkL,1) > 3
                        INxPLm(is4,:) = nanmean(INxPL);
                    else
                        INxPLm(is4,:) = zeros(1,ACxlen);
                    end
                else
                    INxPLm(is4,:) = zeros(1,ACxlen);
                end
            end
            tINxPL = nanmean(INxPLm);
            INxPLpf(is3,:) = tINxPL;
            PFsignlevels(1 : 2,:,is3) = squeeze(nanmean(signrangesL,3));
            %right trials
            INxPRm = zeros(length(pureRpf),ACxlen);
            for is4 = 1 : length(pureRpf)
                PFRtimebounds = Rightimes{is4}(Rightpos{is4} >= stemregions{is2}(1,is3)*2 & Rightpos{is4} <= stemregions{is2}(2,is3)*2);
                splitarmspkR = RPFspikes{is4}(RPFspikespos{is4} >= stemregions{is2}(1,is3)*2 & RPFspikespos{is4} <= stemregions{is2}(2,is3)*2);
                INRpfspk = INRspikes{is4}(INRspikepos{is4} >= stemregions{is2}(1,is3)*2 & INRspikepos{is4} <= stemregions{is2}(2,is3)*2);
                if length(splitarmspkR) >= 3 & length(INRpfspk) >= 30
                    %real
                    INxPR = zeros(length(splitarmspkR),ACxlen);
                    %CC calculation
                    for is5 = 1 : length(splitarmspkR)
                        for is6 = 1 : ACxlen
                            INxPR(is5,is6) = length(find(INRpfspk >= splitarmspkR(is5) - maxxlag + (is6-1) * lag & INRpfspk < splitarmspkR(is5) - maxxlag + is6 * lag));
                        end
                    end
                    %shuffled
                    shINxPR = zeros(length(splitarmspkR),ACxlen,size(shINRspikes{is4},1));
                    for sis1 = 1 : size(shINRspikes{is4},2)
                        shINRtrain = shINRspikes{is4}{sis1}(shINRspikes{is4}{sis1} >= PFRtimebounds(1) & shINRspikes{is4}{sis1} <= PFRtimebounds(end));
                        %CC calculation
                        for sis2 = 1 : length(splitarmspkR)
                            for sis3 = 1 : ACxlen
                                shINxPR(sis2,sis3,sis1) = length(find(shINRtrain >= splitarmspkR(sis2) - maxxlag + (sis3-1) * lag & shINRtrain < splitarmspkR(sis2) - maxxlag + sis3 * lag));
                            end
                        end
                    end
                    tsrR = nanmean(shINxPR,1);
                    tsrR = squeeze(tsrR)';
                    signrangesR(:,:,is4) = quantile(tsrR,[0.01 0.99],1);
                    if isempty(splitarmspkR(isnan(splitarmspkR) == 1)) == 1 & size(splitarmspkR,1) > 3
                        INxPRm(is4,:) = nanmean(INxPR);
                    else
                        INxPRm(is4,:) = zeros(1,ACxlen);
                    end
                else
                    INxPRm(is4,:) = zeros(1,ACxlen);
                end
            end
            tINxPR = nanmean(INxPRm);
            INxPRpf(is3,:) = tINxPR;
            PFsignlevels(3 : 4,:,is3) = squeeze(nanmean(signrangesR,3));
        end
        INxPCppc{is2} = cat(3,INxPLpf,INxPRpf);
        PFsgnl{is2} = PFsignlevels;
        else
            INxPCppc{is2} = NaN;
            PFsgnl{is2} = NaN;
        end
    end
    INxPmaster{is1} = INxPCppc; %IN-PC crosscorrelations
    signifmaster{is1} = PFsgnl; %shuffling-based significance levels
end

PFA_PC_IN_crosscorrelations = struct('IN_PC_crosscorrelations',INxPmaster,'shuffling_based_significance_levels',signifmaster);
filename = [sessionname{sessnum} '_PC_IN_crosscorrelations.mat'];
save(fullfile(resdir,sessionname{sessnum},filename),'PFA_PC_IN_crosscorrelations')