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

%%
spd = getSpd('ACh_spd/Linear track and open area','ACh_OXT004_015.csv');

Ns_t = length(spd.time);
load('ACh_spd/Linear track and open area/ACh_OXT004_015_230411_150537_Fluo.mat')
Nf_t = length(FluoData.Ach_dFF);

NN = min(Ns_t,Nf_t);
Ndown = 10;

Ach_plot = downsample(FluoData.Ach_dFF(4.5e4:NN)*100,Ndown);
Oxt_plot = downsample(FluoData.Oxt_dFF(4.5e4:NN)*100,Ndown);
spd_plot = downsample(spd.data(4.5e4:NN)*100,Ndown);
figure
set(gcf,'Position',[100 100 600 1200])
subplot(3,1,1)
scatter(spd_plot,Ach_plot,20, 'r','filled')
xlim([-2,30])
ylim([-1,1])
hold on
H = lsline;
P = polyfit(spd_plot,Ach_plot,1);
yfit = polyval(P,spd_plot);
hold on;
%plot(spd_plot,yfit,'k-.');
eqn = string(" Linear: y = " + P(1)) + "x + " + string(P(2));
text(min(spd_plot),max(Ach_plot),eqn,"HorizontalAlignment","left","VerticalAlignment","bottom")
box off
mdl = fitlm(spd_plot,Ach_plot);

subplot(3,1,2)
scatter(spd_plot,Oxt_plot,20,'g','filled')
xlim([-2,30])
ylim([-1,1])
hold on
lsline
P = polyfit(spd_plot,Oxt_plot,1);
yfit = polyval(P,spd_plot);
hold on;
%plot(spd_plot,yfit,'k-.');
eqn = string(" Linear: y = " + P(1)) + "x + " + string(P(2));
text(min(spd_plot),max(Oxt_plot),eqn,"HorizontalAlignment","left","VerticalAlignment","bottom")
box off

subplot(3,1,3)
scatter(Ach_plot,Oxt_plot,20,'k','filled')
xlim([-1,1])
ylim([-1,1])
hold on
H = lsline;
P = polyfit(Ach_plot,Oxt_plot,1);
yfit = polyval(P,Ach_plot);
eqn = string(" Linear: y = " + P(1)) + "x + " + string(P(2));
text(min(Ach_plot),max(Oxt_plot),eqn,"HorizontalAlignment","left","VerticalAlignment","bottom")
box off

%% Averaged AcH
clear all

% bad sessions are 'LAMP37_009';
% mediocore session 'LAMP2_042'
animalNames = {
    'LAMP37_012';
    'LAMP37_011';
    'ACh_OXT004_017';
    'ACh_OXT004_016';
    'ACh_OXT004_015';
    'LAMP2_042';};

zAchdFF_all = [];
spd_all = [];
rAchspd_all =[];
for anID = 1:length(animalNames)
    Names_LIMP = dir(['ACh_spd/',animalNames{anID},'*Fluo.mat']);
    filenamesLIMP = {Names_LIMP.name}
    foldernameLIMP = {Names_LIMP.folder};
    load(fullfile(foldernameLIMP{1},filenamesLIMP{1}));
    AchdFF = (FluoData.Ach_dFF)*100;
    Nf = length(AchdFF);
    spd = getSpd('ACh_spd',[animalNames{anID},'.csv']);
    Ns = length(spd.data);
    NN = min(Nf,Ns);
    zAchdFF = zscore(AchdFF(1e4:NN-2e3));
    zspd = spd.data(1e4:NN-2e3)*100;
    zAchdFF_all = [zAchdFF_all; zAchdFF];
    spd_all = [spd_all;zspd];
    for i = 1:15
        idx = find(zspd>(i-1)*2 & zspd<(i*2));
        Achmean(i,anID) = mean(zAchdFF(idx),'omitnan');
    end
    [r,p] = corrcoef(zAchdFF,zspd);
    rAchspd_all(anID) = r(1,2);
    pAchspd_all(anID) = p(1,2);

    
end
%
fig1 = figure();
options.handle     =  fig1;
set(gcf,'Position',[100 100 400 700])
options.color_area = [0 0 0]./255;    
options.color_line = [0 0 0]./255;
options.alpha      = 0.5;
options.line_width = 2;
options.error      = 'sem';
options.x_axis = 1:2:29;
Ns_t = length(spd_all);
Nf_t = length(zAchdFF_all);

NN = min(Ns_t,Nf_t);
Ndown = 50;
subplot(2,1,1)
options.color_area = [0 0 0]./255;    
options.color_line = [0 0 0]./255;
hold on
Achmean = fillmissing(Achmean,"linear");
scatter(downsample(spd_all(1:NN),Ndown),downsample(zAchdFF_all(1:NN),Ndown),10, [1,0,0],'filled')
plot_areaerrorbar(Achmean',options)
xlabel('Speed')
ylabel('AcH')
xlim([0,5])
ylim([-3,3])
box off
disp(['Mean crosscorrelation is :    ', num2str(mean(rAchspd_all))])
disp(['Std of crosscorrelation is :    ', num2str(std(rAchspd_all)/sqrt(length(rAchspd_all)))])
% Averaged OXT
%clear all
% mdl = fitlm(downsample(spd_all(1:NN),Ndown),downsample(zAchdFF_all(1:NN),Ndown));
%     figure;plot(mdl)

% bad sessions are 'ACh_OXT004_017'   'ACh_OXT007_011';
% 'ACh_OXT005_014';  'ACh_OXT005_013'; 'ACh_OXT005_012';
animalNames = {
    'ACh_OXT004_017';
    'ACh_OXT004_016';
    'ACh_OXT004_015';
    'ACh_OXT007_010';
    'ACh_OXT007_012';};
%animalNames = {'ACh_OXT007_012'};
zOxtdFF_all = [];
spd_all = [];
rOxtspd_all = [];
pOxtspd_all = [];
for anID = 1:length(animalNames)
    Names_LIMP = dir(['Oxt_spd/',animalNames{anID},'*Fluo.mat']);
    filenamesLIMP = {Names_LIMP.name}
    foldernameLIMP = {Names_LIMP.folder};
    load(fullfile(foldernameLIMP{1},filenamesLIMP{1}));
    OxtdFF = (FluoData.Oxt_dFF)*100;
    Nf = length(OxtdFF);
    spd = getSpd('Oxt_spd',[animalNames{anID},'.csv']);
    Ns = length(spd.data);
    NN = min(Nf,Ns);
    zOxtdFF = zscore(OxtdFF(1e4:NN-2e3));
    zspd = spd.data(1e4:NN-2e3)*100;
    zOxtdFF_all = [zOxtdFF_all; zOxtdFF];
    spd_all = [spd_all;zspd];
    for i = 1:15
        idx = find(zspd>(i-1)*2 & zspd<(i*2));
        Oxtmean(i,anID) = mean(zOxtdFF(idx),'omitnan');
    end
    r = corrcoef(zOxtdFF,zspd);
    rOxtspd_all(anID) = r(1,2);
    pOxtspd_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 = 1:2:29;
Ns_t = length(spd_all);
Nf_t = length(zOxtdFF_all);

NN = min(Ns_t,Nf_t);
Ndown = 30;
subplot(2,1,2)
options.color_area = [0 0 0]./255;    
options.color_line = [0 0 0]./255;
hold on
Oxtmean = fillmissing(Oxtmean,"linear");
scatter(downsample(spd_all(1:NN),Ndown),downsample(zOxtdFF_all(1:NN),Ndown),10, [0,0.6,0],'filled')
plot_areaerrorbar(Oxtmean',options)
xlabel('Speed')
ylabel('OXT')
xlim([0,5])
ylim([-3,3])
box off
disp(['Mean crosscorrelation is :    ', num2str(mean(rOxtspd_all))])
disp(['Std of crosscorrelation is :    ', num2str(std(rOxtspd_all)/sqrt(length(rOxtspd_all)))])
% Averaged OXT
exportgraphics(gcf,'SpeedvsOXTACH.pdf','BackgroundColor','none','ContentType','vector')

%%
figure
plot(Oxt_plot/std(Oxt_plot)-3,'Color', [0,0.6,0])
hold on
plot(Ach_plot/std(Ach_plot)+7,'Color', [1,0,0])
plot(spd_plot/std(spd_plot),'Color', 'k')
box off

%%
ids = find(spd.data(1e3:NN)>0.1)+1e3;
Oxt_th = FluoData.Oxt_dFF(ids);
Ach_th = FluoData.Ach_dFF(ids);
spd_th = spd.data(ids);

figure
set(gcf,'Position',[100 100 600 1200])
subplot(3,1,1)
scatter(spd_th*100,Ach_th*100,1, 'r')

xlim([-2,40])
ylim([-1,1])
subplot(3,1,2)
scatter(spd_th*100,Oxt_th*100,1)

xlim([-2,40])
ylim([-1,1])
subplot(3,1,3)
scatter(Ach_th*100,Oxt_th*100,1)
xlim([-3,3])
ylim([-1,1])