VentralPallidumAlcohol2019/AnalysisAndFigures/h_PooledDecodeSucroseAlcohol.m

%decodes trial identity from spiking across bins for each individual neuron

clear all
load('RAWpsvpALL.mat')
load('RAWvpppsALL.mat')
load('RAWvppaALL.mat')
load('RAWvpppaALL.mat')
if ~exist('RPPSall'),load('RPPSall.mat'); end
if ~exist('RPSall'),load('RPSall.mat'); end
if ~exist('RvppaALL'),load('RvppaALL.mat'); end
if ~exist('RvpppaALL'),load('RvpppaALL.mat'); end

for PavInst=1; %DS is 1, Pav is 2
    if PavInst==1
        RAW=RAWvpppsALL;
        selection=RPPSall.Structure==10 & RPPSall.CSPlusRatio>=.5 & (RPPSall.CSPlusRatio./(RPPSall.CSPlusRatio+RPPSall.CSMinusRatio))>=.7;
    else if PavInst==2
            RAW=RAWpsvpALL;
            selection=RPSall.Structure==10 & RPSall.CSPlusRatio>=.5 & (RPSall.CSPlusRatio./(RPSall.CSPlusRatio+RPSall.CSMinusRatio))>=.7;
        else if PavInst==3
                RAW=RAWvpppaALL;
                selection=RvpppaALL.Structure==10 & RvpppaALL.CSPlusRatio>=.5 & (RvpppaALL.CSPlusRatio./(RvpppaALL.CSPlusRatio+RvpppaALL.CSMinusRatio))>=.7;
            else if PavInst==4
                    RAW=RAWvppaALL;
                    selection=RvppaALL.Structure==10 & RvppaALL.CSPlusRatio>=.5 & (RvppaALL.CSPlusRatio./(RvppaALL.CSPlusRatio+RvppaALL.CSMinusRatio))>=.7;
                end
            end
        end
    end   
    
    
    
    
    
    TotalNeurons = 0;
    
    
    NumNeurons=[1 5 10 25 50 100]; %matrix of how many neurons to use on each iteration
    repetitions=20; %how many times to run the analysis
    folds = 5; %number of times cross-validated:
    Dura=[0 0.3]; %size of bin used for decoding:
    bins = 1; %number of bins
    binint = 0.3; %spacing of bins
    binstart = 0; %start time of first bin relative to event
    shuffs = 1; %number of shuffled models created
    
    
    
    for i=1:length(NumNeurons)
        PoolMdlAcc{i,1}=NaN(folds,bins);
        PoolMdlAccShuff{i,1}=NaN(folds,bins);
    end
    
    %total number of neurons in all sessions and events per session
    %for i=1:length(RAW)
    TotalNeurons=sum(selection);
    %     Ev1perSess(i)=length(RAW(i).Erast{Ev1,1});
    %     Ev2perSess(i)=length(RAW(i).Erast{Ev2,1});
    %end
    
    %figure out how many trials of each event can be used
    % Ev1s=min(Ev1perSess);
    % Ev2s=min(Ev2perSess);
    
    %Need to keep number of trials used in VP and NAc constant, so have to pick
    %minimum number across all sessions in both regions. Right now with
    %existing sessions it's 20 Ev1s and 20 Ev2s
    Ev1s=27;
    Ev2s=27;
    
    for v = 1:length(NumNeurons)
        for u=1:repetitions
            %pick which neurons to use
            SetupSel=cat(1,ones(NumNeurons(v),1),zeros(TotalNeurons-NumNeurons(v),1));
            Sel=(SetupSel(randperm(length(SetupSel)))==1);
            
            %setup the event identity matrix for decoding
            DecodeRs(1:Ev1s,1)=1;
            DecodeRs((Ev1s+1):(Ev1s+Ev2s),1)=2;
            
            for l=1:bins
                DecodeSpikes=NaN(1,1);
                AllSpikes=NaN(1,1);
                LowHz=zeros(1,NumNeurons(v));
                NN = 1;
                AllNN=1;
                for i=1:length(RAW) %loops through sessions
                    for j= 1:size(RAW(i).Nrast,1) %Number of neurons per session
                        if selection(AllNN)==1
                            Ev1Spikes=NaN(1,1);
                            Ev2Spikes=NaN(1,1);
                            
                            %events being compared
                            if (PavInst==1 | PavInst==2)
                                Ev1=strcmp('CSPlus1', RAW(i).Einfo(:,2));
                                Ev2=strcmp('CSMinus1', RAW(i).Einfo(:,2));
                            else
                                Ev1=strcmp('CSPlus', RAW(i).Einfo(:,2));
                                Ev2=strcmp('CSMinus', RAW(i).Einfo(:,2));
                            end
                            
                            
                            %get number of spikes for this neuron in this bin for all
                            %Ev1 trials
                            [PSR1,N1]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{Ev1},[(Dura(1)+(binstart - binint)+l*binint) (Dura(2)+(binstart - binint)+l*binint)],{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
                            for m=1:length(PSR1)
                                Ev1Spikes(m,1)=sum(PSR1{1,m}>(binstart));
                            end
                            
                            %pick which trials get used for decoding
                            SetupTrials=cat(1,ones(Ev1s,1),zeros(length(PSR1)-Ev1s,1));
                            Trials=(SetupTrials(randperm(length(SetupTrials)))==1);
                            
                            %put spikes from those trials in a matrix
                            AllSpikes(Trials,NN)=Ev1Spikes(Trials,1);
                            
                            %get all the spikes from reward 1p2 trials
                            [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{Ev2},[(Dura(1)+(binstart - binint)+l*binint) (Dura(2)+(binstart - binint)+l*binint)],{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
                            for n=1:length(PSR2)
                                Ev2Spikes(n,1)=sum(PSR2{1,n}>(binstart));
                            end
                            
                            %pick which trials get used for decoding
                            SetupTrials2=cat(1,ones(Ev2s,1),zeros(length(PSR2)-Ev2s,1));
                            Trials2=(SetupTrials2(randperm(length(SetupTrials2)))==1);
                            
                            %put spikes from those trials in a matrix
                            AllSpikes((Ev1s+1):(Ev1s+Ev2s),NN)=Ev2Spikes(Trials2,1);
                            NN=NN+1; %neuron counter
                            AllNN=AllNN+1;
                        else AllNN=AllNN+1;
                        end
                    end
                end
                
                
                
                %Setup spikes matrix for decoding
                DecodeSpikes=AllSpikes(:,Sel);
                
                %find neurons with too few spikes
                
                for t=1:NumNeurons(v)
                    if sum(DecodeSpikes(1:Ev1s,t))<3 || sum(DecodeSpikes((1+Ev1s):(Ev1s+Ev2s),t))<3
                        LowHz(1,t)=1;
                    end
                end
                
                %remove neurons with too few spikes
%                 if sum(LowHz)>0
%                     Sel2=LowHz<1;
%                     DecodeSpikes=DecodeSpikes(:,Sel2);
%                 end
                
                %set up validation result matrices
                CVacc = NaN(folds,1);
                CVaccSh = NaN(folds,1);
                
                %normal model
                Partitions = cvpartition(DecodeRs,'KFold',folds);
                for r = 1:folds
                    
                    SVMModel = fitcdiscr(DecodeSpikes(Partitions.training(r),:),DecodeRs(Partitions.training(r)),'discrimType','diaglinear');
                    [prediction,score] = predict(SVMModel,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
                        SVMModelSh = fitcdiscr(DecodeSpikes(PartitionsSh.training(s),:),DecodeRsSh(PartitionsSh.training(s)),'discrimType','diaglinear');
                        [predictionSh, score] = predict(SVMModelSh,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
                
                PoolMdlAcc{v,1}(u,l) = nanmean(CVacc);
                PoolMdlAccShuff{v,1}(u,l) = nanmean(AccShuff);
                
                fprintf(['Condition #' num2str(v) ', Rep #' num2str(u) ', Bin #' num2str(l) '\n']);
            end %bins
        end
    end
    
    
    
    if PavInst==1
        PoolMdlAccDSSucrose=PoolMdlAcc;
        PoolMdlAccDSSucroseShuff=PoolMdlAccShuff;
        save('PoolMdlAccDSSucrose.mat','PoolMdlAccDSSucrose')
        save('PoolMdlAccDSSucroseShuff.mat','PoolMdlAccDSSucroseShuff')
    else if PavInst==2
            PoolMdlAccPavSucrose=PoolMdlAcc;
            PoolMdlAccPavSucroseShuff=PoolMdlAccShuff;
            save('PoolMdlAccPavSucrose.mat','PoolMdlAccPavSucrose')
            save('PoolMdlAccPavSucroseShuff.mat','PoolMdlAccPavSucroseShuff')
        else if PavInst==3
                PoolMdlAccDSAlcohol=PoolMdlAcc;
                PoolMdlAccDSAlcoholShuff=PoolMdlAccShuff;
                save('PoolMdlAccDSAlcohol.mat','PoolMdlAccDSAlcohol')
                save('PoolMdlAccDSAlcoholShuff.mat','PoolMdlAccDSAlcoholShuff')
            else if PavInst==4
                    PoolMdlAccPavAlcohol=PoolMdlAcc;
                    PoolMdlAccPavAlcoholShuff=PoolMdlAccShuff;
                    save('PoolMdlAccPavAlcohol.mat','PoolMdlAccPavAlcohol')
                    save('PoolMdlAccPavAlcoholShuff.mat','PoolMdlAccPavAlcoholShuff')
                    
                end
            end
        end
    end
    
end