fname = 'cl-long-bb_01_2012-02-09_chr2_lollipop_10803_nowaveform';
dirname = 'C:\Users\USER\Documents\mingfaifong\opticalstim ctrl cnqx apv\matfiles\';
timestamps = [3*3600 15*3600];
fid_label = fname;

load([dirname fname]);

%%
clearvars -except type dirID fidlist time timestamps channel excludechannels thresASDR fid_label unit meta
savedir = 'C:\Users\USER\Documents\mingfaifong\opticalstim ctrl cnqx apv\burst analysis\';
excludechannels = [1 8 33 57 64];
channel = channel(ismember(channel,excludechannels)==0);
time = time(ismember(channel,excludechannels)==0);

%% 
dt = 0.01; %1s

pre.time = time(time<timestamps(1));
pre.channel = channel(time<timestamps(1));
pre.unit= unit(time<timestamps(1));
pre.bin = 0:dt:timestamps(1);
pre.asdr = histc(pre.time,pre.bin)/dt;

drug.time = time(time>=timestamps(1) & time<timestamps(2));
drug.channel = channel(time>=timestamps(1) & time<timestamps(2));
drug.unit= unit(time>=timestamps(1) & time<timestamps(2));
drug.bin = timestamps(1):dt:timestamps(2);
drug.asdr = histc(drug.time,drug.bin)/dt;

type = input('treatment: ');

if strcmp(type,'CNQXstim1X')  %this is for stimulated cultures
    drugendtime = (1445/60)*3600;
    drugstarttime = 300; % 5min after CNQX is added
elseif strcmp(type,'CNQX') || strcmp(type,'vehicle') ||  strcmp(type,'APV') || strcmp(type,'TTX')
    if timestamps(2)<24*3600
        drugendtime = timestamps(2);
    else
        drugendtime = 24*3600;
    end
    drugstarttime = 0;
elseif strcmp(type,'12-hour clamp')
    drugstarttime = timestamps(1);
    drugendtime = timestamps(2);
else
    disp('unknown type');
    keyboard;
end

close all;
figure(1);
stairs(pre.bin+dt/2,pre.asdr,'-k');
hold on;
stairs(drug.bin+dt/2,drug.asdr,'-k');
xlim([0 drugendtime])

thresASDR = input('ASDR threshold: ');

%%
% thresASDR = 1000;
minIBI = 2;

x = find(pre.asdr>thresASDR);
pre.burst.ind = [];
pre.burst.ind = x(1);
c = 1;
for i = 1:length(x)
    if pre.bin(x(i))-minIBI>pre.bin(pre.burst.ind(c))
        c = c+1;
        pre.burst.ind(c) = x(i);
    end
end;
pre.burst.time = pre.bin(pre.burst.ind);

x = find(drug.asdr>thresASDR);
drug.burst.ind = [];
drug.burst.ind = x(1);
c = 1;
for i = 2:length(x)
    if drug.bin(x(i))-minIBI>drug.bin(drug.burst.ind(c))
        c = c+1;
        drug.burst.ind(c) = x(i);
    end
end;
drug.burst.time = drug.bin(drug.burst.ind);

plot(pre.burst.time,pre.asdr(pre.burst.ind),'*r');
plot(drug.burst.time,drug.asdr(drug.burst.ind),'*b');

set(gcf,'Position',[4 549 2547 373])
xlim([0 drugendtime])
set(gca,'FontSize',12,'XTick',[-6:3:24]*3600,'XTickLabel',[-6:3:24],'TickDir','out','TickLength',[0.0025 0.025]);
box off;
xlabel('Time (hours)','FontSize',16);
ylabel('Array-Wide Firing Rate (Hz)','FontSize',16);

title([fid_label ' ASDR (bin size=' num2str(dt) 's) w/ detected bursts (thres=' num2str(thresASDR) 'Hz)']);

export_fig([savedir fid_label '_ASDRwithBursts'],'-png','-transparent');
saveas(gca,[savedir fid_label '_ASDRwithBursts'],'fig');

%% get burst snippet times

snip_min = -2;
snip_max = 8;
snippetWidth = [snip_min:dt:snip_max];

figure(5); 
subplot(2,2,1); hold off;
subplot(2,2,2); hold off;
burstcolor = jet(length(pre.burst.ind));
pre.burst.snippet.FR = [];
pre.burst.snippet.spiketimes = [];
pre.burst.snippet.sniptime = 0;

for i = 1:length(pre.burst.ind)
    if pre.burst.ind(i)+snip_min/dt < 1
        pre.burst.snippet.FR(i,:) = nan(size(snippetWidth));
        pre.burst.snippet.timestamps(i,:) = [pre.bin(1),pre.bin(pre.burst.ind(i)+snip_max/dt)];
    elseif pre.burst.ind(i)+snip_max/dt > length(pre.asdr)
        pre.burst.snippet.FR(i,:) = nan(size(snippetWidth));
        pre.burst.snippet.timestamps(i,:) = [pre.bin(pre.burst.ind(i)+snip_min/dt),pre.bin(end)];
    else
        pre.burst.snippet.FR(i,:) = pre.asdr(pre.burst.ind(i)+snip_min/dt:pre.burst.ind(i)+snip_max/dt);
        pre.burst.snippet.timestamps(i,:) = [pre.bin(pre.burst.ind(i)+snip_min/dt),pre.bin(pre.burst.ind(i)+snip_max/dt)];
    end
    pre.burst.snippet.spiketimes = vertcat(pre.burst.snippet.spiketimes,...
        pre.time(pre.time>pre.burst.snippet.timestamps(i,1)...
        & pre.time<pre.burst.snippet.timestamps(i,2)));
    if i>1 && pre.burst.snippet.timestamps(i,1)<pre.burst.snippet.timestamps(i-1,2)
        extratime = pre.burst.snippet.timestamps(i,2)-pre.burst.snippet.timestamps(i-1,2);
    else
        extratime = 10;
    end
    pre.burst.snippet.sniptime = pre.burst.snippet.sniptime + extratime;
    subplot(2,2,1);
    plot(snippetWidth,pre.burst.snippet.FR(i,:),'Color',burstcolor(i,:));
    hold on;
    subplot(2,2,2);
    plot(snippetWidth,pre.burst.snippet.FR(i,:),'Color',burstcolor(i,:));
    hold on;
end

drug.burst.snippet.FR = [];
drug.burst.snippet.spiketimes = [];
drug.burst.snippet.sniptime = 0;
subplot(2,2,3); hold off;
subplot(2,2,4); hold off;
burstcolor = jet(length(drug.burst.ind));

for i = 1:length(drug.burst.ind)
    if drug.burst.ind(i)+snip_min/dt < 1
        drug.burst.snippet.FR(i,:) = nan(size(snippetWidth));
        drug.burst.snippet.timestamps(i,:) = [drug.bin(1),drug.bin(drug.burst.ind(i)+snip_max/dt)];
    elseif drug.burst.ind(i)+snip_max/dt > length(drug.asdr)
        drug.burst.snippet.FR(i,:) = nan(size(snippetWidth));
        drug.burst.snippet.timestamps(i,:) = [drug.bin(drug.burst.ind(i)+snip_min/dt),drug.bin(end)];
    else
        drug.burst.snippet.FR(i,:) = drug.asdr(drug.burst.ind(i)+snip_min/dt:drug.burst.ind(i)+snip_max/dt);
        drug.burst.snippet.timestamps(i,:) = [drug.bin(drug.burst.ind(i)+snip_min/dt),drug.bin(drug.burst.ind(i)+snip_max/dt)];
    end
    drug.burst.snippet.spiketimes = vertcat(drug.burst.snippet.spiketimes,...
        drug.time(drug.time>drug.burst.snippet.timestamps(i,1)...
        & drug.time<drug.burst.snippet.timestamps(i,2)));
    if i>1 && drug.burst.snippet.timestamps(i,1)<drug.burst.snippet.timestamps(i-1,2)
        extratime = drug.burst.snippet.timestamps(i,2)-drug.burst.snippet.timestamps(i-1,2);
    else
        extratime = 10;
    end
    drug.burst.snippet.sniptime = drug.burst.snippet.sniptime + extratime;
    subplot(2,2,3);
    plot(snippetWidth,drug.burst.snippet.FR(i,:),'Color',burstcolor(i,:));
    hold on;
    subplot(2,2,4);
    plot(snippetWidth,drug.burst.snippet.FR(i,:),'Color',burstcolor(i,:));
    hold on;
end

set(gca,'FontSize',12,'TickDir','out','TickLength',[0.0025 0.025],'XTick',[snip_min:1:snip_max],'yscale','log');
box off;
grid on;
xlim([snip_min snip_max]);

subplot(2,2,2);
set(gca,'FontSize',12,'TickDir','out','TickLength',[0.0025 0.025],'XTick',[snip_min:1:snip_max],'yscale','log');
box off;
grid on;
xlim([snip_min snip_max]);

subplot(2,2,1);
set(gca,'FontSize',12,'TickDir','out','TickLength',[0.0025 0.025],'XTick',[snip_min:1:snip_max]);
box off;
grid on;
xlim([snip_min snip_max]);

subplot(2,2,3);
set(gca,'FontSize',12,'TickDir','out','TickLength',[0.0025 0.025],'XTick',[snip_min:1:snip_max]);
xlabel('Time (sec)','FontSize',16);
ylabel('ASDR (Hz)','FontSize',16);
box off;
grid on;
xlim([snip_min snip_max]);

title([fid_label ' ' type ', Burst Overlay']);
subplot(2,2,1);
title([fid_label ' pre-drug, Burst Overlay']);

set(gcf,'Position',[1281 1 1280 957]);

export_fig([savedir fid_label '_BurstOverlay'],'-png','-transparent');
saveas(gca,[savedir fid_label '_BurstOverlay'],'fig');

%% interburst interval histogram
pre.IBI = diff(pre.burst.time);
drug.IBI = diff(drug.burst.time);
xmax = ceil(max([pre.IBI,drug.IBI])/100)*100;
% xmax = 21*60;
xint = ceil(xmax/10/60); %minutes

dtIBI = 20;
binIBI = 0:dtIBI:max([pre.IBI,drug.IBI]);

figure(2);

subplot(2,1,1);
hist(pre.IBI,binIBI);
h = findobj(gca,'Type','patch');
set(h,'FaceColor','r','EdgeColor','k','LineWidth',1)
xlim([0 xmax]);
ylabel('Counts','FontSize',16);
set(gca,'FontSize',12,'XTick',[0:60*xint:xmax],'XTickLabel',[0:1*xint:xmax/60],'TickDir','out','TickLength',[0.0025 0.025]);
box off;
title([fid_label ' IBI histogram, pre-drug (3 hours)']);
% xlim([0 10*60])

subplot(2,1,2);
hist(drug.IBI,binIBI);
h = findobj(gca,'Type','patch');
set(h,'FaceColor','b','EdgeColor','k','LineWidth',1)
xlim([0 xmax]);
set(gca,'FontSize',12,'XTick',[0:60*xint:xmax],'XTickLabel',[0:1*xint:xmax/60],'TickDir','out','TickLength',[0.0025 0.025]);
box off;
xlabel('Interburst Interval (min)','FontSize',16);
ylabel('Counts','FontSize',16);
title([fid_label ' IBI histogram, with drug (24 hours)']);
% xlim([0 10*60])

set(gcf,'Position',[360 388 589 534]);

export_fig([savedir fid_label '_HistogramIBI'],'-png','-transparent');
saveas(gca,[savedir fid_label '_HistogramIBI'],'fig');

%% interburst interval over time
figure(3);
yint = xint;
ymax = xmax;

hold off;
plot(pre.burst.time(2:end),pre.IBI,'ro');
hold on;
plot(drug.burst.time(2:end),drug.IBI,'bo');
xlim([0 drugendtime])
ylim([0 ymax]);
set(gca,'FontSize',12,'XTick',[0 drugendtime],'XTickLabel',[0 drugendtime]/3600, ...
    'YTick',[0:60*yint:ymax],'YTickLabel',[0:yint:ymax/60], ...
    'TickDir','out','TickLength',[0.0025 0.025]);
box off;
grid on;
axis square;
xlabel('Time (hours)','FontSize',16);
ylabel('Interburst Interval (min)','FontSize',16);
set(gcf,'Position',[360 388 589 534]);
% ylim([0 70*60]);

title([fid_label ', IBIs over time']);

export_fig([savedir fid_label '_IBIoverTime'],'-png','-transparent');
saveas(gca,[savedir fid_label '_IBIoverTime'],'fig');


%% interburst interval FR
for i = 1:size(pre.burst.snippet.timestamps,1)
    if i == 1
        pre.IBI.bounds(i,1) = 0;
    else
        pre.IBI.bounds(i,1) = pre.burst.snippet.timestamps(i,2);
    end
    if i == size(pre.burst.snippet.timestamps,1)
        pre.IBI.bounds(i,2) = drugstarttime;        
    else
        pre.IBI.bounds(i,2) = pre.burst.snippet.timestamps(i+1,1);
    end
end

for i = 1:size(drug.burst.snippet.timestamps,1)
    if i == 1
        drug.IBI.bounds(i,1) = drugstarttime;
    else
        drug.IBI.bounds(i,1) = drug.burst.snippet.timestamps(i,2);
    end
    if i == size(drug.burst.snippet.timestamps,1)
        drug.IBI.bounds(i,2) = drugendtime;        
    else
        drug.IBI.bounds(i,2) = drug.burst.snippet.timestamps(i+1,1);
    end
end

% if strcmp(type,'CNQXstim1X')
%     bin_sec = [-6*3600:3*3600:0,300:3*3600:24*3600+300];
% else
    bin_sec = [0:3*3600:drugendtime];
%     bin_sec = [-6*3600:3*3600:24*3600];
% end

IBI.bounds = [pre.IBI.bounds; drug.IBI.bounds];

for i = 1:length(bin_sec)
    for j = 1:size(IBI.bounds,1)
        if bin_sec(i)>IBI.bounds(j,1) && bin_sec(i)<IBI.bounds(j,2)
            IBI.bounds(j+2:end+1,:) = IBI.bounds(j+1:end,:);
            IBI.bounds(j:j+1,:) = [IBI.bounds(j,1),bin_sec(i);bin_sec(i),IBI.bounds(j,2)];
        end
    end
end

m = 1;
dur_IBI = zeros(1,length(bin_sec)-1);
count_spkIBI = zeros(1,length(bin_sec)-1);
for n = 1:size(IBI.bounds,1)
    if diff(IBI.bounds(n,:))<0
        continue;
    end
    if IBI.bounds(n,1)<bin_sec(m)
        continue;
    elseif IBI.bounds(n,2)>bin_sec(m+1)
        m = m+1;
    end
    dur_IBI(1,m) = dur_IBI(1,m) + diff(IBI.bounds(n,:));
    count_spkIBI(1,m) = count_spkIBI(1,m) + length(time(time>=IBI.bounds(n,1) & time<IBI.bounds(n,2)));
end

if strcmp(type,'CNQXstim1X')
    dur_IBI = [dur_IBI(1:2) dur_IBI(4:end)];
    count_spkIBI = [count_spkIBI(1:2) count_spkIBI(4:end)];
end


%% spk an burst counts
count_spk = [length(pre.time(pre.time>=-6*3600 & pre.time<-3*3600)),...
    length(pre.time(pre.time>=-3*3600 & pre.time<0)),...
    length(drug.time(drug.time>=drugstarttime & drug.time<drugstarttime+3*3600)),...
    length(drug.time(drug.time>=drugstarttime+3*3600 & drug.time<drugstarttime+6*3600)),...
    length(drug.time(drug.time>=drugstarttime+6*3600 & drug.time<drugstarttime+9*3600)),...
    length(drug.time(drug.time>=drugstarttime+9*3600 & drug.time<drugstarttime+12*3600)),...
    length(drug.time(drug.time>=drugstarttime+12*3600 & drug.time<drugstarttime+15*3600)),...
    length(drug.time(drug.time>=drugstarttime+15*3600 & drug.time<drugstarttime+18*3600)),...
    length(drug.time(drug.time>=drugstarttime+18*3600 & drug.time<drugstarttime+21*3600)),...
    length(drug.time(drug.time>=drugstarttime+21*3600 & drug.time<drugstarttime+24*3600))];

count_burst = [length(pre.burst.time(pre.burst.time>=-6*3600 & pre.burst.time<-3*3600)),...
    length(pre.burst.time(pre.burst.time>=-3*3600 & pre.burst.time<0)),...
    length(drug.burst.time(drug.burst.time>=drugstarttime & drug.burst.time<drugstarttime+3*3600)),...
    length(drug.burst.time(drug.burst.time>=drugstarttime+3*3600 & drug.burst.time<drugstarttime+6*3600)),...
    length(drug.burst.time(drug.burst.time>=drugstarttime+6*3600 & drug.burst.time<drugstarttime+9*3600)),...
    length(drug.burst.time(drug.burst.time>=drugstarttime+9*3600 & drug.burst.time<drugstarttime+12*3600)),...
    length(drug.burst.time(drug.burst.time>=drugstarttime+12*3600 & drug.burst.time<drugstarttime+15*3600)),...
    length(drug.burst.time(drug.burst.time>=drugstarttime+15*3600 & drug.burst.time<drugstarttime+18*3600)),...
    length(drug.burst.time(drug.burst.time>=drugstarttime+18*3600 & drug.burst.time<drugstarttime+21*3600)),...
    length(drug.burst.time(drug.burst.time>=drugstarttime+21*3600 & drug.burst.time<drugstarttime+24*3600))];
z = [count_spk,count_burst,count_spkIBI,dur_IBI];

%%
% close all;
save([savedir fid_label '_burst.mat']);