From 5954c87a4c98cbd7dd027c6e955389a58f6155d9 Mon Sep 17 00:00:00 2001 From: Pierre Ablin Date: Mon, 21 Jan 2019 19:09:16 +0100 Subject: [PATCH 1/3] first try --- picard/_core_picard.py | 19 ++++++++++++++++++- picard/_tools.py | 10 ++++++++++ 2 files changed, 28 insertions(+), 1 deletion(-) diff --git a/picard/_core_picard.py b/picard/_core_picard.py index a9d0494..b128051 100644 --- a/picard/_core_picard.py +++ b/picard/_core_picard.py @@ -8,11 +8,13 @@ from scipy.linalg import expm from .densities import Tanh +from ._tools import fuse_mask def core_picard(X, density=Tanh(), ortho=False, extended=False, m=7, max_iter=1000, tol=1e-7, lambda_min=0.01, ls_tries=10, - verbose=False, covariance=None): + verbose=False, covariance=None, ll_reject=None, + reject_every=None): '''Runs the Picard algorithm The algorithm is detailed in:: @@ -83,7 +85,20 @@ def core_picard(X, density=Tanh(), ortho=False, extended=False, m=7, if covariance is None: # Need this for extended covariance = X.dot(X.T) / T C = covariance.copy() + if ll_reject is not None: # create a mask + mask = np.zeros(T, dtype=np.bool_) for n in range(max_iter): + # Rejection + if ll_reject is not None: + if (n + 1) % reject_every == 0: + ll = density.log_lik(Y) + rejected = np.sum(ll > ll_reject, axis=0) > 0 + # Exclude rejected data + Y = Y[:, ~rejected] + # Fuse the mask to keep track of the total mask + mask = fuse_mask(mask, rejected) + # Compute the new loss + current_loss = _loss(Y, W, density, signs, ortho, extended) # Compute the score function psiY, psidY = density.score_and_der(Y) # Compute the relative gradient and the Hessian off-diagonal @@ -169,6 +184,8 @@ def core_picard(X, density=Tanh(), ortho=False, extended=False, m=7, n_iterations=n) if extended: infos['signs'] = signs + if ll_reject is not None: + infos['rejected'] = mask return Y, W, infos diff --git a/picard/_tools.py b/picard/_tools.py index de40304..043d40f 100644 --- a/picard/_tools.py +++ b/picard/_tools.py @@ -89,3 +89,13 @@ def _ica_par(X, fun, max_iter, w_init, verbose): if verbose: print('Running Picard...') return W + + +def fuse_mask(mask1, mask2): + idx_m2 = 0 + N = len(mask1) + for i in range(N): + if not mask1[i]: + mask1[i] = mask2[idx_m2] + idx_m2 += 1 + return mask1 From fc743b83036e6eaa91e8a4d89b025396bc561c59 Mon Sep 17 00:00:00 2001 From: Pierre Ablin Date: Wed, 23 Jan 2019 11:37:49 +0100 Subject: [PATCH 2/3] fix --- picard/_core_picard.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/picard/_core_picard.py b/picard/_core_picard.py index b128051..3fdd9a4 100644 --- a/picard/_core_picard.py +++ b/picard/_core_picard.py @@ -80,25 +80,25 @@ def core_picard(X, density=Tanh(), ortho=False, extended=False, m=7, current_loss = _loss(Y, W, density, signs, ortho, extended) requested_tolerance = False sign_change = False + rejection = ll_reject is not None gradient_norm = 1. if extended: if covariance is None: # Need this for extended covariance = X.dot(X.T) / T C = covariance.copy() - if ll_reject is not None: # create a mask + if rejection: # create a mask mask = np.zeros(T, dtype=np.bool_) for n in range(max_iter): # Rejection - if ll_reject is not None: - if (n + 1) % reject_every == 0: - ll = density.log_lik(Y) - rejected = np.sum(ll > ll_reject, axis=0) > 0 - # Exclude rejected data - Y = Y[:, ~rejected] - # Fuse the mask to keep track of the total mask - mask = fuse_mask(mask, rejected) - # Compute the new loss - current_loss = _loss(Y, W, density, signs, ortho, extended) + if rejection and (n + 1) % reject_every == 0: + ll = density.log_lik(Y) + rejected = np.sum(ll > ll_reject, axis=0) > 0 + # Exclude rejected data + Y = Y[:, ~rejected] + # Fuse the mask to keep track of the total mask + mask = fuse_mask(mask, rejected) + # Compute the new loss + current_loss = _loss(Y, W, density, signs, ortho, extended) # Compute the score function psiY, psidY = density.score_and_der(Y) # Compute the relative gradient and the Hessian off-diagonal @@ -184,7 +184,7 @@ def core_picard(X, density=Tanh(), ortho=False, extended=False, m=7, n_iterations=n) if extended: infos['signs'] = signs - if ll_reject is not None: + if rejection: infos['rejected'] = mask return Y, W, infos From 26732985ecd2d462db69bfe28028e398290fcb26 Mon Sep 17 00:00:00 2001 From: Pierre Ablin Date: Mon, 28 Jan 2019 18:22:48 +0100 Subject: [PATCH 3/3] added test and solver doc --- picard/solver.py | 41 +++++++++++++++++++++++++++++++++---- picard/tests/test_solver.py | 25 ++++++++++++++++++++++ 2 files changed, 62 insertions(+), 4 deletions(-) diff --git a/picard/solver.py b/picard/solver.py index fc4bee9..05627eb 100644 --- a/picard/solver.py +++ b/picard/solver.py @@ -15,9 +15,9 @@ def picard(X, fun='tanh', n_components=None, ortho=True, extended=None, whiten=True, return_X_mean=False, return_n_iter=False, - centering=True, max_iter=100, tol=1e-07, m=7, ls_tries=10, - lambda_min=0.01, check_fun=True, w_init=None, fastica_it=None, - random_state=None, verbose=False): + centering=True, ll_reject=None, reject_every=5, max_iter=100, + tol=1e-07, m=7, ls_tries=10, lambda_min=0.01, check_fun=True, + w_init=None, fastica_it=None, random_state=None, verbose=False): """Perform Independent Component Analysis. Parameters @@ -67,6 +67,14 @@ def picard(X, fun='tanh', n_components=None, ortho=True, extended=None, centering : bool, optional If True, X is mean corrected. + ll_reject : None or float, optional + If None, no sample rejection is performed. Otherwise, samples with a + negative log-likelihood > reject are rejected every `reject_every` + iteration. + + reject_every : int, optional + Samples are rejected every `reject_every` iteration + max_iter : int, optional Maximum number of iterations to perform. @@ -129,9 +137,15 @@ def picard(X, fun='tanh', n_components=None, ortho=True, extended=None, n_iter : int Number of iterations taken to converge. This is returned only when return_n_iter is set to `True`. + + rejected : boolean array, shape (n_samples, ) + A mask indicating which samples have been rejected (1 = rejected). + Returned if `ll_reject` is not None. """ random_state = check_random_state(random_state) + rejection = (ll_reject is not None) + # Behaves like Fastica if ortho, standard infomax if not ortho if extended is None: extended = ortho @@ -200,7 +214,8 @@ def picard(X, fun='tanh', n_components=None, ortho=True, extended=None, kwargs = {'density': fun, 'm': m, 'max_iter': max_iter, 'tol': tol, 'lambda_min': lambda_min, 'ls_tries': ls_tries, 'verbose': verbose, 'ortho': ortho, 'extended': extended, - 'covariance': covariance} + 'covariance': covariance, 'll_reject': ll_reject, + 'reject_every': reject_every} with np.errstate(all='raise'): # So code breaks if overflow/Nans Y, W, infos = core_picard(X1, **kwargs) @@ -217,6 +232,24 @@ def picard(X, fun='tanh', n_components=None, ortho=True, extended=None, if not whiten: K = None n_iter = infos['n_iterations'] + if rejection: + rejected = infos['rejected'] + if return_X_mean: + if centering: + if return_n_iter: + return K, W, Y, X_mean, n_iter, rejected + else: + return K, W, Y, X_mean, rejected + else: + if return_n_iter: + return K, W, Y, np.zeros(p), n_iter, rejected + else: + return K, W, Y, np.zeros(p), rejected + else: + if return_n_iter: + return K, W, Y, n_iter, rejected + else: + return K, W, Y, rejected if return_X_mean: if centering: if return_n_iter: diff --git a/picard/tests/test_solver.py b/picard/tests/test_solver.py index 89c1cd7..725ef0e 100644 --- a/picard/tests/test_solver.py +++ b/picard/tests/test_solver.py @@ -96,6 +96,31 @@ def test_pre_fastica(): err_msg=err_msg) +def test_rejection(): + N, T = 3, 1000 + rng = np.random.RandomState(42) + fun = 'tanh' + extended = False + A = rng.randn(N, N) + S = rng.laplace(size=(N, T)) + X = np.dot(A, S) + # Add some strong artifacts + X[1, 10] = 1e2 + X[2, 20] = 1e2 + X[0, 20] = -1e2 + X[2, 100] = 1e2 + K, W, Y, rejected = picard(X, fun=fun, ortho=False, extended=extended, + random_state=0, ll_reject=15) + for loc in [10, 20, 100]: + assert rejected[loc] + assert np.sum(rejected) == 3 + assert_equal(W.shape, A.shape) + assert_equal(K.shape, A.shape) + WA = W.dot(K).dot(A) + WA = permute(WA) # Permute and scale + assert_allclose(WA, np.eye(N), rtol=0, atol=1e-1) + + def test_picard(): N, T = 3, 1000 rng = np.random.RandomState(42)