clear; clc;
warning('off');

sesslist = {'\Non-stimulation\LAMP3_016_190619_180636';
            '\Stimulation\LAMP3_014_190618_163825';};
        
%% parameters
ip = 2; nchan = 64;
session = sesslist{ip,1};
outsr = 100;
datapath = ['F:','\Dropbox\AchProject\Fig3 data\'];

%% load achfile and posfile
% ach 
achfile = dir([datapath,session,filesep,'LAMP*.mat']);
load([achfile.folder,filesep,achfile.name],[achfile.name(1:end-4),'_AChSig1']);
achdata = eval([achfile.name(1:end-4),'_AChSig1']);
load([achfile.folder,filesep,achfile.name],[achfile.name(1:end-4),'_OpTrTTL']);
achopt = eval([achfile.name(1:end-4),'_OpTrTTL']);
load([achfile.folder,filesep,achfile.name],[achfile.name(1:end-4),'_OpGeSti']);
opgeopt = eval([achfile.name(1:end-4),'_OpGeSti']);

achtime(1) = achopt.times(find(achopt.values >= 1,1));
achtime(2) = achopt.times(find(achopt.values >= 1,1,'last'));
syncachts = achdata.times - achtime(1);

% position
load([datapath,session,filesep,'scvajusted.mat'],'pos');

fig1 = figure('Color','w');
scatter(pos.rx,pos.ry,'.');
xlabel('x'); ylabel('y');
title('position');
savefig(fig1,[datapath,session,filesep,'position']);

%% processing ach signal
edgepts = 90*outsr;
achsig = downsample(achdata.values,round(1/(achdata.interval*outsr)));
achts = downsample(syncachts,round(1/(achdata.interval*outsr)));

endpts = 180*outsr;
achsig = achsig(edgepts:end-endpts);
achts = achts(edgepts:end-endpts);
opgesig = interp1(opgeopt.times-achtime(1),opgeopt.values,achts);

p = polyfit(achts,achsig,5);
Ybaseline = polyval(p,achts);
dFoverF = smoothdata((achsig - Ybaseline)./Ybaseline,'movmean',1*outsr);

fig2 = figure('Color','w');
subplot(1,2,1)
plot(achts,achsig); hold on;
plot(achts,Ybaseline,'LineWidth',2);
plot(achts, dFoverF);
title('Ach dFoverF');
subplot(1,2,2)
plot(achts,opgesig);
title('Optogenetics Stimulation');
savefig(fig2,[datapath,session,filesep,'dFoverF']);

%%
posx = interp1(pos.t,pos.rx,achts);
posy = interp1(pos.t,pos.ry,achts);
spd = [0;sqrt(diff(posx,1).^2+diff(posy,1).^2)]*outsr;
spd(spd>1) = nan;
spd = smoothdata(spd,'movmean',0.5*outsr);

sti_on = find(opgesig>1);

fig1 = figure('Color','w');
scatter(posx,posy,'g.'); hold on;
scatter(posx(sti_on),posy(sti_on),'b.'); hold on;
xlabel('x'); ylabel('y');
title('position');
savefig(fig1,[datapath,session,filesep,'position']);

%%
epochs_ids = zeros(0);
jumps = find(diff(sti_on)>3*outsr);
epochs_ids(:,1) = [1;sti_on(jumps'+1)];
epochs_ids(:,2) = [sti_on(jumps');sti_on(end)];
epochs_ids(diff(epochs_ids,1,2)<15*outsr,:) = [];

figure;
plot(achts,opgesig); hold on;
achplot = zeros(0);
speedplot = zeros(0);
for ie = 1:size(epochs_ids,1)
    plot([1,1]*achts(epochs_ids(ie,1)),[0,10],'b-'); hold on;
    plot([1,1]*achts(epochs_ids(ie,2)),[0,10],'r-');
    epochs_ids(ie,2) = epochs_ids(ie,1) + min(diff(epochs_ids,1,2));
    plotwindow = epochs_ids(ie,:) + [-1,1]*10*outsr;
    achplot(ie,:) = dFoverF(plotwindow(1):plotwindow(2));
    speedplot(ie,:) = spd(plotwindow(1):plotwindow(2));
end

figure('Color','white');
yyaxis left;
options.color_area = [128 193 219]./255;    % Blue theme
options.color_line = [ 52 148 186]./255;
% options.color_area = [243 169 114]./255;    % Orange theme
% options.color_line = [236 112  22]./255;
options.alpha      = 0.5;
options.line_width = 2;
options.error      = 'sem';
options.x_axis     = [1:size(achplot,2)]./outsr-10;
h = plot_areaerrorbar(achplot,options);
hold on;
pch = patch([0,20,20,0],[-5,-5,5,5],'b');
pch.FaceAlpha = 0.1; pch.EdgeColor = 'none';
legend([h(1)],sprintf('n=%d',size(achplot,1)));
xlim([-10,30]); ylim([-0.05,0.05]); 
yticks([-0.05,0,0.05]); yticklabels([-5,0,5]);
xlabel('Time (s)'); ylabel('Ach (\Delta F/F %)');
yyaxis right;
% options.color_area = [128 193 219]./255;    % Blue theme
% options.color_line = [ 52 148 186]./255;
options.color_area = [243 169 114]./255;    % Orange theme
options.color_line = [236 112  22]./255;
options.alpha      = 0.5;
options.line_width = 2;
options.error      = 'sem';
options.x_axis     = [1:size(speedplot,2)]./outsr-10;
h = plot_areaerrorbar(speedplot,options);
hold on; spdlim = 0.5;
pch = patch([0,20,20,0],[-5,-5,5,5],'b');
pch.FaceAlpha = 0.1; pch.EdgeColor = 'none';
legend([h(1)],sprintf('n=%d',size(speedplot,1)));
xlim([-10,30]); ylim([0,1]*spdlim); 
yticks([0,0.5,1]*spdlim); yticklabels([0,0.5,1]*spdlim*100)
ylabel('Speed (cm/s)');
savefig(gcf,'perievent_sti','compact');

save([datapath,session,filesep,'achplot.mat'],'achplot','speedplot');

%%
load([datapath,session,filesep,'achplot.mat'],'achplot','speedplot');

boxdata = [];
boxdata(:,1) = mean(achplot(:,1:1000),2);
boxdata(:,2) = mean(achplot(:,1001:3000),2);
boxdata(:,3) = mean(achplot(:,3001:4000),2);

ca1 = [128 128 128]./255;    % grey theme
cl1 = [172 172 172]./255;

ca2 = [128 193 219]./255;    % Blue theme
cl2 = [ 52 148 186]./255;

figure('Color','white');
iosr.statistics.boxPlot([1,3],boxdata(:,[1,3]),'theme','colorall',...
            'showScatter',true,'scatterMarker','.','symbolMarker','none',...
            'scatterColor',cl1 ,'scatterSize', 128, 'scatterAlpha', 0.5,...
            'boxWidth',0.4,'boxColor',ca1,'showViolin',false); hold on;
iosr.statistics.boxPlot([2],boxdata(:,[2]),'theme','colorall',...
            'showScatter',true,'scatterMarker','.','symbolMarker','none',...
            'scatterColor',cl2 ,'scatterSize', 128, 'scatterAlpha', 0.5,...
            'boxWidth',0.4,'boxColor',ca2,'showViolin',false); hold on;
xticks([1,2,3]); xticklabels({'pre-delay ','delay','post-delay'})
set(gca,'xcolor',[1 1 1]); xlim([0.5,3.5]);
yticks([-0.06:0.02:0.06]); yticklabels([-0.06:0.02:0.06].*100);
ylabel('Ach (\Delta F/F %)');  ylim([-0.07,0.07]);
 
sigstar([1,2],0.0001);
sigstar([2,3],0.0001);