%decodes trial identity from spiking across bins for each individual neuron
clear all;
load ('R_2R.mat');
load ('RAW.mat');

folds = 5; %number of times cross-validated:
shuffs = 1; %number of shuffled models created

%load parameters
BinDura=R_2R.Param.BinDura;
bins=R_2R.Param.bins;
binint=R_2R.Param.binint;
binstart=R_2R.Param.binstart;



%setup variables
UnitDec=[];
NN = 0;

for i=1:length(RAW) %loops through sessions
    if strcmp('NA',RAW(i).Type(1:2)) | strcmp('VP',RAW(i).Type(1:2)) %check if it's a suc vs mal session
        %events being compared
        RD1=strcmp('RD1', RAW(i).Einfo(:,2));
        RD2=strcmp('RD2', RAW(i).Einfo(:,2));
        for j= 1:size(RAW(i).Nrast,1) %Number of neurons per session
            NN=NN+1; %neuron counter
            for l=1:bins

                LowHz=zeros(1,1);
                DecodeSpikes=NaN(1,1);
                DecodeRs=NaN(1,1);

                [PSR1,N1]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{RD1},[(BinDura(1)+(binstart - binint)+l*binint) (BinDura(2)+(binstart - binint)+l*binint)],{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
                for m=1:length(PSR1)
                    DecodeSpikes(m,1)=sum(PSR1{1,m}>(binstart));
                    DecodeRs(m,1)=1;
                end

                %get all the spikes from reward 1p2 trials
                [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{RD2},[(BinDura(1)+(binstart - binint)+l*binint) (BinDura(2)+(binstart - binint)+l*binint)],{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
                for n=1:length(PSR2)
                    DecodeSpikes(n+m,1)=sum(PSR2{1,n}>(binstart));
                    DecodeRs(n+m,1)=2;
                end

                if sum(DecodeSpikes(1:N1,1))<7 || sum(DecodeSpikes((1+N1):(N1+N2),1))<7 %prevents errors with LDA on conditions with no variance
                    LowHz(1,1)=1;
                end


                %if one neuron or more has too few spikes, so LDA has error because not
                %enough variance, this removes those neurons from the analysis
                if sum(LowHz)==0

                    %creating models
                    CVacc = NaN(folds,1);
                    CVaccSh = NaN(folds,1);
                    %normal model
                    Partitions = cvpartition(DecodeRs,'KFold',folds);
                    for r = 1:folds

                        LDAModel = fitcdiscr(DecodeSpikes(Partitions.training(r),:),DecodeRs(Partitions.training(r)));
                        prediction = predict(LDAModel,DecodeSpikes(Partitions.test(r),:));
                        actual = DecodeRs(Partitions.test(r));
                        correct = prediction - actual;
                        CVacc(r) = sum(correct==0) / length(correct);

                    end

                    %shuffled model
                    for q=1:shuffs
                        DecodeRsSh=DecodeRs(randperm(length(DecodeRs)));
                        PartitionsSh = cvpartition(DecodeRsSh,'KFold',folds);
                        for s = 1:folds
                            LDAModelSh = fitcdiscr(DecodeSpikes(PartitionsSh.training(s),:),DecodeRsSh(PartitionsSh.training(s)));
                            predictionSh = predict(LDAModelSh,DecodeSpikes(PartitionsSh.test(s),:));
                            actualSh = DecodeRsSh(PartitionsSh.test(s));
                            correctSh = predictionSh - actualSh;
                            CVaccSh(s) = sum(correctSh==0) / length(correctSh);
                        end
                        AccShuff(q,1) = nanmean(CVaccSh);
                    end
                    UnitDec.True(NN,l) = nanmean(CVacc);
                    UnitDec.Shuff(NN,l) = nanmean(AccShuff);
                else 
                    UnitDec.True(NN,l) = NaN;
                    UnitDec.Shuff(NN,l) = NaN;
                end
            end

            fprintf('Neuron ID # %d\n',NN);
        end %neurons: FOR j= 1:size(RAW(i).Nrast,1)
    end %checking if the right session type
end %sessions: FOR i=1:length(RAW)

save('UnitDec_2R.mat','UnitDec');