%decodes trial identity from spiking across bins for each individual neuron
clear all
load('RAWvppsALL.mat')
load('RAWvpppsALL.mat')
load('RAWvppasALL.mat')
load('RAWvpppasALL.mat')

for PavInst=1:4; %DS is 1, Pav is 2
    if PavInst==1
        RAW=RAWvpppsALL;
    else if PavInst==2
            RAW=RAWvppsALL;
        else if PavInst==3
                RAW=RAWvpppasALL;
            else if PavInst==4
                    RAW=RAWvppasALL;
                end
            end
        end
    end
    folds = 5; %number of times cross-validated:
    Dura=[0 0.3]; %size of bin used for decoding:
    bins = 13; %number of bins
    binint = 0.3; %spacing of bins
    binstart = -.9; %start time of first bin relative to event
    shuffs = 1; %number of shuffled models created
    trials = 30; %max number of trials to use
    
    
    NN = 0;
    
    
    %events being compared
    
    
    UnitMdlAcc=NaN(100,bins);
    UnitMdlAccShuff=NaN(100,bins);
    
    for i=1:length(RAW) %loops through sessions
        if (PavInst==1 | PavInst==2)
            RD1=strcmp('CSPlus1', RAW(i).Einfo(:,2));
            RD2=strcmp('CSMinus1', RAW(i).Einfo(:,2));
        else
            RD1=strcmp('CSPlus', RAW(i).Einfo(:,2));
            RD2=strcmp('CSMinus', RAW(i).Einfo(:,2));
        end
        
        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);
                
                %find all spikes surrounding ev1
                [PSR1,N1]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{RD1},[(Dura(1)+(binstart - binint)+l*binint) (Dura(2)+(binstart - binint)+l*binint)],{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
                
                %only do at most the number of trials listed above
                if N1>trials
                    trials1=trials;
                else
                    trials1=N1;
                end
                
                %get the number of spikes on each trial
                for m=1:trials1
                    DecodeSpikes(m,1)=sum(PSR1{1,m}>(binstart));
                    DecodeRs(m,1)=1;
                end
               
                %find all spikes surrounding ev2
                [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{RD2},[(Dura(1)+(binstart - binint)+l*binint) (Dura(2)+(binstart - binint)+l*binint)],{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)

                %only do at most the number of trials listed above
                if N2>trials
                    trials2=trials;
                else
                    trials2=N2;
                end
                
                for n=1:trials2
                    DecodeSpikes(n+m,1)=sum(PSR2{1,n}>(binstart));
                    DecodeRs(n+m,1)=2;
                end
                
                if sum(DecodeSpikes(1:trials1,1))<7 || sum(DecodeSpikes((1+trials1):(trials1+trials2),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
                    UnitMdlAcc(NN,l) = nanmean(CVacc);
                    UnitMdlAccShuff(NN,l) = nanmean(AccShuff);
                else
                    UnitMdlAcc(NN,l) = NaN;
                    UnitMdlAccShuff(NN,l) = NaN;
                end
            end
            
            fprintf('Neuron ID # %d\n',NN);
        end %neurons: FOR j= 1:size(RAW(i).Nrast,1)
    end %sessions: FOR i=1:length(RAW)
    
    if PavInst==1
        UnitMdlAccDSSucrose=UnitMdlAcc;
        UnitMdlAccDSShuffSucrose=UnitMdlAccShuff;
        save('UnitMdlDSSucrose.mat','UnitMdlAccDSSucrose')
        save('UnitMdlDSShuffSucrose.mat','UnitMdlAccDSShuffSucrose')
    else if PavInst==2
            UnitMdlAccPavSucrose=UnitMdlAcc;
            UnitMdlAccPavShuffSucrose=UnitMdlAccShuff;
            save('UnitMdlPavSucrose.mat','UnitMdlAccPavSucrose')
            save('UnitMdlPavShuffSucrose.mat','UnitMdlAccPavShuffSucrose')
        else if PavInst==3
                UnitMdlAccDSAlcohol=UnitMdlAcc;
                UnitMdlAccDSShuffAlcohol=UnitMdlAccShuff;
                save('UnitMdlDSAlcohol.mat','UnitMdlAccDSAlcohol')
                save('UnitMdlDSShuffAlcohol.mat','UnitMdlAccDSShuffAlcohol')                
            else if PavInst==4
                    UnitMdlAccPavAlcohol=UnitMdlAcc;
                    UnitMdlAccPavShuffAlcohol=UnitMdlAccShuff;
                    save('UnitMdlPavAlcohol.mat','UnitMdlAccPavAlcohol')
                    save('UnitMdlPavShuffAlcohol.mat','UnitMdlAccPavShuffAlcohol')
                    
                end
            end
        end
    end
end