Add tools for testing and make an example with CPI (#265)

* remove not necessary function for testing

* add a check on the number of featires for X

* add assertion on the model

* improve the check_fit

* update data_generation

* add function for the generation of data for the tests

* improve test for CPI

* fix tests

* improve geenration of model

* change multivariate function but need to fix the function reid

* change the generation of data

* Add TODO

It requires some reflexion about the method.
The threshold based on the coefficient for estimating the weight seems
strange.
The equation signal to ratio seems also weird

* add tests for senario

* fix tests

* Update src/hidimstat/_utils/scenario.py

Co-authored-by: bthirion <[email protected]>

* Fix docstring in the PR

* fix noise_mag

* improve test of knockoff

* small improvement

* fix estimation of variance

* clean noise_std

* change name of the sigma

* add the possibility of continous support

* dfix test by increasing the snr and find a good seed

* fix test for noise_std

* fix tests

* fix tests

* fix errror

* fix dcrt example

* fix bug in exampel by modification of example

* fix size of tests

* Apply suggestions from code review

Co-authored-by: Joseph Paillard <[email protected]>

* fix time and seed for test

* Fix number of threadpool to 1

* increase coverage

* add a warning

* add test for warning

* Modify snr for better approch

* Apply suggestions from code review

Co-authored-by: bthirion <[email protected]>

* Update src/hidimstat/_utils/scenario.py

Co-authored-by: bthirion <[email protected]>

* modify senteces of the tests

* fix situation of 0 support

* change name of noise including spatial information

* change mane of the noise fonction

* Change parameter continous

* fix range of paraemter for rho_noise_time

* fix message of the error

* fix message of error

* fix test

* fix error in name parameters

* fix name in the error

* fix tests

* fix generation of data parameters

* fix the missmathc of feature

* fix error when the support of noise was zero

* fix bug in the assertion

* fix assertion

* Update src/hidimstat/_utils/scenario.py

Co-authored-by: bthirion <[email protected]>

* transform beta in boolean array

* rename noise serial

* Improve comment of shuffle

* done

* remove assertion on number of jobs

* Improve docstring

* fix bug in tests

* fix tests cpi

* Remove unessesary tests

* fix test noise_std

* replace n_times by n_targets

* fix change name

* Chage name for index

* refactor data generation with spatial

* fix tests for senario

* chnage the example

* update regression test

* update test

* update noise_std test

* fix test_senario

* chnage to snr

* Increase the robustness ofreid test

I modified the parameters for having a more robust tests

* update docstring

* Put back cov in the tests

* Modify warning message

* fix tests

* simplify the assert

* chnage name of the file

* fix docstring

* change sigma by signal_noise_ratio

* Update the example

* change error message

* Improve test of the noise

* fix format of pyproject

* Improve the different type of indexing

* fix a bug in the condition

* Update test/test_conditional_permutation_importance.py

Co-authored-by: Joseph Paillard <[email protected]>

* rename the main function

* update the condition for the number of features

* remove the count of index

* improve the test and the coverage

* update the modification

* Apply suggestions from code review

Co-authored-by: bthirion <[email protected]>

* finish merge

* fix example

* remove the categoriacl in function of cpi test

* fix rename

* rename some tests

* Update examples/plot_dcrt_example.py

fix error in the computation of type-1 error.

Co-authored-by: Joseph Paillard <[email protected]>

* fix format

---------

Co-authored-by: bthirion <[email protected]>
Co-authored-by: Joseph Paillard <[email protected]>

Author: 3 people

examples/plot_2D_simulation_example.py

@@ -69,7 +69,7 @@
     ensemble_clustered_inference_pvalue,
 )
 from hidimstat.statistical_tools.p_values import zscore_from_pval
-from hidimstat._utils.scenario import multivariate_simulation
+from hidimstat._utils.scenario import multivariate_simulation_spatial
 
 #############################################################################
 # Specific plotting functions
@@ -167,12 +167,12 @@ def plot(maps, titles):
 shape = (40, 40)
 n_features = shape[1] * shape[0]
 roi_size = 4  # size of the edge of the four predictive regions
-sigma = 2.0  # noise standard deviation
+signal_noise_ratio = 10.0  # noise standard deviation
 smooth_X = 1.0  # level of spatial smoothing introduced by the Gaussian filter
 
 # generating the data
-X_init, y, beta, epsilon, _, _ = multivariate_simulation(
-    n_samples, shape, roi_size, sigma, smooth_X, seed=1
+X_init, y, beta, epsilon = multivariate_simulation_spatial(
+    n_samples, shape, roi_size, signal_noise_ratio, smooth_X, seed=1
 )
 
 ##############################################################################

examples/plot_dcrt_example.py

@@ -20,7 +20,7 @@
 from sklearn.linear_model import LassoCV
 
 from hidimstat import D0CRT
-from hidimstat._utils.scenario import multivariate_1D_simulation
+from hidimstat._utils.scenario import multivariate_simulation
 
 #############################################################################
 # Processing the computations
@@ -38,14 +38,20 @@
     # Number of relevant variables
     n_signal = 2
     # Signal-to-noise ratio
-    snr = 4
+    signal_noise_ratio = 4
     # Correlation coefficient
     rho = 0.8
     # Nominal false positive rate
     alpha = 5e-2
 
-    X, y, _, __ = multivariate_1D_simulation(
-        n_samples=n, n_features=p, support_size=n_signal, rho=rho, seed=sim_ind
+    X, y, beta_true, noise = multivariate_simulation(
+        n_samples=n,
+        n_features=p,
+        support_size=n_signal,
+        rho=rho,
+        signal_noise_ratio=signal_noise_ratio,
+        shuffle=True,
+        seed=sim_ind,
     )
 
     # Applying a reLu function on the outcome y to get non-linear relationships
@@ -58,8 +64,9 @@
     results_list.append(
         {
             "model": "Lasso",
-            "type-1 error": sum(pvals_lasso[n_signal:] < alpha) / (p - n_signal),
-            "power": sum(pvals_lasso[:n_signal] < alpha) / (n_signal),
+            "type-1 error": sum(pvals_lasso[np.logical_not(beta_true)] < alpha)
+            / (p - n_signal),
+            "power": sum(pvals_lasso[beta_true] < alpha) / (n_signal),
         }
     )
 
@@ -73,8 +80,9 @@
     results_list.append(
         {
             "model": "RF",
-            "type-1 error": sum(pvals_forest[n_signal:] < alpha) / (p - n_signal),
-            "power": sum(pvals_forest[:n_signal] < alpha) / (n_signal),
+            "type-1 error": sum(pvals_forest[np.logical_not(beta_true)] < alpha)
+            / (n_signal),
+            "power": sum(pvals_forest[beta_true] < alpha) / (n_signal),
         }
     )
 

examples/plot_importance_classification_iris.py

@@ -116,7 +116,7 @@ def run_one_fold(X, y, model, train_index, test_index, vim_name="CFI", groups=No
     GridSearchCV(SVC(kernel="rbf"), {"C": np.logspace(-3, 3, 10)}),
 ]
 cv = KFold(n_splits=5, shuffle=True, random_state=0)
-groups = {ft: i for i, ft in enumerate(dataset.feature_names)}
+groups = {ft: [i] for i, ft in enumerate(dataset.feature_names)}
 out_list = Parallel(n_jobs=5)(
     delayed(run_one_fold)(
         X, y, model, train_index, test_index, vim_name=vim_name, groups=groups

examples/plot_knockoff_aggregation.py

@@ -31,7 +31,7 @@
     model_x_knockoff_pvalue,
 )
 from hidimstat.statistical_tools.multiple_testing import fdp_power
-from hidimstat._utils.scenario import multivariate_1D_simulation_AR
+from hidimstat._utils.scenario import multivariate_simulation
 
 
 #############################################################################
@@ -49,14 +49,14 @@
 # Number of variables
 n_features = 150
 # Correlation parameter
-rho = 0.4
+rho = 0.5
 # Ratio of number of variables with non-zero coefficients over total
 # coefficients
 sparsity = 0.2
 # Desired controlled False Discovery Rate (FDR) level
 fdr = 0.1
 # signal noise ration
-snr = 10
+signal_noise_ratio = 10
 # number of repetitions for the bootstraps
 n_bootstraps = 25
 # seed for the random generator
@@ -73,11 +73,26 @@
 #######################################################################
 # Define the function for running the three procedures on the same data
 # ---------------------------------------------------------------------
-def single_run(n_samples, n_features, rho, sparsity, snr, fdr, n_bootstraps, seed=None):
+def single_run(
+    n_samples,
+    n_features,
+    rho,
+    sparsity,
+    signal_noise_ratio,
+    fdr,
+    n_bootstraps,
+    seed=None,
+):
     # Generate data
-    X, y, _, non_zero_index = multivariate_1D_simulation_AR(
-        n_samples, n_features, rho=rho, sparsity=sparsity, seed=seed, snr=snr
+    X, y, beta_true, noise = multivariate_simulation(
+        n_samples,
+        n_features,
+        rho=rho,
+        support_size=int(n_features * sparsity),
+        signal_noise_ratio=signal_noise_ratio,
+        seed=seed,
     )
+    non_zero_index = np.where(beta_true)[0]
 
     # Use model-X Knockoffs [1]
     selected, test_scores, threshold, X_tildes = model_x_knockoff(
@@ -165,7 +180,14 @@ def effect_number_samples(n_samples):
     parallel = Parallel(n_jobs, verbose=joblib_verbose)
     results = parallel(
         delayed(single_run)(
-            n_samples, n_features, rho, sparsity, snr, fdr, n_bootstraps, seed=seed
+            n_samples,
+            n_features,
+            rho,
+            sparsity,
+            signal_noise_ratio,
+            fdr,
+            n_bootstraps,
+            seed=seed,
         )
         for seed in seed_list
     )

pyproject.toml

@@ -110,3 +110,9 @@ markers = ["slow: marks tests as slow (deselect with '-m \"not slow\"')"]
 # pytest-timeout
 timeout = 60           # on individual test should not take more than 10s
 session_timeout = 1200 # all the tests should be run in 5 min
+
+[tool.pytest_env]
+OPENBLAS_NUM_THREADS = 1
+BLIS_NUM_THREADS = 1
+MKL_NUM_THREADS = 1
+OMP_NUM_THREADS = 1

src/hidimstat/_utils/exception.py

@@ -0,0 +1,12 @@
+class InternalError(BaseException):
+    """
+    Create an error for internal error of the library
+
+    Parameters
+    ----------
+    message: str
+        Message of explanation of the error
+    """
+
+    def __init__(self, message):
+        self.message = message