diff --git a/clar/data/artificial.py b/clar/data/artificial.py index 1dafff9..65ed596 100644 --- a/clar/data/artificial.py +++ b/clar/data/artificial.py @@ -107,8 +107,7 @@ def get_S_star( if noise_type == "Gaussian_iid": S_star = np.eye(n_channels) elif noise_type == "Gaussian_multivariate": - vect = rho_noise ** np.arange(n_channels) - S_star = toeplitz(vect, vect) + S_star = toeplitz(rho_noise ** np.arange(n_channels)) else: raise ValueError("Unknown noise type %s" % noise_type) return S_star @@ -146,7 +145,7 @@ def get_dictionary( elif dictionary_type == 'Gaussian': X = rng.randn(n_channels, n_sources) else: - raise NotImplementedError("No dictionary '{}' in maxsparse" + raise NotImplementedError("No dictionary '{}' in clar" .format(dictionary_type)) normalize(X) return X @@ -181,23 +180,6 @@ def get_toeplitz_dictionary( The dictionary. """ rng = check_random_state(seed) - vect = rho ** np.arange(n_sources) - covar = toeplitz(vect, vect) + covar = toeplitz(rho ** np.arange(n_sources)) X = rng.multivariate_normal(np.zeros(n_sources), covar, n_channels) return X - - -def decimate(M, n_channels, axis, seed): - if n_channels in (M.shape[0], -1): - return M - - n_channels_max = M.shape[0] - rng = check_random_state(seed) - - to_choose = rng.choice(np.arange(n_channels_max), n_channels) - to_choose.sort() - if axis.__contains__(1): - M = M[:, to_choose] - if axis.__contains__(0): - M = M[to_choose, :] - return M diff --git a/clar/duality_gap.py b/clar/duality_gap.py index b4329d1..8a429f5 100644 --- a/clar/duality_gap.py +++ b/clar/duality_gap.py @@ -1,6 +1,5 @@ import numpy as np from numpy.linalg import norm, slogdet -from numba import njit from clar.utils import l_2_inf from clar.utils import l_2_1 @@ -24,7 +23,6 @@ def get_p_obj_mrce( return p_obj -@njit def get_p_obj_mtl(R, B, alpha): n_sensors, n_times = R.shape p_obj = (R ** 2).sum() / (2 * n_times * n_sensors) \ @@ -32,7 +30,6 @@ def get_p_obj_mtl(R, B, alpha): return p_obj -@njit def get_d_obj_mtl(Y, Theta, alpha): n_sensors, n_times = Y.shape d_obj = alpha * (Theta * Y).sum() - \ @@ -40,7 +37,6 @@ def get_d_obj_mtl(Y, Theta, alpha): return d_obj -@njit def get_feasible_theta_mtl(R, X, alpha): n_sensors, n_times = R.shape scaling_factor = l_2_inf(X.T @ R) @@ -48,7 +44,6 @@ def get_feasible_theta_mtl(R, X, alpha): return R / scaling_factor -@njit def get_p_obj_me(R_all_epochs, B, S_inv_R, S_trace, alpha): n_epochs, n_channels, n_times = R_all_epochs.shape p_obj = (R_all_epochs * S_inv_R).sum() @@ -58,7 +53,6 @@ def get_p_obj_me(R_all_epochs, B, S_inv_R, S_trace, alpha): return p_obj -@njit def get_d_obj_me(all_epochs, Theta, sigma_min, alpha): n_epochs, n_channels, n_times = all_epochs.shape d_obj = alpha * (all_epochs * Theta).sum() / n_epochs @@ -69,7 +63,6 @@ def get_d_obj_me(all_epochs, Theta, sigma_min, alpha): return d_obj -@njit def get_d_obj(Y, Theta, sigma_min, alpha): n_channels, n_times = Y.shape d_obj = alpha * (Y * Theta).sum() @@ -78,7 +71,6 @@ def get_d_obj(Y, Theta, sigma_min, alpha): return d_obj -@njit def get_p_obj(R, B, S_trace, alpha, S_inv_R): n_channels, n_times = R.shape p_obj = (R * S_inv_R).sum() / (2. * n_channels * n_times) @@ -87,7 +79,6 @@ def get_p_obj(R, B, S_trace, alpha, S_inv_R): return p_obj -@njit def get_feasible_theta(X, alpha, S_inv_R): n_channels, n_times = S_inv_R.shape scaling_factor = max( @@ -97,7 +88,6 @@ def get_feasible_theta(X, alpha, S_inv_R): return S_inv_R / scaling_factor -@njit def get_feasible_theta_me(X, alpha, S_inv_R): n_epochs, n_channels, n_times = S_inv_R.shape S_inv_R_mean = np.zeros((n_channels, n_times), dtype=np.float64) @@ -118,7 +108,6 @@ def get_feasible_theta_me(X, alpha, S_inv_R): return S_inv_R / scaling_factor -@njit def get_duality_gap_mtl(X, Y, B, alpha): R = Y - X @ B p_obj = get_p_obj_mtl(R, B, alpha) @@ -127,7 +116,6 @@ def get_duality_gap_mtl(X, Y, B, alpha): return p_obj, d_obj -@njit def get_duality_gap( R, X, Y, B, S_trace, S_inv_R, sigma_min, alpha): p_obj = get_p_obj(R, B, S_trace, alpha, S_inv_R) @@ -136,7 +124,6 @@ def get_duality_gap( return p_obj, d_obj -@njit def get_duality_gap_me( X, all_epochs, B, S_trace, S_inv, sigma_min, alpha): diff --git a/clar/solvers.py b/clar/solvers.py index be15b3a..1c2893d 100644 --- a/clar/solvers.py +++ b/clar/solvers.py @@ -37,16 +37,8 @@ def get_path( # save the results mask = np.abs(B_hat).sum(axis=1) != 0 str_pourcentage_alpha = '%0.10f' % pourcentage_alpha - if pb_name == "MTLME": - n_sources = X.shape[1] - n_epochs, _, n_times = measurement.shape - B_reshaped = B_hat.reshape((n_sources, n_epochs, n_times)) - B_reshaped = B_reshaped.mean(axis=1) - dict_masks[str_pourcentage_alpha] = mask - dict_dense_Bs[str_pourcentage_alpha] = B_reshaped[mask, :] - else: - dict_masks[str_pourcentage_alpha] = mask - dict_dense_Bs[str_pourcentage_alpha] = B_hat[mask, :] + dict_masks[str_pourcentage_alpha] = mask + dict_dense_Bs[str_pourcentage_alpha] = B_hat[mask, :] assert len(dict_dense_Bs.keys()) == len(list_pourcentage_alpha) return dict_masks, dict_dense_Bs @@ -89,7 +81,7 @@ def solver( S is updated every S times. pb_name: str choose the problem you want to solve between - "MTL", "MTLME", "SGCL", "CLAR" and "mrce" + "MTL", "SGCL", "CLAR" and "mrce" use_accel: bool States if you want to use accelratio while computing the dual. heur_stop: bool @@ -110,11 +102,7 @@ def solver( print("--------- %s -----------------" % pb_name) if B0 is None: - if pb_name != "MTLME": - B = np.zeros((n_sources, n_times), dtype=float) - else: - n_epochs, _, n_times = all_epochs.shape - B = np.zeros((n_sources, n_times * n_epochs), dtype=float) + B = np.zeros((n_sources, n_times), dtype=float) else: B = B0.copy().astype(np.float64) @@ -125,15 +113,6 @@ def solver( observations = all_epochs[None, :, :] elif pb_name in ("CLAR", "mrce"): observations = all_epochs - elif pb_name == "MTLME": - if all_epochs.ndim != 3: - raise ValueError( - "Wrong number of dimensions, expected 2, " - "got %d " % all_epochs.ndim) - observations = all_epochs.transpose((1, 0, 2)) - observations = observations.reshape(observations.shape[0], -1) - observations = observations.reshape((1, *observations.shape)) - n_epochs, _, n_times = all_epochs.shape else: raise ValueError("Unknown solver %s" % pb_name) @@ -166,11 +145,9 @@ def solver_( Y += all_epochs[l, :, :] Y2 /= n_epochs Y /= n_epochs - elif pb_name in ("MTL", "SGCL", "MTLME"): + elif pb_name in ("MTL", "SGCL"): Y = all_epochs[0] Y2 = None - elif pb_name == "MTLME": - Y = all_epochs if use_accel: K = 6 @@ -192,7 +169,7 @@ def solver_( primal_first, _ = get_duality_gap( Y, X, Y, B_first, S_trace_first, S_inv_R, sigma_min, alpha) - elif pb_name in("MTL", "MTLME"): + elif pb_name in "MTL": primal_first, _ = get_duality_gap_mtl( X, Y, B_first, alpha) elif pb_name == "mrce": @@ -231,7 +208,7 @@ def solver_( S_trace, S_inv = clp_sqrt(ZZT, sigma_min) S_inv_R = np.asfortranarray(S_inv @ R) S_inv_X = S_inv @ X - elif pb_name in ("MTL", "MTLME"): + elif pb_name in "MTL": # this else case is for MTL # dummy variables for njit to work: S_trace = n_sensors @@ -289,7 +266,7 @@ def solver_( if verbose: print("gap_acc: %.2e" % (p_obj - d_obj_acc)) gaps_acc.append(p_obj - d_obj_acc) - elif pb_name in ("MTL", "MTLME"): + elif pb_name in "MTL": p_obj, d_obj = get_duality_gap_mtl( X, Y, B, alpha) elif pb_name == "mrce": @@ -333,7 +310,7 @@ def update_S(Y, X, B, Y2, sigma_min, pb_name): Z = Y - X @ B ZZT = Z @ Z.T / n_times S_trace, S_inv = clp_sqrt(ZZT, sigma_min) - elif pb_name in ("MTL", "MTLME"): + elif pb_name in "MTL": # this else case is for MTL # dummy variables for njit to work: S_trace = n_sensors @@ -349,7 +326,7 @@ def update_B( n_sensors, n_times = Y.shape n_sources = X.shape[1] - is_not_MTL = pb_name not in ("MTL", "MTLME") + is_not_MTL = pb_name not in "MTL" active_set = np.ones(n_sources) diff --git a/clar/tests/test_clar.py b/clar/tests/test_clar.py deleted file mode 100644 index 70e2e81..0000000 --- a/clar/tests/test_clar.py +++ /dev/null @@ -1,44 +0,0 @@ -import numpy as np - -from clar.solvers import solver -from clar.data.artificial import get_data_me -from clar.utils import get_sigma_min, get_alpha_max - - -def test_clar1(): - test_clar( - n_channels=20, n_times=30, n_epochs=50, - n_sources=50, n_active=3, n_iter=10**4, tol=1e-7) - - -def test_clar( - noise_type="Gaussian_iid", rho_noise=0.3, - SNR=0.5, n_channels=20, n_times=30, n_sources=10, n_epochs=50, - n_active=3, gap_freq=50, active_set_freq=1, S_freq=10, - n_iter=10**4, alpha_under_alpha_max=0.2, tol=1e-7): - X, all_epochs, _, _ = get_data_me( - dictionary_type="Gaussian", noise_type=noise_type, - n_epochs=n_epochs, n_channels=n_channels, n_times=n_times, - n_sources=n_sources, n_active=n_active, rho_noise=rho_noise, - SNR=SNR - ) - Y = np.mean(all_epochs, axis=0) - sigma_min = get_sigma_min(Y) - alpha_max = get_alpha_max(X, all_epochs, sigma_min, pb_name="CLAR") - alpha = alpha_max * \ - alpha_under_alpha_max - - print("alpha = %.2e" % alpha) - print("alpha = %.2e" % alpha) - print("sigma_min = %.2e" % sigma_min) - gaps_me = solver( - X, all_epochs, alpha, sigma_min, B0=None, - n_iter=n_iter, - gap_freq=gap_freq, active_set_freq=active_set_freq, - S_freq=S_freq, tol=tol, pb_name="CLAR")[-1] - gap_me = gaps_me[-1] - assert gap_me < tol - - -if __name__ == '__main__': - test_clar1() diff --git a/clar/tests/test_mrce.py b/clar/tests/test_mrce.py deleted file mode 100644 index 44b51ca..0000000 --- a/clar/tests/test_mrce.py +++ /dev/null @@ -1,58 +0,0 @@ -import numpy as np -from clar.solvers import solver, update_sigma_glasso -from clar.utils import ( - get_alpha_max, get_sigma_min) -from clar.data.artificial import get_data_me - - -def test_update_sigma_glasso(): - rho_noise = 0.8 - SNR = 1 - n_epochs, n_channels, n_sources, n_times = 5, 20, 10, 30 - - all_epochs = get_data_me( - dictionary_type="Gaussian", noise_type="Gaussian_multivariate", - n_epochs=n_epochs, n_channels=n_channels, n_times=n_times, - n_sources=n_sources, n_active=3, rho_noise=rho_noise, - SNR=SNR)[1] - - emp_cov = np.zeros((n_channels, n_channels)) - for i in range(n_epochs): - emp_cov += all_epochs[i, :, :] @ all_epochs[i, :, :].T - emp_cov /= (n_times * n_epochs) - - alpha_prec = 0.001 - - update_sigma_glasso( - emp_cov, alpha_prec, enet_tol=1e-4, max_iter=100) - - -def test_mrce(): - rho_noise = 0.6 - SNR = 1 - n_epochs, n_channels, n_sources, n_times = 5, 20, 10, 30 - pb_name = "mrce" - tol = 1e-4 - X, all_epochs = get_data_me( - dictionary_type="Gaussian", noise_type="Gaussian_multivariate", - n_epochs=n_epochs, n_channels=n_channels, n_times=n_times, - n_sources=n_sources, n_active=3, rho_noise=rho_noise, - SNR=SNR)[:2] - - alpha_Sigma_inv = 0.01 - - Y = np.mean(all_epochs, axis=0) - sigma_min = get_sigma_min(Y) - alpha_max = get_alpha_max( - X, all_epochs, sigma_min, pb_name=pb_name, - alpha_Sigma_inv=alpha_Sigma_inv) - - alpha = alpha_max * 0.9 - - solver( - X, all_epochs, alpha, sigma_min, B0=None, - tol=tol, pb_name=pb_name, n_iter=1000, alpha_Sigma_inv=alpha_Sigma_inv) - - -if __name__ == '__main__': - test_mrce() diff --git a/clar/tests/test_mtl.py b/clar/tests/test_mtl.py deleted file mode 100644 index 758f113..0000000 --- a/clar/tests/test_mtl.py +++ /dev/null @@ -1,71 +0,0 @@ -import numpy as np -from clar.solvers import solver -from clar.utils import ( - get_alpha_max, get_sigma_min) -from clar.data.artificial import get_data_me - - -def test_mtl(): - rho_noise = 0.3 - SNR = 1 - n_epochs, n_channels, n_sources, n_times = 5, 20, 10, 30 - pb_name = "MTL" - tol = 1e-7 - - X, all_epochs, _, _ = get_data_me( - dictionary_type="Gaussian", noise_type="Gaussian_iid", - n_epochs=n_epochs, n_channels=n_channels, n_times=n_times, - n_sources=n_sources, n_active=3, rho_noise=rho_noise, - SNR=SNR) - - Y = np.mean(all_epochs, axis=0) - sigma_min = get_sigma_min(Y) - alpha_max = get_alpha_max(X, Y, sigma_min, pb_name=pb_name) - - alpha_div = 5 - alpha = alpha_max / alpha_div - - B_mtl, _, E, gaps = solver( - X, Y, alpha, sigma_min, B0=None, - tol=tol, pb_name=pb_name, n_iter=10000) - gap = gaps[-1] - np.testing.assert_array_less(gap, tol) - - _, _, E, gaps = solver( - X, Y, alpha, sigma_min, B0=B_mtl, - tol=tol, pb_name=pb_name, n_iter=10000) - np.testing.assert_equal(len(E), 2) - gap = gaps[-1] - np.testing.assert_array_less(gap, tol * E[0]) - - -def test_mtl_me(): - rho_noise = 0.3 - SNR = 1 - n_epochs, n_channels, n_sources, n_times = 5, 20, 10, 30 - pb_name = "MTLME" - tol = 1e-7 - - X, all_epochs, _, _ = get_data_me( - dictionary_type="Gaussian", noise_type="Gaussian_iid", - n_epochs=n_epochs, n_channels=n_channels, n_times=n_times, - n_sources=n_sources, n_active=3, rho_noise=rho_noise, - SNR=SNR) - - Y = np.mean(all_epochs, axis=0) - sigma_min = get_sigma_min(Y) - - alpha_max = get_alpha_max(X, all_epochs, sigma_min, pb_name=pb_name) - - alpha_div = 1.1 - alpha = alpha_max / alpha_div - - gap = solver( - X, all_epochs, alpha, sigma_min, B0=None, - tol=tol, pb_name=pb_name, n_iter=10000)[-1] - np.testing.assert_array_less(gap, tol) - - -if __name__ == '__main__': - test_mtl() - test_mtl_me() diff --git a/clar/tests/test_sgcl.py b/clar/tests/test_sgcl.py deleted file mode 100644 index 811bfae..0000000 --- a/clar/tests/test_sgcl.py +++ /dev/null @@ -1,49 +0,0 @@ -import numpy as np - -from clar.solvers import solver -from clar.utils import get_alpha_max, get_sigma_min -from clar.data.artificial import get_data_me - - -def test_sgcl1(): - tests_sgcl( - n_channels=30, n_times=30, - n_sources=10, n_active=1, n_iter=10**7, tol=10**-4) - - -def tests_sgcl( - noise_type="Gaussian_iid", rho_noise=0.3, - SNR=0.5, n_channels=20, n_times=30, n_sources=10, n_active=3, - gap_freq=100, active_set_freq=1, S_freq=10, n_iter=10**6, - p_alpha_max=0.9, tol=10**-4): - - X, all_epochs, _, _ = get_data_me( - dictionary_type="Gaussian", noise_type=noise_type, - n_channels=n_channels, n_times=n_times, n_sources=n_sources, - n_active=n_active, rho_noise=rho_noise, SNR=SNR, - n_epochs=50) - - Y = all_epochs.mean(axis=0) - sigma_min = get_sigma_min(Y) - alpha_max = get_alpha_max(X, Y, sigma_min, "SGCL") - alpha = alpha_max * p_alpha_max - print("alpha = %.2e" % alpha) - print("sigma_min = %.2e" % sigma_min) - - all_epochs = np.zeros((1, *Y.shape)) - all_epochs[0] = Y - - _, _, E, (gaps, gaps_accel) = solver( - X, Y, alpha, sigma_min, B0=None, n_iter=n_iter, - gap_freq=gap_freq, active_set_freq=active_set_freq, - S_freq=S_freq, pb_name="SGCL", use_accel=True, tol=tol, - verbose=True) - - log_gap = np.log10(gaps[-1]) - log_gap_accel = np.log10(gaps_accel[-1]) - assert log_gap_accel < np.log10(tol) * E[0] or \ - log_gap < np.log10(tol) * E[0] - - -if __name__ == '__main__': - test_sgcl1() diff --git a/clar/tests/test_solvers.py b/clar/tests/test_solvers.py new file mode 100644 index 0000000..85daddb --- /dev/null +++ b/clar/tests/test_solvers.py @@ -0,0 +1,140 @@ +import numpy as np +import pytest + +from clar.solvers import solver, get_path +from clar.utils import get_alpha_max, get_sigma_min +from clar.data.artificial import get_data_me + + +def test_mtl(): + rho_noise = 0.3 + SNR = 1 + n_epochs, n_channels, n_sources, n_times = 5, 20, 10, 30 + pb_name = "MTL" + tol = 1e-10 + + X, all_epochs, _, _ = get_data_me( + dictionary_type="Gaussian", noise_type="Gaussian_iid", + n_epochs=n_epochs, n_channels=n_channels, n_times=n_times, + n_sources=n_sources, n_active=3, rho_noise=rho_noise, + SNR=SNR) + + Y = np.mean(all_epochs, axis=0) + sigma_min = get_sigma_min(Y) + alpha_max = get_alpha_max(X, Y, sigma_min, pb_name=pb_name) + + alpha_div = 5 + alpha = alpha_max / alpha_div + + B_mtl, _, E, gaps = solver( + X, Y, alpha, sigma_min, B0=None, + tol=tol, pb_name=pb_name, n_iter=10000) + np.testing.assert_array_less(gaps[-1], tol) + + _, _, E, gaps = solver( + X, Y, alpha, sigma_min, B0=B_mtl, + tol=tol, pb_name=pb_name, n_iter=10000) + np.testing.assert_equal(len(E), 2) + np.testing.assert_array_less(gaps[-1], tol * E[0]) + + +@pytest.mark.parametrize("n_sources", [10, 15, 20]) +def tests_sgcl(n_sources): + rho_noise = 0.3 + SNR = 0.5 + n_channels = 20 + n_times = 30 + n_active = 3 + n_iter = 10**6 + tol = 10**-4 + + X, all_epochs, _, _ = get_data_me( + n_channels=n_channels, n_times=n_times, n_sources=n_sources, + n_active=n_active, rho_noise=rho_noise, SNR=SNR, + n_epochs=50) + + Y = all_epochs.mean(axis=0) + sigma_min = get_sigma_min(Y) + alpha_max = get_alpha_max(X, Y, sigma_min, "SGCL") + alpha = alpha_max * 0.1 + + all_epochs = np.zeros((1, *Y.shape)) + all_epochs[0] = Y + + _, _, E, (gaps, gaps_accel) = solver( + X, Y, alpha, sigma_min, B0=None, n_iter=n_iter, + pb_name="SGCL", use_accel=True, tol=tol, + verbose=True) + + np.testing.assert_array_less( + np.minimum(gaps[-1], gaps_accel[-1]), tol * E[0]) + + +@pytest.mark.parametrize("n_sources", [10, 15, 20]) +def test_clar(n_sources): + rho_noise = 0.3 + SNR = 0.5 + n_channels = 20 + n_times = 30 + n_epochs = 50 + n_active = 3 + n_iter = 10**4 + p_alpha = 0.2 + tol = 1e-7 + X, all_epochs, _, _ = get_data_me( + dictionary_type="Toeplitz", + n_epochs=n_epochs, n_channels=n_channels, n_times=n_times, + n_sources=n_sources, n_active=n_active, rho_noise=rho_noise, + SNR=SNR + ) + Y = np.mean(all_epochs, axis=0) + sigma_min = get_sigma_min(Y) + alpha_max = get_alpha_max(X, all_epochs, sigma_min, pb_name="CLAR") + alpha = alpha_max * p_alpha + + gaps_me = solver( + X, all_epochs, alpha, sigma_min, B0=None, + n_iter=n_iter, tol=tol, pb_name="CLAR")[-1] + gap_me = gaps_me[-1] + np.testing.assert_array_less(gap_me, tol) + + p_alphas = np.geomspace(1, 0.1, 10) + dict_masks = get_path( + X, all_epochs, p_alphas, alpha_max, sigma_min, + n_iter=10**4, tol=10**-4)[0] + + old_size_supp = 0 + for supp in dict_masks.values(): + size_supp = supp.sum() + np.testing.assert_array_less(old_size_supp - size_supp, 1) + old_size_supp = size_supp + + +def test_mrce(): + rho_noise = 0.6 + SNR = 1 + n_epochs, n_channels, n_sources, n_times = 5, 20, 10, 30 + pb_name = "mrce" + tol = 1e-4 + X, all_epochs = get_data_me( + dictionary_type="Gaussian", noise_type="Gaussian_multivariate", + n_epochs=n_epochs, n_channels=n_channels, n_times=n_times, + n_sources=n_sources, n_active=3, rho_noise=rho_noise, + SNR=SNR)[:2] + + alpha_Sigma_inv = 0.01 + + Y = np.mean(all_epochs, axis=0) + sigma_min = get_sigma_min(Y) + alpha_max = get_alpha_max( + X, all_epochs, sigma_min, pb_name=pb_name, + alpha_Sigma_inv=alpha_Sigma_inv) + + alpha = alpha_max * 0.1 + + Es = solver( + X, all_epochs, alpha, sigma_min, B0=None, + tol=tol, pb_name=pb_name, n_iter=1000, + alpha_Sigma_inv=alpha_Sigma_inv)[-2] + + np.testing.assert_array_less(Es[-1] - Es[-2], 1e-10) diff --git a/clar/utils.py b/clar/utils.py index 43c66fd..6c179e2 100644 --- a/clar/utils.py +++ b/clar/utils.py @@ -1,8 +1,6 @@ import numpy as np from numpy.linalg import norm - from numba import njit - from sklearn.covariance import graphical_lasso @@ -20,27 +18,6 @@ def sqrtm(ZZT): return eigvecs @ (eigvals * eigvecs.T) -def get_S_Sinv(ZZT, sigma_min=1e-6): - """Take the square root and inverse of square root of a - symmetric definite matrix. - - Output: (float, np.array, shape (n_sensors, n_sensors)) - (trace of Sigma updated, inverse of Sigma updated) - """ - eigvals, eigvecs = np.linalg.eigh(ZZT) - eigvals = np.maximum(0, eigvals) - eigvals = np.sqrt(eigvals) - div_eigvals = 1 / eigvals - mask = (eigvals < sigma_min * eigvals.max()) - - print(eigvals) - print('Number of eigvals clipped: %d' % mask.sum()) - div_eigvals[mask] = 0 - eigvals = np.expand_dims(eigvals, axis=1) - div_eigvals = np.expand_dims(div_eigvals, axis=1) - return eigvecs @ (eigvals * eigvecs.T), eigvecs @ (div_eigvals * eigvecs.T) - - @njit def BST(u, tau): """ @@ -153,12 +130,6 @@ def l_2_1(A): # return norm(A, axis=1, ord=2).sum() -def get_alpha_max_mtl(X, Y): - n_sensors, n_times = Y.shape - alpha_max = l_2_inf(X.T @ Y) / (n_times * n_sensors) - return alpha_max - - def get_emp_cov(R): if R.ndim != 3: raise ValueError( @@ -223,80 +194,12 @@ def get_alpha_max(X, observation, sigma_min, pb_name, alpha_Sigma_inv=None): Sigma_inv = graphical_lasso( emp_cov, alpha_Sigma_inv, max_iter=10 ** 6)[-1] alpha_max = l_2_inf(X.T @ Sigma_inv @ Y) / (n_channels * n_times) - elif pb_name == "glasso": - assert observation.ndim == 2 - assert alpha_Sigma_inv is not None - emp_cov = observation @ observation.T / n_times - Sigma_inv = graphical_lasso(emp_cov, alpha_Sigma_inv)[-1] - alpha_max = l_2_inf(X.T @ Sigma_inv @ Y) / (n_channels * n_times) - elif pb_name == "mrce": - assert observation.ndim == 2 - assert alpha_Sigma_inv is not None - emp_cov = observation @ observation.T / n_times - Sigma_inv = graphical_lasso( - emp_cov, alpha_Sigma_inv, max_iter=10 ** 6)[-1] - alpha_max = np.abs(X.T @ Sigma_inv @ Y).max() / (n_channels * n_times) else: raise NotImplementedError( "No solver '{}' in sgcl".format(pb_name)) return alpha_max -def get_alpha_max_sgcl(X, Y, sigma_min): - """Function to compute the maximal alpha before obtaining all zeros. - """ - n_sensors, n_times = Y.shape - _, Sigma_max_inv = clp_sqrt( - Y @ Y.T / n_times, sigma_min) - result = l_2_inf(X.T @ Sigma_max_inv @ Y) - result /= (n_sensors * n_times) - return result - - -def get_alpha_max_me(X, all_epochs, sigma_min): - """Function to compute the maximal alpha before obtaining all zeros. - """ - n_epochs, n_sensors, n_times = all_epochs.shape - Y = all_epochs.mean(axis=0) - - cov_Yl = 0 - for l in range(n_epochs): - cov_Yl += all_epochs[l, :, :] @ all_epochs[l, :, :].T - cov_Yl /= (n_epochs * n_times) - - _, Sigma_max_inv = clp_sqrt( - cov_Yl, sigma_min) - result = l_2_inf(X.T @ Sigma_max_inv @ Y) - result /= (n_sensors * n_times) - return result - - def get_sigma_min(Y): sigma_min = norm(Y, ord='fro') / (np.sqrt(Y.shape[1] * Y.shape[0]) * 1000) return sigma_min - - -def get_relative_log_res( - X, Y, B_star, Sigma_inv_star, B_hat, Sigma_inv_hat, me=False): - if me: - res = get_norm_res_me(X, Y, B_hat, Sigma_inv_hat) / \ - get_norm_res_me(X, Y, B_star, Sigma_inv_star) - else: - res = get_norm_res(X, Y, B_hat, Sigma_inv_hat) / \ - get_norm_res(X, Y, B_star, Sigma_inv_star) - return np.log10(res) - - -def get_norm_res(X, Y, B, Sigma_inv, ord='fro'): - R = Y - X @ B - res = norm(R.T @ (Sigma_inv @ R), ord=ord) - return res - - -def get_norm_res_me(X, all_epochs, B, Sigma_inv, ord='fro'): - R = all_epochs - X @ B - n_epochs = R.shape[0] - res = 0 - for l in range(n_epochs): - res += norm(R[l, :, :].T @ Sigma_inv @ R[l, :, :], ord=ord) - return res