clean up and commenting

cont_ca_sampler.m

@@ -1,18 +1,17 @@
 function SAMPLES = cont_ca_sampler(Y,params)
 
-% add the capability of dealing with missing numbers
+% MCMC continuous time sampler of spiketimes given a fluorescence trace Y
 
-% Continuous time sampler
+% Inputs:
 % Y                     data
-% P                     intialization parameters (discrete time constant P.g required)
 % params                parameters structure
 % params.g              discrete time constant(s) (estimated if not provided)
 % params.sn             initializer for noise (estimated if not provided)
 % params.b              initializer for baseline (estimated if not provided)
 % params.c1             initializer for initial concentration (estimated if not provided)
 % params.Nsamples       number of samples after burn in (default 500)
 % params.B              number of burn in samples (default 200)
-% params.marg           flag for marginalized sampler (default 0)
+% params.marg           flag for marginalized sampler for baseline and initial concentration (default 0)
 % params.upd_gam        flag for updating gamma (default 1)
 % params.gam_step       number of samples after which gamma is updated (default 50)
 % params.std_move       standard deviation of shifting kernel (default 3*Dt)
@@ -45,6 +44,9 @@
 
 % Author: Eftychios A. Pnevmatikakis and Josh Merel
 
+% Reference: Pnevmatikakis et al., Bayesian spike inference from calcium
+%                imaging data, Asilomar 2013.
+
 Y = Y(:);
 T = length(Y);
 isanY = ~isnan(Y);                          % Deal with possible missing entries

plot_continuous_samples.m

@@ -1,5 +1,18 @@
 function plot_continuous_samples(SAMPLES,Y)
 
+% plot results of MCMC sampler
+% The mean calcium sample, spike sampler raster plot and samples for
+% amplitude, number of spikes, discrete time constants, noise variance, 
+% baseline and initial concentration are generated, together with their 
+% autocorrelation functions. If the marginalized flag was used, then the 
+% posterior pdfs of baseline and initial concentration are plotted.
+
+% Inputs:
+% SAMPLES:  structure with SAMPLES obtained from cont_ca_sampler.m
+% Y:        inpurt fluorescence trace
+
+% Author: Eftychios A. Pnevmatikakis, 2016, Simons Foundation
+
 T = length(Y);
 N = length(SAMPLES.ns);
 show_gamma = 1;
@@ -33,7 +46,6 @@ function plot_continuous_samples(SAMPLES,Y)
 end
 
 
-
 if length(SAMPLES.Cb) == 2
     marg = 1;       % marginalized sampler
 else

sampling_demo_ar2.m

@@ -1,6 +1,8 @@
-clear all;
+% creating a demo for the MCMC sampler
+
+clearvars;
 addpath utilities
-dt = 1e-1;
+dt = 1e-1;   % Data is generated with provided dt but sampled with Dt = 1
 T = 7000;
 ld = 0.05;   % rate spikes per second
 
@@ -10,22 +12,23 @@
 hmax = tau_decay/(tau_decay+tau_rise)*(tau_rise/(tau_decay+tau_rise))^(tau_rise/tau_decay);
 [g,h1] = tau_c2d(tau_rise,tau_decay,dt);
 
-b = hmax/4;
-cin = [.2*b,.15*b];
-c = [cin,filter(h1,[1,-g],s(3:end),filtic(h1,[1,-g],cin))] + b;
-c_true = c(round(1/dt):round(1/dt):round(T/dt));
-sg = hmax/4;        % noise level
-y = c_true + sg*randn(1,length(c_true));
+b = hmax/4;                                                         % baseline value
+cin = [.2*b,.15*b];                                                 % initial concentration
+c = [cin,filter(h1,[1,-g],s(3:end),filtic(h1,[1,-g],cin))] + b;     % true calcium at original resolution
+c_true = c(round(1/dt):round(1/dt):round(T/dt));                    % true calcium at sampled resolution
+sg = hmax/4;                                                        % noise level
+y = c_true + sg*randn(1,length(c_true));                            % noisy observations
 
 figure;subplot(2,1,1); plot(dt:dt:T,c); hold all; scatter(1:T,y,'r*'); 
                        legend('True Calcium','Observed Values');
-%%  constrained foopsi
-[g2,h2] = tau_c2d(tau_rise,tau_decay,1);
-P.f = 1;
-P.g = g2;
-P.sn = sg;
-params.marg = 0;
+
+%%  Run constrained deconvolution approach (constrained foopsi)
+
+[g2,h2] = tau_c2d(tau_rise,tau_decay,1);   % generate discrete time constants
+
 [ca_foopsi,cb,c1,~,~,spikes_foopsi] = constrained_foopsi(y,[],[],g2);
+
+% do some plotting
 spiketimes{1} =  find(s)*dt;
 spikeRaster = samples_cell2mat(spiketimes,T,1);
 f = find(spikeRaster);
@@ -44,7 +47,8 @@
     title('Foopsi Spikes','FontWeight','bold','Fontsize',14); xlabel('Timestep','FontWeight','bold','Fontsize',16);
     legend('Foopsi Spikes','Ground Truth');
     drawnow;
-%% MCMC 
+
+%% run continuous time MCMC sampler
 
 params.p = 2;
 params.g = g2;
@@ -54,5 +58,7 @@
 params.c1 = c1;
 params.sn = sg;
 params.marg = 0;
-SAMPLES2 = cont_ca_sampler(y,params);    %% MCMC   
-plot_continuous_samples(SAMPLES2,y(:));
+
+SAMPLES = cont_ca_sampler(y,params);    %% MCMC   
+plot_continuous_samples(SAMPLES,y(:));
+M = plot_marginals(SAMPLES.ss,T,spikeRaster);

utilities/addSpike.m

@@ -1,12 +1,33 @@
-function [newSpikeTrain, newCalcium, newLL] = addSpike(oldSpikeTrain,oldCalcium,oldLL,filter,tau,obsCalcium,timeToAdd,indx,Dt,A)
+function [newSpikeTrain, newCalcium, newLL] = addSpike(oldSpikeTrain,oldCalcium,oldLL,filters,tau,obsCalcium,timeToAdd,indx,Dt,A)
+
+% Add a given spike to the existing spike train.
+
+% Inputs:
+% oldSpikeTrain:        current spike train
+% oldCalcium:           current noiseless calcium trace
+% oldLL:                current value of the log-likelihood function
+% filters:              exponential rise and decay kernels for calcium transient
+% tau:                  continuous time rise and decay time constants
+% obsCalcium:           observed fluorescence trace
+% timetoAdd:            time of the spike to be added
+% indx:                 place where the new spike is added in the existing spike train vector
+% Dt:                   time-bin width
+% A:                    spike amplitude
+
+% Outputs:
+% newSpikeTrain:        new vector of spike times
+% newCalcium:           new noiseless calcium trace
+% newLL:                new value of the log-likelihood function
+
+% Author: Eftychios A. Pnevmatikakis and Josh Merel
 
     tau_h = tau(1);
     tau_d = tau(2);
 
-    ef_h = filter{1,1};
-    ef_d = filter{1,2};
-    ef_nh = filter{2,1};
-    ef_nd = filter{2,2};
+    ef_h = filters{1,1};
+    ef_d = filters{1,2};
+    ef_nh = filters{2,1};
+    ef_nd = filters{2,2};
 
     newSpikeTrain = [oldSpikeTrain(1:indx-1) timeToAdd oldSpikeTrain(indx:end)]; %possibly inefficient, change if problematic (only likely to be a problem for large numbers of spikes)
 

utilities/get_initial_sample.m

@@ -1,7 +1,16 @@
 function SAM = get_initial_sample(Y,params)
 
 % obtain initial sample by performing sparse noise-constrained deconvolution
-% Author: Eftychios A. Pnevmatikakis
+
+% Inputs:
+% Y:        observed fluorescence trace
+% params:   parameter struct (results of constrained foopsi can be passed here to avoid unnecessary computation)
+
+% Output:
+% SAM:      sample structure with values for spikes in continuous time,
+%               amplitude, baseline, noise, initial concentration and time constants
+
+% Author: Eftychios A. Pnevmatikakis, 2016, Simons Foundation
 
 
 if isfield(params,'p'); options.p = params.p; else options.p = 1; end