fixed fano slop eproject

projects/project_fano_slope/solution/counthist.m

@@ -0,0 +1,11 @@
+function [counts, cbins] = counthist(spikes, tmin, tmax, T, cmax)
+tbins = tmin+T/2:T:tmax;
+cbins = 0.5:cmax;
+counts = zeros(1, length(cbins));
+for k = 1:length(spikes)
+    times = spikes{k};
+    n = hist(times((times>=tmin)&(times<=tmax)), tbins);
+    counts = counts + hist(n, cbins);
+end
+counts = counts / sum(counts);
+end

projects/project_fano_slope/solution/discriminability.m

@@ -0,0 +1,23 @@
+function [d, thresholds, true1s, false1s, true2s, false2s, pratio] = discriminability(spikes1, spikes2, tmax, T, cmax)
+[c1, b1] = counthist(spikes1, 0.0, tmax, T, cmax);
+[c2, b2] = counthist(spikes2, 0.0, tmax, T, cmax);
+thresholds = 0:cmax;
+true1s = zeros(length(thresholds), 1);
+true2s = zeros(length(thresholds), 1);
+false1s = zeros(length(thresholds), 1);
+false2s = zeros(length(thresholds), 1);
+for k = 1:length(thresholds)
+    th = thresholds(k);
+    t1 = sum(c1(b1<=th));
+    f1 = sum(c1(b1>th));
+    t2 = sum(c2(b2>=th));
+    f2 = sum(c2(b2<th));
+    true1s(k) = t1;
+    true2s(k) = t2;
+    false1s(k) = f1;
+    false2s(k) = f2;
+end
+%pratio = (true1s + true2s)./(false1s+false2s);
+pratio = (true1s + true2s)/2;
+d = max(pratio);
+end

projects/project_fano_slope/solution/fanoslope.m

@@ -0,0 +1,92 @@
+%% 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;
+ds = 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);
+    ds(k) = d;
+end
+figure()
+plot(gains, ds)
+
+

projects/project_fano_slope/solution/ficurve.m

@@ -0,0 +1,8 @@
+function rates = ficurve(trials, inputs, tmax, gain)
+% compute f-I curve.
+rates = zeros(length(inputs), 1);
+for k=1:length(inputs)
+	spikes = lifboltzmannspikes(trials, inputs(k), tmax, gain);
+	rates(k) = firingrate(spikes, 0.0, tmax);
+end
+end

projects/project_fano_slope/solution/firingrate.m

@@ -0,0 +1,9 @@
+function rate = firingrate(spikes, tmin, tmax)
+% mean firing rate between tmin and tmax.
+rates = zeros(length(spikes), 1);
+for i = 1:length(spikes)
+	times= spikes{i};
+	rates(i) = length(times((times>=tmin)&(times<=tmax)))/(tmax-tmin);
+end
+rate = mean(rates);
+end

projects/project_fano_slope/solution/lifboltzmannspikes.m

@@ -0,0 +1,35 @@
+function spikes = lifboltzmannspikes(trials, input, tmax, gain)
+% Generate spike times of a leaky integrate-and-fire neuron.
+% trials: the number of trials to be generated.
+% input: the stimulus intensity.
+% tmax: the duration of a trial.
+% gain: gain of the neuron, i.e. the slope factor of the boltzmann input.
+
+    tau = 0.01;
+    D = 1e-2;
+    imax = 25;
+    ithresh = 10;
+    slope = gain;
+    vreset = 0.0;
+    vthresh = 10.0;
+    dt = 1e-4;
+    
+    n = ceil(tmax/dt);
+    inb = imax./(1.0+exp(-slope.*(input - ithresh)));
+    spikes = cell(trials, 1);
+    for k=1:trials
+        times = [];
+        j = 1;
+        v = vreset;
+        noise = sqrt(2.0*D)*randn(n, 1)/sqrt(dt);
+        for i=1:n
+            v = v + (- v + noise(i) + inb)*dt/tau;
+            if v >= vthresh
+                v = vreset;
+                times(j) = i*dt;
+                j = j + 1;
+            end
+        end
+        spikes{k} = times;
+    end
+end

projects/project_fano_slope/solution/savefigpdf.m

@@ -0,0 +1,28 @@
+function savefigpdf(fig, name, width, height)
+% Saves figure fig in pdf file name.pdf with appropriately set page size
+% and fonts
+
+% default width:
+if nargin < 3
+    width = 11.7;
+end
+% default height:
+if nargin < 4
+    height = 9.0;
+end
+
+% paper:
+set(fig, 'PaperUnits', 'centimeters');
+set(fig, 'PaperSize', [width height]);
+set(fig, 'PaperPosition', [0.0 0.0 width height]);
+set(fig, 'Color', 'white')
+
+% font:
+set(findall(fig, 'type', 'axes'), 'FontSize', 12)
+set(findall(fig, 'type', 'text'), 'FontSize', 12)
+
+% save:
+saveas(fig, name, 'pdf')
+
+end
+

projects/project_fano_slope/solution/spikeraster.m

@@ -0,0 +1,30 @@
+function spikeraster(spikes, tmin, tmax)
+% Display a spike raster of the spike times given in spikes.
+%
+% spikeraster(spikes, tmax)
+%   spikes: a cell array of vectors of spike times in seconds
+%   tmin: plot spike raster starting at tmin seconds
+%   tmax: plot spike raster upto tmax seconds
+
+ntrials = length(spikes);
+for k = 1:ntrials
+    times = spikes{k};
+    times = times((times>=tmin) & (times<=tmax));
+    if tmax < 1.5
+        times = 1000.0*times; % conversion to ms
+    end
+    for i = 1:length( times )
+      line([times(i) times(i)],[k-0.4 k+0.4], 'Color', 'k'); 
+    end
+end
+if (tmax-tmin) < 1.5
+    xlabel('Time [ms]');
+    xlim([1000.0*tmin 1000.0*tmax]);
+else
+    xlabel('Time [s]');
+    xlim([tmin tmax]);
+end
+ylabel('Trials');
+ylim([0.3 ntrials+0.7 ]);
+end
+