From e4c06940929f1e504aec0a70f54dce1bd4640ddb Mon Sep 17 00:00:00 2001
From: epnev <eap2111@columbia.edu>
Date: Thu, 28 Jan 2016 14:34:16 -0500
Subject: [PATCH] combine update_temporal_compoents & parallel
combine update_temporal_compoents & update_temporal_components_parallel
into a single function. The default option is to use the parallel
implementation if the parallel processing toolbox is present
---
@Sources2D/Sources2D.m | 8 +-
CNMFSetParms.m | 118 ++++++-----
demo_script.m | 4 +-
update_temporal_components.m | 295 ++++++++++++++++++--------
update_temporal_components_parallel.m | 233 +-------------------
5 files changed, 275 insertions(+), 383 deletions(-)
diff --git a/@Sources2D/Sources2D.m b/@Sources2D/Sources2D.m
index a0e54e8..ef48481 100644
--- a/@Sources2D/Sources2D.m
+++ b/@Sources2D/Sources2D.m
@@ -56,13 +56,7 @@ function updateTemporal(obj, Y)
[obj.C, obj.f, obj.P, obj.S] = update_temporal_components(...
Y, obj.A, obj.b, obj.C, obj.f, obj.P, obj.options);
end
-
- %% update temporal components in parallel
- function updateTemporalParallel(obj, Y)
- [obj.C, obj.f, obj.P, obj.S] = update_temporal_components_parallel(...
- Y, obj.A, obj.b, obj.C, obj.f, obj.P, obj.options);
- end
-
+
%% merge found components
function [nr, merged_ROIs] = merge(obj, Y)
[obj.A, obj.C, nr, merged_ROIs, obj.P, obj.S] = merge_components(...
diff --git a/CNMFSetParms.m b/CNMFSetParms.m
index 854c953..bd5d4c4 100644
--- a/CNMFSetParms.m
+++ b/CNMFSetParms.m
@@ -8,76 +8,77 @@
Names = [
% dataset info
- 'd1 ' % number of rows
- 'd2 ' % number of cols
+ 'd1 ' % number of rows
+ 'd2 ' % number of cols
% INITIALIZATION (initialize_components.m)
- 'ssub ' % spatial downsampling factor (default: 1)
- 'tsub ' % temporal downsampling factor (default: 1)
- 'init_method ' % initialization method ('greedy','sparse_NMF') (default: 'greedy')
+ 'ssub ' % spatial downsampling factor (default: 1)
+ 'tsub ' % temporal downsampling factor (default: 1)
+ 'init_method ' % initialization method ('greedy','sparse_NMF') (default: 'greedy')
% greedyROI parameters (greedyROI2d.m)
- 'gSig ' % half size of neurons to be found (default: [5,5])
- 'gSiz ' % half size of bounding box for each neuron (default: 2*gSig+1)
- 'nb ' % number of background components (default: 1)
- 'nIter ' % maximum number of rank-1 NMF iterations during refining
- 'med_app ' % number of timesteps to be interleaved for fast (approximate) median calculation (default: 1)
- 'save_memory ' % process data sequentially to save memory (default: 0)
- 'chunkSiz ' % filter this number of timesteps each time (default: 100)
- 'windowSiz ' % size of window over which is computed sequentially (default: 32 x 32)
+ 'gSig ' % half size of neurons to be found (default: [5,5])
+ 'gSiz ' % half size of bounding box for each neuron (default: 2*gSig+1)
+ 'nb ' % number of background components (default: 1)
+ 'nIter ' % maximum number of rank-1 NMF iterations during refining
+ 'med_app ' % number of timesteps to be interleaved for fast (approximate) median calculation (default: 1)
+ 'save_memory ' % process data sequentially to save memory (default: 0)
+ 'chunkSiz ' % filter this number of timesteps each time (default: 100)
+ 'windowSiz ' % size of window over which is computed sequentially (default: 32 x 32)
% sparse_NMF parameters (sparse_NMF_initialization.m)
- 'snmf_max_iter ' % max # of sparse NMF iterations
- 'err_thr ' % relative change threshold for stopping sparse_NMF
- 'eta ' % frobenious norm factor *max(Y(:))^2
- 'beta ' % sparsity factor
+ 'snmf_max_iter ' % max # of sparse NMF iterations
+ 'err_thr ' % relative change threshold for stopping sparse_NMF
+ 'eta ' % frobenious norm factor *max(Y(:))^2
+ 'beta ' % sparsity factor
% HALS parameters (HALS_2d.m)
- 'bSiz ' % expand kernel for HALS growing (default: 3)
- 'maxIter ' % maximum number of HALS iterations (default: 5)
+ 'bSiz ' % expand kernel for HALS growing (default: 3)
+ 'maxIter ' % maximum number of HALS iterations (default: 5)
% Noise and AR coefficients calculation (preprocess_data.m)
- 'noise_range ' % frequency range over which to estimate the noise (default: [0.25,0.5])
- 'noise_method ' % method for which to estimate the noise level (default: 'logmexp')
- 'flag_g ' % compute global AR coefficients (default: false)
- 'lags ' % number of extra lags when computing the AR coefficients (default: 5)
- 'include_noise ' % include early lags when computing AR coefs (default: 0)
- 'pixels ' % pixels to include when computing the AR coefs (default: 1:numel(Y)/size(Y,ndims(Y)))
- 'split_data ' % split data into patches for memory reasons (default: 0)
- 'block_size ' % block size for estimating noise std in patches (default: [64,64])
+ 'noise_range ' % frequency range over which to estimate the noise (default: [0.25,0.5])
+ 'noise_method ' % method for which to estimate the noise level (default: 'logmexp')
+ 'flag_g ' % compute global AR coefficients (default: false)
+ 'lags ' % number of extra lags when computing the AR coefficients (default: 5)
+ 'include_noise ' % include early lags when computing AR coefs (default: 0)
+ 'pixels ' % pixels to include when computing the AR coefs (default: 1:numel(Y)/size(Y,ndims(Y)))
+ 'split_data ' % split data into patches for memory reasons (default: 0)
+ 'block_size ' % block size for estimating noise std in patches (default: [64,64])
% UPDATING SPATIAL COMPONENTS (unpdate_spatial_components.m)
- 'search_method ' % method for determining footprint of spatial components 'ellipse' or 'dilate' (default: 'ellipse')
- 'use_parallel ' % update pixels in parallel (default: 1 if present)
+ 'search_method ' % method for determining footprint of spatial components 'ellipse' or 'dilate' (default: 'ellipse')
+ 'use_parallel ' % update pixels in parallel (default: 1 if present)
% determine_search_location.m
- 'min_size ' % minimum size of ellipse axis (default: 3)
- 'max_size ' % maximum size of ellipse axis (default: 8)
- 'dist ' % expansion factor of ellipse (default: 3)
- 'se ' % morphological element for dilation (default: strel('disk',4,0))
+ 'min_size ' % minimum size of ellipse axis (default: 3)
+ 'max_size ' % maximum size of ellipse axis (default: 8)
+ 'dist ' % expansion factor of ellipse (default: 3)
+ 'se ' % morphological element for dilation (default: strel('disk',4,0))
% threshold_components.m
- 'nrgthr ' % energy threshold (default: 0.9999)
- 'clos_op ' % morphological element for closing (default: strel('square',3))
- 'medw ' % size of median filter (default: [3,3])
+ 'nrgthr ' % energy threshold (default: 0.9999)
+ 'clos_op ' % morphological element for closing (default: strel('square',3))
+ 'medw ' % size of median filter (default: [3,3])
% UPDATING TEMPORAL COMPONENTS (update_temporal_components.m)
- 'deconv_method ' % method for spike deconvolution (default: 'constrained_foopsi')
- 'restimate_g ' % flag for updating the time constants for each component (default: 1)
- 'temporal_iter ' % number of block-coordinate descent iterations (default 2)
+ 'deconv_method ' % method for spike deconvolution (default: 'constrained_foopsi')
+ 'restimate_g ' % flag for updating the time constants for each component (default: 1)
+ 'temporal_iter ' % number of block-coordinate descent iterations (default: 2)
+ 'temporal_parallel ' % flag for parallel updating of temporal components (default: true if present)
% CONSTRAINED DECONVOLUTION (constrained_foopsi.m)
- 'method ' % methods for performing spike inference ('dual','cvx','spgl1','lars') (default:'cvx')
- 'bas_nonneg ' % flag for setting the baseline lower bound. if 1, then b >= 0 else b >= min(y) (default 0)
- 'fudge_factor ' % scaling constant to reduce bias in the time constant estimation (default 1 - no scaling)
- 'resparse ' % number of times that the solution is resparsened (default: 0)
+ 'method ' % methods for performing spike inference ('dual','cvx','spgl1','lars') (default:'cvx')
+ 'bas_nonneg ' % flag for setting the baseline lower bound. if 1, then b >= 0 else b >= min(y) (default 0)
+ 'fudge_factor ' % scaling constant to reduce bias in the time constant estimation (default 1 - no scaling)
+ 'resparse ' % number of times that the solution is resparsened (default: 0)
% MERGING (merge_ROIs.m)
- 'merge_thr ' % merging threshold (default: 0.85)
- 'fast_merge ' % flag for using fast merging (default 1)
+ 'merge_thr ' % merging threshold (default: 0.85)
+ 'fast_merge ' % flag for using fast merging (default 1)
% VIDEO (make_patch_video.m)
- 'ind ' % indeces of components to be shown (deafult: 1:4)
- 'skip_frame ' % skip frames when showing the video (default: 1 (no skipping))
- 'sx ' % half size of representative patches (default: 16)
- 'make_avi ' % flag for saving avi video (default: 0)
- 'show_background' % flag for displaying the background in the denoised panel (default: 1)
- 'show_contours ' % flag for showing the contour plots of the patches in the FoV (default: 0)
- 'cmap ' % colormap for plotting (default: 'default')
- 'name ' % name of saved video file (default: based on current date)
+ 'ind ' % indeces of components to be shown (deafult: 1:4)
+ 'skip_frame ' % skip frames when showing the video (default: 1 (no skipping))
+ 'sx ' % half size of representative patches (default: 16)
+ 'make_avi ' % flag for saving avi video (default: 0)
+ 'show_background ' % flag for displaying the background in the denoised panel (default: 1)
+ 'show_contours ' % flag for showing the contour plots of the patches in the FoV (default: 0)
+ 'cmap ' % colormap for plotting (default: 'default')
+ 'name ' % name of saved video file (default: based on current date)
% PLOT COMPONENTS (view_patches.m)
- 'plot_df ' % flag for displaying DF/F estimates (default: 1)
- 'make_gif ' % save animation (default: 0)
- 'save_avi ' % save video (default: 0)
- 'pause_time ' % time to pause between each component (default: Inf, user has to click)
+ 'plot_df ' % flag for displaying DF/F estimates (default: 1)
+ 'make_gif ' % save animation (default: 0)
+ 'save_avi ' % save video (default: 0)
+ 'pause_time ' % time to pause between each component (default: Inf, user has to click)
];
[m,n] = size(Names);
@@ -207,6 +208,7 @@
{'constrained_foopsi'}
{1}
{2}
+ {~isempty(which('parpool'))}
% CONSTRAINED DECONVOLUTION (constrained_foopsi.m)
{'cvx'}
{0}
@@ -216,7 +218,7 @@
{0.85}
{1}
% VIDEO (make_patch_video.m)
- {[1:4]}
+ {1:4}
{1}
{16}
{0}
diff --git a/demo_script.m b/demo_script.m
index 3baf9f4..6b5f384 100644
--- a/demo_script.m
+++ b/demo_script.m
@@ -8,6 +8,7 @@
num2read=2000; % user input: how many frames to read (optional, default until the end)
Y = bigread2(nam,sframe,num2read);
+Y = Y - min(Y(:));
if ~isa(Y,'double'); Y = double(Y); end % convert to double
[d1,d2,T] = size(Y); % dimensions of dataset
@@ -25,7 +26,7 @@
'search_method','ellipse','dist',3,... % search locations when updating spatial components
'deconv_method','constrained_foopsi',... % activity deconvolution method
'temporal_iter',2,... % number of block-coordinate descent steps
- 'fudge_factor',0.98,... % bias correction for AR coefficients
+ 'fudge_factor',0.98,... % bias correction for AR coefficients
'merge_thr',merge_thr... % merging threshold
);
%% Data pre-processing
@@ -50,7 +51,6 @@
[A,b] = update_spatial_components(Yr,Cin,fin,Ain,P,options);
%% update temporal components
-% consider using update_temporal_components_parallel for increased speed
[C,f,P,S] = update_temporal_components(Yr,A,b,Cin,fin,P,options);
%% merge found components
diff --git a/update_temporal_components.m b/update_temporal_components.m
index f257e83..e686d9e 100644
--- a/update_temporal_components.m
+++ b/update_temporal_components.m
@@ -19,6 +19,9 @@
% subject to: G*c_j >= 0
% ||Y(i,:) - A*C - b*f|| <= sn(i)*sqrt(T)
+% The update can happen either in parallel (default) or serial by tuning options.temporal_parallel.
+% In the case of parallel implementation the methods 'MCEM_foopsi' and 'noise_constrained' are not supported
+
% INPUTS:
% Y: raw data ( d X T matrix)
% A: spatial footprints (d x nr matrix)
@@ -53,6 +56,7 @@
if ~isfield(options,'temporal_iter') || isempty(options.temporal_iter); ITER = defoptions.temporal_iter; else ITER = options.temporal_iter; end % number of block-coordinate descent iterations
if ~isfield(options,'bas_nonneg'); options.bas_nonneg = defoptions.bas_nonneg; end
if ~isfield(options,'fudge_factor'); options.fudge_factor = defoptions.fudge_factor; end
+if ~isfield(options,'temporal_parallel'); options.temporal_parallel = defoptions.temporal_parallel; end
if isfield(P,'interp'); Y_interp = P.interp; else Y_interp = sparse(d,T); end % missing data
if isfield(P,'unsaturatedPix'); unsaturatedPix = P.unsaturatedPix; else unsaturatedPix = 1:d; end % saturated pixels
@@ -62,8 +66,14 @@
if (strcmpi(method,'noise_constrained') || strcmpi(method,'project')) && ~isfield(P,'g')
options.flag_g = 1;
- if ~isfield(P,'p') || isempty(P.p); p = 2; end;
- [P,Y] = preprocess_data(Y,p,options);
+ if ~isfield(P,'p') || isempty(P.p); P.p = 2; end;
+ p = P.p;
+ P = arpfit(Yr,p,options,P.sn);
+ if ~iscell(P.g)
+ G = make_G_matrix(T,P.g);
+ end
+else
+ G = speye(T);
end
ff = find(sum(A)==0);
@@ -117,99 +127,208 @@
P.c1 = cell(K,1);
P.neuron_sn = cell(K,1);
end
+ if strcmpi(method,'MCMC')
+ params.B = 300;
+ params.Nsamples = 400;
+ params.p = P.p;
+ else
+ params = [];
+ end
end
-for iter = 1:ITER
+p = P.p;
+if options.temporal_parallel
+ for iter = 1:ITER
+ [O,lo] = update_order(A(:,1:K));
+ for jo = 1:length(O)
+ Ytemp = YrA(:,O{jo}(:)) + (diag(nA(O{jo}))*Cin(O{jo},:))';
+ Ctemp = zeros(length(O{jo}),T);
+ Stemp = zeros(length(O{jo}),T);
+ btemp = zeros(length(O{jo}),1);
+ sntemp = btemp;
+ c1temp = btemp;
+ gtemp = cell(length(O{jo}),1);
+ nT = nA(O{jo});
+ % FN added the part below in order to save SAMPLES as a field of P
+ if strcmpi(method,'MCMC')
+ clear samples_mcmc
+ samples_mcmc(length(O{jo})) = struct();
+ [samples_mcmc.Cb] = deal(zeros(params.Nsamples,1));
+ [samples_mcmc.Cin] = deal(zeros(params.Nsamples,1));
+ [samples_mcmc.sn2] = deal(zeros(params.Nsamples,1));
+ [samples_mcmc.ns] = deal(zeros(params.Nsamples,1));
+ [samples_mcmc.ss] = deal(cell(params.Nsamples,1));
+ [samples_mcmc.ld] = deal(zeros(params.Nsamples,1));
+ [samples_mcmc.Am] = deal(zeros(params.Nsamples,1));
+ [samples_mcmc.g] = deal(zeros(params.Nsamples,1));
+ [samples_mcmc.params] = deal(struct('lam_', [], 'spiketimes_', [], 'A_', [], 'b_', [], 'C_in', [], 'sg', [], 'g', []));
+ end
+ parfor jj = 1:length(O{jo})
+ if p == 0 % p = 0 (no dynamics assumed)
+ cc = max(Ytemp(:,jj)/nT(jj),0);
+ Ctemp(jj,:) = full(cc');
+ Stemp(jj,:) = C(jj,:);
+ else
+ switch method
+ case 'project'
+ maxy = max(Ytemp(:,jj)/nT(jj));
+ cc = plain_foopsi(Ytemp(:,jj)/nT(jj)/maxy,G);
+ Ctemp(jj,:) = full(cc')*maxy;
+ Stemp(jj,:) = Ctemp(jj,:)*G';
+ case 'constrained_foopsi'
+ %if restimate_g
+ [cc,cb,c1,gn,sn,spk] = constrained_foopsi(Ytemp(:,jj)/nT(jj),[],[],[],[],options);
+ %else
+ % [cc,cb,c1,gn,sn,spk] = constrained_foopsi(Ytemp(:,jj)/nA(jj),[],[],P.g,[],options);
+ %end
+ gd = max(roots([1,-gn'])); % decay time constant for initial concentration
+ gd_vec = gd.^((0:T-1));
+ Ctemp(jj,:) = full(cc(:)' + cb + c1*gd_vec);
+ Stemp(jj,:) = spk(:)';
+ Ytemp(:,jj) = Ytemp(:,jj) - nT(jj)*Ctemp(jj,:)';
+ btemp(jj) = cb;
+ c1temp(jj) = c1;
+ sntemp(jj) = sn;
+ gtemp{jj} = gn(:)';
+ case 'MCMC'
+ SAMPLES = cont_ca_sampler(Ytemp(:,jj)/nT(jj),params);
+ Ctemp(jj,:) = make_mean_sample(SAMPLES,Ytemp(:,jj)/nT(jj));
+ Stemp(jj,:) = mean(samples_cell2mat(SAMPLES.ss,T));
+ btemp(jj) = mean(SAMPLES.Cb);
+ c1temp(jj) = mean(SAMPLES.Cin);
+ sntemp(jj) = sqrt(mean(SAMPLES.sn2));
+ gtemp{jj} = mean(exp(-1./SAMPLES.g))';
+ samples_mcmc(jj) = SAMPLES; % FN added.
+ end
+ end
+ end
+ if p > 0
+ if strcmpi(method,'constrained_foopsi') || strcmpi(method,'MCMC');
+ P.b(O{jo}) = num2cell(btemp);
+ P.c1(O{jo}) = num2cell(c1temp);
+ P.neuron_sn(O{jo}) = num2cell(sntemp);
+ for jj = 1:length(O{jo})
+ P.gn(O{jo}(jj)) = gtemp(jj);
+ end
+ YrA(:,O{jo}(:)) = Ytemp;
+ C(O{jo}(:),:) = Ctemp;
+ S(O{jo}(:),:) = Stemp;
+
+ if strcmpi(method,'MCMC');
+ P.samples_mcmc(O{jo}) = samples_mcmc; % FN added, a useful parameter to have.
+ end
+ end
+ end
+ fprintf('%i out of %i components updated \n',sum(lo(1:jo)),K);
+ end
+ ii = K + 1;
+ YrA(:,ii) = YrA(:,ii) + nA(ii)*Cin(ii,:)';
+ cc = max(YrA(:,ii)/nA(ii),0);
+ C(ii,:) = full(cc');
+ YrA(:,ii) = YrA(:,ii) - nA(ii)*C(ii,:)';
+ % if mod(jj,10) == 0
+ % fprintf('%i out of total %i temporal components updated \n',jj,K);
+ % end
+ %disp(norm(Fin(1:nr,:) - F,'fro')/norm(F,'fro'));
+ if norm(Cin - C,'fro')/norm(C,'fro') <= 1e-3
+ % stop if the overall temporal component does not change by much
+ break;
+ else
+ Cin = C;
+ end
+ end
+else
+ for iter = 1:ITER
perm = randperm(K+size(b,2));
- for jj = 1:K
- ii = perm(jj);
- if ii<=K
- if P.p == 0 % p = 0 (no dynamics assumed)
+ for jj = 1:K
+ ii = perm(jj);
+ if ii<=K
+ if P.p == 0 % p = 0 (no dynamics assumed)
+ YrA(:,ii) = YrA(:,ii) + nA(ii)*Cin(ii,:)';
+ maxy = max(YrA(:,ii)/nA(ii));
+ cc = max(YrA(:,ii)/nA(ii)/maxy,0);
+ C(ii,:) = full(cc')*maxy;
+ YrA(:,ii) = YrA(:,ii) - nA(ii)*C(ii,:)';
+ S(ii,:) = C(ii,:);
+ else
+ switch method
+ case 'project'
+ YrA(:,ii) = YrA(:,ii) + nA(ii)*Cin(ii,:)';
+ maxy = max(YrA(:,ii)/nA(ii));
+ cc = plain_foopsi(YrA(:,ii)/nA(ii)/maxy,G);
+ C(ii,:) = full(cc')*maxy;
+ YrA(:,ii) = YrA(:,ii) - nA(ii)*C(ii,:)';
+ S(ii,:) = C(ii,:)*G';
+ case 'constrained_foopsi'
+ YrA(:,ii) = YrA(:,ii) + nA(ii)*Cin(ii,:)';
+ if restimate_g
+ [cc,cb,c1,gn,sn,spk] = constrained_foopsi(YrA(:,ii)/nA(ii),[],[],[],[],options);
+ P.gn{ii} = gn;
+ else
+ [cc,cb,c1,gn,sn,spk] = constrained_foopsi(YrA(:,ii)/nA(ii),[],[],P.g,[],options);
+ end
+ gd = max(roots([1,-gn'])); % decay time constant for initial concentration
+ gd_vec = gd.^((0:T-1));
+ C(ii,:) = full(cc(:)' + cb + c1*gd_vec);
+ S(ii,:) = spk(:)';
+ YrA(:,ii) = YrA(:,ii) - nA(ii)*C(ii,:)';
+ P.b{ii} = cb;
+ P.c1{ii} = c1;
+ P.neuron_sn{ii} = sn;
+ case 'MCEM_foopsi'
+ options.p = length(P.g);
+ YrA(:,ii) = YrA(:,ii) + nA(ii)*Cin(ii,:)';
+ [cc,cb,c1,gn,sn,spk] = MCEM_foopsi(YrA(:,ii)/nA(ii),[],[],P.g,[],options);
+ gd = max(roots([1,-gn.g(:)']));
+ gd_vec = gd.^((0:T-1));
+ C(ii,:) = full(cc(:)' + cb + c1*gd_vec);
+ S(ii,:) = spk(:)';
+ YrA(:,ii) = YrA(:,ii) - nA(ii)*C(ii,:)';
+ P.b{ii} = cb;
+ P.c1{ii} = c1;
+ P.neuron_sn{ii} = sn;
+ P.gn{ii} = gn.g;
+ case 'MCMC'
+ params.B = 300;
+ params.Nsamples = 400;
+ params.p = P.p; %length(P.g);
+ YrA(:,ii) = YrA(:,ii) + nA(ii)*Cin(ii,:)';
+ SAMPLES = cont_ca_sampler(YrA(:,ii)/nA(ii),params);
+ C(ii,:) = make_mean_sample(SAMPLES,YrA(:,ii)/nA(ii));
+ S(ii,:) = mean(samples_cell2mat(SAMPLES.ss,T));
+ YrA(:,ii) = YrA(:,ii) - nA(ii)*C(ii,:)';
+ P.b{ii} = mean(SAMPLES.Cb);
+ P.c1{ii} = mean(SAMPLES.Cin);
+ P.neuron_sn{ii} = sqrt(mean(SAMPLES.sn2));
+ P.gn{ii} = mean(exp(-1./SAMPLES.g));
+ P.samples_mcmc(ii) = SAMPLES; % FN added, a useful parameter to have.
+ case 'noise_constrained'
+ Y_res = Y_res + A(:,ii)*Cin(ii,:);
+ [~,srt] = sort(A(:,ii),'descend');
+ ff = srt(1:mc);
+ [cc,LD(:,ii)] = lagrangian_foopsi_temporal(Y_res(ff,:),A(ff,ii),T*P.sn(unsaturatedPix(ff)).^2,G,LD(:,ii));
+ C(ii,:) = full(cc');
+ Y_res = Y_res - A(:,ii)*cc';
+ S(ii,:) = C(ii,:)*G';
+ end
+ end
+ else
YrA(:,ii) = YrA(:,ii) + nA(ii)*Cin(ii,:)';
- maxy = max(YrA(:,ii)/nA(ii));
- cc = max(YrA(:,ii)/nA(ii)/maxy,0);
- C(ii,:) = full(cc')*maxy;
+ cc = max(YrA(:,ii)/nA(ii),0);
+ C(ii,:) = full(cc');
YrA(:,ii) = YrA(:,ii) - nA(ii)*C(ii,:)';
- S(ii,:) = C(ii,:);
- else
- switch method
- case 'project'
- YrA(:,ii) = YrA(:,ii) + nA(ii)*Cin(ii,:)';
- maxy = max(YrA(:,ii)/nA(ii));
- cc = plain_foopsi(YrA(:,ii)/nA(ii)/maxy,G);
- C(ii,:) = full(cc')*maxy;
- YrA(:,ii) = YrA(:,ii) - nA(ii)*C(ii,:)';
- S(ii,:) = C(ii,:)*G';
- case 'constrained_foopsi'
- YrA(:,ii) = YrA(:,ii) + nA(ii)*Cin(ii,:)';
- if restimate_g
- [cc,cb,c1,gn,sn,spk] = constrained_foopsi(YrA(:,ii)/nA(ii),[],[],[],[],options);
- P.gn{ii} = gn;
- else
- [cc,cb,c1,gn,sn,spk] = constrained_foopsi(YrA(:,ii)/nA(ii),[],[],P.g,[],options);
- end
- gd = max(roots([1,-gn'])); % decay time constant for initial concentration
- gd_vec = gd.^((0:T-1));
- C(ii,:) = full(cc(:)' + cb + c1*gd_vec);
- S(ii,:) = spk(:)';
- YrA(:,ii) = YrA(:,ii) - nA(ii)*C(ii,:)';
- P.b{ii} = cb;
- P.c1{ii} = c1;
- P.neuron_sn{ii} = sn;
- case 'MCEM_foopsi'
- options.p = length(P.g);
- YrA(:,ii) = YrA(:,ii) + nA(ii)*Cin(ii,:)';
- [cc,cb,c1,gn,sn,spk] = MCEM_foopsi(YrA(:,ii)/nA(ii),[],[],P.g,[],options);
- gd = max(roots([1,-gn.g(:)']));
- gd_vec = gd.^((0:T-1));
- C(ii,:) = full(cc(:)' + cb + c1*gd_vec);
- S(ii,:) = spk(:)';
- YrA(:,ii) = YrA(:,ii) - nA(ii)*C(ii,:)';
- P.b{ii} = cb;
- P.c1{ii} = c1;
- P.neuron_sn{ii} = sn;
- P.gn{ii} = gn.g;
- case 'MCMC'
- params.B = 300;
- params.Nsamples = 400;
- params.p = P.p; %length(P.g);
- YrA(:,ii) = YrA(:,ii) + nA(ii)*Cin(ii,:)';
- SAMPLES = cont_ca_sampler(YrA(:,ii)/nA(ii),params);
- C(ii,:) = make_mean_sample(SAMPLES,YrA(:,ii)/nA(ii));
- S(ii,:) = mean(samples_cell2mat(SAMPLES.ss,T));
- YrA(:,ii) = YrA(:,ii) - nA(ii)*C(ii,:)';
- P.b{ii} = mean(SAMPLES.Cb);
- P.c1{ii} = mean(SAMPLES.Cin);
- P.neuron_sn{ii} = sqrt(mean(SAMPLES.sn2));
- P.gn{ii} = mean(exp(-1./SAMPLES.g));
- P.samples_mcmc(ii) = SAMPLES; % FN added, a useful parameter to have.
- case 'noise_constrained'
- Y_res = Y_res + A(:,ii)*Cin(ii,:);
- [~,srt] = sort(A(:,ii),'descend');
- ff = srt(1:mc);
- [cc,LD(:,ii)] = lagrangian_foopsi_temporal(Y_res(ff,:),A(ff,ii),T*P.sn(unsaturatedPix(ff)).^2,G,LD(:,ii));
- C(ii,:) = full(cc');
- Y_res = Y_res - A(:,ii)*cc';
- S(ii,:) = C(ii,:)*G';
- end
end
- else
- YrA(:,ii) = YrA(:,ii) + nA(ii)*Cin(ii,:)';
- cc = max(YrA(:,ii)/nA(ii),0);
- C(ii,:) = full(cc');
- YrA(:,ii) = YrA(:,ii) - nA(ii)*C(ii,:)';
+ if mod(jj,10) == 0
+ fprintf('%i out of total %i temporal components updated \n',jj,K);
+ end
end
- if mod(jj,10) == 0
- fprintf('%i out of total %i temporal components updated \n',jj,K);
+ %disp(norm(Fin(1:nr,:) - F,'fro')/norm(F,'fro'));
+ if norm(Cin - C,'fro')/norm(C,'fro') <= 1e-3
+ % stop if the overall temporal component does not change by much
+ break;
+ else
+ Cin = C;
end
end
- %disp(norm(Fin(1:nr,:) - F,'fro')/norm(F,'fro'));
- if norm(Cin - C,'fro')/norm(C,'fro') <= 1e-3
- % stop if the overall temporal component does not change by much
- break;
- else
- Cin = C;
- end
end
-
f = C(K+1:end,:);
-C = C(1:K,:);
+C = C(1:K,:);
\ No newline at end of file
diff --git a/update_temporal_components_parallel.m b/update_temporal_components_parallel.m
index bcacf38..77d84a8 100644
--- a/update_temporal_components_parallel.m
+++ b/update_temporal_components_parallel.m
@@ -2,236 +2,13 @@
% update temporal components and background given spatial components in
% parallel, by forming of sequence of vertex covers.
-% A variety of different methods can be used and are separated into 2 classes:
-
-% 1-d approaches, where for each component a 1-d trace is computed by removing
-% the effect of all the other components and then averaging with the corresponding
-% spatial footprint. Then each trace is denoised. This corresponds to a block-coordinate approach
-% 4 different 1-d approaches are included, and any custom method
-% can be easily incorporated:
-% 'project': The trace is projected to satisfy the constraints by the (known) calcium indicator dynamics
-% 'constrained_foopsi': The noise constrained deconvolution approach is used. Time constants can be re-estimated (default)
-% 'MCEM_foopsi': Alternating between constrained_foopsi and a MH approach for re-estimating the time constants
-% 'MCMC': A fully Bayesian method (slowest, but usually most accurate)
-
-% multi-dimensional approaches: (slowest)
-% 'noise_constrained':
-% C(j,:) = argmin_{c_j} sum(G*c_j),
-% subject to: G*c_j >= 0
-% ||Y(i,:) - A*C - b*f|| <= sn(i)*sqrt(T)
-
-% INPUTS:
-% Y: raw data ( d X T matrix)
-% A: spatial footprints (d x nr matrix)
-% b: spatial background (d x 1 vector)
-% Cin: current estimate of temporal components (nr X T matrix)
-% fin: current estimate of temporal background (1 x T vector)
-% P: struct for neuron parameters
-% options: struct for algorithm parameters
-
-% LD: Lagrange multipliers (needed only for 'noise_constrained' method).
-%
-% OUTPUTS:
-% C: temporal components (nr X T matrix)
-% f: temporal background (1 x T vector)
-% P: struct for neuron parameters
-% S: deconvolved activity
+% Note that the parallel implementation is now implemented by default from
+% update_temporal_components. update_temporal_components_parallel will be
+% removed in the next release.
% Written by:
% Eftychios A. Pnevmatikakis, Simons Foundation, 2015
-[d,T] = size(Y);
-if isempty(P) || nargin < 6
- active_pixels = find(sum(A,2)); % pixels where the greedy method found activity
- unsaturated_pixels = find_unsaturatedPixels(Y); % pixels that do not exhibit saturation
- options.pixels = intersect(active_pixels,unsaturated_pixels); % base estimates only on unsaturated, active pixels
-end
-
-defoptions = CNMFSetParms;
-if nargin < 7 || isempty(options); options = []; end
-if ~isfield(options,'deconv_method') || isempty(options.deconv_method); method = defoptions.deconv_method; else method = options.deconv_method; end % choose method
-if ~isfield(options,'restimate_g') || isempty(options.restimate_g); restimate_g = defoptions.restimate_g; else restimate_g = options.restimate_g; end % re-estimate time constant (only with constrained foopsi)
-if ~isfield(options,'temporal_iter') || isempty(options.temporal_iter); ITER = defoptions.temporal_iter; else ITER = options.temporal_iter; end % number of block-coordinate descent iterations
-if ~isfield(options,'bas_nonneg'); options.bas_nonneg = defoptions.bas_nonneg; end
-if ~isfield(options,'fudge_factor'); options.fudge_factor = defoptions.fudge_factor; end
-
-if isfield(P,'interp'); Y_interp = P.interp; else Y_interp = sparse(d,T); end % missing data
-if isfield(P,'unsaturatedPix'); unsaturatedPix = P.unsaturatedPix; else unsaturatedPix = 1:d; end % saturated pixels
-
-mis_data = find(Y_interp); % interpolate any missing data before deconvolution
-Y(mis_data) = Y_interp(mis_data);
-
-if (strcmpi(method,'noise_constrained') || strcmpi(method,'project')) && ~isfield(P,'g')
- options.flag_g = 1;
- if ~isfield(P,'p') || isempty(P.p); P.p = 2; end;
- p = P.p;
- P = arpfit(Yr,p,options,P.sn);
- if ~iscell(P.g)
- G = make_G_matrix(T,P.g);
- end
-else
- G = speye(T);
-end
-
-ff = find(sum(A)==0);
-if ~isempty(ff)
- A(:,ff) = [];
- if exist('Cin','var')
- if ~isempty(Cin)
- Cin(ff,:) = [];
- end
- end
-end
-
-if isempty(Cin) || nargin < 4
- Cin = max((A'*A)\(A'*Y),0);
- ITER = max(ITER,3);
-end
-
-if isempty(b) || isempty(fin) || nargin < 5
- if isempty(b) || nargin < 3
- [b,fin] = nnmf(max(Y - A*Cin,0),1);
- else
- fin = max(b'*Y/norm(b)^2,0);
- end
-end
-
-saturatedPix = setdiff(1:d,unsaturatedPix); % remove any saturated pixels
-Ysat = Y(saturatedPix,:);
-Asat = A(saturatedPix,:);
-bsat = b(saturatedPix,:);
-Y = Y(unsaturatedPix,:);
-A = A(unsaturatedPix,:);
-b = b(unsaturatedPix,:);
-d = length(unsaturatedPix);
-
-K = size(A,2);
-A = [A,b];
-S = zeros(size(Cin));
-Cin = [Cin;fin];
-C = Cin;
-
-if strcmpi(method,'noise_constrained')
- Y_res = Y - A*Cin;
- mc = min(d,15); % number of constraints to be considered
- LD = 10*ones(mc,K);
-else
- nA = sum(A.^2);
- YrA = Y'*A - Cin'*(A'*A);
- if strcmpi(method,'constrained_foopsi') || strcmpi(method,'MCEM_foopsi')
- P.gn = cell(K,1);
- P.b = cell(K,1);
- P.c1 = cell(K,1);
- P.neuron_sn = cell(K,1);
- end
- if strcmpi(method,'MCMC')
- params.B = 300;
- params.Nsamples = 400;
- params.p = P.p;
- else
- params = [];
- end
-end
-p = P.p;
-for iter = 1:ITER
- [O,lo] = update_order(A(:,1:K));
- for jo = 1:length(O)
- Ytemp = YrA(:,O{jo}(:)) + (diag(nA(O{jo}))*Cin(O{jo},:))';
- Ctemp = zeros(length(O{jo}),T);
- Stemp = zeros(length(O{jo}),T);
- btemp = zeros(length(O{jo}),1);
- sntemp = btemp;
- c1temp = btemp;
- gtemp = cell(length(O{jo}),1);
- nT = nA(O{jo});
- % FN added the part below in order to save SAMPLES as a field of P
- if strcmpi(method,'MCMC')
- clear samples_mcmc
- samples_mcmc(length(O{jo})) = struct();
- [samples_mcmc.Cb] = deal(zeros(params.Nsamples,1));
- [samples_mcmc.Cin] = deal(zeros(params.Nsamples,1));
- [samples_mcmc.sn2] = deal(zeros(params.Nsamples,1));
- [samples_mcmc.ns] = deal(zeros(params.Nsamples,1));
- [samples_mcmc.ss] = deal(cell(params.Nsamples,1));
- [samples_mcmc.ld] = deal(zeros(params.Nsamples,1));
- [samples_mcmc.Am] = deal(zeros(params.Nsamples,1));
- [samples_mcmc.g] = deal(zeros(params.Nsamples,1));
- [samples_mcmc.params] = deal(struct('lam_', [], 'spiketimes_', [], 'A_', [], 'b_', [], 'C_in', [], 'sg', [], 'g', []));
- end
- parfor jj = 1:length(O{jo})
- if p == 0 % p = 0 (no dynamics assumed)
- cc = max(Ytemp(:,jj)/nT(jj),0);
- Ctemp(jj,:) = full(cc');
- Stemp(jj,:) = C(jj,:);
- else
- switch method
- case 'project'
- maxy = max(Ytemp(:,jj)/nT(jj));
- cc = plain_foopsi(Ytemp(:,jj)/nT(jj)/maxy,G);
- Ctemp(jj,:) = full(cc')*maxy;
- Stemp(jj,:) = Ctemp(jj,:)*G';
- case 'constrained_foopsi'
- %if restimate_g
- [cc,cb,c1,gn,sn,spk] = constrained_foopsi(Ytemp(:,jj)/nT(jj),[],[],[],[],options);
- %else
- % [cc,cb,c1,gn,sn,spk] = constrained_foopsi(Ytemp(:,jj)/nA(jj),[],[],P.g,[],options);
- %end
- gd = max(roots([1,-gn'])); % decay time constant for initial concentration
- gd_vec = gd.^((0:T-1));
- Ctemp(jj,:) = full(cc(:)' + cb + c1*gd_vec);
- Stemp(jj,:) = spk(:)';
- Ytemp(:,jj) = Ytemp(:,jj) - nT(jj)*Ctemp(jj,:)';
- btemp(jj) = cb;
- c1temp(jj) = c1;
- sntemp(jj) = sn;
- gtemp{jj} = gn(:)';
- case 'MCMC'
- SAMPLES = cont_ca_sampler(Ytemp(:,jj)/nT(jj),params);
- Ctemp(jj,:) = make_mean_sample(SAMPLES,Ytemp(:,jj)/nT(jj));
- Stemp(jj,:) = mean(samples_cell2mat(SAMPLES.ss,T));
- btemp(jj) = mean(SAMPLES.Cb);
- c1temp(jj) = mean(SAMPLES.Cin);
- sntemp(jj) = sqrt(mean(SAMPLES.sn2));
- gtemp{jj} = mean(exp(-1./SAMPLES.g))';
- samples_mcmc(jj) = SAMPLES; % FN added.
- end
- end
- end
- if p > 0
- if strcmpi(method,'constrained_foopsi') || strcmpi(method,'MCMC');
- P.b(O{jo}) = num2cell(btemp);
- P.c1(O{jo}) = num2cell(c1temp);
- P.neuron_sn(O{jo}) = num2cell(sntemp);
- for jj = 1:length(O{jo})
- P.gn(O{jo}(jj)) = gtemp(jj);
- end
- YrA(:,O{jo}(:)) = Ytemp;
- C(O{jo}(:),:) = Ctemp;
- S(O{jo}(:),:) = Stemp;
-
- if strcmpi(method,'MCMC');
- P.samples_mcmc(O{jo}) = samples_mcmc; % FN added, a useful parameter to have.
- end
- end
- end
- fprintf('%i out of %i components updated \n',sum(lo(1:jo)),K);
- end
- ii = K + 1;
- YrA(:,ii) = YrA(:,ii) + nA(ii)*Cin(ii,:)';
- cc = max(YrA(:,ii)/nA(ii),0);
- C(ii,:) = full(cc');
- YrA(:,ii) = YrA(:,ii) - nA(ii)*C(ii,:)';
-% if mod(jj,10) == 0
-% fprintf('%i out of total %i temporal components updated \n',jj,K);
-% end
- %disp(norm(Fin(1:nr,:) - F,'fro')/norm(F,'fro'));
- if norm(Cin - C,'fro')/norm(C,'fro') <= 1e-3
- % stop if the overall temporal component does not change by much
- break;
- else
- Cin = C;
- end
-end
+if ~isempty(options); options.temporal_parallel = 1; end
-f = C(K+1:end,:);
-C = C(1:K,:);
+[C,f,P,S] = update_temporal_components(Y,A,b,Cin,fin,P,options);
\ No newline at end of file