Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

nb=0 #1172

Open
wants to merge 3 commits into
base: dev
Choose a base branch
from
Open

nb=0 #1172

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
315ac76
enable nb=0 (#1068); update deprecated selem in skimage.dilation (#1130)
j-friedrich Sep 19, 2023
b8f0153
fix bug when using patches and nb=0
j-friedrich Sep 19, 2023
132fca4
fix bug when using tsub or ssub and nb=0
j-friedrich Sep 20, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions caiman/base/rois.py
View file
Open in desktop
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since the selem keyword was deprecated at scikit learn, and we still have it, it might be nice to define what it means in this function: e.g., in the docs for it maybe just mention it stands for 'structuring element'

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you mean skimage? I remember playing with this years ago, it means something like the "where stuff is multiplier" =S

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes oops scikitimage

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this explains selem, the structuring element https://scikit-image.org/docs/stable/auto_examples/applications/plot_morphology.html

Original file line number Diff line number Diff line change
Expand Up @@ -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')

Expand Down
27 changes: 18 additions & 9 deletions caiman/source_extraction/cnmf/initialization.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down Expand Up @@ -970,11 +974,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, 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

#%%
Expand Down Expand Up @@ -1140,7 +1146,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)
Expand All @@ -1166,13 +1172,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
Expand Down Expand Up @@ -1385,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')
Expand Down
32 changes: 17 additions & 15 deletions caiman/source_extraction/cnmf/spatial.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We may want to consider making this threshold overridable or at least less magical

f = Y[not_px, :].mean(0)
else:
print('estimating f')
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There are a few print statements in here -- in final draft convert to logger info or debug statements?

f = 0
for xxx in np.where(not_px)[0]:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we give this variable a better name?

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))
Expand Down
3 changes: 2 additions & 1 deletion caiman/source_extraction/cnmf/temporal.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you add a comment explaining how this works?

Cin = np.vstack((Cin, fin))
C = Cin.copy()
nA = np.ravel(A.power(2).sum(axis=0)) + np.finfo(np.float32).eps

Expand Down