%% generate figure 6 and figure 6 - supplement figure 1
% fourrier analysis - hierarchical clustering - Random forest approach on
% early training dataset
% statistics for the early training dataset


clc;clear;close all;
load ('Rearlytraining_light.mat');
load ('Celltype_earlyTraining.mat');
load('TRN_DT5_earlytrain.mat')

Xevent=[-0.25 0.25];
Yaxis=[-5 10];Xaxis=[1 357];
Ishow=find(Tm>=Xevent(1) & Tm<=Xevent(2));
time=[1:1:357];
BrainRegion=[10 20];
% ----------- COLORS ---------
myblue=[0 113/255 188/255];
mypurple=[200/255 0 255/255];
mydarkblue=[0 0/255 255/255];
TickSize=[0.015 0.02];
Clim=[-5 5];
Eventnames={'';'LI';'';'LP1';'';'LP2';'';'LP3';'';'LP4';'';'LP5';'';'PE'};

DSconcat_ACQ_allsession=[];Coord_allsession=[];Ninfo_allsession=[];DSmeanz_ACQ_allsession=[];
TRN_allsession=[];EXCINH_allsession=[];sessionID_allsession=[];
for i=1:length(Ses)
    
    %normalization based on max
    for k=1:length(Erefnames)
        Ev(k).PSTH_norm1=normalize(Ses(i).Ev(k).PSTHraw,Ses(i).Ev(k).PSTHrawBL,1);
    end
    
    DSconcat_ACQ=[];DSmeanz_ACQ=[];
    %concat with normalized data
    DSconcat_ACQ = Ev(7).PSTH_norm1(:,Ishow);
    DSmeanz_ACQ = Ses(i).Ev(7).Meanz(:,1);
    for k=1:5
        DSconcat_ACQ = cat(2,DSconcat_ACQ,Ev(k).PSTH_norm1(:,Ishow));
        DSmeanz_ACQ = cat(2,DSmeanz_ACQ,Ses(i).Ev(k).MeanzPRE(:,1));
        DSmeanz_ACQ = cat(2,DSmeanz_ACQ,Ses(i).Ev(k).Meanz(:,1));
    end
    DSconcat_ACQ = cat(2,DSconcat_ACQ,Ev(6).PSTH_norm1(:,Ishow));
    DSmeanz_ACQ = cat(2,DSmeanz_ACQ,Ses(i).Ev(6).MeanzPRE(:,1));
    
    DSconcat_ACQ_allsession=cat(1,DSconcat_ACQ_allsession,DSconcat_ACQ);
    Coord_allsession=cat(1,Coord_allsession,Ses(i).Coord(:,4));
    Ninfo_allsession=cat(1,Ninfo_allsession,Ses(i).Ninfo);
    DSmeanz_ACQ_allsession=cat(1,DSmeanz_ACQ_allsession,DSmeanz_ACQ);
    TRN_allsession=cat(1,TRN_allsession,sesTRN(i).TRN(:,1));
    sessionID_allsession=cat(1,sessionID_allsession,repmat(i,length(Ses(i).TRN(:,1)),1));
    if i>3
        EXCINH_allsession=cat(1,EXCINH_allsession,sesTRN(i).TRN(:,2));
    end
end

s2={'YVG04','YVG05','YVG08','YVJ11','YVJ15'};
for j=1:5
    for i=1:length(Ninfo_allsession)
        RatID(i,j)=strncmp(Ninfo_allsession{i,1},s2(j),5);
    end
end
RatID_ACQ=sum(RatID,2);
% concat with Zscore
DSconcatZ_ACQ_allsession=[];DSevent_ACQ_allsession=[];
for i=1:length(Ses)
    DSconcat_ACQ=[];DSevent_ACQ=[];
    
    DSconcat_ACQ=Ses(i).Ev(7).PSTHz(:,Ishow);
    DSevent_ACQ=Ses(i).Ev(7).Meanz(:,1);
    for k=1:6 %loop through event
        DSconcat_ACQ=cat(2,DSconcat_ACQ,Ses(i).Ev(k).PSTHz(:,Ishow));
        DSevent_ACQ = cat(2,DSevent_ACQ,Ses(i).Ev(k).Meanz(:,1));
    end
    DSconcatZ_ACQ_allsession=cat(1,DSconcatZ_ACQ_allsession,DSconcat_ACQ);
    DSevent_ACQ_allsession=cat(1,DSevent_ACQ_allsession,DSevent_ACQ);
end

selection=Coord_allsession~=0 & Celltype(:,2)==1 & ~isnan(mean(DSconcat_ACQ_allsession,2));
DSconcat_ACQshort=DSconcat_ACQ_allsession(selection,:);

DSevent_ACQshort=DSevent_ACQ_allsession(selection,:);
Coord_short=Coord_allsession(selection,:);
Ninfo_short=Ninfo_allsession(selection,:);

%% Clustering on phasic and non phasic neurons with Fourier analysis 

close all;
f_low=1;f_high=4;

[X_all,f_all,P_all,Nt_all]=myfft(DSconcat_ACQshort',0.01,1);
total_power=sum(P_all);
P_all_norm = bsxfun(@rdivide,P_all, total_power);
%figure;plot(f_all,P_all_norm);set(gca,'ylim',[0 0.6]);xlabel('Hz');ylabel('Normalized power');
%title('Normalized FFT of all neurons');

a=find(f_all>f_low);b=find(f_all<f_high);
interval_to_sum=[a(1) b(end)];
lower_window_sum_all=sum(P_all_norm(1:a(1),:));
upper_window_sum_all=sum(P_all_norm(a(1):b(end),:));

ratio_all=lower_window_sum_all./upper_window_sum_all;
close all

X_ratio=[[lower_window_sum_all]' [upper_window_sum_all]']; %#ok<NBRAK>
eva_Ch_ratio = evalclusters(X_ratio,'linkage','CalinskiHarabasz','KList',[1:10]);
phasicID=eva_Ch_ratio.OptimalY;

Z = linkage(X_ratio,'ward'); % X is the data set, neuron in rows, features in columns
c = cluster(Z,'Maxclust',3); % determine the number cluster here)
D = pdist(X_ratio);
leafOrder = optimalleaforder(Z,D);
cutoff = median([Z(end-2,3) Z(end-1,3)]);
% 
% figure
% subplot(9,4,[29 33])
% dendrogram(Z,'ColorThreshold',cutoff, 'Reorder',leafOrder);
%% Figure Fourier transform

c_ratio=eva_Ch_ratio.OptimalY;

for i=1:length(c_ratio)
    if c_ratio(i,1)==1
        color(i,:)=[255/255 0/255 0/255]; %phasic red
    elseif c_ratio(i,1)==2
        color(i,:)=[0 255/255 0/255]; %sustained green
    end
end

figure
subplot(5,3,1)
scatter(lower_window_sum_all,upper_window_sum_all,[],color)
hold on
title(strcat('hierarchical clustering'))
xlabel('power in low freq') % x-axis label
ylabel('power in int freq')% y-axis label
legend('phasic','non-phasic')
hold off

%% visualization classes

for i=1:2 
    ind=find(phasicID==i);
    t=num2str(i);
    Clim=[-5 5];
    DSconcat_ACQ_class=[];SORTimg=[];MaxValDLS=[];MeanValDLS=[]; MaxTimeDLS=[];TMP=[];NNid=[];
    DSconcat_ACQ_class=DSconcat_ACQshort(ind,:);
    
    for NN=1:size(DSconcat_ACQ_class,1)
        MeanValDLS(NN,1)=mean(DSconcat_ACQ_class(NN,52:306),2);
        [MaxValDLS(NN,1),MaxInd]=max(DSconcat_ACQ_class(NN,:));
        MaxTimeDLS(NN,1)=time(time(1)+MaxInd-1);
        if i==2
            MaxTimePhasic=MaxTimeDLS;
        end
    end
    if i==1
        TMP=MeanValDLS;
        subplot(5,3,[4 7 10 13])
    elseif i==2 
        TMP=MeanValDLS;
        subplot(5,3,[5 8 11 14])
    end
    TMP(isnan(TMP))=0; %To place the neurons with no onset/duration/peak at the top of the color-coded map
    [~,SORTimg]=sort(TMP);
    %NNid=cell2mat(Ninfo(phasic_neuron_selection,3));NNid2=NNid(SORTimg);
    DSconcat_ACQ_class=DSconcat_ACQ_class(SORTimg,:);
    
    
    imagesc(time,[1 size(DSconcat_ACQ_class,1)],DSconcat_ACQ_class);%colorbar;axis tight;
    colormap(jet);
    hold on
    plot([51 51],[size(DSconcat_ACQ_class,1)+0.5 0.5],'k')
    plot([102 102],[size(DSconcat_ACQ_class,1)+0.5 0.5],'k')
    plot([153 153],[size(DSconcat_ACQ_class,1)+0.5 0.5],'k')
    plot([204 204],[size(DSconcat_ACQ_class,1)+0.5 0.5],'k')
    plot([255 255],[size(DSconcat_ACQ_class,1)+0.5 0.5],'k')
    plot([306 306],[size(DSconcat_ACQ_class,1)+0.5 0.5],'k')
    set(gca,'XTick',[0:25.5:357]);
    set(gca,'xticklabel',Eventnames)
    hold off
    if i==1
        ylabel('Non-Phasic neurons')
    else
        ylabel('Phasic neurons')
    end
end

ClassACQ(1:length(Celltype),1)=0;
ClassACQ(selection,1)=phasicID;


%% Run random forest for phasic neurons

load('ClassHC.mat' );
XTrainingDATA=ClassPHAS(:,2:4);
YTrainingDATA=ClassPHAS(:,1);

LIwin=[26:51];%post lever insertion window
LPwin=[52:280];%pre-1stLP to pre-lastLP window
PEwin=[307:332];% pre-PE window

binLI=mean(DSconcat_ACQ_allsession(:,LIwin),2);
binLP=mean(DSconcat_ACQ_allsession(:,LPwin),2);
binPE=mean(DSconcat_ACQ_allsession(:,PEwin),2);
binnedDS=cat(2,binLI,binLP,binPE);   
XTestingDATA=binnedDS(ClassACQ(:,1)==2,:);

model = classRF_train(XTrainingDATA, YTrainingDATA); % train Random Forest model

% prepare the data for FOR loop
X = {XTrainingDATA, XTestingDATA};
Y = {YTrainingDATA};
classifier = {'Raw Data','Random Forest'};
mc = ['r','g']; msym = ['o','o']; msiz = [3,15];


for i = 1 : 2 % iterate through training and test data
    x = X{i}; y = Y;
    for j = 1 : 2 % iterate through different classifiers
        
        switch(j) % y_hat = running models ? model predictions : y
            case 1 % raw data
                y_hat = y;
            case 2 % Random Forest
                y_hat = classRF_predict(x ,model); % Random Forest predictions
                if i==2
                    Class_Acq(:,1)=y_hat;
                end
        end
    end
end

ClassACQ(1:length(Celltype),2)=0;
ClassACQ(ClassACQ(:,1)==2,2)=Class_Acq;

%% Random forest for non phasic neurons

XTrainingDATA_NP=ClassNONPHASIC(:,2:4);
YTrainingDATA_NP=ClassNONPHASIC(:,1);
XTestingDATA_NP=binnedDS(ClassACQ(:,1)==1,:);

model = classRF_train(XTrainingDATA_NP, YTrainingDATA_NP); % train Random Forest model

% prepare the data for FOR loop
X = {XTrainingDATA_NP, XTestingDATA_NP};
Y = {YTrainingDATA_NP};
classifier = {'Raw Data','Random Forest'};
mc = ['r','g']; msym = ['o','o']; msiz = [3,15];


for i = 1 : 2 % iterate through training and test data
    x = X{i}; y = Y;
    for j = 1 : 2 % iterate through different classifiers
        
        switch(j) % y_hat = running models ? model predictions : y
            case 1 % raw data
                y_hat = y;
            case 2 % Random Forest
                y_hat = classRF_predict(x ,model); % Random Forest predictions
                if i==2
                    Class_Acq_NP(:,1)=y_hat;
                end
        end
    end
end

ClassACQ(ClassACQ(:,1)==1,2)=Class_Acq_NP+3;

%% selection of neurons per session per class
% 1-3 is for phasic and 4-5 for non phasic

for i=1:max(Celltype(:,1))
    InDSsel=Ses(i).Coord(:,4)~=0;
    Ses(i).Class(InDSsel,1)=ClassACQ(Celltype(:,1)==i,2);
    Ses(i).Celltype(InDSsel,1)=Celltype(Celltype(:,1)==i,2);
end

for k=1:max(ClassACQ(:,2))
    DLSconcat_ACQ_3ds(k).DLS=[];
    DMSconcat_ACQ_3ds(k).DMS=[];
    for i=4:13
        DLSselection=Ses(i).Coord(:,4)==10 & Ses(i).Class(:,1)==k & Ses(i).Celltype(:,1)==1;
        DMSselection=Ses(i).Coord(:,4)==20 & Ses(i).Class(:,1)==k & Ses(i).Celltype(:,1)==1;
        
        DLSconcat_ACQ=[];DMSconcat_ACQ=[];
        DLSconcat_ACQ = cat(2,DLSconcat_ACQ,Ses(i).Ev(7).PSTHz(DLSselection,Ishow));
        DMSconcat_ACQ = cat(2,DMSconcat_ACQ,Ses(i).Ev(7).PSTHz(DMSselection,Ishow));
        for j=1:6
            DLSconcat_ACQ = cat(2,DLSconcat_ACQ,Ses(i).Ev(j).PSTHz(DLSselection,Ishow));
            DMSconcat_ACQ = cat(2,DMSconcat_ACQ,Ses(i).Ev(j).PSTHz(DMSselection,Ishow));
        end
        MeanDLS(k).img(i-3,:)=nanmean(DLSconcat_ACQ,1);
        MeanDMS(k).img(i-3,:)=nanmean(DMSconcat_ACQ,1);
        if i>10
            DLSconcat_ACQ_3ds(k).DLS=cat(1,DLSconcat_ACQ_3ds(k).DLS,DLSconcat_ACQ);
            DMSconcat_ACQ_3ds(k).DMS=cat(1,DMSconcat_ACQ_3ds(k).DMS,DMSconcat_ACQ);
            
            DLSsem(k,:)=nanste(DLSconcat_ACQ(:,time),1);
            DMSsem(k,:)=nanste(DMSconcat_ACQ(:,time),1);
        end
    end
end

%% Proportion of neurons in phasic vs non pHasic classes

for j=1:2
    for i=4:13
        DLSselectionP=Ses(i).Coord(:,4)==BrainRegion(1) & (Ses(i).Class(:,1)<4) & Ses(i).Celltype(:,1)==1 ;
        DMSselectionP=Ses(i).Coord(:,4)==BrainRegion(2) & (Ses(i).Class(:,1)<4) & Ses(i).Celltype(:,1)==1 ;

        if j==1;DSselection=logical(DLSselectionP); else DSselection=logical(DMSselectionP);end
        
        DS_Nb_tot_phasic(i-3,j)=sum(DSselection);
        PercentPhasic(i-3,1+(j-1)*2)=100*sum(DSselection)/sum(Ses(i).Coord(:,4)==BrainRegion(j) & Ses(i).Celltype(:,1)==1);
        PercentPhasic(i-3,2+(j-1)*2)=100-PercentPhasic(i-3,1+(j-1)*2);
    end   
end

subplot(5,3,2) % plot DLS phasic vs non phasic
bp1=barh(PercentPhasic(:,1:2),1,'stacked');
axis([0 100 0.5 10.5])
xlabel('% DLS neurons')
ylabel('sessions')
subplot(5,3,3) % plot DMS phasic vs non phasic
bp2=barh(PercentPhasic(:,3:4),1,'stacked');
axis([0 100 0.5 10.5])
xlabel('% DMS neurons')
ylabel('sessions')

%% heat maps classes of neurons
Eventnames={'';'LI';'';'LP1';'';'LP2';'';'LP3';'';'LP4';'';'LP5';'';'PE'};
titlegraphDLS={'DLS - Start','DLS - Stop', 'DLS - middle','DLS - INH','DLS - EXC'};
titlegraphDMS={'DMS - Start','DMS - Stop', 'DMS - middle','DMS - INH','DMS - EXC'};
%titlePSTH={'Inhibited neurons', 'Excited neurons'};

BrainRegion=[10 20];
figure
nplot=0;

for k=1:max(ClassACQ(:,2)) % loop thru non phasic classes 
    if k<4 || k==5
        nplot=nplot+1;
    elseif k==4
        nplot=nplot+10;
    end
    subplot(6,4,2+(nplot-1))
    imagesc(time,[1 size(MeanDLS(k).img,1)],MeanDLS(k).img,Clim); %colorbar;axis tight;
    colormap(jet);
    %colormap(plasma);
    set(gca,'XTick',[0:25.5:357]);
    set(gca,'xticklabel',Eventnames)
    title(titlegraphDLS(k));
    ylabel('sessions');
    hold on
    plot([51 51],[size(MeanDLS(k).img,1)+0.5 0.5],'k')
    plot([102 102],[size(MeanDLS(k).img,1)+0.5 0.5],'k')
    plot([153 153],[size(MeanDLS(k).img,1)+0.5 0.5],'k')
    plot([204 204],[size(MeanDLS(k).img,1)+0.5 0.5],'k')
    plot([255 255],[size(MeanDLS(k).img,1)+0.5 0.5],'k')
    plot([306 306],[size(MeanDLS(k).img,1)+0.5 0.5],'k')
    hold off
    
    
    subplot(6,4,6+(nplot-1))
    imagesc(time,[1 size(MeanDMS(k).img,1)],MeanDMS(k).img,Clim); %colorbar;axis tight;
    colormap(jet);
    set(gca,'XTick',[0:25.5:357]);
    set(gca,'xticklabel',Eventnames)
    title(titlegraphDMS(k));
    ylabel('sessions');
    hold on
    plot([51 51],[size(MeanDMS(k).img,1)+0.5 0.5],'k')
    plot([102 102],[size(MeanDMS(k).img,1)+0.5 0.5],'k')
    plot([153 153],[size(MeanDMS(k).img,1)+0.5 0.5],'k')
    plot([204 204],[size(MeanDMS(k).img,1)+0.5 0.5],'k')
    plot([255 255],[size(MeanDMS(k).img,1)+0.5 0.5],'k')
    plot([306 306],[size(MeanDMS(k).img,1)+0.5 0.5],'k')
    hold off
    
    subplot(6,4,10+(nplot-1))
    
    DLSpsth=nanmean(DLSconcat_ACQ_3ds(k).DLS,1);
    DLS_sem=nanste(DLSconcat_ACQ_3ds(k).DLS,1);
    DLSup=DLSpsth+DLS_sem;
    DLSdown=DLSpsth-DLS_sem;
    
    DMSpsth=nanmean(DMSconcat_ACQ_3ds(k).DMS,1);
    DMS_sem=nanste(DMSconcat_ACQ_3ds(k).DMS,1);
    DMSup=DMSpsth+DMS_sem;
    DMSdown=DMSpsth-DMS_sem;
    
    hold on;
    patch([time,time(end:-1:1)],[DLSup,DLSdown(end:-1:1)],myblue,'EdgeColor','none');alpha(0.5);
    plot(time,DLSpsth,'Color',myblue,'linewidth',1.5);
    patch([time,time(end:-1:1)],[DMSup,DMSdown(end:-1:1)],'r','EdgeColor','none');alpha(0.5);
    plot(time,DMSpsth,'r','linewidth',1.5);
    
    yaxis=[round(min(DMSdown)-0.5,0) round(max(DMSup)+0.5,0)];
    plot([0 357], [0 0],'LineStyle',':','Color','k');
    plot([51 51],yaxis,'LineStyle',':','Color','k')
    plot([102 102],yaxis,'LineStyle',':','Color','k')
    plot([153 153],yaxis,'LineStyle',':','Color','k')
    plot([204 204],yaxis,'LineStyle',':','Color','k')
    plot([255 255],yaxis,'LineStyle',':','Color','k')
    plot([306 306],yaxis,'LineStyle',':','Color','k')
    axis([Xaxis,yaxis]);
    set(gca,'XTick',0:25.5:357);
    set(gca,'xticklabel',Eventnames)
    ylabel('z-score');
    hold off
end


%% Proportion of neurons
for k=1:3
    for i=4:13   
        DLSselection=Ses(i).Coord(:,4)==BrainRegion(1) & (Ses(i).Class(:,1)==k) & Ses(i).Celltype(:,1)==1;
        DMSselection=Ses(i).Coord(:,4)==BrainRegion(2) & (Ses(i).Class(:,1)==k) & Ses(i).Celltype(:,1)==1;
        
        DLS_prop_Phasic(i-3,k)=100*sum(DLSselection)/(sum(Ses(i).Coord(:,4)==BrainRegion(1) & Ses(i).Celltype(:,1)==1 & Ses(i).Class(:,1)<4 & Ses(i).Class(:,1)~=0));
        DMS_prop_Phasic(i-3,k)=100*sum(DMSselection)/(sum(Ses(i).Coord(:,4)==BrainRegion(2) & Ses(i).Celltype(:,1)==1 & Ses(i).Class(:,1)<4 & Ses(i).Class(:,1)~=0));
        
        DLSselectionNumber(i-3,k)=sum(DLSselection);
        DMSselectionNumber(i-3,k)=sum(DMSselection);
    end
end

subplot(6,4,1)
b1=barh(DLS_prop_Phasic,1,'stacked');
axis([0 100 0.5 10.5])
xlabel('% DLS neurons')
ylabel('sessions')

subplot(6,4,5)
b2=barh(DMS_prop_Phasic,1,'stacked');
axis([0 100 0.5 10.5])
xlabel('% DMS neurons')
ylabel('sessions')


for k=1:2
    for i=4:13   
        DLSselection=Ses(i).Coord(:,4)==BrainRegion(1) & (Ses(i).Class(:,1)==k+3) & Ses(i).Celltype(:,1)==1;
        DMSselection=Ses(i).Coord(:,4)==BrainRegion(2) & (Ses(i).Class(:,1)==k+3) & Ses(i).Celltype(:,1)==1;
        
        DLS_prop_nonPhasic(i-3,k)=100*sum(DLSselection)/(sum(Ses(i).Coord(:,4)==BrainRegion(1) & Ses(i).Celltype(:,1)==1 & Ses(i).Class(:,1)>3 & Ses(i).Class(:,1)~=0));
        DMS_prop_nonPhasic(i-3,k)=100*sum(DMSselection)/(sum(Ses(i).Coord(:,4)==BrainRegion(2) & Ses(i).Celltype(:,1)==1 & Ses(i).Class(:,1)>3 & Ses(i).Class(:,1)~=0));
        
        DLSselectionNumber(i-3,k)=sum(DLSselection);
        DMSselectionNumber(i-3,k)=sum(DMSselection);
    end
end

subplot(6,4,13)
b1=barh(DLS_prop_nonPhasic,1,'stacked');
axis([0 100 0.5 10.5])
xlabel('% DLS neurons')
ylabel('sessions')

subplot(6,4,17)
b2=barh(DMS_prop_nonPhasic,1,'stacked');
axis([0 100 0.5 10.5])
xlabel('% DMS neurons')
ylabel('sessions')

%% %% make table for statistics

selectDT5Session=sessionID_allsession(:,1)>3;

table_ACQ=cat(2,sessionID_allsession(selectDT5Session,1),Coord_allsession(selectDT5Session,1),Celltype(selectDT5Session,2),ClassACQ(selectDT5Session,1),ClassACQ(selectDT5Session,2),DSmeanz_ACQ_allsession(selectDT5Session,:),TRN_allsession(selectDT5Session,1),EXCINH_allsession(:,1), RatID_ACQ(selectDT5Session,1));
save('table_ACQ.mat','table_ACQ')


nbevent=[1:12];

selectionSTAT=table_ACQ(:,2)~=0 & table_ACQ(:,3)==1 & table_ACQ(:,18)~=0;
sessionCat=categorical(table_ACQ(selectionSTAT,1));
tableSTAT=array2table(table_ACQ(selectionSTAT,1:20),'VariableNames',{'session','region','celltype','phasicID','class','win1','win2','win3','win4','win5','win6','win7','win8','win9','win10','win11','win12','TRN','EXCINH','habit'});
tableSTAT.session=categorical(tableSTAT.session);
tableSTAT.region=categorical(tableSTAT.region);
tableSTAT.habit=categorical(tableSTAT.habit);
 
 rm = fitrm(tableSTAT,'win1-win12~region*session','WithinDesign',nbevent');
 ranovatbl_event = ranova(rm);   
 Btw_ranovatbl_event = anova(rm);
 SphericityAssumption=mauchly(rm);
 posthoc = multcompare(rm,'session', 'by', 'region');
 
 
 rmh = fitrm(tableSTAT,'win1-win12~region*session*habit','WithinDesign',nbevent');
 ranovatbl_event_habit = ranova(rmh);   
 Btw_ranovatbl_event_habit = anova(rmh);
 SphericityAssumption_habit=mauchly(rmh);
 posthoc_habit = multcompare(rmh,'region', 'by', 'habit');

 
%  for i=1:max(table_ACQ(:,5))
for i=1:3 %loop thru phasic neurons, problem, not enough neurons in non phasic classes
    selectionSTAT1=table_ACQ(:,2)~=0 & table_ACQ(:,3)==1 & table_ACQ(:,5)==i; 
    
    tableSTAT1=array2table(table_ACQ(selectionSTAT1,1:20),'VariableNames',{'session','region','celltype','phasicID','class','win1','win2','win3','win4','win5','win6','win7','win8','win9','win10','win11','win12','TRN','EXCINH','habit'});
    tableSTAT1.session=categorical(tableSTAT1.session);
    tableSTAT1.region=categorical(tableSTAT1.region);
    tableSTAT1.habit=categorical(tableSTAT1.habit);
    
    rm = fitrm(tableSTAT1,'win1-win12~region*session','WithinDesign',nbevent');
    stat(i).ranovatbl_event = ranova(rm);
    stat(i).Btw_ranovatbl_event = anova(rm);
    
    [~,stat(i).LIevent,~]=anovan(table_ACQ(selectionSTAT1,6),{table_ACQ(selectionSTAT1,2)},'varnames',{'region'}, 'display','off');
    [~,stat(i).LRevent,~]=anovan(table_ACQ(selectionSTAT1,16),{table_ACQ(selectionSTAT1,2)},'varnames',{'region'}, 'display','off');
    [~,stat(i).PEevent,~]=anovan(table_ACQ(selectionSTAT1,17),{table_ACQ(selectionSTAT1,2)},'varnames',{'region'}, 'display','off');

    rmh = fitrm(tableSTAT1,'win1-win12~region*session*habit','WithinDesign',nbevent');
    stat(i).ranovatbl_event_habit = ranova(rmh);
    stat(i).Btw_ranovatbl_habit = anova(rmh); 
    stat(i).posthoc_habit = multcompare(rmh,'region', 'by', 'habit');
end


 for i=1:10
     selectionSTAT2=table_ACQ(:,2)~=0 & table_ACQ(:,3)==1 & table_ACQ(:,1)==i+3 & table_ACQ(:,4)==1; 
    [chi(i).tbl1,chi(i).chi2,chi(i).p1,chi(i).labels] = crosstab(table_ACQ(selectionSTAT2,2),table_ACQ(selectionSTAT2,5));
 end
 
 save('Rearlytraining_light.mat','Erefnames','Ses','Tm')