fid_prefix = 'orion363';
load(['C:\Users\USER\Documents\mingfaifong\datafiles\' fid_prefix '.mat']);

dt = 0.001;
t_smooth = 0.002; %smooth data over 5ms
ISI_max = 5*60; % maximum detected ISI is 5 minutes
t_burst_min = 0.005; %minimum burst length is 10ms

time = time(logical(clean))-events.time(1);
channel = channel(logical(clean));
tmp = [time,channel];
tmp = sortrows(tmp,1);
time = tmp(:,1);
channel = tmp(:,2);
clear tmp;

if strfind(char(events.comment(1)),'CNQX')
    treatment = 'CNQX';
elseif strfind(char(events.comment(1)),'TTX')
    treatment = 'TTX';
else
    treatment = 'vehicle';
end

title1 = [fid_prefix ' before ' treatment];
title2 = [fid_prefix ' after ' treatment];
title3 = [fid_prefix ' after wash'];

pre.time = time(time>-6*3600 & time<=0);
pre.channel = channel(time>-6*3600 & time<=0);
pre.ISI = [NaN; diff(pre.time)];
pre.bins = 0.5*dt:dt:max(pre.ISI(pre.ISI<ISI_max));
pre.hist = hist(pre.ISI,pre.bins);
pre.smooth = smooth(pre.hist.*pre.bins,round(t_smooth/dt));
[xmax imax xmin imin] = extrema(pre.smooth);
pre.burstISIcut = imin(dt*imin>=t_burst_min);
pre.burstISIcut = min(pre.burstISIcut);
pre.isburst = zeros(size(pre.ISI));
pre.isburst(pre.ISI<pre.bins(pre.burstISIcut)) = 1;
pre.FR.overall = length(pre.isburst)/(6*3600);
pre.FR.burst = sum(pre.isburst)/(6*3600);
pre.FR.IBI = sum(abs(pre.isburst-1))/(6*3600);
clear xmax imax xmin imin;

drug.time = time(time>0 & time<=48*3600);
drug.channel = channel(time>0 & time<=48*3600);
drug.ISI = [NaN; diff(drug.time)];
drug.bins = 0.5*dt:dt:max(drug.ISI(drug.ISI<ISI_max));
drug.hist = hist(drug.ISI,drug.bins);
drug.smooth = smooth(drug.hist.*drug.bins,round(t_smooth/dt));
[xmax imax xmin imin] = extrema(drug.smooth);
drug.burstISIcut = imin(dt*imin>=t_burst_min);
drug.burstISIcut = min(drug.burstISIcut);
drug.isburst = zeros(size(drug.ISI));
drug.isburst(drug.ISI<drug.bins(drug.burstISIcut)) = 1;
drug.FR.overall = length(drug.isburst)/(48*3600);
drug.FR.burst = sum(drug.isburst)/(48*3600);
drug.FR.IBI = sum(abs(drug.isburst-1))/(48*3600);
clear xmax imax xmin imin;

if length(events.comment)>2 && strcmp(events.comment(3),'wash')
    wash.time = time(time>=events.time(3)-events.time(1));
    wash.time = wash.time-max(wash.time)+52*3600; %set the wash to be at 52 hours
    wash.time = wash.time(wash.time<58*3600);
    wash.channel = channel(time>=events.time(3)-events.time(1));
    wash.ISI = [NaN; diff(wash.time)];
    wash.bins = 0.5*dt:dt:max(wash.ISI(wash.ISI<ISI_max));
    wash.hist = hist(wash.ISI,wash.bins);
    wash.smooth = smooth(wash.hist.*wash.bins,round(t_smooth/dt));
    [xmax imax xmin imin] = extrema(wash.smooth);
    wash.burstISIcut = imin(dt*imin>=t_burst_min);
    wash.burstISIcut = min(wash.burstISIcut);
    wash.isburst = zeros(size(wash.ISI));
    wash.isburst(wash.ISI<wash.bins(wash.burstISIcut)) = 1;
    wash.FR.overall = length(wash.isburst)/(6*3600);
    wash.FR.burst = sum(wash.isburst)/(6*3600);
    wash.FR.IBI = sum(abs(wash.isburst-1))/(6*3600);
    clear xmax imax xmin imin;
    n_subplots = 3;
else
    n_subplots = 2;
end

close all;
figure('units','normalized','outerposition',[1 0.5 1 0.5])

subplot(1,n_subplots,1);
bar(pre.bins,pre.hist.*pre.bins,1);
hold on;
plot(pre.bins,pre.smooth,'--r','LineWidth',2);
plot(pre.bins(pre.burstISIcut),pre.smooth(pre.burstISIcut),'o','MarkerSize',10,'Color',[0.5 0.5 0.5]);
text(pre.bins(pre.burstISIcut),pre.smooth(pre.burstISIcut)*1.25,[num2str(round(100000*pre.bins(pre.burstISIcut))/100) 'ms'],'HorizontalAlignment','center');
axis([0 pre.bins(pre.burstISIcut)*4 0 pre.smooth(pre.burstISIcut)*4]);
grid on;
set(gca,'TickDir','out');
xlabel('Interspike Interval (seconds)'); ylabel('Counts * ISI');
title([title1,10,'ISI Distribution',10,'(normalized to ISI length)']);

disp([10,fid_prefix,' firing rate data',10,'--------------------------------------']);
disp(['over 6 hours before ' treatment,':',10,...
    ' + overall: ',num2str(pre.FR.overall),'Hz',10,...
    ' + in bursts: ',num2str(pre.FR.burst),'Hz',10,...
    ' + IBI: ',num2str(pre.FR.IBI),'Hz']);



subplot(1,n_subplots,2);
if strmatch('TTX',treatment)
else
    bar(drug.bins,drug.hist.*drug.bins,1);
    hold on;
    plot(drug.bins,drug.smooth,'--r','LineWidth',2);
    plot(drug.bins(drug.burstISIcut),drug.smooth(drug.burstISIcut),'o','MarkerSize',10,'Color',[0.5 0.5 0.5]);
    text(drug.bins(drug.burstISIcut),drug.smooth(drug.burstISIcut)*1.25,[num2str(round(100000*drug.bins(drug.burstISIcut))/100) 'ms'],'HorizontalAlignment','center');
    axis([0 drug.bins(drug.burstISIcut)*4 0 drug.smooth(drug.burstISIcut)*4]);
    
    disp([10,'over 48 hours after ' treatment,':',10,...
    ' + overall: ',num2str(drug.FR.overall),'Hz',10,...
    ' + in bursts: ',num2str(drug.FR.burst),'Hz',10,...
    ' + IBI: ',num2str(drug.FR.IBI),'Hz']);
    disp([10,'percent change from pre-drug condition after ' treatment,':',10,...
    ' + overall: ',num2str(drug.FR.overall/pre.FR.overall*100),'%',10,...
    ' + in bursts: ',num2str(drug.FR.burst/pre.FR.burst*100),'%',10,...
    ' + IBI: ',num2str(drug.FR.IBI/pre.FR.IBI*100),'%']);
end
grid on;
set(gca,'TickDir','out');
xlabel('Interspike Interval (seconds)'); ylabel('Counts * ISI');
title([title2,10,'ISI Distribution',10,'(normalized to ISI length)']);



if n_subplots==3
    subplot(1,n_subplots,3);
    bar(wash.bins,wash.hist.*wash.bins,1);
    hold on;
    plot(wash.bins,wash.smooth,'--r','LineWidth',2);
    plot(wash.bins(wash.burstISIcut),wash.smooth(wash.burstISIcut),'o','MarkerSize',10,'Color',[0.5 0.5 0.5]);
    axis([0 wash.bins(wash.burstISIcut)*4 0 wash.smooth(wash.burstISIcut)*4]);
    text(wash.bins(wash.burstISIcut),wash.smooth(wash.burstISIcut)*1.25,[num2str(round(100000*wash.bins(wash.burstISIcut))/100) 'ms'],'HorizontalAlignment','center');
    grid on;
    set(gca,'TickDir','out');
    xlabel('Interspike Interval (seconds)'); ylabel('Counts * ISI');
    title([title2,10,'ISI Distribution',10,'(normalized to ISI length)']);
    
    disp([10,'for 6 hours after wash:',10,...
    ' + overall: ',num2str(wash.FR.overall),'Hz',10,...
    ' + in bursts: ',num2str(wash.FR.burst),'Hz',10,...
    ' + IBI: ',num2str(wash.FR.IBI),'Hz']);

    disp([10,'percent change from pre-drug condition after wash:',10,...
    ' + overall: ',num2str(wash.FR.overall/pre.FR.overall*100),'%',10,...
    ' + in bursts: ',num2str(wash.FR.burst/pre.FR.burst*100),'%',10,...
    ' + IBI: ',num2str(wash.FR.IBI/pre.FR.IBI*100),'%']);
end

disp([10,'========================================================']);
%%
print('-depsc',['C:\Users\USER\Documents\mingfaifong\datafiles\figures\' fid_prefix '_ISICounts_bin1ms.eps']);
screen2png(gcf,['C:\Users\USER\Documents\mingfaifong\datafiles\figures\' fid_prefix '_ISICounts_bin1ms.png']);
saveas(gcf,['C:\Users\USER\Documents\mingfaifong\datafiles\figures\' fid_prefix '_ISICounts_bin1ms.fig'],'fig');
