clear; clc;
warning('off');

sesslist = {'LAMP3_014_190618_163825';
            'LAMP3_021_190629_155831';
            'LAMP3_023_190630_171532';
            'LAMP3_025_190701_130903';};
        
%% parameters
ip = 4; nchan = 64;
session = sesslist{ip,1};
outsr = 100;
datapath = 'C:\Users\Liang\Dropbox\AchProject\Fig3 data\Stimulation\';

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

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

if ismember(ip,[1,2])
    endpts = 180*outsr;
    achsig = achsig(edgepts:end-endpts);
    achts = achts(edgepts:end-endpts);
elseif ismember(ip,[3,4])
    endpts = 700*outsr;
    achsig = achsig(edgepts:endpts);
    achts = achts(edgepts:endpts);
end

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

%% ach map
posx = interp1(pos.t,pos.rx,achts);
posy = interp1(pos.t,pos.ry,achts);
edges = -0.4:0.01:0.4;
[Np,Xedges,Yedges,binX,binY] = histcounts2(posx,posy,edges,edges);

AChmapc = zeros(length(Xedges),length(Yedges));
AChmapcs = zeros(length(Xedges),length(Yedges));
for ix = 1:length(Xedges)
    for iy = 1:length(Yedges)
        aprphase = binX == ix & binY == iy;
        AChmapc(iy,ix) = nanmean(dFoverF(aprphase));
        AChmapcs(iy,ix) = nanstd(dFoverF(aprphase));
    end
end
save([datapath,session,filesep,'Achmapc.mat'],'AChmapc','AChmapcs');

AChmapc_z = (AChmapc-nanmean(AChmapc,'all'))./nanstd(AChmapc,[],'all');
AChmapc_z(isnan(AChmapc_z)) = -3;
AChmapc_z_sm = imgaussfilt(AChmapc_z,1,'FilterDomain','spatial');
alphamap = ~isnan(AChmapc);

fig3 = figure('Color','w');
colormap('jet');

im = imagesc(AChmapc_z_sm); set(gca,'YDir','normal');
xticks([]); yticks([]); clim([-1,1]*3); 
set(im,'AlphaData',~isnan(alphamap));
cb = colorbar(); cb.Label.String = 'dFoverF (z-scored)'; 

savefig(fig3,[datapath,session,filesep,'Ach_map']);

clear fig1 fig2 fig3;