From 515539e3bb27e6be9d3d2e0d59234ebccf5ca6b2 Mon Sep 17 00:00:00 2001 From: Anuja Negi Date: Wed, 20 Nov 2024 11:57:13 +0100 Subject: [PATCH] alpha max fix for iterative L1 --- bsi_zoo/estimators.py | 19 ++++- bsi_zoo/run_benchmark.py | 176 ++++++++++++++++++++------------------- 2 files changed, 106 insertions(+), 89 deletions(-) diff --git a/bsi_zoo/estimators.py b/bsi_zoo/estimators.py index 25447c3..a26807d 100644 --- a/bsi_zoo/estimators.py +++ b/bsi_zoo/estimators.py @@ -253,7 +253,7 @@ def _gamma_map_opt( Parameters ---------- - M : array, shape=(n_sensors, n_times) + : array, shape=(n_sensors, n_times) Observation. G : array, shape=(n_sensors, n_sources) Forward operator. @@ -530,6 +530,14 @@ def gprime(w): return x +def norm_l2inf(A, n_orient, copy=True): + from math import sqrt + """L2-inf norm.""" + if A.size == 0: + return 0.0 + if copy: + A = A.copy() + return sqrt(np.max(groups_norm2(A, n_orient))) def iterative_L1(L, y, alpha=0.2, n_orient=1, max_iter=1000, max_iter_reweighting=10): """Iterative Type-I estimator with L1 regularizer. @@ -578,9 +586,16 @@ def gprime(w): grp_norms = np.sqrt(groups_norm2(w.copy(), n_orient)) return np.repeat(grp_norms, n_orient).ravel() + eps - alpha_max = abs(L.T.dot(y)).max() / len(L) + if n_orient==1: + alpha_max = abs(L.T.dot(y)).max() / len(L) + else: + n_dip_per_pos = 3 + alpha_max = norm_l2inf(np.dot(L.T, y), n_dip_per_pos) + alpha = alpha * alpha_max + # y->M + # L->gain x = _solve_reweighted_lasso( L, y, alpha, n_orient, weights, max_iter, max_iter_reweighting, gprime ) diff --git a/bsi_zoo/run_benchmark.py b/bsi_zoo/run_benchmark.py index b73855f..a8bd2d9 100644 --- a/bsi_zoo/run_benchmark.py +++ b/bsi_zoo/run_benchmark.py @@ -19,8 +19,10 @@ from bsi_zoo.config import get_leadfield_path n_jobs = 20 -nruns = 10 -spatial_cv = [False, True] +nruns = 1 +# spatial_cv = [False, True] +spatial_cv = [False] + # subjects = ["CC120166", "CC120264", "CC120313", "CC120309"] metrics = [ @@ -38,81 +40,81 @@ for do_spatial_cv in spatial_cv: for subject in subjects: - """Fixed orientation parameters for the benchmark""" - - orientation_type = "fixed" - data_args_I = { - # "n_sensors": [50], - "n_times": [10], - # "n_sources": [200], - "nnz": nnzs, - "cov_type": ["diag"], - "path_to_leadfield": [get_leadfield_path(subject, type=orientation_type)], - "orientation_type": [orientation_type], - "alpha": alpha_SNR, # this is actually SNR - } - - data_args_II = { - # "n_sensors": [50], - "n_times": [10], - # "n_sources": [200], - "nnz": nnzs, - "cov_type": ["full"], - "path_to_leadfield": [get_leadfield_path(subject, type=orientation_type)], - "orientation_type": [orientation_type], - "alpha": alpha_SNR, # this is actually SNR - } - - estimators = [ - (fake_solver, data_args_I, {"alpha": estimator_alphas_I}, {}), - (eloreta, data_args_I, {"alpha": estimator_alphas_II}, {}), - (iterative_L1, data_args_I, {"alpha": estimator_alphas_I}, {}), - (iterative_L2, data_args_I, {"alpha": estimator_alphas_I}, {}), - (iterative_sqrt, data_args_I, {"alpha": estimator_alphas_I}, {}), - (iterative_L1_typeII, data_args_II, {"alpha": estimator_alphas_I}, {}), - (iterative_L2_typeII, data_args_II, {"alpha": estimator_alphas_I}, {}), - # (gamma_map, data_args_II, {"alpha": estimator_alphas_I}, {"update_mode": 1}), - (gamma_map, data_args_II, {"alpha": estimator_alphas_II}, {"update_mode": 2}), - # (gamma_map, data_args_II, {"alpha": estimator_alphas_I}, {"update_mode": 3}), - ] - - df_results = [] - for estimator, data_args, estimator_args, estimator_extra_params in estimators: - benchmark = Benchmark( - estimator, - subject, - metrics, - data_args, - estimator_args, - random_state=42, - memory=memory, - n_jobs=n_jobs, - do_spatial_cv=do_spatial_cv, - estimator_extra_params=estimator_extra_params, - ) - results = benchmark.run(nruns=nruns) - df_results.append(results) - # save results - data_path = Path("bsi_zoo/data/updated_alpha_grid") - data_path.mkdir(exist_ok=True) - if do_spatial_cv: - FILE_NAME = f"{estimator}_{subject}_{data_args['orientation_type'][0]}_spatialCV_{time.strftime('%b-%d-%Y_%H%M', time.localtime())}.pkl" - else: - FILE_NAME = f"{estimator}_{subject}_{data_args['orientation_type'][0]}_{time.strftime('%b-%d-%Y_%H%M', time.localtime())}.pkl" - results.to_pickle(data_path / FILE_NAME) - - - df_results = pd.concat(df_results, axis=0) - - data_path = Path("bsi_zoo/data/ramen") - data_path.mkdir(exist_ok=True) - if do_spatial_cv: - FILE_NAME = f"benchmark_data_{subject}_{data_args['orientation_type'][0]}_spatialCV_{time.strftime('%b-%d-%Y_%H%M', time.localtime())}.pkl" - else: - FILE_NAME = f"benchmark_data_{subject}_{data_args['orientation_type'][0]}_{time.strftime('%b-%d-%Y_%H%M', time.localtime())}.pkl" - df_results.to_pickle(data_path / FILE_NAME) - - print(df_results) + # """Fixed orientation parameters for the benchmark""" + + # orientation_type = "fixed" + # data_args_I = { + # # "n_sensors": [50], + # "n_times": [10], + # # "n_sources": [200], + # "nnz": nnzs, + # "cov_type": ["diag"], + # "path_to_leadfield": [get_leadfield_path(subject, type=orientation_type)], + # "orientation_type": [orientation_type], + # "alpha": alpha_SNR, # this is actually SNR + # } + + # data_args_II = { + # # "n_sensors": [50], + # "n_times": [10], + # # "n_sources": [200], + # "nnz": nnzs, + # "cov_type": ["full"], + # "path_to_leadfield": [get_leadfield_path(subject, type=orientation_type)], + # "orientation_type": [orientation_type], + # "alpha": alpha_SNR, # this is actually SNR + # } + + # estimators = [ + # (fake_solver, data_args_I, {"alpha": estimator_alphas_I}, {}), + # (eloreta, data_args_I, {"alpha": estimator_alphas_II}, {}), + # (iterative_L1, data_args_I, {"alpha": estimator_alphas_I}, {}), + # (iterative_L2, data_args_I, {"alpha": estimator_alphas_I}, {}), + # (iterative_sqrt, data_args_I, {"alpha": estimator_alphas_I}, {}), + # (iterative_L1_typeII, data_args_II, {"alpha": estimator_alphas_I}, {}), + # (iterative_L2_typeII, data_args_II, {"alpha": estimator_alphas_I}, {}), + # # (gamma_map, data_args_II, {"alpha": estimator_alphas_I}, {"update_mode": 1}), + # (gamma_map, data_args_II, {"alpha": estimator_alphas_II}, {"update_mode": 2}), + # # (gamma_map, data_args_II, {"alpha": estimator_alphas_I}, {"update_mode": 3}), + # ] + + # df_results = [] + # for estimator, data_args, estimator_args, estimator_extra_params in estimators: + # benchmark = Benchmark( + # estimator, + # subject, + # metrics, + # data_args, + # estimator_args, + # random_state=42, + # memory=memory, + # n_jobs=n_jobs, + # do_spatial_cv=do_spatial_cv, + # estimator_extra_params=estimator_extra_params, + # ) + # results = benchmark.run(nruns=nruns) + # df_results.append(results) + # # save results + # data_path = Path("bsi_zoo/data/updated_alpha_grid") + # data_path.mkdir(exist_ok=True) + # if do_spatial_cv: + # FILE_NAME = f"{estimator}_{subject}_{data_args['orientation_type'][0]}_spatialCV_{time.strftime('%b-%d-%Y_%H%M', time.localtime())}.pkl" + # else: + # FILE_NAME = f"{estimator}_{subject}_{data_args['orientation_type'][0]}_{time.strftime('%b-%d-%Y_%H%M', time.localtime())}.pkl" + # results.to_pickle(data_path / FILE_NAME) + + + # df_results = pd.concat(df_results, axis=0) + + # data_path = Path("bsi_zoo/data/ramen") + # data_path.mkdir(exist_ok=True) + # if do_spatial_cv: + # FILE_NAME = f"benchmark_data_{subject}_{data_args['orientation_type'][0]}_spatialCV_{time.strftime('%b-%d-%Y_%H%M', time.localtime())}.pkl" + # else: + # FILE_NAME = f"benchmark_data_{subject}_{data_args['orientation_type'][0]}_{time.strftime('%b-%d-%Y_%H%M', time.localtime())}.pkl" + # df_results.to_pickle(data_path / FILE_NAME) + + # print(df_results) """ Free orientation parameters for the benchmark """ @@ -142,22 +144,22 @@ if spatial_cv: # currently no support for type II methods estimators = [ - (fake_solver, data_args_I, {"alpha": estimator_alphas_I}, {}), + # (fake_solver, data_args_I, {"alpha": estimator_alphas_I}, {}), (iterative_L1, data_args_I, {"alpha": estimator_alphas_I}, {}), (iterative_L2, data_args_I, {"alpha": estimator_alphas_I}, {}), (iterative_sqrt, data_args_I, {"alpha": estimator_alphas_I}, {}), ] else: estimators = [ - (fake_solver, data_args_I, {"alpha": estimator_alphas_I}, {}), - (eloreta, data_args_I, {"alpha": estimator_alphas_II}, {}), + # (fake_solver, data_args_I, {"alpha": estimator_alphas_I}, {}), + # (eloreta, data_args_I, {"alpha": estimator_alphas_II}, {}), (iterative_L1, data_args_I, {"alpha": estimator_alphas_I}, {}), - (iterative_L2, data_args_I, {"alpha": estimator_alphas_I}, {}), - (iterative_sqrt, data_args_I, {"alpha": estimator_alphas_I}, {}), - (iterative_L1_typeII, data_args_II, {"alpha": estimator_alphas_I}, {}), - (iterative_L2_typeII, data_args_II, {"alpha": estimator_alphas_I}, {}), + # (iterative_L2, data_args_I, {"alpha": estimator_alphas_I}, {}), + # (iterative_sqrt, data_args_I, {"alpha": estimator_alphas_I}, {}), + # (iterative_L1_typeII, data_args_II, {"alpha": estimator_alphas_I}, {}), + # (iterative_L2_typeII, data_args_II, {"alpha": estimator_alphas_I}, {}), # (gamma_map, data_args_II, {"alpha": estimator_alphas_I}, {"update_mode": 1}), - (gamma_map, data_args_II, {"alpha": estimator_alphas_II}, {"update_mode": 2}), + # (gamma_map, data_args_II, {"alpha": estimator_alphas_II}, {"update_mode": 2}), # (gamma_map, data_args_II, {"alpha": estimator_alphas_I}, {"update_mode": 3}), ] @@ -178,7 +180,7 @@ results = benchmark.run(nruns=nruns) df_results.append(results) # save results - data_path = Path("bsi_zoo/data/free2") + data_path = Path("bsi_zoo/data/free3") data_path.mkdir(exist_ok=True) if do_spatial_cv: @@ -189,7 +191,7 @@ df_results = pd.concat(df_results, axis=0) - data_path = Path("bsi_zoo/data/free2") + data_path = Path("bsi_zoo/data/free3") data_path.mkdir(exist_ok=True) if do_spatial_cv: FILE_NAME = f"benchmark_data_{subject}_{data_args['orientation_type'][0]}_spatialCV_{time.strftime('%b-%d-%Y_%H%M', time.localtime())}.pkl"