Merge branch 'master' of whale:scientificComputing

projects/project_onset_fi/solution/boltzmannFit.m

@@ -0,0 +1,83 @@
+function [best_params] = boltzmannFit(contrasts, onset_fi)
+% use a gradient descent approach to fit the fi curve owith a Boltzmann
+% sigmoidal function
+%
+% Arguments:
+% constrasts: the x-values (must be sorted ascendingly)
+% onset_fi: the measured neuronal responses
+%
+% Returns the parameters leading to the best fit.
+
+% we will make some educated guesses for the start parameters; minimum:
+% minimum fi response; saturation: maximum fi response; inflection: zero
+% contrast; slope: response range/contrast range
+start_params = [min(onset_fi), max(onset_fi), 0, ...
+    (max(onset_fi) - min(onset_fi))/(max(contrasts) - min(contrasts))];
+
+best_params = gradientDescent(start_params, contrasts, onset_fi);
+
+
+%% subfunctions that could be separate m-files, but are only used in this context
+
+    function bp = gradientDescent(p, x, y)
+    % performs the gradient descent. Input arguments: 
+    % p: the list of starting parameters for the search
+    % x: the x-values
+    % y: the measured y-values
+    %
+    % returns the parameters leading to the best fit
+        
+       delta = 0.001;
+       gradient = getGradient(p, x, y, delta);
+       eps = [0.001 0.001, 0.00001, 0.00001]; % key point; use smaller changes for the more delicate parameters
+       while norm(gradient) > 0.5
+           p = p - eps .* gradient;
+           gradient = getGradient(p, x, y, delta);
+       end  
+       bp = p;
+    end
+
+    function g = getGradient(p, x, y, delta)
+    % calculate the gradient in the error surface for a small step
+    % along all parameters
+    % p: the parameter list
+    % x: the x-values
+    % y: the measured y-values
+    % delta: the step used or the gradient
+    %
+    % returns the gradient
+        g = zeros(size(p));
+        
+        for i = 1:length(p)
+            np = p;
+            np(i) = np(i) + delta;
+            g(i) = (costFunction(np, x, y) - costFunction(p, x, y)) / delta;
+        end
+    end
+
+    function cost = costFunction(p, x, y)
+        % cost function that relates the model to the data. Applies a
+        % Boltzman model and calculates the error in a mean square error
+        % sense.
+        % p: the set of parameters used for the model
+        % x: the x-values
+        % y: the measured y-values
+        %
+        % returns the costs related to that parameterization
+        
+        y_est = boltzmannModel(p, x);
+        cost = msqError(y, y_est);
+    end
+
+    function e = msqError(y, y_est)
+    % calculates the deviation between y and y_est in the mean square
+    % error sense
+    % y and y_est must have the same size
+    %
+    % returns the error
+        e = mean((y - y_est).^2);
+    end
+
+    
+
+end

projects/project_onset_fi/solution/boltzmannModel.m

@@ -0,0 +1,10 @@
+function y = boltzmannModel(p, x)
+% computes the Boltmannfunction given the parameters given in p
+% the following order is assumed
+% p(1) = minimum y-value
+% p(2) = maximum y-value, i.e. the saturation
+% p(3) = the position of the inflection point
+% p(4) = the slope at the inflection
+
+y = ((p(1) - p(2))) ./ (1 + exp((x-p(3))/p(4))) + p(2);
+end

projects/project_onset_fi/solution/convolutionRate.m

@@ -0,0 +1,27 @@
+function [time, rate] = convolution_rate(spikes, sigma, dt, t_max)
+% PSTH computed with convolution method. 
+%
+% [time, rate] = convolution_rate(spikes, sigma, dt, t_max)
+%
+% Arguments: 
+%   spikes: a vector containing the spike times.
+%   sigma : the standard deviation of the Gaussian kernel in seconds.
+%   dt    : the temporal resolution in seconds.
+%   t_max : the trial duration in seconds.
+%
+% Returns:    two vectors containing the time and the rate.
+
+time = 0:dt:t_max - dt;
+rate = zeros(size(time));
+spike_indices = round(spikes / dt);
+rate(spike_indices) = 1;
+kernel = gaussKernel(sigma, dt);
+
+rate = conv(rate, kernel, 'same');
+end
+
+
+function y = gaussKernel(s, step)
+    x = -4 * s:step:4 * s;
+    y = exp(-0.5 .* (x ./ s) .^ 2) ./ sqrt(2 * pi) / s;
+end

projects/project_onset_fi/solution/getFICurve.m

@@ -0,0 +1,12 @@
+function [fi, fi_std, sort_contrasts] = getFICurve(firing_rates, time, duration, contrasts)
+
+[sort_contrasts, sort_idx] = sort(contrasts);
+fi = zeros(length(sort_contrasts), 1);
+fi_std = zeros(length(sort_contrasts), 1); 
+
+for i = 1:length(sort_contrasts)
+    responses = firing_rates{sort_idx(i)};
+    onset_responses = mean(responses((time > 0) & (time <= duration), :),1);
+    fi(i) = mean(onset_responses);
+    fi_std(i) = std(onset_responses); 
+end

projects/project_onset_fi/solution/getFiringRates.m

@@ -0,0 +1,10 @@
+function [time, rates] = getFiringRates(all_spikes, min_t, max_t, dt)
+
+
+time = min_t:dt:max_t-dt;
+rates = zeros(length(time), length(all_spikes));
+for i = 1:length(all_spikes)
+    spike_times = all_spikes{i}/1000 - min_t;
+    spike_times(spike_times == 0.0) = dt;
+    [t, rates(:, i)] = convolutionRate(spike_times, 0.01, dt, max_t-min_t);
+end

projects/project_onset_fi/solution/main.m

@@ -0,0 +1,36 @@
+%% solution for project onset-FI
+
+clear all 
+close all
+
+dt = 1./20000;
+t_min = -0.2;
+t_max = 1.825;
+data_folder = strcat('..', filesep, 'data');
+files = dir(strcat(data_folder, filesep, '*.mat'));
+
+contrasts = zeros(length(files), 1);
+rates = {};
+for i = 1:length(files)
+    data = load(strcat(data_folder, filesep, files(i).name));
+    contrasts(i) = data.contrast;
+    spikes_times = data.spike_times;
+    [t, rates{i}] = getFiringRates(spikes_times, t_min, t_max, dt);
+end
+
+%% create a plot of the average response for the different contrasts
+plotAverageResponse(rates, t, contrasts, 'averageResponse')
+
+%% extract the onset and the steady-state firing rate
+[fi, fi_std, contr] = getFICurve(rates, t, 0.015, contrasts);
+
+%% plot FI curve 
+plotFICurve(fi, fi_std, contr, 'ficurve')
+
+%% fit FICurve with Boltzmann function
+best_params = boltzmannFit(contr, fi);
+
+% plot the fi curve with the fit
+plotFICurveFit(fi, fi_std, contr, best_params, 'ficurvefit')
+
+

projects/project_onset_fi/solution/plotAverageResponse.m

@@ -0,0 +1,47 @@
+function [] = plotAverageResponse(rates, time, contrasts, filename)
+[sort_contrasts, sort_idx] = sort(contrasts);
+
+figure
+set(gcf, 'paperunits', 'centimeters', 'papersize', [20 20], ...
+    'paperposition', [0.0 0.0 20, 20], 'color', 'white')
+axes = cell(size(sort_contrasts));
+for i = 1:length(sort_contrasts)
+    subplot(ceil(length(sort_contrasts)/2), 2, i);
+    responses = rates{sort_idx(i)};
+    avg_response = mean(responses, 2)';
+    std_response = std(responses, [], 2)';
+    f = fill([time fliplr(time)], [avg_response + std_response fliplr(avg_response - std_response)], ...
+            'b', 'EdgeColor','none', 'displayname', 'std');
+    alpha(f, 0.25)
+    hold on
+    plot(time, avg_response, 'b', 'displayname', 'average response')
+    text(1.2, 500, sprintf('contrast: %i %', sort_contrasts(i)), ...
+        'fontsize', 9, 'FontWeight', 'bold')
+    xlim([time(1) time(end)])
+    ylim([0, 600])
+    set(gca, 'XTick', time(1):0.2:time(end), 'YTick', 0:200:600)
+    if mod(i,2) == 1
+        ylabel('firing rate [Hz]', 'fontsize', 11);
+    else
+        set(gca, 'yticklabels', []);
+    end
+    if i > length(sort_contrasts) - 2 
+        xlabel('time [s]', 'fontsize', 11)
+    else
+        set(gca, 'xticklabels', []);
+    end
+   
+    box('off')
+    axes{i} = gca();
+    if i == 1
+        l = legend('show');
+        l.FontSize = 9;
+        l.Box = 'off';
+        l.Orientation = 'horizontal';
+        l.Position = [0.25, 0.96, .5, 0.01];
+    end
+end
+
+axes{1}.Position = [axes{3}.Position(1) axes{2}.Position(2) ... 
+    axes{2}.Position(3) axes{3}.Position(4)];
+saveas(gcf, filename, 'pdf')

projects/project_onset_fi/solution/plotFICurve.m

@@ -0,0 +1,16 @@
+function plotFICurve(fi, fi_std, contrasts, filename)
+
+figure()
+set(gcf, 'paperunits', 'centimeters', 'papersize', [10 10], ...
+    'paperposition', [0.0 0.0 10, 10], 'color', 'white')
+
+errorbar(contrasts, fi, fi_std, 'o', 'displayname', 'onset response')
+xlim([-20, 20])
+ylim([0, 600])
+xlabel('stimulus contrast [%]', 'fontsize', 11)
+ylabel('firing rate [Hz]', 'fontsize', 11)
+box('off')
+
+if ~isempty(filename) 
+    saveas(gcf, filename, 'pdf')
+end

projects/project_onset_fi/solution/plotFICurveFit.m

@@ -0,0 +1,8 @@
+function plotFICurveFit(fi, fi_error, contrasts, fit_parameter, filename)
+
+plotFICurve(fi, fi_error, contrasts, '')
+fit = boltzmannModel(fit_parameter, min(contrasts):0.1:max(contrasts));
+hold on 
+plot(min(contrasts):0.1:max(contrasts), fit, 'displayname', 'Boltzmann fit')
+legend('Location', 'northwest')
+saveas(gcf, filename, 'pdf')