From 315ac76c80aa3209de04042838a5ed65dbcbc669 Mon Sep 17 00:00:00 2001 From: Johannes Friedrich Date: Tue, 19 Sep 2023 13:44:49 -0700 Subject: [PATCH 1/3] enable nb=0 (#1068); update deprecated selem in skimage.dilation (#1130) --- caiman/base/rois.py | 4 +-- .../source_extraction/cnmf/initialization.py | 21 +++++++----- caiman/source_extraction/cnmf/spatial.py | 32 ++++++++++--------- 3 files changed, 32 insertions(+), 25 deletions(-) diff --git a/caiman/base/rois.py b/caiman/base/rois.py index e268aeeff..a3671795f 100644 --- a/caiman/base/rois.py +++ b/caiman/base/rois.py @@ -131,9 +131,9 @@ def extract_binary_masks_from_structural_channel(Y, for i in range(areas[1]): temp = (areas[0] == i + 1) if expand_method == 'dilation': - temp = dilation(temp, selem=selem) + temp = dilation(temp, footprint=selem) elif expand_method == 'closing': - temp = closing(temp, selem=selem) + temp = closing(temp, footprint=selem) A[:, i] = temp.flatten('F') diff --git a/caiman/source_extraction/cnmf/initialization.py b/caiman/source_extraction/cnmf/initialization.py index df734561e..eeed626c3 100644 --- a/caiman/source_extraction/cnmf/initialization.py +++ b/caiman/source_extraction/cnmf/initialization.py @@ -970,11 +970,13 @@ def onclick(event): res = np.reshape(Y, (np.prod(d[0:-1]), d[-1]), order='F') + med.flatten(order='F')[:, None] -# model = NMF(n_components=nb, init='random', random_state=0) - model = NMF(n_components=nb, init='nndsvdar') - b_in = model.fit_transform(np.maximum(res, 0)).astype(np.float32) - f_in = model.components_.astype(np.float32) - + if nb > 0: + model = NMF(n_components=nb, init='nndsvdar') + b_in = model.fit_transform(np.maximum(res, 0)).astype(np.float32) + f_in = model.components_.astype(np.float32) + else: + b_in = np.empty((A.shape[0], 0), dtype=np.float32) + f_in = np.empty((0, C.shape[1]), dtype=np.float32) return A, C, np.array(center, dtype='uint16'), b_in, f_in #%% @@ -1140,7 +1142,7 @@ def hals(Y, A, C, b, f, bSiz=3, maxIter=5): else: ind_A = A>1e-10 - ind_A = spr.csc_matrix(ind_A) # indicator of nonnero pixels + ind_A = spr.csc_matrix(ind_A) # indicator of nonzero pixels def HALS4activity(Yr, A, C, iters=2): U = A.T.dot(Yr) @@ -1166,13 +1168,16 @@ def HALS4shape(Yr, A, C, iters=2): return A Ab = np.c_[A, b] - Cf = np.r_[C, f.reshape(nb, -1)] + Cf = np.r_[C, f] for _ in range(maxIter): Cf = HALS4activity(np.reshape( Y, (np.prod(dims), T), order='F'), Ab, Cf) Ab = HALS4shape(np.reshape(Y, (np.prod(dims), T), order='F'), Ab, Cf) - return Ab[:, :-nb], Cf[:-nb], Ab[:, -nb:], Cf[-nb:].reshape(nb, -1) + if nb == 0: + return Ab, Cf, b, f + else: + return Ab[:, :-nb], Cf[:-nb], Ab[:, -nb:], Cf[-nb:].reshape(nb, -1) @profile diff --git a/caiman/source_extraction/cnmf/spatial.py b/caiman/source_extraction/cnmf/spatial.py index dc8a0d7ef..29809fa66 100644 --- a/caiman/source_extraction/cnmf/spatial.py +++ b/caiman/source_extraction/cnmf/spatial.py @@ -1050,22 +1050,24 @@ def computing_indicator(Y, A_in, b, C, f, nb, method, dims, min_size, max_size, px = (np.sum(dist_indicator, axis=1) > 0) not_px = ~px - if nb>1: - f = NMF(nb, init='nndsvda').fit(np.maximum(Y[not_px, :], 0)).components_ - else: - if Y.shape[-1] < 30000: - f = Y[not_px, :].mean(0) + if nb > 0: + if nb > 1: + f = NMF(nb, init='nndsvda').fit(np.maximum(Y[not_px, :], 0)).components_ else: - print('estimating f') - f = 0 - for xxx in np.where(not_px)[0]: - f += Y[xxx] - f /= not_px.sum() - - f = np.atleast_2d(f) - - Y_resf = np.dot(Y, f.T) - b = np.maximum(Y_resf, 0) / (np.linalg.norm(f)**2) + if Y.shape[-1] < 30000: + f = Y[not_px, :].mean(0) + else: + print('estimating f') + f = 0 + for xxx in np.where(not_px)[0]: + f += Y[xxx] + f /= not_px.sum() + f = np.atleast_2d(f) + Y_resf = np.dot(Y, f.T) + b = np.maximum(Y_resf, 0) / (np.linalg.norm(f)**2) + else: + f = np.empty((0, Y.shape[-1]), dtype='float32') + b = np.empty((Y.shape[0], 0), dtype='float32') C = np.maximum(csr_matrix(dist_indicator_av.T).dot( Y) - dist_indicator_av.T.dot(b).dot(f), 0) A_in = scipy.sparse.coo_matrix(A_in.astype(np.float32)) From b8f0153a30a0b74b0f6b41ddba3c21cc12310eeb Mon Sep 17 00:00:00 2001 From: Johannes Friedrich Date: Tue, 19 Sep 2023 14:05:44 -0700 Subject: [PATCH 2/3] fix bug when using patches and nb=0 --- caiman/source_extraction/cnmf/temporal.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/caiman/source_extraction/cnmf/temporal.py b/caiman/source_extraction/cnmf/temporal.py index 1f03be530..a1c9973b2 100644 --- a/caiman/source_extraction/cnmf/temporal.py +++ b/caiman/source_extraction/cnmf/temporal.py @@ -197,7 +197,8 @@ def update_temporal_components(Y, A, b, Cin, fin, bl=None, c1=None, g=None, sn=N A = scipy.sparse.hstack((A, b)).tocsc() S = np.zeros(np.shape(Cin)) - Cin = np.vstack((Cin, fin)) + if fin is not None: + Cin = np.vstack((Cin, fin)) C = Cin.copy() nA = np.ravel(A.power(2).sum(axis=0)) + np.finfo(np.float32).eps From 132fca4e9fa2ea91928f730e2316d3dc7fb7fd03 Mon Sep 17 00:00:00 2001 From: Johannes Friedrich Date: Wed, 20 Sep 2023 09:12:36 -0700 Subject: [PATCH 3/3] fix bug when using tsub or ssub and nb=0 --- caiman/source_extraction/cnmf/initialization.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/caiman/source_extraction/cnmf/initialization.py b/caiman/source_extraction/cnmf/initialization.py index eeed626c3..8aa26cbbd 100644 --- a/caiman/source_extraction/cnmf/initialization.py +++ b/caiman/source_extraction/cnmf/initialization.py @@ -444,6 +444,10 @@ def initialize_components(Y, K=30, gSig=[5, 5], gSiz=None, ssub=1, tsub=1, nIter f_in = resize(np.atleast_2d(f_in), [b_in.shape[-1], T]) + elif nb == 0: + b_in = np.empty((np.prod(d), 0), dtype=np.float32, order='F') + f_in = np.empty((0, T), dtype=np.float32) + if Ain.size > 0: Cin = resize(Cin, [K, T]) center = np.asarray( @@ -975,7 +979,7 @@ def onclick(event): b_in = model.fit_transform(np.maximum(res, 0)).astype(np.float32) f_in = model.components_.astype(np.float32) else: - b_in = np.empty((A.shape[0], 0), dtype=np.float32) + b_in = np.empty((A.shape[0], 0), dtype=np.float32, order='F') f_in = np.empty((0, C.shape[1]), dtype=np.float32) return A, C, np.array(center, dtype='uint16'), b_in, f_in @@ -1390,7 +1394,7 @@ def compute_B(b0, W, B): # actually computes -B to efficiently compute Y-B in p b_in, s_in, f_in = spr.linalg.svds(B, k=nb) f_in *= s_in[:, np.newaxis] else: - b_in = np.empty((A.shape[0], 0)) + b_in = np.empty((A.shape[0], 0), order='F') f_in = np.empty((0, T)) if nb == 0: logging.info('Returning background as b0 and W')