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

%% f-I curves:
figure()
gains = [0.1, 0.2, 0.5, 1.0];
for j = 1:length(gains)
	gain = gains(j);
	inputs = 0.0:1.0:60.0;
	rates = ficurve(trials, inputs, tmax, gain);
	plot(inputs, rates);
	hold on;
end
hold off;

%% generate and plot spiketrains for two inputs:
I1 = 10.0;
I2 = 11.0;
gain = 0.1;
spikes1 = lifboltzmannspikes(trials, I1, tmax, gain);
spikes2 = lifboltzmannspikes(trials, I2, tmax, gain);
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:
cmax = 20;
T = 0.1;
figure()
for k = 1:length(gains)
    gain = gains(k);
	spikes1 = lifboltzmannspikes(trials, I1, tmax, gain);
	spikes2 = lifboltzmannspikes(trials, I2, tmax, gain);
    [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('gain=%g', gain))
    hold off;
end

%% discrimination measure:
T = 0.1;
gain = 0.1;
cmax = 15;
spikes1 = lifboltzmannspikes(trials, I1, tmax, gain);
spikes2 = lifboltzmannspikes(trials, I2, tmax, gain);
[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:
T = 0.1;
gains = 0.01:0.01:1.0;
cmax = 100;
dstt = zeros(length(gains), 1);
dsff = zeros(length(gains), 1);
for k = 1:length(gains)
    gain = gains(k);
	spikes1 = lifboltzmannspikes(trials, I1, tmax, gain);
	spikes2 = lifboltzmannspikes(trials, I2, tmax, gain);
    [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);
    dstt(k) = d;
    dsff(k) = min(false1s + false2s);
end
figure()
plot(gains, dstt);
hold on;
plot(gains, dsff);
hold off;