set(0, 'DefaultLineLineWidth', 2);
set(0,'defaultTextFontSize', 12);
set(0,'defaultAxesFontSize',12);
set(0,'DefaultFigureColormap',jet)
set(0,'defaultfigurecolor',[1 1 1])
set(0,'DefaultAxesTitleFontWeight', 'normal')

clc
clear all

%% Only using Figure 8 maze ACh data which is AWAKE state.  Ignore first  200 s. 
clear all
animalNames = {
    'LIMP003_008';
    'LAMP37_011';
    'LAMP3_010';
    };

maxInds = [1e5,1e5,1000];
zAChdFF_all = [];
thetaACh_all = [];
rAchtheta_all = [];
pAchtheta_all = [];
for anID = 1:length(animalNames)
    Names_ACh = dir(['ACh_theta/',animalNames{anID},'*ACh_theta.mat']);
    filenamesACh = {Names_ACh.name};
    foldernameACh = {Names_ACh.folder};
    load(fullfile(foldernameACh{1},filenamesACh{1}));
    topti = lfp_fbtime;
    fs = THLFPs.samplingRate;
    lfp = double(THLFPs.data(ceil(topti(1)/2e4*fs):end));
    lfp(abs(lfp)>5000) = nan;
    lfp = fillmissing(lfp,'linear');
    [spectrogrammt,tm,fm] = MTSpectrogram(double(lfp),'window',2.0);
    bands = SpectrogramBands(spectrogrammt,fm);

    fs_dff  = 100;
    dt = mean(diff(tm));
    Ndown = round(fs_dff*dt);
    dFF_ACh_plot = downsample_yz(dFoverF,Ndown);

    N11 = min(length(dFF_ACh_plot),length(bands.theta));
    N1 = min(N11, maxInds(anID))
    zAChdFF = (zscore(dFF_ACh_plot(300:N1)))';
    ztheta = zscore(bands.theta(300:N1));
    while  ~isempty(find(ztheta>3, 1))
        idx_turn = find(ztheta>3);
        ztheta(idx_turn(1)-2:idx_turn(1)+2) = [];
        zAChdFF(idx_turn(1)-2:idx_turn(1)+2) = [];
    end
    zAChdFF_all = [zAChdFF_all; zAChdFF];
    thetaACh_all = [thetaACh_all;ztheta];
    for i = 1:9
        idx = find(ztheta>(i-4)*0.5 & ztheta<(i-3)*0.5);
        length(idx);
        AChmean(i,anID) = mean(zAChdFF(idx),'omitnan');
    end
    figure
    set(gcf,'Position',[100 100 400 400])
    scatter(zAChdFF,ztheta,20,[1,0,0],'filled')
    xlim([-2,6])
    ylim([-3,3])
    box off
    [r p] = corrcoef(zAChdFF,ztheta);
    rAchtheta_all(anID) = r(1,2);
    pAchtheta_all(anID) = p(1,2)
end
%%
fig1 = figure();
options.handle     =  fig1;
options.color_area = [0,0,0];    
options.color_line = [0,0,0];
options.alpha      = 0.4;
options.line_width = 3;
options.error      = 'sem';
options.x_axis = -1.25:0.5:2.75;
set(gcf,'Position',[100 100 400 700])
subplot(2,1,1)
Ndown =1;
hold on
scatter(downsample(thetaACh_all,Ndown),downsample(zAChdFF_all,Ndown),10, [1,0,0],'filled')
plot_areaerrorbar(AChmean',options)
xlabel('Theta power (zscored)')
ylabel('ACh')
xlim([-1.5,3])
ylim([-3,3])
box off
disp(['Mean crosscorrelation is :    ', num2str(mean(rAchtheta_all))])
disp(['Std of crosscorrelation is :    ', num2str(std(rAchtheta_all))])

corrcoef(zAChdFF_all,thetaACh_all)
%% Theta power is better than theta score for OXT relationship
%%
clear all
animalNames = {
    'OXYT014_001';
    'OXYT010_004';
    'OXYT012_001';
    };
statenames = {'WAKE','','NREM','','REM'};
zOxtdFF_all = [];
thetaOxt_all = [];
rOxttheta_all = [];
pOxttheta_all = [];
for anID = 1:length(animalNames)
    Names_OXT = dir(['OXT_theta/',animalNames{anID},'*OXT_theta.mat']);
    filenamesOXT = {Names_OXT.name};
    foldernameOXT = {Names_OXT.folder};
    load(fullfile(foldernameOXT{1},filenamesOXT{1}));
    topti = lfp_fbtime;
    fs = THLFPs.samplingRate;
    lfp = double(THLFPs.data(ceil(topti(1)/2e4*fs):floor(topti(2)*fs/2e4)));
    lfp(abs(lfp)>5000) = nan;
    lfp = fillmissing(lfp,'linear');
    [spectrogrammt,tm,fm] = MTSpectrogram(double(lfp),'window',2.0);
    bands = SpectrogramBands(spectrogrammt,fm);

    fs_dff  = 100;
    dt = mean(diff(tm));
    Ndown = round(fs_dff*dt);
    dFF_oxt_plot = downsample_yz(dFoverF,Ndown);

    sState = bz_INTtoIDX(SleepState.ints, 'statenames', statenames);
    N1 = min(length(dFF_oxt_plot),length(bands.theta));
    idx_wake = find(sState.states(round(topti(1)/2e4):end) ==1);
    idx_wake(idx_wake>N1)= [];
    
    zOxtdFF = (zscore(dFF_oxt_plot(idx_wake)))';
    ztheta = zscore(bands.theta(idx_wake));
    zOxtdFF_all = [zOxtdFF_all; zOxtdFF];
    thetaOxt_all = [thetaOxt_all;ztheta];
    for i = 1:12
        idx = find(ztheta>(i-6)*0.5 & ztheta<(i-5)*0.5);
        length(idx)
        Oxtmean(i,anID) = mean(zOxtdFF(idx),'omitnan');
    end
%     figure
%     set(gcf,'Position',[100 100 400 400])
%     scatter(zscore(dFF_oxt_plot(idx_wake)),zscore(bands.theta(idx_wake)),20,[0,0.6,0],'filled')
%     xlim([-2,6])
%     %ylim([-3,3])
%     box off
     [r,p] = corrcoef(zOxtdFF,ztheta);
     rOxttheta_all(anID) = r(1,2);
     pOxttheta_all(anID) = p(1,2);
end
%%
options.handle     =  fig1;
options.color_area = [0 0 0]./255;    
options.color_line = [0 0 0]./255;
options.alpha      = 0.4;
options.line_width = 3;
options.error      = 'sem';
options.x_axis = -2.25:0.5:3.25;
set(gcf,'Position',[100 100 400 700])
subplot(2,1,2)
Ndown =2;
options.color_area = [0 0 0]./255;    
options.color_line = [0 0 0]./255;
scatter(downsample(thetaOxt_all,Ndown),downsample(zOxtdFF_all,Ndown),10, [0,0.6,0],'filled')
hold on
plot_areaerrorbar(Oxtmean',options)
xlabel('Theta power (zscored)')
ylabel('OXT')
xlim([-2.5,3.5])
ylim([-3,3])
box off
hold off

disp(['Mean crosscorrelation is :    ', num2str(mean(rOxttheta_all))])
disp(['Std of crosscorrelation is :    ', num2str(std(rOxttheta_all))])
exportgraphics(gcf,'ThetaPowervsOXTACH.pdf','BackgroundColor','none','ContentType','vector')
%%
clear all
animalNames = {
    'ACh_OXT004_017';
    };
animalNames = {'ACh_OXT004_015'};
zAChdFF_all = [];
theta_all = [];
for anID = 1:length(animalNames)
    Names_ACh = dir(['ACh_theta/',animalNames{anID},'*LFPs.mat']);
    filenamesACh = {Names_ACh.name};
    foldernameACh = {Names_ACh.folder};
    load(fullfile(foldernameACh{1},filenamesACh{1}));
    topti = lfp_fbtime;
    fs = THLFPs.samplingRate;
    lfp = double(THLFPs.data(ceil(topti(1)/2e4*fs):end));
    lfp(abs(lfp)>5000) = nan;
    lfp = fillmissing(lfp,'linear');
    [spectrogrammt,tm,fm] = MTSpectrogram(double(lfp),'window',2.0);
    bands = SpectrogramBands(spectrogrammt,fm);

    Names_ACh = dir(['ACh_theta/',animalNames{anID},'*Fluo.mat']);
    filenamesACh = {Names_ACh.name};
    foldernameACh = {Names_ACh.folder};
    load(fullfile(foldernameACh{1},filenamesACh{1}));
    fs_dff  = 100;
    dt = mean(diff(tm));
    Ndown = round(fs_dff*dt);
    dFF_ACh_plot = downsample_yz(FluoData.Ach_dFF,Ndown);

    N1 = min(length(dFF_ACh_plot),length(bands.theta));
    
    zAChdFF = (zscore(dFF_ACh_plot(50:N1)))';
    ztheta = zscore(bands.theta(50:N1));
%     while  ~isempty(find(ztheta>3))
%         idx_turn = find(ztheta>3);
%         ztheta(idx_turn(1)-2:idx_turn(1)+2) = [];
%         zAChdFF(idx_turn(1)-2:idx_turn(1)+2) = [];
%     end
    zAChdFF_all = [zAChdFF_all; zAChdFF];
    theta_all = [theta_all;ztheta];
    for i = 1:12
        idx = find(ztheta>(i-6)*0.5 & ztheta<(i-5)*0.5);
        length(idx);
        AChmean(i,anID) = mean(zAChdFF(idx),'omitnan');
    end
    figure
    set(gcf,'Position',[100 100 400 400])
    scatter(zAChdFF,ztheta,20,[1,0,0],'filled')
    xlim([-2,6])
    ylim([-3,3])
    box off
    corrcoef(zAChdFF,ztheta)


end
%%
figure
plot(zAChdFF);
hold on
plot(ztheta,'r')


%% Simultaneous recording of ACh and OXT  using linear track , data is not useful. OXT and theta has no relation
clear all
animalNames = {
    'ACh_OXT004_017';
    };
animalNames = {'ACh_OXT004_016'};
fs =1250;
for anID = 1:length(animalNames)
    Names_OXT = dir(['OXT_theta/',animalNames{anID},'*LFPs.mat']);
    filenamesOXT = {Names_OXT.name};
    foldernameOXT = {Names_OXT.folder};
    load(fullfile(foldernameOXT{1},filenamesOXT{1}));
    topti = lfp_fbtime;
    lfp = THLFPs.data(ceil(topti(1)*fs/2e4):end);
    lfp(abs(lfp)>5000) = nan;
    lfp = fillmissing(lfp,'linear');
    [spectrogrammt,tm,fm] = MTSpectrogram(double(lfp),'window',2.0);
    bands = SpectrogramBands(spectrogrammt,fm);


    Names_OXT = dir(['OXT_theta/',animalNames{anID},'*Fluo.mat']);
    filenamesOXT = {Names_OXT.name};
    foldernameOXT = {Names_OXT.folder};
    load(fullfile(foldernameOXT{1},filenamesOXT{1}));

    fs_dff  = 100;
    dt = mean(diff(tm));
    Ndown = round(fs_dff*dt);
    dFF_oxt_plot = downsample_yz(FluoData.Oxt_dFF,Ndown);
    N1 = min(length(dFF_oxt_plot),length(bands.theta));
    figure
    set(gcf,'Position',[100 100 400 400])
    scatter(zscore(dFF_oxt_plot(1:N1)),zscore(bands.ratios.hippocampus(1:N1)),20,[0,0.6,0],'filled')
    xlim([-2,6])
    %ylim([-3,3])
    box off
    corrcoef(zscore(dFF_oxt_plot(1:N1)),zscore(zscore(bands.ratios.hippocampus(1:N1))))

end

%%
figure
plot(zscore(dFF_oxt_plot));
hold on
plot(find(sState.states(ceil(topti(1)/2e4):end) ==1),3,'k|','MarkerSize',1);
plot(find(sState.states(ceil(topti(1)/2e4):end) ==3),3.25,'b|','MarkerSize',1);
plot(find(sState.states(ceil(topti(1)/2e4):end) ==5),3.5,'m|','MarkerSize',1);

plot(zscore(bands.theta),'r')

%%


%%
load("OXT_theta/OXYT014_001_dFoverF.mat");
%
fs_dff  = 100;
dt = mean(diff(tm));
Ndown = round(fs_dff*dt);
dFF_oxt_plot = downsample_yz(dFoverF,Ndown);
N1 = min(length(dFF_oxt_plot),length(bands.theta));
figure
set(gcf,'Position',[100 100 600 800])
subplot(2,1,2)
scatter(zscore(dFF_oxt_plot(1:N1)),zscore(bands.ratio(1:N1)),20,'b','filled')
box off

corrcoef(zscore(dFF_oxt_plot(1:N1)),zscore(bands.ratio(1:N1)))
% plot(dFF_oxt_plot(1:N))
% hold on
figure(2)
set(gcf,'Position',[100 100 300 500])
[c_oxt,lags_oxt] = xcorr(dFF_oxt_plot(1:N1),bands.ratio(1:N1),1000,'coeff');
plot(lags_oxt*dt,c_oxt,'Color',[0,0.6,0])
xlim([-15,15])
%%
load("ACh_theta/ACh_OXT004_017_230412_151435_LFPs.mat");
lfp2 = THLFPs.data_aligned;
Nf = 1000;
frerange = linspace(0,250,Nf);
lfp2(abs(lfp2)>5000) = nan;
lfp2 = fillmissing(lfp2,'linear');
[spectrogrammt2,tm2,fm2] = MTSpectrogram(double(lfp2),'window',1.0);
bands2 = SpectrogramBands(spectrogrammt2,fm2);
%%
load("ACh_theta/ACh_OXT004_017_230412_151435_dFoverF.mat");
fs_dff  = 100;
dt = mean(diff(tm2));
Ndown = round(fs_dff*dt);
dFF_ach_plot = downsample_yz(dFoverF,Ndown);
N2 = min(length(dFF_ach_plot),length(bands2.theta));
figure(2)
hold on
set(gcf,'Position',[100 100 300 500])
 %plot(bands2.theta(1:N)/max(bands2.theta(1:N))/20)
 [c_ach,lags_ach] = xcorr(dFF_ach_plot(1:N2),bands2.ratio(1:N2),1000,'coeff');
 plot(lags_ach*dt,c_ach,'Color',[1,0,0])
 xlim([-15,15])

 %%
figure(2)
set(gcf,'Position',[100 100 400 700])
hold on
subplot(3,1,1)
scatter(zscore(bands2.ratio(1:N2)),zscore(dFF_ach_plot(1:N2)),20,'r','filled')
box off
corrcoef(zscore(dFF_ach_plot(1:N2)),zscore(bands2.ratio(1:N2)))
hold on
lsline
ylim([-4,4])
xlim([-1,6])
mdl = fitlm(zscore(bands2.ratio(1:N2)),zscore(dFF_ach_plot(1:N2)));
mdl.Rsquared.Ordinary
hold on
subplot(3,1,2)
scatter(zscore(bands.ratio(1:N1)),zscore(dFF_oxt_plot(1:N1)),20,[0,0.6,0],'filled')
box off
corrcoef(zscore(dFF_oxt_plot(1:N1)),zscore(bands.ratio(1:N1)))
hold on
lsline
ylim([-5,3])
xlim([-1,6])

mdl = fitlm(zscore(dFF_oxt_plot(1:N1)),zscore(bands.ratio(1:N1)));
mdl.Rsquared.Ordinary
subplot(3,1,3)
[c_oxt,lags_oxt] = xcorr(zscore(bands.ratio(1:N1)),zscore(dFF_oxt_plot(1:N1)),100,'normalized');
plot(lags_oxt*dt,c_oxt,'Color',[0,0.6,0])
hold on
[c_ach,lags_ach] = xcorr(zscore(bands2.ratio(1:N2)),zscore(dFF_ach_plot(1:N2)),100,'normalized');
plot(lags_ach*dt,c_ach,'Color',[1,0,0])
xlim([-20,20])
ylim([-0.4,0.65])
xline(0,'-.k')
box off


%% Averaged
zAchdFF = zscore(dFF_ach_plot(1:N2));
zOxtdFF = zscore(dFF_oxt_plot(1:N1));

bratio1 = zscore(bands.ratio(1:N1));
bratio2 = zscore(bands2.ratio(1:N2));
for i = 1:6
    idx = find(bratio1>(i-2) & bratio1<i);
    Oxtmean(i) = mean(zOxtdFF(idx));
    idx = find(bratio2>(i-2) & bratio2<i);
    Achmean(i) = mean(zAchdFF(idx));
end
%%
figure
[c,lags] = xcorr(dFF_ach_plot(1:N),bands2.theta(1:N));
plot(lags,c)
xlim([-100,100])
