mind-inria
diff --git a/‎docs/src/api.rst‎
Lines changed: 1 addition & 3 deletions b/‎docs/src/api.rst‎
Lines changed: 1 addition & 3 deletions
diff --git a/‎examples/plot_2D_simulation_example.py‎
Lines changed: 53 additions & 53 deletions b/‎examples/plot_2D_simulation_example.py‎
Lines changed: 53 additions & 53 deletions
diff --git a/‎examples/plot_fmri_data_example.py‎
Lines changed: 41 additions & 26 deletions b/‎examples/plot_fmri_data_example.py‎
Lines changed: 41 additions & 26 deletions
@@ -29,6 +29,7 @@ Feature Importance Classes
    PFI
    D0CRT
    ModelXKnockoff
+   DesparsifiedLasso
 
 Feature Importance functions
 ============================
@@ -39,9 +40,6 @@ Feature Importance functions
 
    clustered_inference
    clustered_inference_pvalue
-   desparsified_lasso
-   desparsified_lasso_pvalue
-   desparsified_group_lasso_pvalue
    ensemble_clustered_inference
    ensemble_clustered_inference_pvalue
 
 
@@ -46,22 +46,6 @@
 .. footbibliography::
 
 """
-import matplotlib.pyplot as plt
-import numpy as np
-from sklearn.cluster import FeatureAgglomeration
-from sklearn.feature_extraction import image
-from sklearn.preprocessing import StandardScaler
-
-from hidimstat._utils.scenario import multivariate_simulation_spatial
-from hidimstat.desparsified_lasso import desparsified_lasso, desparsified_lasso_pvalue
-from hidimstat.ensemble_clustered_inference import (
-    clustered_inference,
-    clustered_inference_pvalue,
-    ensemble_clustered_inference,
-    ensemble_clustered_inference_pvalue,
-)
-from hidimstat.statistical_tools.p_values import zscore_from_pval
-
 # %%
 # Generating the data
 # -------------------
@@ -71,6 +55,9 @@
 # example.
 
 # simulation parameters
+
+from hidimstat._utils.scenario import multivariate_simulation_spatial
+
 n_samples = 100
 shape = (40, 40)
 n_features = shape[1] * shape[0]
@@ -83,6 +70,7 @@
     n_samples, shape, roi_size, signal_noise_ratio, smooth_X, seed=0
 )
 
+
 # %%
 # Choosing inference parameters
 # -----------------------------
@@ -110,7 +98,8 @@
 delta = 6
 
 # number of worker
-n_jobs = 3
+n_jobs = 4
+
 
 # %%
 # Computing z-score thresholds for support estimation
@@ -125,12 +114,15 @@
 # consists in dividing by the number of clusters.
 
 
+from hidimstat.statistical_tools.p_values import zscore_from_pval
+
 # computing the z-score thresholds for feature selection
 correction_no_cluster = 1.0 / n_features
 correction_cluster = 1.0 / n_clusters
 thr_c = zscore_from_pval((fwer_target / 2) * correction_cluster)
 thr_nc = zscore_from_pval((fwer_target / 2) * correction_no_cluster)
 
+
 # %%
 # Inference with several algorithms
 # ---------------------------------
@@ -139,6 +131,9 @@
 # the theoretical tolerance region.
 
 
+import numpy as np
+
+
 # The following function builds a 2D map with four active regions that are
 # enfolded by thin tolerance regions.
 def weight_map_2D_extended(shape, roi_size, delta):
@@ -174,6 +169,7 @@ def weight_map_2D_extended(shape, roi_size, delta):
 # compute true support with visible spatial tolerance
 beta_extended = weight_map_2D_extended(shape, roi_size, delta)
 
+
 # %%
 # Now, we compute the support estimated by a high-dimensional statistical
 # inference method that does not leverage the data structure.
@@ -183,56 +179,53 @@ def weight_map_2D_extended(shape, roi_size, delta):
 # and referred to as Desparsified Lasso.
 
 
+from hidimstat import DesparsifiedLasso
+
 # compute desparsified lasso
-beta_hat, sigma_hat, precision_diagonal = desparsified_lasso(
-    X_init,
-    y,
-    n_jobs=n_jobs,
-    random_state=0,
-)
-pval, pval_corr, one_minus_pval, one_minus_pval_corr, cb_min, cb_max = (
-    desparsified_lasso_pvalue(
-        X_init.shape[0],
-        beta_hat,
-        sigma_hat,
-        precision_diagonal,
-    )
-)
+desparsified_lasso = DesparsifiedLasso(n_jobs=n_jobs, random_state=0)
+desparsified_lasso.fit_importance(X_init, y)
 
 # compute estimated support (first method)
-zscore = zscore_from_pval(pval, one_minus_pval)
+zscore = zscore_from_pval(
+    desparsified_lasso.pvalues_, desparsified_lasso.one_minus_pvalues_
+)
 selected_dl = zscore > thr_nc  # use the "no clustering threshold"
 
 # compute estimated support (second method)
 selected_dl = np.logical_or(
-    pval_corr < fwer_target / 2, one_minus_pval_corr < fwer_target / 2
+    desparsified_lasso.pvalues_corr_ < fwer_target / 2,
+    desparsified_lasso.one_minus_pvalues_corr_ < fwer_target / 2,
 )
 
+
 # %%
 # Now, we compute the support estimated using a clustered inference algorithm
 # (c.f. :footcite:t:`chevalier2022spatially`) called Clustered Desparsified Lasso
 # (CluDL) since it uses the Desparsified Lasso technique after clustering the data.
 
 # Define the FeatureAgglomeration object that performs the clustering.
 # This object is necessary to run the current algorithm and the following one.
+
+from sklearn.cluster import FeatureAgglomeration
+from sklearn.feature_extraction import image
+from sklearn.preprocessing import StandardScaler
+
+from hidimstat.ensemble_clustered_inference import (
+    clustered_inference,
+    clustered_inference_pvalue,
+)
+
 connectivity = image.grid_to_graph(n_x=shape[0], n_y=shape[1])
 ward = FeatureAgglomeration(
     n_clusters=n_clusters, connectivity=connectivity, linkage="ward"
 )
 
 # clustered desparsified lasso (CluDL)
-ward_, beta_hat, theta_hat, omega_diag = clustered_inference(
+ward_, desparsified_lasso_ = clustered_inference(
     X_init, y, ward, scaler_sampling=StandardScaler(), random_state=0
 )
 beta_hat, pval, pval_corr, one_minus_pval, one_minus_pval_corr = (
-    clustered_inference_pvalue(
-        n_samples,
-        False,
-        ward_,
-        beta_hat,
-        theta_hat,
-        omega_diag,
-    )
+    clustered_inference_pvalue(n_samples, False, ward_, desparsified_lasso_)
 )
 
 # compute estimated support (first method)
@@ -244,40 +237,46 @@ def weight_map_2D_extended(shape, roi_size, delta):
     pval_corr < fwer_target / 2, one_minus_pval_corr < fwer_target / 2
 )
 
+
 # %%
 # Finally, we compute the support estimated by an ensembled clustered
 # inference algorithm (c.f. :footcite:t:`chevalier2022spatially`). This algorithm is called
 # Ensemble of Clustered Desparsified Lasso (EnCluDL) since it runs several
 # CluDL algorithms with different clustering choices. The different CluDL
 # solutions are then aggregated into one.
 
+from hidimstat.ensemble_clustered_inference import (
+    ensemble_clustered_inference,
+    ensemble_clustered_inference_pvalue,
+)
+
 # ensemble of clustered desparsified lasso (EnCluDL)
-list_ward, list_beta_hat, list_theta_hat, list_omega_diag = (
-    ensemble_clustered_inference(
-        X_init,
-        y,
-        ward,
-        scaler_sampling=StandardScaler(),
-        random_state=0,
-    )
+list_ward, list_desparsified_lasso = ensemble_clustered_inference(
+    X_init,
+    y,
+    ward,
+    scaler_sampling=StandardScaler(),
+    random_state=0,
+    n_jobs=n_jobs,
 )
 beta_hat, selected_ecdl = ensemble_clustered_inference_pvalue(
     n_samples,
     False,
     list_ward,
-    list_beta_hat,
-    list_theta_hat,
-    list_omega_diag,
+    list_desparsified_lasso,
     fdr=fwer_target,
 )
 
+
 # %%
 # Results
 # -------
 #
 # Now we plot the true support, the theoretical tolerance regions and
 # the estimated supports for every method.
 
+import matplotlib.pyplot as plt
+
 
 # To generate a plot that exhibits
 # the true support and the estimated supports for every method,
@@ -342,6 +341,7 @@ def plot(maps, titles):
 
 plot(maps, titles)
 
+
 # %%
 # Analysis of the results
 # -----------------------
 
@@ -38,12 +38,15 @@
 from nilearn.image import mean_img
 from nilearn.maskers import NiftiMasker
 from nilearn.plotting import plot_stat_map, show
+from sklearn.base import clone
 from sklearn.cluster import FeatureAgglomeration
 from sklearn.feature_extraction import image
+from sklearn.linear_model import LassoCV
+from sklearn.model_selection import KFold
 from sklearn.preprocessing import StandardScaler
 from sklearn.utils import Bunch
 
-from hidimstat.desparsified_lasso import desparsified_lasso, desparsified_lasso_pvalue
+from hidimstat.desparsified_lasso import DesparsifiedLasso
 from hidimstat.ensemble_clustered_inference import (
     clustered_inference,
     clustered_inference_pvalue,
@@ -144,19 +147,34 @@ def preprocess_haxby(subject=2, memory=None):
 new_hard_limit = limit_5G if hard < 0 else min(limit_5G, hard)
 resource.setrlimit(resource.RLIMIT_AS, (new_soft_limit, new_hard_limit))
 
+# Default estimator
+estimator = LassoCV(
+    eps=1e-2,
+    fit_intercept=False,
+    cv=KFold(n_splits=5, shuffle=True, random_state=0),
+    tol=1e-2,
+    max_iter=6000,
+    random_state=1,
+    n_jobs=1,
+)
+
+
 # First, we try to recover the discriminative pattern by computing
 # p-values from desparsified lasso.
 # Due to the size of the X, it's not possible to use this method with a limit
 # of 5 G for memory. To handle this problem, the following methods use some
 # feature aggregation methods.
 #
 try:
-    beta_hat, sigma_hat, precision_diagonal = desparsified_lasso(
-        X, y, noise_method="median", max_iteration=1000, random_state=0, n_jobs=n_jobs
-    )
-    pval_dl, _, one_minus_pval_dl, _, cb_min, cb_max = desparsified_lasso_pvalue(
-        X.shape[0], beta_hat, sigma_hat, precision_diagonal
+    desparsified_lasso = DesparsifiedLasso(
+        noise_method="median",
+        estimator=clone(estimator),
+        random_state=0,
+        n_jobs=n_jobs,
     )
+    desparsified_lasso.fit_importance(X, y)
+    pval_dl = desparsified_lasso.pvalues_
+    one_minus_pval_dl = desparsified_lasso.one_minus_pvalues_
 except MemoryError as err:
     pval_dl = None
     one_minus_pval_dl = None
@@ -165,17 +183,18 @@ def preprocess_haxby(subject=2, memory=None):
 # %%
 # Now, the clustered inference algorithm which combines parcellation
 # and high-dimensional inference (c.f. References).
-ward_, beta_hat, theta_hat, omega_diag = clustered_inference(
+ward_, cl_desparsified_lasso = clustered_inference(
     X,
     y,
     ward,
     scaler_sampling=StandardScaler(),
-    tolerance=1e-2,
+    estimator=clone(estimator),
+    tolerance_reid=1e-2,
     random_state=1,
     n_jobs=n_jobs,
 )
 beta_hat, pval_cdl, _, one_minus_pval_cdl, _ = clustered_inference_pvalue(
-    X.shape[0], None, ward_, beta_hat, theta_hat, omega_diag
+    X.shape[0], None, ward_, cl_desparsified_lasso
 )
 
 # %%
@@ -185,28 +204,24 @@ def preprocess_haxby(subject=2, memory=None):
 # which means that 5 different parcellations are considered and
 # then 5 statistical maps are produced and aggregated into one.
 # However you might benefit from clustering randomization taking
-# `n_bootstraps=25` or `n_bootstraps=100`, also we set `n_jobs`.
-list_ward, list_beta_hat, list_theta_hat, list_omega_diag = (
-    ensemble_clustered_inference(
-        X,
-        y,
-        ward,
-        groups=groups,
-        scaler_sampling=StandardScaler(),
-        n_bootstraps=5,
-        max_iteration=6000,
-        tolerance=1e-2,
-        random_state=2,
-        n_jobs=n_jobs,
-    )
+# `n_bootstraps=25` or `n_bootstraps=100`, also we set `n_jobs=n_jobs`.
+list_ward, list_cl_desparsified_lasso = ensemble_clustered_inference(
+    X,
+    y,
+    ward,
+    groups=groups,
+    scaler_sampling=StandardScaler(),
+    n_bootstraps=5,
+    estimator=clone(estimator),
+    tolerance_reid=1e-2,
+    random_state=2,
+    n_jobs=n_jobs,
 )
 beta_hat, selected = ensemble_clustered_inference_pvalue(
     X.shape[0],
     False,
     list_ward,
-    list_beta_hat,
-    list_theta_hat,
-    list_omega_diag,
+    list_cl_desparsified_lasso,
     fdr=0.1,
 )