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

projects/project_input_resistance/input_resistance.tex

@@ -0,0 +1,77 @@
+\documentclass[addpoints,11pt]{exam}
+\usepackage{url}
+\usepackage{color}
+\usepackage{hyperref}
+
+\pagestyle{headandfoot}
+\runningheadrule
+\firstpageheadrule
+\firstpageheader{Scientific Computing}{Project Assignment}{WS 2016/17}
+%\runningheader{Homework 01}{Page \thepage\ of \numpages}{23. October 2014}
+\firstpagefooter{}{}{{\bf Supervisor:} Jan Grewe}
+\runningfooter{}{}{}
+\pointsinmargin
+\bracketedpoints
+
+%\printanswers
+%\shadedsolutions
+
+
+\begin{document}
+%%%%%%%%%%%%%%%%%%%%% Submission instructions %%%%%%%%%%%%%%%%%%%%%%%%%
+\sffamily
+% \begin{flushright}
+% \gradetable[h][questions]
+% \end{flushright}
+
+\begin{center}
+  \input{../disclaimer.tex}
+\end{center}
+
+%%%%%%%%%%%%%% Questions %%%%%%%%%%%%%%%%%%%%%%%%%
+\section*{Estimating cellular properties of different cell types.}
+You will analyze data from intracellular \textit{in vitro} recordings
+of pyramidal neurons from two different maps of the electrosensory
+lateral line lobe (ELL) of the weakly electric fish
+\textit{Apteronotus leptorhynchus}. The membrane resistance and the
+membrane capacitance are fundamental properties of a neuron that have
+a great influence on the coding properties of the cell. They are
+typically estimated by injecting pulses of hyperpolarizing current
+into the cell. From the respective responses we can calculate the
+membrane resistance by applying Ohm's law ($U = R \cdot I$). To
+estimate the membrane capacitance we need to fit an exponential
+function of the form $y = a \cdot e^{(-x/\tau)}$ to the response to get the
+membrane time-constant $\tau$. With the knowledge of $R$ and $\tau$ we
+can estimate the capacitance $C$ from the simple relation $\tau = R
+\cdot C$.
+
+\begin{questions}
+  \question{} The accompanying dataset (input\_resistance.zip)
+  contains datasets from cells originating from two different parts of
+  the ELL, the medial segment (MS) and the centro-medial segment
+  (CMS). Each mat-file contains four variables. (i) \textit{V} the
+  average membrane potential of 20 repeated current injections, (ii)
+  \textit{V\_std} the across-trial standard deviation of the
+  responses, (iii) \textit{t} a vector representing the recording
+  time (in ms), and (iv) \textit{I} a vector containing the time-course of the
+  injected current.
+  
+  \begin{parts}
+    \part{} Create plots of the raw data. Plot the average response as
+    a function of time. This plot should also show the across-trial
+    variability. Also plot the time-course of the injected
+    current. \\[0.5ex]
+    \part{} Estimate the input resistances of each cell.\\[0.5ex]
+    \part{} Fit an exponential to the initial few milliseconds of the
+    current-on response. Use a gradient-descent approach to do
+    this.\\ It is very important to understand the exponential decay
+    function. If you are unsure, play with the function and understand
+    how the parameters influence the decay. (Hint: It might be
+    necessary to transform the data a bit.)\\[0.5ex]
+    \part{} Estimate the membrane capacitance of each cell. Compare
+    $R$, $I$ and $\tau$ between cells of the two segments.\\[0.5ex]
+    \part{} Optional: use a double exponential and see, if the fit improves.
+  \end{parts}
+\end{questions}
+
+\end{document}

projects/project_input_resistance/solution/exponentialDecay.m

@@ -0,0 +1,5 @@
+function y = exponentialDecay(p, x)
+
+  y = p(1) .* exp(-x./p(2));
+
+

projects/project_input_resistance/solution/gradientDescent.m

@@ -0,0 +1,14 @@
+function [position, errors] = gradientDescent(x, y, position)
+
+gradient = [];
+errors = [];
+count = 1;
+eps = 0.02;
+close all
+while isempty(gradient) || norm(gradient) > 0.025
+    gradient = lsqGradient(position, x,y);
+    % disp(gradient)
+    errors(count) = lsqError(position, x, y);
+    position = position - eps .* gradient';
+    count = count + 1;
+end

projects/project_input_resistance/solution/lsqError.m

@@ -0,0 +1,12 @@
+function error = lsqError(parameter, x, y)
+% Objective function for fitting a exponential equation to data.
+%
+% Arguments: parameter, vector containing the parameters for the exp. decay
+%            x, vector of the input values
+%            y, vector of the corresponding measured output values
+%
+% Returns: the estimation error in terms of the mean sqaure error
+
+  y_est = exponentialDecay(parameter, x);
+  error = meanSquareError(y, y_est);
+end

projects/project_input_resistance/solution/lsqGradient.m

@@ -0,0 +1,9 @@
+function gradient = lsqGradient(parameter, x, y)
+
+  h = 1e-1;   % stepsize for derivatives
+  gradient = zeros(length(parameter),1);
+  for i = 1:length(parameter)
+      new_p = parameter;
+      new_p(i) = parameter(i) + h;
+      gradient(i) = (lsqError(new_p, x, y) - lsqError(parameter, x, y))/ h;
+  end

projects/project_input_resistance/solution/meanSquareError.m

@@ -0,0 +1,10 @@
+function error = meanSquareError(y, y_est)
+% Mean squared error between observed and predicted values.
+%
+% Arguments:    y, vector of observed values.
+%               y_est, vector of predicted values.
+%
+% Returns:      the error in the mean-squared-deviation sense.
+
+  error = mean((y - y_est).^2);
+end

projects/project_input_resistance/solution/project_input_resistance.m

@@ -0,0 +1,43 @@
+clear
+maps = {'CMS', 'MS'};
+
+taus = zeros(length(dir('data/resistance/resistance*.mat')), 1);
+cell_types = cell(size(taus));
+count = 1;
+size(taus)
+for j = 1:length(maps)
+    files = dir(strcat('data/resistance/resistance_', maps{j}, '*.mat'));
+    for i = 1:length(files)
+        load(strcat('data/resistance/', files(i).name));
+        disp(files(i).name)
+        x = t(t>0.5 & t<15);
+        x = x - x(1);
+        y = V(t>0.5 & t<15);
+        y = y - y(end);
+        initial_params = [1, 2.5];
+        [params, errors] = gradientDescent(x, y, initial_params);
+        disp(count)
+        taus(count) = params(2);
+        cell_types{count} = maps{j};
+        count = count + 1;
+        plot_fit(x, y, errors, params);
+    end
+end
+
+
+function plot_fit(x, y, errors, parameter)
+    figure()
+    subplot(2,1,1)
+    hold on
+    scatter(x,y, 'displayname', 'data')
+    f_x = exponentialDecay(parameter, x);
+    plot(x, f_x, 'displayname', 'fit')
+    xlabel('Input')
+    ylabel('Output')
+    grid on
+    legend show
+    subplot(2,1,2)
+    plot(errors)
+    xlabel('optimization steps')
+    ylabel('error')
+end