clear all;
close all;

% fid_label = 'trouble22467'; type = 'CNQX';
% load('C:\Users\USER\Documents\mingfaifong\2013-01-21_trouble22467_12-14DIV_MEA_CNQX\trouble22467_12-14DIV.mat'); thresASDR = 2000;

% fid_label = 'trouble22468'; type = 'CNQXstim1X';
% load('C:\Users\USER\Documents\mingfaifong\2013-01-21_trouble22468_12-14DIV_MEA_CNQX-ChR2-stim-CL\trouble22468_12-14DIV.mat'); thresASDR = 5000;

fid_label = 'usain22471'; type = 'CNQXstim1X';
load('C:\Users\USER\Documents\mingfaifong\2013-01-24_usain22471_12-14DIV_MEA_CNQX-ChR2-stim-CL\usain22471_12-14DIV_01.mat'); thresASDR = 5000;

% fid_label = 'usain22493'; type = 'CNQX';
% load('C:\Users\USER\Documents\mingfaifong\2013-01-24_usain22493_12-14DIV_MEA_CNQX\usain22493_12-14DIV.mat'); thresASDR = 5000;

% fid_label = 'wingnut21332'; type = 'CNQXstim1X';
% load('C:\Users\USER\Documents\mingfaifong\2013-02-03_wingnut21332_10-13DIV_MEA_CNQX-ChR2-stim-CL\wingnut21332_10-12DIV_01.mat'); thresASDR = 2000;

% fid_label = 'wingnut22468'; type = 'CNQX';
% load('C:\Users\USER\Documents\mingfaifong\2013-02-03_wingnut22468_10-12DIV_MEA_CNQX\wingnut22468_10-12DIV.mat'); thresASDR = 2000;

% fid_label = 'xinhai22471'; type = 'CNQXstim1X';
% load('C:\Users\USER\Documents\mingfaifong\2013-02-08_xinhai22471_11-13DIV_MEA_CNQX-ChR2-stim-CL\xinhai22471_11-13DIV.mat'); thresASDR = 5000;

% fid_label = 'xinhai22493' ;type = 'CNQX';
% load('C:\Users\USER\Documents\mingfaifong\2013-02-08_xinhai22493_11-13DIV_MEA_CNQX\xinhai22493_11-13DIV.mat'); thresASDR = 500;

% fid_label = 'zarathustra21332'; type = 'CNQXstim2X';
% load('C:\Users\USER\Documents\mingfaifong\2013-02-19_zarathustra21332_12-14DIV_MEA_CNQX-ChR2-stim-CL\zarathustra21332_12-14DIV.mat'); thresASDR = 5000;

% fid_label = 'atlas22493'; type = 'CNQX';
% load('C:\Users\USER\Documents\mingfaifong\2013-03-01_atlas22493_10-12DIV_MEA_CNQX\atlas22493_10-12DIV.mat'); thresASDR = 5000;

% fid_label = 'barbie22468'; type = 'CNQXstim2X';
% load('C:\Users\USER\Documents\mingfaifong\2013-03-08_barbie22468_11-13DIV_MEA_CNQX-ChR2-stim-2x-CL\barbie22468_11-13DIV.mat'); thresASDR = 5000;

% fid_label = 'dojo22471'; type = 'CNQXstim2X';
% load('C:\Users\USER\Documents\mingfaifong\2013-03-19_dojo22471_11-13DIV_MEA_CNQX-ChR2-stim-2x-CL\dojo22471_11-13DIV.mat'); thresASDR = 2000;

% fid_label = 'dojo22493'; type = 'CNQX';
% load('C:\Users\USER\Documents\mingfaifong\2013-03-19_dojo22493_11-13DIV_MEA_CNQX\dojo22493_11-13DIV.mat'); thresASDR = 1000;

% fid_label = 'jujuI5372'; type = 'CNQX';
% load('C:\Users\USER\Documents\mingfaifong\2012-10-28_jujuI5372_11-13DIV_MEA_CNQX\jujuI5372_11-13DIV.mat'); thresASDR = 300;
% 
% fid_label = 'jujuI5388'; type = 'CNQXstim1X';
% load('C:\Users\USER\Documents\mingfaifong\2012-10-28_jujuI5388_11-13DIV_MEA_CNQX-ChR2-stim-CL\jujuI5388_11-13DIV.mat'); thresASDR = 10000;


%%
dt = 0.01; %1s

pre.time = time(time>=-3*3600 & time<0);
pre.channel = channel(time>=-3*3600 & time<0);
pre.unit= unit(time>=-3*3600 & time<0);
pre.bin = -3*3600:dt:0;
pre.asdr = histc(pre.time,pre.bin)/dt;

if strcmp(type,'CNQX')
    if timestamps(2)<24*3600
        drugendtime = timestamps(2);
    else
        drugendtime = 24*3600;
    end
    drugstarttime = 0;
else %this is for stimulated cultures
    drugendtime = (1445/60)*3600;
    drugstarttime = 300; % 5min after CNQX is added
end

drug.time = time(time>=drugstarttime & time<drugendtime);
drug.channel = channel(time>=drugstarttime & time<drugendtime);
drug.unit = unit(time>=drugstarttime & time<drugendtime);
drug.bin = drugstarttime:dt:drugendtime;
drug.asdr = histc(drug.time,drug.bin)/dt;

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

% asdr.bin = 0:20:ceil(max([pre.asdr;drug.asdr])/10)*10;
% pre.asdrsize = histc(pre.asdr,asdr.bin);
% drug.asdrsize = histc(drug.asdr,asdr.bin);
% figure(2);
% bar(asdr.bin,[pre.asdrsize,drug.asdrsize],'stacked');
% set(gca,'yscale','log');

%%
% 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);
plot(pre.burst.time,pre.asdr(pre.burst.ind),'*r');

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(drug.burst.time,drug.asdr(drug.burst.ind),'*b');

set(gcf,'Position',[4 549 2547 373])
xlim([-3*3600 drugendtime])
set(gca,'FontSize',12,'XTick',[-3:3:24]*3600,'XTickLabel',[-3: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(['C:\Users\USER\Documents\mingfaifong\2013-05-28_burstanalysis\' fid_label '_ASDRwithBursts'],'-png','-transparent');
saveas(gca,['C:\Users\USER\Documents\mingfaifong\2013-05-28_burstanalysis\' fid_label '_ASDRwithBursts'],'fig');

%% get burst snippet

% snip_min = -2;
% snip_max = 8;
% snippetWidth = [snip_min:dt:snip_max];
% burstcolor = jet(length(pre.burst.ind));
% 
% figure(5); 
% subplot(2,2,1); hold off;
% subplot(2,2,2); hold off;
% pre.burst.snippet = [];
% for i = 1:length(pre.burst.ind)
%     pre.burst.snippet(i,:) = pre.asdr(pre.burst.ind(i)+snip_min/dt:pre.burst.ind(i)+snip_max/dt);
%     subplot(2,2,1);
%     plot(snippetWidth,pre.burst.snippet(i,:),'Color',burstcolor(i,:));
%     hold on;
%     subplot(2,2,2);
%     plot(snippetWidth,pre.burst.snippet(i,:),'Color',burstcolor(i,:));
%     hold on;
% end
% 
% subplot(2,2,3); hold off;
% subplot(2,2,4); hold off;
% burstcolor = jet(length(drug.burst.ind));
% drug.burst.snippet = [];
% for i = 1:length(drug.burst.ind)
%     drug.burst.snippet(i,:) = drug.asdr(drug.burst.ind(i)+snip_min/dt:drug.burst.ind(i)+snip_max/dt);
%     subplot(2,2,3);
%     plot(snippetWidth,drug.burst.snippet(i,:),'Color',burstcolor(i,:));
%     hold on;
%     subplot(2,2,4);
%     plot(snippetWidth,drug.burst.snippet(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 ' in CNQX, Burst Overlay']);
% subplot(2,2,1);
% title([fid_label ' pre-drug, Burst Overlay']);


%% 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
pre.burst.snippet.sniptime/3600
x1 = pre.time(~ismember(pre.time,pre.burst.snippet.spiketimes));
avg_interburstFR1 = length(x1)/(3*3600-pre.burst.snippet.sniptime)
avg_overall1 = length(pre.time)/(3*3600)

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
drug.burst.snippet.sniptime/3600
x2 = drug.time(~ismember(drug.time,drug.burst.snippet.spiketimes));
avg_interburstFR2 = length(x2)/(24*3600-drug.burst.snippet.sniptime)
avg_overall2 = length(drug.time)/(24*3600)


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 ' in CNQX, Burst Overlay']);
subplot(2,2,1);
title([fid_label ' pre-drug, Burst Overlay']);

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

export_fig(['C:\Users\USER\Documents\mingfaifong\2013-05-28_burstanalysis\' fid_label '_BurstOverlay'],'-png','-transparent');
saveas(gca,['C:\Users\USER\Documents\mingfaifong\2013-05-28_burstanalysis\' 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 = 70*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(['C:\Users\USER\Documents\mingfaifong\2013-05-28_burstanalysis\' fid_label '_HistogramIBI'],'-png','-transparent');
saveas(gca,['C:\Users\USER\Documents\mingfaifong\2013-05-28_burstanalysis\' 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([-3 24]*3600)
ylim([0 ymax]);
set(gca,'FontSize',12,'XTick',[-3:3:24]*3600,'XTickLabel',[-3:3:24], ...
    '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(['C:\Users\USER\Documents\mingfaifong\2013-05-28_burstanalysis\' fid_label '_IBIoverTime'],'-png','-transparent');
saveas(gca,['C:\Users\USER\Documents\mingfaifong\2013-05-28_burstanalysis\' fid_label '_IBIoverTime'],'fig');


%%
% close all;
save(['C:\Users\USER\Documents\mingfaifong\2013-05-28_burstanalysis\' fid_label '_burst.mat']);