Merge branch 'master' of whale.am28.uni-tuebingen.de:scientificComputing

projects/Makefile

@@ -1,4 +1,4 @@
-all: projects evalutation
+all: projects evaluation
 
 evaluation: evaluation.pdf
 evaluation.pdf: evaluation.tex

projects/evaluation.tex

@@ -15,7 +15,7 @@
 \begin{document}
 
 \sffamily
-\section*{Scientific computing WS16/17}
+\section*{Scientific computing WS17/18}
 
 \begin{tabular}{|p{0.15\textwidth}|p{0.07\textwidth}|p{0.07\textwidth}|p{0.07\textwidth}|p{0.07\textwidth}|p{0.07\textwidth}|p{0.07\textwidth}|p{0.07\textwidth}|p{0.07\textwidth}|p{0.07\textwidth}|}
 \hline

projects/project_lif/lif.tex

@@ -58,7 +58,7 @@ time = [0.0:dt:tmax]; % t_i
     \part Response of the passive membrane to a step input. 
 
     Set $V_0=0$. Construct a vector for the input $E(t)$ such that
-    $E(t)=0$ for $t\le 20$\,ms and $t\ge 70$\,ms and $E(t)=10$\,mV for
+    $E(t)=0$ for $t\le 20$\,ms or $t\ge 70$\,ms, and $E(t)=10$\,mV for
     $20$\,ms $<t<70$\,ms. Plot $E(t)$ and the resulting $V(t)$ for
     $t_{max}=120$\,ms.
 
@@ -92,7 +92,7 @@ time = [0.0:dt:tmax]; % t_i
     spike'' only means that we note down the time of the threshold
     crossing as a time where an action potential occurred. The
     waveform of the action potential is not modeled. Here we use a
-    voltage threshold of one.
+    voltage threshold of 1\,mV.
 
     Write a function that implements this leaky integrate-and-fire
     neuron by expanding the function for the passive neuron
@@ -114,7 +114,6 @@ time = [0.0:dt:tmax]; % t_i
       \label{firingrate}
       r = \frac{n-1}{t_n - t_1}
     \end{equation}
-
     What do you observe? Does the firing rate encode the frequency of
     the stimulus?
   \end{parts}

projects/project_noiseficurves/lifspikes.m

@@ -1,8 +1,8 @@
-function spikes = lifspikes(trials, input, tmax, D)
+function spikes = lifspikes(trials, current, tmax, D)
 % Generate spike times of a leaky integrate-and-fire neuron
 % trials: the number of trials to be generated
-% input: the stimulus either as a single value or as a vector
-% tmax: duration of a trial
+% current: the stimulus either as a single value or as a vector
+% tmax: duration of a trial if input is a single number
 % D: the strength of additive white noise
     
     tau = 0.01;
@@ -13,7 +13,11 @@ function spikes = lifspikes(trials, input, tmax, D)
     vthresh = 10.0;
     dt = 1e-4;
     
-    n = ceil(tmax/dt);
+    n = length( current );
+    if n <= 1
+        n = ceil(tmax/dt);
+        current = zeros( n, 1 ) + current;
+    end
     spikes = cell(trials, 1);
     for k=1:trials
         times = [];
@@ -21,7 +25,7 @@ function spikes = lifspikes(trials, input, tmax, D)
         v = vreset + (vthresh-vreset)*rand();
         noise = sqrt(2.0*D)*randn(n, 1)/sqrt(dt);
         for i=1:n
-            v = v + (- v + noise(i) + input)*dt/tau;
+            v = v + (- v + noise(i) + current(i))*dt/tau;
             if v >= vthresh
                 v = vreset;
                 times(j) = i*dt;

projects/project_noiseficurves/noiseficurves.tex

@@ -9,9 +9,6 @@
 
 \input{../instructions.tex}
 
-
-%%%%%%%%%%%%%% Questions %%%%%%%%%%%%%%%%%%%%%%%%%
-%\section{REPLACE BY SUBTHRESHOLD RESONANCE PROJECT!}
 \begin{questions}
   \question You are recording the activity of a neuron in response to
   constant stimuli of intensity $I$ (think of that, for example,
@@ -19,24 +16,32 @@
 
   Measure the tuning curve (also called the intensity-response curve) of the
   neuron.  That is, what is the mean firing rate of the neuron's response
-  as a function of the input $I$?
+  as a function of the constant input current $I$?
 
   How does the intensity-response curve of a neuron depend on the
   level of the intrinsic noise of the neuron?
 
+  How can intrinsic noise be usefull for encoding stimuli?
+
   The neuron is implemented in the file \texttt{lifspikes.m}.  Call it
-  with the following parameters:
+  with the following parameters:\\[-7ex]
     \begin{lstlisting}
 trials = 10;
 tmax = 50.0;
-input = 10.0;  % the input I
-Dnoise = 1.0;  % noise strength
-spikes = lifspikes(trials, input, tmax, Dnoise);
+current = 10.0;  % the constant input current I
+Dnoise = 1.0;    % noise strength
+spikes = lifspikes(trials, current, tmax, Dnoise);
     \end{lstlisting}
     The returned \texttt{spikes} is a cell array with \texttt{trials}
     elements, each being a vector of spike times (in seconds) computed
-    for a duration of \texttt{tmax} seconds.  The input is set via the
-    \texttt{input} variable, the noise strength via \texttt{Dnoise}.
+    for a duration of \texttt{tmax} seconds.  The input current is set
+    via the \texttt{current} variable, the strength of the intrinsic
+    noise via \texttt{Dnoise}. If \texttt{current} is a single number,
+    then an input current of that intensity is simulated for
+    \texttt{tmax} seconds. Alternatively, \texttt{current} can be a
+    vector containing an input current that changes in time. In this
+    case, \texttt{tmax} is ignored, and you have to provide a value
+    for the input current for every 0.0001\,seconds.
 
     Think of calling the \texttt{lifspikes()} function as a simple way
     of doing an electrophysiological experiment. You are presenting a
@@ -52,20 +57,17 @@ spikes = lifspikes(trials, input, tmax, Dnoise);
     and plot neuron's $f$-$I$ curve, i.e. the mean firing rate (number
     of spikes within the recording time \texttt{tmax} divided by
     \texttt{tmax} and averaged over trials) as a function of the input
-    for inputs ranging from 0 to 20.
+    current for inputs ranging from 0 to 20.
 
     How are different stimulus intensities encoded by the firing rate
     of this neuron?
 
     \part Compute the $f$-$I$ curves of neurons with various noise
-    strengths \texttt{Dnoise}. Use $D_{noise} = 1e-3$, $1e-2$, and
-    $1e-1$. 
+    strengths \texttt{Dnoise}. Use for example $D_{noise} = 1e-3$,
+    $1e-2$, and $1e-1$.
 
     How does the intrinsic noise influence the response curve?
 
-    How is the encoding of stimuli influenced by increasing intrinsic
-    noise?
-
     What are possible sources of this intrinsic noise?
 
     \part Show spike raster plots and interspike interval histograms
@@ -74,14 +76,21 @@ spikes = lifspikes(trials, input, tmax, Dnoise);
     responses of the four different neurons to the same input, or by
     the same resulting mean firing rate.
 
-    \part How does the coefficient of variation $CV_{isi}$ (standard
-    deviation divided by mean) of the interspike intervalls depend on
-    the input and the noise level?
+    \part Let's now use as an input to the neuron a 1\,s long sine
+    wave $I(t) = I_0 + A \sin(2\pi f t)$ with offset current $I_0$,
+    amplitude $A$, and frequency $f$. Set $I_0=5$, $A=4$, and
+    $f=5$\,Hz.
+
+    Do you get a response of the noiseless ($D_{noise}=0$) neuron?
+
+    What happens if you increase the noise strength?
+
+    What happens at really large noise strengths?
+
+    Generate some example plots that illustrate your findings.
 
-    \part Based o your results, discuss how intrinsic noise might
-    improve and how it might deteriote the encoding of different
-    stimulus intensities.
-    
+    Explain the encoding of the sine wave based on your findings
+    regarding the $f$-$I$ curves.
 
  \end{parts}
 

projects/project_noiseficurves/solution/lifspikes.m

@@ -1,8 +1,8 @@
-function spikes = lifspikes(trials, input, tmax, D)
+function spikes = lifspikes(trials, current, tmax, D)
 % Generate spike times of a leaky integrate-and-fire neuron
 % trials: the number of trials to be generated
-% input: the stimulus either as a single value or as a vector
-% tmax: duration of a trial
+% current: the stimulus either as a single value or as a vector
+% tmax: duration of a trial if input is a single number
 % D: the strength of additive white noise
     
     tau = 0.01;
@@ -13,7 +13,11 @@ function spikes = lifspikes(trials, input, tmax, D)
     vthresh = 10.0;
     dt = 1e-4;
     
-    n = ceil(tmax/dt);
+    n = length( current );
+    if n <= 1
+        n = ceil(tmax/dt);
+        current = zeros( n, 1 ) + current;
+    end
     spikes = cell(trials, 1);
     for k=1:trials
         times = [];
@@ -21,7 +25,7 @@ function spikes = lifspikes(trials, input, tmax, D)
         v = vreset + (vthresh-vreset)*rand();
         noise = sqrt(2.0*D)*randn(n, 1)/sqrt(dt);
         for i=1:n
-            v = v + (- v + noise(i) + input)*dt/tau;
+            v = v + (- v + noise(i) + current(i))*dt/tau;
             if v >= vthresh
                 v = vreset;
                 times(j) = i*dt;

projects/project_noiseficurves/solution/noiseficurves.m

@@ -7,20 +7,34 @@ figure()
 Ds = [0, 0.001, 0.01, 0.1];
 for j = 1:length(Ds)
 	D = Ds(j);
-	inputs = 0.0:0.5:20.0;
-	rates = ficurve(trials, inputs, tmax, D);
-	plot(inputs, rates);
+	currents = 0.0:0.5:20.0;
+	rates = ficurve(trials, currents, tmax, D);
+	plot(currents, rates);
 	hold on;
 end
 hold off;
 
 %% spike raster and CVs
-input = 12.0;
+figure()
+current = 12.0;
 for j = 1:length(Ds)
 	D = Ds(j);
-	spikes = lifspikes(trials, input, tmax, D);
+	spikes = lifspikes(trials, current, tmax, D);
 	subplot(4, 2, 2*j-1);
 	spikeraster(spikes, 0.0, 1.0);
 	subplot(4, 2, 2*j);
 	isih(spikes, [0:0.001:0.04]);
 end
+
+%% subthreshold resonance:
+time = [0.0:0.0001:1.0];
+current = 5.0 + 4.0*sin(2.0*pi*5.0*time);
+D = 0.1;
+spikes = lifspikes(trials, current, tmax, D);
+subplot(2, 1, 1);
+spikeraster(spikes, 0.0, 1.0);
+subplot(2, 1, 2);
+plot(time, current);
+
+
+

projects/project_populationvector/code/generatedata.py

@@ -69,15 +69,15 @@ plt.show()
 
 # tuning curves:
 nunits = 6
-unitphases = np.linspace(0.0, 1.0, nunits) + 0.05*np.random.randn(nunits)/float(nunits)
+unitphases = np.linspace(0.04, 0.96, nunits) + 0.05*np.random.randn(nunits)/float(nunits)
 unitgains = 15.0 + 5.0*(2.0*np.random.rand(nunits)-1.0)
 nangles = 12
 angles = 180.0*np.arange(nangles)/nangles
 for unit, (phase, gain) in enumerate(zip(unitphases, unitgains)):
-    print '%.1f %.0f' % (gain, phase*180.0)
+    print 'gain=%.1f phase=%.0f' % (gain, phase*180.0)
     allspikes = []
     for k, angle in enumerate(angles):
-        spikes = lifadaptspikes(0.5*(1.0-np.cos(2.0*np.pi*(angle/180.0-phase))), gain)
+        spikes = lifadaptspikes(0.5*(1.0+np.cos(2.0*np.pi*(angle/180.0-phase))), gain)
         allspikes.append(spikes)
     spikesobj = np.zeros((len(allspikes), len(allspikes[0])), dtype=np.object)
     for k in range(len(allspikes)):
@@ -89,10 +89,10 @@ for unit, (phase, gain) in enumerate(zip(unitphases, unitgains)):
 nangles = 50
 angles = 180.0*np.random.rand(nangles)
 for k, angle in enumerate(angles):
-    print '%.0f' % angle
+    print 'angle = %.0f' % angle
     allspikes = []
     for unit, (phase, gain) in enumerate(zip(unitphases, unitgains)):
-        spikes = lifadaptspikes(0.5*(1.0-np.cos(2.0*np.pi*(angle/180.0-phase))), gain)
+        spikes = lifadaptspikes(0.5*(1.0+np.cos(2.0*np.pi*(angle/180.0-phase))), gain)
         allspikes.append(spikes)
     spikesobj = np.zeros((len(allspikes), len(allspikes[0])), dtype=np.object)
     for i in range(len(allspikes)):

projects/project_populationvector/populationvector.tex

@@ -35,8 +35,6 @@
   \texttt{spikes} variables of the \texttt{population*.mat} files.
   The \texttt{angle} variable holds the angle of the presented bar.
 
-%NOTE: the orientation is angle plus 90 degree!!!!!!
-
 \continue
   \begin{parts}
     \part Illustrate the spiking activity of the V1 cells in response
@@ -50,13 +48,11 @@
     of the neurons.
 
     \part Fit the function \[ r(\varphi) =
-    g(1-\cos(\varphi-\varphi_0))/2 \] to the measured tuning curves in
+    g(1+\cos(\varphi-\varphi_0))/2 \] to the measured tuning curves in
     order to estimated the orientation angle at which the neurons
     respond strongest. In this function $\varphi_0$ is the position of
-    the peak (really? How exactly is $\varphi_0$ related to the
-    position of the peak? Do you find a better function where
-    $\varphi_0$ is identical with the peak position?) and $g$ is a
-    gain factor that sets the maximum firing rate.
+    the peak and $g$ is a gain factor that sets the maximum firing
+    rate.
 
     \part How can the orientation angle of the presented bar be read
     out from one trial of the population activity of the 6 neurons?
@@ -70,8 +66,8 @@
     An alternative read out is maximum likelihood (see script).
 
     Load one of the \texttt{population*.mat} files, illustrate the
-    data, and estimate the orientation angle of the bar by the two
-    different methods.
+    data, and estimate the orientation angle of the bar from single
+    trial data by the two different methods.
 
     \part Compare, illustrate and discuss the performance of your two
     decoding methods by using all of the recorded responses (all
@@ -81,16 +77,16 @@
   \end{parts}
 \end{questions}
 
-%NOTE: change data generation such that the phase variable is indeed
-%the maximum response and not the minumum!
-
 \end{document}
 
+
 gains and angles of the 6 neurons:
 
-14.6 0
-17.1 36
-17.6 72
-14.1 107
-10.7 144
-11.4 181
+gain=10.7 phase=5
+gain=18.0 phase=38
+gain=11.3 phase=71
+gain=14.1 phase=108
+gain=19.0 phase=138
+gain=16.4 phase=174
+
+

projects/project_populationvector/solution/populationvector.m

@@ -1,7 +1,9 @@
 close all
-files = dir('../unit*.mat');
+datapath = '../';
+% datapath = '../code/';
+files = dir(strcat(datapath, 'unit*.mat'));
 for file = files'
-	a = load(strcat('../', file.name));
+	a = load(strcat(datapath, file.name));
 	spikes = a.spikes;
 	angles = a.angles;
 	figure()
@@ -14,13 +16,13 @@ end
 
 %% tuning curves:
 close all
-cosine = @(p,xdata)0.5*p(1).*(1.0-cos(2.0*pi*(xdata/180.0-p(2))));
-files = dir('../unit*.mat');
+cosine = @(p,xdata)0.5*p(1).*(1.0+cos(2.0*pi*(xdata/180.0-p(2))));
+files = dir(strcat(datapath, 'unit*.mat'));
 phases = zeros(length(files), 1);
 figure()
 for j = 1:length(files)
 	file = files(j);
-	a = load(strcat('../', file.name));
+	a = load(strcat(datapath, file.name));
 	spikes = a.spikes;
 	angles = a.angles;
 	rates = zeros(size(spikes, 1), 1);
@@ -32,10 +34,13 @@ for j = 1:length(files)
 	p0 = [mr, angles(maxi)/180.0-0.5];
 	%p = p0;
 	p = lsqcurvefit(cosine, p0, angles, rates');
-	phase = p(2)*180.0 + 90.0;
+	phase = p(2)*180.0;
     if phase > 180.0
         phase = phase - 180.0;
     end
+    if phase < 0.0
+        phase = phase + 180.0;
+    end
     phases(j) = phase;
 	subplot(2, 3, j);
 	plot(angles, rates, 'b');
@@ -49,40 +54,45 @@ end
 
 
 %% read out:
-a = load('../population04.mat');
+a = load(strcat(datapath, 'population04.mat'));
 spikes = a.spikes;
 angle = a.angle;
-unitphases = a.phases*180.0 + 90.0;
+unitphases = a.phases*180.0;
 unitphases(unitphases>180.0) = unitphases(unitphases>180.0) - 180.0;
 figure();
 subplot(1, 3, 1);
 angleestimates1 = zeros(size(spikes, 2), 1);
 angleestimates2 = zeros(size(spikes, 2), 1);
+[x, inx] = sort(phases);
+% loop over trials:
 for j = 1:size(spikes, 2)
     rates = zeros(size(spikes, 1), 1);
     for k = 1:size(spikes, 1)
         r = firingrate(spikes(k, j), 0.0, 0.2);
         rates(k) = r;
     end
-    [x, inx] = sort(phases);
     plot(phases(inx), rates(inx), '-o');
     hold on;
     angleestimates1(j) = popvecangle(phases, rates);
     [m, i] = max(rates);
     angleestimates2(j) = phases(i);
 end
+xlabel('preferred angle')
+ylabel('firing rate')
 hold off;
 subplot(1, 3, 2);
 hist(angleestimates1);
+xlabel('stimulus angle')
 subplot(1, 3, 3);
 hist(angleestimates2);
+xlabel('stimulus angle')
 angle
 mean(angleestimates1)
 mean(angleestimates2)
 
 
 %% read out robustness:
-files = dir('../population*.mat');
+files = dir(strcat(datapath, 'population*.mat'));
 angles = zeros(length(files), 1);
 e1m = zeros(length(files), 1);
 e1s = zeros(length(files), 1);
@@ -90,7 +100,7 @@ e2m = zeros(length(files), 1);
 e2s = zeros(length(files), 1);
 for i = 1:length(files)
 	file = files(i);
-	a = load(strcat('../', file.name));
+	a = load(strcat(datapath, file.name));
     spikes = a.spikes;
     angle = a.angle;
     angleestimates1 = zeros(size(spikes, 2), 1);
@@ -114,5 +124,9 @@ end
 figure();
 subplot(1, 2, 1);
 scatter(angles, e1m);
+xlabel('stimuluis angle')
+ylabel('estimated angle')
 subplot(1, 2, 2);
 scatter(angles, e2m);
+xlabel('stimuluis angle')
+ylabel('estimated angle')