clear; clc;
warning('off');

sesslist = {'LAMP3_004_190612_151313';
            'LAMP3_006_190613_194519';
            'LAMP3_011_190617_142002';
            'LAMP3_015_190618_172714'};

%% parameters
ip = 3; nchan = 64;
session = sesslist{ip,1};
outsr = 100;
datapath = 'C:\Users\Liang\Dropbox\AchProject\Fig2 data\Ripple_ACh\';

%% load achfile and lfpfile
% ach 
achfile = dir([datapath,session,filesep,session(1:9),'.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'));

% lfp
load([datapath,session,filesep,'LFPs.mat'],'LFPs');
info = read_Intan_RHD2000_file([datapath,session],'*.rhd');
digitalin = read_Intan_digitalin_file([datapath,session,filesep,'digitalin.dat'],info.num_board_dig_in_channels)';

lfptime(1) = find(digitalin >= 1,1) / info.sample_rate;
lfptime(2) = find(digitalin >= 1,1,'last') / info.sample_rate;

ripfile = dir([datapath,session,filesep,session,'.ripples.events.mat']);
load([ripfile.folder,filesep,ripfile.name])
ripples.peaks = ripples.peaks - lfptime(1);

% sync timestamp
syncachts = achdata.times - achtime(1);
synclfpts = LFPs.timestamps - lfptime(1);

%% processing ach signal
edgepts = 90*outsr; endpts = 180*outsr;
achsig = downsample(achdata.values,round(1/(achdata.interval*outsr)));
achts = downsample(syncachts,round(1/(achdata.interval*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);

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

%% LPFs power
[spectrogram,chrt,freq] = MTSpectrogram(double(LFPs.data),'window',1);
bands = SpectrogramBands(spectrogram,freq);
ripplepower = interp1(chrt-lfptime(1),pow2db(bands.ripples),achts);
ripplepower_sm = smoothdata(ripplepower,'movmean',1*outsr);

fig2 = figure;
plot(achts,ripplepower_sm);
title('Ripple power');
savefig(fig2,[datapath,session,filesep,'Ripple_Power']);


syncts = achts;
%% scatter plot
win = 5;
lfpids = length(syncts);
nwins = floor(length(syncts)/(win*outsr));
ripplepower_sm_bined = nanmean(reshape(ripplepower_sm(1:nwins*win*outsr),win*outsr,[]),1)';
dFoverF_bined = nanmean(reshape(dFoverF(1:nwins*win*outsr),win*outsr,[]),1)';

fig3 = figure;
scatter(dFoverF_bined,ripplepower_sm_bined,'filled');
title('dFoverF vs. Ripple power');
savefig(fig3,[datapath,session,filesep,'dFoverF_vs_Ripple_Power']);

scatdata = [dFoverF_bined,ripplepower_sm_bined];
save([datapath,session,filesep,'scatdata.mat'],'scatdata');

clear fig1 fig2 fig3;


%% scatter plot
win = 5;
nwins = floor(length(syncts)/(win*outsr));
achtsdown = nanmean(reshape(achts(1:nwins*win*outsr),win*outsr,[]),1)';
dFoverF_bined = nanmean(reshape(dFoverF(1:nwins*win*outsr),win*outsr,[]),1)';

achtsdown_edges = [achtsdown-0.5*win;achtsdown(end)+0.5*win];
ripple_bined = histcounts(ripples.peaks,achtsdown_edges);
ripple_bined_sm = smoothdata(ripple_bined,'gaussian',6)'./win;

fig4 = figure;
subplot(2,1,1)
scatter(dFoverF_bined*100,ripple_bined_sm,'filled');
xlabel('dFoverF (%)'); ylabel('Ripple rate (Hz)')
title('dFoverF vs. Ripple rate');
xlim([-15,15]);
subplot(2,1,2)
achedge = linspace(-15,15,7);
achc = achedge(1:end-1) + 0.5*diff(achedge(1:2));
[N,~,binX] = histcounts(dFoverF_bined*100,achedge);
histmean = zeros(0);
ca2 = [128 193 219]./255;    % Blue theme
cl2 = [ 52 148 186]./255;
for ib = 1:length(achedge)-1
    histmean{ib,1} = ripple_bined_sm(binX == ib);
    iosr.statistics.boxPlot(ib,histmean{ib,1},'theme','colorall',...
        'showScatter',true,'scatterMarker','.','symbolMarker','none',...
        'scatterColor',cl2 ,'scatterSize', 128, 'scatterAlpha', 0.5,...
        'boxWidth',0.4,'boxColor',ca2,'showViolin',false); hold on;
end
xlim([1,length(achedge)-1]+0.25*0.5*diff(achedge(1:2))*[-1,1]);
xticks(1:length(achedge)-1); xticklabels(achc);
xlabel('dFoverF (%)'); ylabel('Ripple rate (Hz)')
title('dFoverF vs. Ripple rate');
savefig(fig4,[datapath,session,filesep,'dFoverF_vs_Ripple_rate']);

scatdata_ripple = [dFoverF_bined,ripple_bined_sm];
save([datapath,session,filesep,'scatdata_ripple.mat'],'scatdata_ripple');

%% peri-event
ripple_peri_event = zeros(0);
ripples.peaks(ripples.peaks<achts(1)+30 | ripples.peaks>achts(end)-30) = [];
for irip = 1:length(ripples.peaks)
    plotwindow = [-30,30] + ripples.peaks(irip);
    achids = find(achts > plotwindow(1) & achts <= plotwindow(2));
    if isempty(achids), continue; end
    ripple_peri_event(irip,:) = dFoverF(achids);
end

fig5 = figure;

imagesc([1:size(ripple_peri_event,2)./outsr]'-30,[1:size(ripple_peri_event,1)]',...
        zscore(ripple_peri_event,[],2));
xlabel('Times from ripple peaks (s)'); ylabel('Ripple events #');

savefig(fig5,[datapath,session,filesep,'ripple_peri_event']);
save([datapath,session,filesep,'ripple_peri_event.mat'],'ripple_peri_event');