%% general settings for the model neuron:
trials = 10;
tmax = 50.0;

%% generate and plot spiketrains for two inputs:
I1 = 14.0;
I2 = 15.0;
spikes1 = lifspikes(trials, I1, tmax);
spikes2 = lifspikes(trials, I2, tmax);
subplot(1, 2, 1);
tmin = 10.0;
spikeraster(spikes1, tmin, tmin+2.0);
title(sprintf('I_1=%g', I1))
subplot(1, 2, 2);
spikeraster(spikes2, tmin, tmin+2.0);
title(sprintf('I_2=%g', I2))
%savefigpdf(gcf(), 'spikeraster.pdf')

%% spike count histograms:
Ts = [0.01 0.1 0.3 1.0];
cmax = 100;
figure()
for k = 1:length(Ts)
    T = Ts(k);
    [c1, b1] = counthist(spikes1, 0.0, tmax, T, cmax);
    [c2, b2] = counthist(spikes2, 0.0, tmax, T, cmax);
    subplot(2, 2, k)
    bar(b1, c1, 'r');
    hold on;
    bar(b2, c2, 'b');
    xlim([0 cmax])
    title(sprintf('T=%gms', 1000.0*T))
    hold off;
end

%% discrimination measure:
T = 0.1;
cmax = 15;
[d, thresholds, true1s, false1s, true2s, false2s, pratio] = discriminability(spikes1, spikes2, tmax, T, cmax);
figure()
subplot(1, 3, 1);
plot(thresholds, true1s, 'b');
hold on;
plot(thresholds, true2s, 'b');
plot(thresholds, false1s, 'r');
plot(thresholds, false2s, 'r');
xlim([0 cmax])
hold off;
% Ratio:
subplot(1, 3, 2);
fprintf('discriminability = %g\n', d);
plot(thresholds, pratio);
% ROC:
subplot(1, 3, 3);
plot(false2s, true1s);

%% discriminability:
Ts = 0.01:0.01:1.0;
cmax = 100;
ds = zeros(length(Ts), 1);
for k = 1:length(Ts)
    T = Ts(k);
    [c1, b1] = counthist(spikes1, 0.0, tmax, T, cmax);
    [c2, b2] = counthist(spikes2, 0.0, tmax, T, cmax);
    [d, thresholds, true1s, false1s, true2s, false2s, pratio] = discriminability(spikes1, spikes2, tmax, T, cmax);
    ds(k) = d;
end
figure()
plot(Ts, ds)