clear; clc;
warning('off');

sesslist = {'\Non-stimulation\LAMP3_016_190619_180636';
            '\Stimulation\LAMP3_014_190618_163825';};
        
%% parameters
ip = 1; 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']);

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);

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

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

%%
plot(pos.rx,pos.ry,'k.'); hold on;
delayids = find(pos.rx>-0.085 & pos.rx<0.115 & pos.ry < -0.085);
xlim([-0.4,0.4]); xticks([-0.4,0,0.4]); xlabel('x (m)');
ylim([-0.4,0.4]); yticks([-0.4,0,0.4]); ylabel('y (m)');
title('Adjusted csv');

delayepochs(:,2) = pos.t(pos.rx(1:end-1)>-0.085 & pos.rx(1:end-1)<0.115 & pos.ry(1:end-1) < -0.085 ...
                        & pos.rx(2:end)>-0.085 & pos.rx(2:end)<0.115 & pos.ry(2:end) > -0.085);
delayepochs(:,1) = delayepochs(:,2) - 20;
for ie = 1:size(delayepochs,1)
    delayepochs_ids = pos.t > delayepochs(ie,1) & pos.t < delayepochs(ie,2);
    plot(pos.rx(delayepochs_ids),pos.ry(delayepochs_ids),'g.');
end

%%
achplot = zeros(0);
for ie = 4:size(delayepochs,1)-4
    plotwindow = delayepochs(ie,:) + [-1,1]*10;
    delayepochs_ids = achts > plotwindow(1) & achts < plotwindow(2);
    achplot(ie,:) = dFoverF(delayepochs_ids);
end

figure;
% 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('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);

achplot = zeros(0);
speedplot = zeros(0);
for ie = 4:size(delayepochs,1)-4
    plotwindow = delayepochs(ie,:) + [-1,1]*10;
    delayepochs_ids = achts > plotwindow(1) & achts < plotwindow(2);
    achplot(ie,:) = dFoverF(delayepochs_ids);
    speedplot(ie,:) = spd(delayepochs_ids);
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_non','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,'theme','colorall',...
            'showScatter',true,'scatterMarker','.','symbolMarker','none',...
            'scatterColor',cl1 ,'scatterSize', 128, 'scatterAlpha', 0.5,...
            'boxWidth',0.4,'boxColor',ca1,'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.06,0.06]);

sigstar([1,2],0.55);
sigstar([2,3],0.25);
