% Jon Newman
% Multi-parameter activity assay for cultures transfected with
% optical exciters and silencers

% File name
close all
figure
fid = 'max-power-preserving-triangle_2013-12-16_hoss-r3_smiley_22568_01';
s = loadspk([fid '.spk']);
d = loaddig([fid '.dig']);

%% Decode dig file
b = dec2bin(d.state,8);
u = bin2dec(b(:,1:end-3));
c = bin2dec(b(:,end-2:end-1));
t = bin2dec(b(:,end));

stm.time = d.time(logical(t));
stm.channel = u(logical(t));
% stm.type = u(logical(t));

uu = unique(u(c == 0));
uu = uu(uu~=0);

%%
ss = SqueakySpk(fid,s.meta.fs_Hz,1,s,stm);
ss.RemoveChannel([1 8 33 57 64 35]);

%%
 ss.PeriStimHistogram(1,[-25 100])


xlimit = [0 2500];

%% Find start and end time for all u
tu_on = [];
tu_off = [];
for i = 1:length(uu)
    tu_on = [tu_on d.time(find(u == uu(i),1, 'first'))];
    tu_off = [tu_off d.time(find(u == uu(i),1, 'last'))];
end


%% Calculate ASDR for whole record
% figure
% hold on
% subplot(311)
% ss.RandScat;
% xlim(xlimit)
% subplot(312)
ss.ASDR(0.1,[],[],[],[],0);
% xlim(xlimit)
% subplot(313)
% stairs(d.time,u,'k')
% xlim(xlimit)

%% For each unique u in uu find the average firing rate 

% num detected units
nu = length(unique(ss.unit));
fru = [];
for i = 1:length(uu)
    fru = [fru mean(ss.asdr.asdr(ss.asdr.bin >= tu_on(i) & ss.asdr.bin <= tu_off(i)))/nu];
end

% fr_nom = mean(ss.asdr.asdr(ss.asdr.bin >= 0 & ss.asdr.bin <= min(tu_on)))/nu;
fr_nom = mean(ss.asdr.asdr(ss.asdr.bin >= max(tu_on)+120))/nu;


uFR = [];
sFR = [];
db = 0.1;
bins = 0:db:0.9;

% bin results
for b = bins
    uFR = [uFR,mean(fru(uu>= b & uu < b + db))];
    sFR = [sFR,std(fru(uu>= b & uu < b + db))];
end

%% Plot result
fo = fit(uu,fru','poly2');
us = 0:0.01:1;
est = fo(us);

% figure
semilogy(uu, fru, 'bo')
hold on
semilogy(us,est,'k')


%% Halo
uu = u(u~=0 & c == 1);

%% Find start and end time for all u
tu_on = d.time(ismember(u,uu));
tu_off = tu_on + 20;

% %% Calculate ASDR for whole record
% figure
% hold on
% subplot(311)
% ss.RandScat;
% xlim(xlimit)
% subplot(312)
ss.ASDR(0.1,[],[],[],[],0);
% xlim(xlimit)
% subplot(313)
% stairs(d.time,u,'k')
% xlim(xlimit)

%% For each unique u in uu find the average firing rate 

% num detected units
nu = length(unique(ss.unit));
fru = [];
for i = 1:length(uu)
    fru = [fru mean(ss.asdr.asdr(ss.asdr.bin >= tu_on(i) & ss.asdr.bin <= tu_off(i)-10))/nu];
end
[uuu ind] = sort(uu);
fru = fru(ind);

% fr_nom = mean(ss.asdr.asdr(ss.asdr.bin >= 0 & ss.asdr.bin <= tu_on(1)))/nu;
fr_nom = mean(ss.asdr.asdr(ss.asdr.bin >= 0 & ss.asdr.bin < 600))/nu;

uFR = [];
sFR = [];
db = 0.1;
bins = 0:db:0.9;

% bin results
for b = bins
    uFR = [uFR,mean(fru(uuu>= b & uuu < b + db))];
    sFR = [sFR,std(fru(uuu>= b & uuu < b + db))];
end

fo = fit(uuu,fru','poly2');
us = 0:0.01:1;
est = fo(us);

%% Plot result
semilogy(uuu,fru, 'yo')
semilogy(us,est,'k')
semilogy([0 1],[fr_nom fr_nom],'k--')
set(gca,'Fontsize',12,'TickDir','Out');
axis([0 1 0.01 20])
axis square
box off
xlabel('u','Fontsize',12)
ylabel('Hz/unit','Fontsize',12)

export_fig('-r400','-transparent',fid)



