function [PCT,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] = placeTMdecoderRH(SoFM,how2DC,selectedTimes,titp,thetaCycles,seltimes,selspktimes,selvects,preselvects,...
    constrialnum,purelefts,purerights,alltrialtimes,alltrials,trialspeeds,spatbins2,spatfullbins,allChoTrials,LChoTrials,RChoTrials)

%template-based decoding of position: templates are generated from spatial
%firing rate vectors, decoding windows correspond to theta cycles. Decoding
%options (how2DC): Full: all FR vectors of the session are used; Partial with ratio:
%only a predetermined proportion of trials are used. Trials are randomly
%drawn from the session; Left: templates are built from left trials; Right:
%right trials contribute to the template.
%Decoding can also be restricted to the stem (SoFM == 's') or the full maze
%can be decoded (SoFM == 'f'); This function is used by
%placefieldanalyzer_RH2.m

% selection of theta peaks based on inst. phase
TaT = selectedTimes;
srnB = 1250;
if ~isempty(titp)
    if iscell(titp)
        titp = cat(2,titp{:});
    end
    dang = diff(titp);
    reversals = find(dang < -6);
    revinsec = reversals / srnB;
    revonmaze = revinsec(revinsec > TaT(1) & revinsec < TaT(end));
%     revonmaze = revonmaze - TaT(1);
else
    revonmaze = thetaCycles;
    revonmaze = revonmaze(revonmaze > TaT(1) & revonmaze < TaT(end));
end

%%
%within theta cycle spike times

PCT = cell(1,size(seltimes,2));
% PCP = cell(1,size(seltimes,2));
for ct1 = 1 : size(seltimes,2)
    percycletime = cell(1,size(revonmaze,2) - 1);
%     percyclephase = cell(1,size(revonmaze,2) - 1);
    for ct2 = 1 : size(revonmaze,2) - 1
        percycletime{ct2} = seltimes{ct1}(seltimes{ct1} > revonmaze(ct2) & seltimes{ct1} < revonmaze(ct2 + 1));
%         percyclephase{ct2} = selphase{ct1}(seltimes{ct1} > revonmaze(ct2) & seltimes{ct1} < revonmaze(ct2 + 1));
    end
    PCT{ct1} = percycletime;
%     PCP{ct1} = percyclephase;
end

%temporal order of per-cycle firing
rearrthetatimes = cell(1,size(revonmaze,2));
firstspiketimes = zeros(size(seltimes,2),size(revonmaze,2));
ofst = zeros(size(seltimes,2),size(revonmaze,2));
thetatimeorders = zeros(size(seltimes,2),size(revonmaze,2));
percyclespknum = zeros(size(seltimes,2),size(revonmaze,2) - 1);
for ct3 = 1 : size(revonmaze,2) - 1
    pco = cell(1,size(seltimes,2));
    for ct4 = 1 : size(seltimes,2)
        pco{ct4} = PCT{ct4}{ct3};
        if isempty(pco{ct4}) == 0
            firstspiketimes(ct4,ct3) = pco{ct4}(1);
        else
            firstspiketimes(ct4,ct3) = NaN;
        end
    end
    [so,thetatimeorders(:,ct3)] = sort(firstspiketimes(:,ct3));
    ofst(:,ct3) = firstspiketimes(thetatimeorders(:,ct3),ct3);
    rearrthetatimes{ct3} = pco(thetatimeorders(:,ct3));
    thetatimeorders(isnan(firstspiketimes(:,ct3)) == 1,ct3) = 0;
    for ct5 = 1 : size(seltimes,2)
        percyclespknum(ct5,ct3) = size(PCT{ct5}{ct3},1);
    end
end

%%
%Template-based decoding: theta cycle bins
%Theta cycles per trial

ttarm = cell(1,constrialnum);
alltempos = cell(1,constrialnum);
if strcmp(SoFM,'s')
    for jj1 = 1 : constrialnum
        ttarm{jj1} = alltrialtimes{jj1}(alltrials{jj1} > stemm(1) & alltrials{jj1} < stemm(2));
        alltempos{jj1} = alltrials{jj1}(alltrials{jj1} > stemm(1) & alltrials{jj1} < stemm(2));
    end
elseif strcmp(SoFM,'f')
    for jj1 = 1 : constrialnum
        ttarm{jj1} = alltrialtimes{jj1};
        alltempos{jj1} = alltrials{jj1};
    end
end

timebins = 0 : 0.5 : 5;
% rvm = revonmaze + alltime(1);
rvm = revonmaze;
trialcycles = cell(1,constrialnum);
trialcyclepos = cell(1,constrialnum);
for tj1 = 1 : constrialnum
    trialcycles{tj1} = rvm(rvm > ttarm{tj1}(1) & rvm < ttarm{tj1}(end));
    if ~isempty(trialcycles{tj1})
        ttcy = trialcycles{tj1}(2 : end) - trialcycles{tj1}(1);
    else
        ttcy = [];
    end
    [cyc1,cyc2,trialcyclepos{tj1}] = histcounts(ttcy,timebins);
end

%position and speed in theta cycles
posincycles = cell(1,constrialnum);
poscycleinbins = cell(1,constrialnum);
thetaspeed = cell(1,constrialnum);
for tj2 = 1 : constrialnum
    pincbins = zeros(1,length(trialcycles{tj2}) - 1);
    thspeed = zeros(1,length(trialcycles{tj2}) - 1);
    for tj3 = 1 : length(trialcycles{tj2}) - 1
        tpcb = mean(alltempos{tj2}(ttarm{tj2} >= trialcycles{tj2}(tj3) & ttarm{tj2} < trialcycles{tj2}(tj3 + 1)));
        if ~isempty(tpcb)
            pincbins(tj3) = mean(alltempos{tj2}(ttarm{tj2} >= trialcycles{tj2}(tj3) & ttarm{tj2} < trialcycles{tj2}(tj3 + 1)));
        else
            pincbins(tj3) = NaN;
        end
        thsp = mean(trialspeeds{tj2}(ttarm{tj2} >= trialcycles{tj2}(tj3) & ttarm{tj2} < trialcycles{tj2}(tj3 + 1)));
        if ~isempty(thsp)
            thspeed(tj3) = mean(trialspeeds{tj2}(ttarm{tj2} >= trialcycles{tj2}(tj3) & ttarm{tj2} < trialcycles{tj2}(tj3 + 1)));
        else
            thspeed(tj3) = NaN;
        end
    end
    posincycles{tj2} = pincbins;
    [tcn1,tcn2,poscycleinbins{tj2}] = histcounts(pincbins,spatfullbins);
%     if strcmp(SoFM,'s')
%         [tcn1,tcn2,poscycleinbins{tj2}] = histcounts(pincbins,spatbins2);
%     elseif strcmp(SoFM,'f')
%         [tcn1,tcn2,poscycleinbins{tj2}] = histcounts(pincbins,spatfullbins);
%     end
    thetaspeed{tj2} = thspeed;
end

%fillup population time vectors

% if strcmp(SoFM,'s')
%     selspktimes = unitPTP_arm(clusteredPCS);
% elseif strcmp(SoFM,'f')
%     selspktimes = unitPTP_full(clusteredPCS);
% end

thetavectors = cell(1,constrialnum);
thetaISI = cell(1,constrialnum);
for jt1 = 1 : constrialnum
    tcvm = zeros(length(selspktimes),length(trialcycles{jt1}) - 1);
    tISI = zeros(length(selspktimes),length(trialcycles{jt1}) - 1);
    for jt2 = 1 : length(trialcycles{jt1}) - 1
        for jt3 = 1 : length(selspktimes)
            tcvm(jt3,jt2) = size(selspktimes{jt3}{jt1}(selspktimes{jt3}{jt1} >= trialcycles{jt1}(jt2) & selspktimes{jt3}{jt1} < trialcycles{jt1}(jt2 + 1)),1);
            tISI(jt3,jt2) = mean(diff(selspktimes{jt3}{jt1}(selspktimes{jt3}{jt1} >= trialcycles{jt1}(jt2) & selspktimes{jt3}{jt1} < trialcycles{jt1}(jt2 + 1))));
            if isnan(tISI(jt3,jt2))
                tISI(jt3,jt2) = 0;
            end
        end
    end
    thetavectors{jt1} = tcvm;
    thetaISI{jt1} = tISI;
end

%%
%fillup population spatial vectors and generation of spatial templates
templateTnum = 10; % no of trials for template building
%only stem decoding is no longer used--------------------------------------
% if strcmp(SoFM,'s')
%     spatialbins = spatbins2;
%     switch how2DCm
%         case 'F'
%             templatesize = constrialnum;
%             trialSV = cell(1,constrialnum);
%             for vj1 = 1 : constrialnum
%                 pertrialSM = zeros(length(selvects),length(spatialbins));
%                 for vj2 = 1 : length(spatialbins)
%                     for vj3 = 1 : length(selvects)
%                         pertrialSM(vj3,vj2) = selvects{vj3}(vj2,vj1);
%                     end
%                 end
%                 trialSV{vj1} = pertrialSM;
%             end
%             spatempmtx = zeros(length(seltimes),length(spatbins2) - 1,templatesize);
%             for tg1 = 1 : templatesize
%                 spatempmtx(:,:,tg1) = trialSV{tg1}(:,1 : end - 1);
%             end
%             spatemp = nanmean(spatempmtx,3);
%         case 'P'
%             rndord = round(rand(constrialnum) * constrialnum);
%             notrial = fix(constrialnum * prop);
%             templatesize = notrial;
%             trstep = round(constrialnum / notrial);
%             seltrials = rndord(1 : trstep : constrialnum);
%             seltrials(seltrials == 0) = 1;
%             sselvects = cell(1,length(preselvects));
%             for vjj = 1 : length(preselvects)
%                 sselvects{vjj} = preselvects{vjj}(:,seltrials);
%             end
%             trialSV = cell(1,notrial);
%             for vj1 = 1 : notrial
%                 pertrialSM = zeros(length(sselvects),length(spatialbins));
%                 for vj2 = 1 : length(spatialbins)
%                     for vj3 = 1 : length(sselvects)
%                         pertrialSM(vj3,vj2) = sselvects{vj3}(vj2,vj1);
%                     end
%                 end
%                 trialSV{vj1} = pertrialSM;
%             end
%             spatempmtx = zeros(length(seltimes),length(spatbins2) - 1,templatesize);
%             for tg1 = 1 : templatesize
%                 spatempmtx(:,:,tg1) = trialSV{tg1}(:,1 : end - 1);
%             end
%             spatemp = nanmean(spatempmtx,3);
%         case 'L'
%             templatesize = length(purelefts);
%             sselvects = cell(1,length(preselvects));
%             for vjj = 1 : length(preselvects)
%                 sselvects{vjj} = preselvects{vjj}(:,purelefts);
%             end
%             trialSV = cell(1,length(purelefts));
%             for vj1 = 1 : length(purelefts)
%                 pertrialSM = zeros(length(sselvects),length(spatialbins));
%                 for vj2 = 1 : length(spatialbins)
%                     for vj3 = 1 : length(sselvects)
%                         pertrialSM(vj3,vj2) = sselvects{vj3}(vj2,vj1);
%                     end
%                 end
%                 trialSV{vj1} = pertrialSM;
%             end
%             spatempmtx = zeros(length(seltimes),length(spatbins2) - 1,templatesize);
%             for tg1 = 1 : templatesize
%                 spatempmtx(:,:,tg1) = trialSV{tg1}(:,1 : end - 1);
%             end
%             spatemp = nanmean(spatempmtx,3);
%         case 'R'
%             templatesize = length(purerights);
%             sselvects = cell(1,length(preselvects));
%             for vjj = 1 : length(preselvects)
%                 sselvects{vjj} = preselvects{vjj}(:,purerights);
%             end
%             trialSV = cell(1,length(purerights));
%             for vj1 = 1 : length(purerights)
%                 pertrialSM = zeros(length(sselvects),length(spatialbins));
%                 for vj2 = 1 : length(spatialbins)
%                     for vj3 = 1 : length(sselvects)
%                         pertrialSM(vj3,vj2) = sselvects{vj3}(vj2,vj1);
%                     end
%                 end
%                 trialSV{vj1} = pertrialSM;
%             end
%             spatempmtx = zeros(length(seltimes),length(spatbins2) - 1,templatesize);
%             for tg1 = 1 : templatesize
%                 spatempmtx(:,:,tg1) = trialSV{tg1}(:,1 : end - 1);
%             end
%             spatemp = nanmean(spatempmtx,3);
%     end
% end
%--------------------------------------------------------------------------

%%
if strcmp(SoFM,'f')
    randAlltrials = NaN;
    randLtrials = NaN;
    randRtrials = NaN;
    switch how2DC
        case 'P'
            templatesize = templateTnum;
            if isnan(allChoTrials)
                randAlltrials = randperm(constrialnum,templatesize);
            else
                randAlltrials = allChoTrials;
            end
            sselvects = selvects;
            spatialbins = spatfullbins;
            trialSV = cell(1,templatesize);
            for vj1 = 1 : templatesize
                pertrialSM = zeros(length(sselvects),length(spatialbins) - 1);
                for vj2 = 1 : length(spatfullbins) - 1
                    for vj3 = 1 : length(sselvects)
                        pertrialSM(vj3,vj2) = sselvects{vj3}(vj2,randAlltrials(vj1));
                    end
                end
                trialSV{vj1} = pertrialSM;
            end
            spatempmtx = zeros(length(sselvects),length(spatfullbins) - 1,templatesize);
            for tg1 = 1 : templatesize
                spatempmtx(:,:,tg1) = trialSV{tg1};
            end
            spatemp = nanmean(spatempmtx,3);
            spatemp(isinf(spatemp)) = NaN;
        case 'L'
            if length(purelefts) <= 5
                templatesize = length(purelefts);
            elseif length(purelefts) > 5 & length(purelefts) <= 10
                templatesize = 5;
            elseif length(purelefts) > 10
                templatesize = 10;
            end
            sselvects = cell(1,length(preselvects));
            if isnan(LChoTrials) & (templatesize == 5 | templatesize == 10)
                randLtrials = randperm(length(purelefts),templatesize);
                for vjj = 1 : length(preselvects)
                    sselvects{vjj} = preselvects{vjj}(:,purelefts(randLtrials));
                end
            elseif isnan(LChoTrials) & templatesize < 5
                randLtrials = purelefts;
                for vjj = 1 : length(preselvects)
                    sselvects{vjj} = preselvects{vjj}(:,randLtrials);
                end
            else
                randLtrials = LChoTrials;
                for vjj = 1 : length(preselvects)
                    sselvects{vjj} = preselvects{vjj}(:,randLtrials);
                end
            end
            trialSV = cell(1,templatesize);
            for vj1 = 1 : templatesize
                pertrialSM = zeros(length(sselvects),length(spatfullbins));
                for vj2 = 1 : length(spatfullbins) - 1
                    for vj3 = 1 : length(sselvects)
                        pertrialSM(vj3,vj2) = sselvects{vj3}(vj2,vj1);
                    end
                end
                trialSV{vj1} = pertrialSM;
            end
            spatempmtx = zeros(length(sselvects),length(spatfullbins) - 1,templatesize);
            for tg1 = 1 : templatesize
                spatempmtx(:,:,tg1) = trialSV{tg1}(:,1 : end - 1);
            end
            spatemp = nanmean(spatempmtx,3);
%             A = spatemp;
%             B = arrayfun(@(x) inf2nan(A),A,'UniformOutput',false);
%             spatemp = B{1};
            spatemp(isinf(spatemp)) = NaN;
        case 'R'
            if length(purerights) <= 5
                templatesize = length(purerights);
            elseif length(purerights) > 5 & length(purerights) <= 10
                templatesize = 5;
            elseif length(purerights) > 10
                templatesize = 10;
            end
            sselvects = cell(1,length(preselvects));
            if isnan(RChoTrials) & (templatesize == 5 | templatesize == 10)
                randRtrials = randperm(length(purerights),templatesize);
                for vjj = 1 : length(preselvects)
                    sselvects{vjj} = preselvects{vjj}(:,purerights(randRtrials));
                end
            elseif isnan(LChoTrials) & templatesize < 5
                randRtrials = purerights;
                for vjj = 1 : length(preselvects)
                    sselvects{vjj} = preselvects{vjj}(:,randRtrials);
                end
            else
                randRtrials = RChoTrials;
                for vjj = 1 : length(preselvects)
                    sselvects{vjj} = preselvects{vjj}(:,randRtrials);
                end
            end
            trialSV = cell(1,templatesize);
            for vj1 = 1 : templatesize
                pertrialSM = zeros(length(sselvects),length(spatfullbins));
                for vj2 = 1 : length(spatfullbins) - 1
                    for vj3 = 1 : length(sselvects)
                        pertrialSM(vj3,vj2) = sselvects{vj3}(vj2,vj1);
                    end
                end
                trialSV{vj1} = pertrialSM;
            end
            spatempmtx = zeros(length(sselvects),length(spatfullbins) - 1,templatesize);
            for tg1 = 1 : templatesize
                spatempmtx(:,:,tg1) = trialSV{tg1}(:,1 : end - 1);
            end
            spatemp = nanmean(spatempmtx,3);
            spatemp(isinf(spatemp)) = NaN;
    end
end

%%
%correlation of the time and spatial vectors and estimation of position
%based on theta cycle vectors: only stem
if strcmp(SoFM,'s')
    %not used (full maze docoding is done instead)
    esthetapos = cell(1,constrialnum);
    decodederror = cell(1,constrialnum);
    quandecerr = zeros(3,constrialnum);
    for it1 = 1 : constrialnum
        thecodedpos = zeros(3,size(thetavectors{it1},2));
        for it2 = 1 : size(thetavectors{it1},2)
            tptx = zeros(1,length(spatbins2) - 1);
            tptp = zeros(1,length(spatbins2) - 1);
            for it3 = 1 : length(spatbins2) - 1
                [tptx(it3),tptp(it3)] = corr(thetavectors{it1}(:,it2),spatemp(:,it3),'Type','spearman');
            end
            tptx(isnan(tptx)) = 0;
            if sum(tptx) == 0
                tpmtpos = NaN;
            else
                tpmtpos = spatbins2(tptx == max(tptx));
            end
            thecodedpos(1,it2) = tpmtpos(1);
            maxcorr = tptx(tptx == max(tptx));
            maxcorrP = abs(log10(tptp(tptx == max(tptx))));
            thecodedpos(2,it2) = maxcorr(1);
            thecodedpos(3,it2) = maxcorrP(1);
        end
        esthetapos{it1} = thecodedpos;
        decodederror{it1} = esthetapos{it1}(1,:) - posincycles{it1};
        nonandecerr = decodederror{it1}(~isnan(decodederror{it1}));
        quandecerr(:,it1) = quantile(nonandecerr,[0.25 0.5 0.75]);
    end
    estposvect = cat(2,esthetapos{:});
    realposvect = cat(2,posincycles{:});
elseif strcmp(SoFM,'f')
%     esthetapos = cell(1,constrialnum);
%     decodederror = cell(1,constrialnum);
%     quandecerr = zeros(3,constrialnum);
    ethP_CS = cell(1,constrialnum);
    de_CS = cell(1,constrialnum);
    qde_CS = zeros(3,constrialnum);
    for it1 = 1 : constrialnum
%         disp(['Decoding trial: ',num2str(it1)])
%         thecodedpos = zeros(3,size(thetavectors{it1},2));
        tcp_CS = zeros(2,size(thetavectors{it1},2));
        for it2 = 1 : size(thetavectors{it1},2)
%             tptx = zeros(1,length(spatfullbins) - 1);
%             tptp = zeros(1,length(spatfullbins) - 1);
            CS = zeros(1,length(spatfullbins) - 1);
            for it3 = 1 : length(spatfullbins) - 1
                %Spearman rank correlation
%                 [tptx(it3),tptp(it3)] = corr(thetavectors{it1}(:,it2),spatemp(:,it3),'Type','spearman');
                %cosine similarity
                CS(it3) = sum(thetavectors{it1}(:,it2) .* spatemp(:,it3)) / (sqrt(sum(thetavectors{it1}(:,it2).^2)) * sqrt(sum(spatemp(:,it3).^2)));
            end
%             tptx(isnan(tptx)) = 0;
%             if sum(tptx) == 0
%                 tpmtpos = NaN;
%             else
%                 tpmtpos = spatfullbins(tptx == max(tptx));
%             end
            CS(isnan(CS)) = 0;
            if sum(CS) == 0
                tcspos = NaN;
            else
                tcspos = spatfullbins(CS == max(CS));
            end
            % Corr.
%             thecodedpos(1,it2) = tpmtpos(1);
%             maxcorr = tptx(tptx == max(tptx));
%             maxcorrP = abs(log10(tptp(tptx == max(tptx))));
%             thecodedpos(2,it2) = maxcorr(1);
%             thecodedpos(3,it2) = maxcorrP(1);
            % CS
            tcp_CS(1,it2) = tcspos(1);
            tcp_CS(2,it2) = max(CS);
        end
%         esthetapos{it1} = thecodedpos;
%         decodederror{it1} = esthetapos{it1}(1,:) - posincycles{it1};
%         nonandecerr = decodederror{it1}(~isnan(decodederror{it1}));
%         quandecerr(1,it1) = quantile(nonandecerr',0.25);
%         quandecerr(2,it1) = quantile(nonandecerr',0.5);
%         quandecerr(3,it1) = quantile(nonandecerr',0.75);
        ethP_CS{it1} = tcp_CS;
        de_CS{it1} = ethP_CS{it1}(1,:) - posincycles{it1};
        nnde_CS = de_CS{it1}(~isnan(de_CS{it1}));
        qde_CS(1,it1) = quantile(nnde_CS',0.25);
        qde_CS(2,it1) = quantile(nnde_CS',0.5);
        qde_CS(3,it1) = quantile(nnde_CS',0.75);
    end
    switch how2DC
        case 'P'
            trials2rem = ones(1,length(ethP_CS));
            trials2rem(randAlltrials) = 0;
            ethP_CS = ethP_CS(trials2rem == 1);
            posincycles = posincycles(trials2rem == 1);
            de_CS = de_CS(trials2rem == 1);
            ttcycPos = trialcyclepos(trials2rem == 1);
            actrials = 1:constrialnum;
            actrials = actrials(trials2rem == 1);
        case 'L'
            trials2rem = ones(1,length(ethP_CS));
            trials2rem(randLtrials) = 0;
            ethP_CS = ethP_CS(trials2rem == 1);
            posincycles = posincycles(trials2rem == 1);
            de_CS = de_CS(trials2rem == 1);
            ttcycPos = trialcyclepos(trials2rem == 1);
            actrials = 1:constrialnum;
            actrials = actrials(trials2rem == 1);
        case 'R'
            trials2rem = ones(1,length(ethP_CS));
            trials2rem(randRtrials) = 0;
            ethP_CS = ethP_CS(trials2rem == 1);
            posincycles = posincycles(trials2rem == 1);
            de_CS = de_CS(trials2rem == 1);
            ttcycPos = trialcyclepos(trials2rem == 1);
            actrials = 1:constrialnum;
            actrials = actrials(trials2rem == 1);
    end
%     estposvect = cat(2,esthetapos{:});
    estPv_CS = cat(2,ethP_CS{:});
    realposvect = cat(2,posincycles{:});
end

timeposdistr = zeros(length(timebins),length(spatfullbins) - 1);
for pd1 = 1 : length(timebins)
    posintimebins = cell(1,length(actrials));
    for pd2 = 1 : length(actrials)
        posintimebins{pd2} = ethP_CS{pd2}(1,ttcycPos{pd2} == pd1);
    end
    ptbv = cat(2,posintimebins{:});
    timeposdistr(pd1,:) = histcounts(ptbv,spatfullbins,'normalization','probability');
end