diff --git a/experiments/SIAM/.gitignore b/experiments/SIAM/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..90b53353a4f477d61fa4a41a5f8dc03061b41d9d
--- /dev/null
+++ b/experiments/SIAM/.gitignore
@@ -0,0 +1,2 @@
+*.eps
+*.pdf
\ No newline at end of file
diff --git a/experiments/SIAM/setup.py b/experiments/SIAM/setup.py
index f853db30f9b0f5e7deb3d839e9d375a71d5fa7d2..4cfe71b2480e3fe5c8436843d79ad9e49fd2e3d2 100755
--- a/experiments/SIAM/setup.py
+++ b/experiments/SIAM/setup.py
@@ -23,6 +23,8 @@ class Setup(object):
 
       self.n_rep = data["n_rep"]
 
+      self.eval_gap_it = data["eval_gap_it"]
+
       self.list_sequence = data["sequences"]
       self.list_ratio_lbd = data["list_ratio_lbd"]
 
diff --git a/experiments/SIAM/setups.json b/experiments/SIAM/setups.json
index 6d86f2fae26bcea84b38d6a0a83b82f7b6866025..a89e488f8717425b7731537b3a6dad1eb256478b 100755
--- a/experiments/SIAM/setups.json
+++ b/experiments/SIAM/setups.json
@@ -1,9 +1,9 @@
 {
-	"setup1a": {
+	"setupSIAM": {
 		"remark": "OSCAR 1, 2, 3",
 
-		"m": 10,
-		"n": 30,
+		"m": 100,
+		"n": 300,
 
 		"dictionaries": [
 			"gaussian 0.",
@@ -14,6 +14,8 @@
 
 		"n_rep": 50,
 
+		"eval_gap_it": 20,
+
 		"sequences": [
 			["oscar-lim", 1, 0.9],
 			["oscar-lim", 1.0, 0.1],
@@ -42,6 +44,8 @@
 
 		"n_rep": 50,
 
+		"eval_gap_it": 20,
+
 		"sequences": [
    			["exp-lim", 0.99999, 0.9],
    			["exp-lim", 0.99999, 0.5],
diff --git a/experiments/SIAM/slopepb.py b/experiments/SIAM/slopepb.py
index f43760f9812824773a0dfbfea16f3514a003a2e0..c301bd6c0a56c29e5b47f2427cd4df2383787a25 100644
--- a/experiments/SIAM/slopepb.py
+++ b/experiments/SIAM/slopepb.py
@@ -21,6 +21,15 @@ class SlopePb(object):
       self.lbd = ratio_lbd * self.lbdmax
 
 
+   def get_dual_scaling_coeff(self, vecr):
+
+      beta_dual = np.sort(np.abs(self.matA.T @ vecr))[::-1]
+      beta_dual = np.cumsum(beta_dual) / \
+         np.cumsum(self.lbd * self.vec_gammas)
+
+      return 1. / max(np.max(beta_dual), 1.)
+
+
    def make_dual_scaling(self, vecr):
 
       beta_dual = np.sort(np.abs(self.matA.T @ vecr))[::-1]
diff --git a/experiments/SIAM/xp_0_balls/make_fig_paper.sh b/experiments/SIAM/xp_0_balls/make_fig_paper.sh
new file mode 100755
index 0000000000000000000000000000000000000000..02554698e3ce2952dc11b52195654208a26351f2
--- /dev/null
+++ b/experiments/SIAM/xp_0_balls/make_fig_paper.sh
@@ -0,0 +1,3 @@
+python xp_a_accuracy_sol.py --id SIAM
+python xp_b_screening.py --id SIAM
+python xp_c_viz_paper.py --id SIAM
\ No newline at end of file
diff --git a/experiments/SIAM/xp_0_balls/process_data.py b/experiments/SIAM/xp_0_balls/process_data.py
new file mode 100644
index 0000000000000000000000000000000000000000..d8a4f63bcd5db7e8f43f704e9275d0e1713c2dd9
--- /dev/null
+++ b/experiments/SIAM/xp_0_balls/process_data.py
@@ -0,0 +1,56 @@
+# -*- coding: utf-8 -*-
+import numpy as np
+
+from xps.SIAM.setup import Setup
+from xps.SIAM.xp_1_balls import xpparams
+
+
+def process(setup, log=True):
+   
+   # ---- load ----
+   solutions_filename = f"results/xp_setup{setup.setup_id}.npz"
+   solutions = np.load(solutions_filename, allow_pickle=True)
+
+   screenings_filename = f"results/xp_setup{setup.setup_id}_screening.npz"
+   screenings = np.load(screenings_filename, allow_pickle=True)
+
+      # mat_results_gap[1, i_dic, i_seq, rep, i_ratio]
+   mat_pvopt = solutions["mat_pvopt"]
+
+   mat_nb_zero_detected = screenings['mat_nb_zero_detected']
+   nb_test = mat_nb_zero_detected.shape[0]
+   list_tests = screenings["list_test"]
+
+   # ---- log ----
+   if log:
+      print("Experiment info")
+      # print(f"- run with version {solutions["version"]}")
+      print(f"- setup{setup.setup_id}")
+      print("")
+
+   # ---- processing ----
+   mat_pc_detected =np.zeros(
+      (nb_test, xpparams.nb_point, setup.nb_dic, setup.nb_sequence, setup.nb_ratio_lbd, setup.n_rep)
+   )
+
+   for i_dic in range(setup.nb_dic):
+      for i_seq in range(setup.nb_sequence):
+         for i_ratio in range(setup.nb_ratio_lbd):
+            for rep in range(setup.n_rep):
+
+               nb_0 = np.sum(mat_pvopt[i_dic, i_seq, i_ratio, rep, :] == 0)
+               if nb_0 > 0:
+                  mat_detected = mat_nb_zero_detected[:, :, i_dic, i_seq, i_ratio, rep]
+                  mat_pc_detected[:, :, i_dic, i_seq, i_ratio, rep] = mat_detected / float(nb_0)
+
+
+   # ---- return ----
+   return {
+      "mat_pc_detected": np.mean(mat_pc_detected, axis=5),
+      "list_tests": list_tests,
+   }
+
+   # {
+   #    "results_mat_nb": results_mat_nb,
+   #    "results_mat_nbz": results_mat_nbz,
+   # }
\ No newline at end of file
diff --git a/experiments/SIAM/xp_0_balls/xp_a_accuracy_sol.py b/experiments/SIAM/xp_0_balls/xp_a_accuracy_sol.py
new file mode 100755
index 0000000000000000000000000000000000000000..795f54c63ff5a06c82c943601181c990126022c4
--- /dev/null
+++ b/experiments/SIAM/xp_0_balls/xp_a_accuracy_sol.py
@@ -0,0 +1,153 @@
+# -*- coding: utf-8 -*-
+import time, argparse, sys
+
+import numpy as np
+
+# Algorithm import
+from src import __version__
+from src.solver.parameters import SlopeParameters, EnumLipchitzOptions
+from src.solver.slope import primal_func, dual_func, slope_gp
+
+# Generative models
+from src.utils import get_lambda_max, gamma_sequence_generator
+from src.dictionaries import generate_dic
+
+# Screening
+from src.screening.gap_test_all import GapTestAll
+
+# XP import
+from experiments.SIAM.slopepb import SlopePb
+from experiments.SIAM.setup import Setup
+from experiments.SIAM.xp_1_balls import xpparams
+
+parser=argparse.ArgumentParser()
+parser.add_argument('--erase', help='erase existing results', action="store_true")
+parser.add_argument('--id', help='setup id', type=str)
+args=parser.parse_args()
+
+
+# -------------------------
+#       Load  Setup
+# -------------------------
+
+setup = Setup(args.id)
+out_file_name = f"results/xp_setup{args.id}.npz"
+
+mat_seed = np.random.randint(
+   0, 2**32-1,
+   size=(setup.nb_dic, setup.nb_sequence, setup.n_rep), 
+)
+
+mat_pvopt = np.zeros(
+   (setup.nb_dic, setup.nb_sequence, setup.nb_ratio_lbd, setup.n_rep, setup.n)
+)
+
+mat_dvopt = np.zeros(
+   (setup.nb_dic, setup.nb_sequence, setup.nb_ratio_lbd, setup.n_rep, setup.m)
+)
+
+
+try:
+   if args.erase:
+      raise FileNotFoundError
+
+   # Try to load existing results
+   load_results = np.load(out_file_name, allow_pickle=True)
+   mat_seed  = load_results["mat_seed"]
+   mat_pvopt = load_results["mat_pvopt"]
+   mat_dvopt = load_results["mat_dvopt"]
+
+except FileNotFoundError:
+   pass
+
+
+# -------------------------
+#            Xp
+# -------------------------
+
+# For each trial
+# 1. Compute high accuracy solution
+# 2. Create ideal safe ball
+# 3. increase radius
+
+nb_xp = setup.nb_dic * setup.n_rep * setup.nb_sequence * setup.nb_ratio_lbd
+i_xp = 0
+for i_dic in range(setup.nb_dic):
+   for i_seq, seq in enumerate(setup.list_sequence):
+      for rep in range(setup.n_rep):
+
+         np.random.seed(mat_seed[i_dic, i_seq, rep])
+
+         # ---- 1. Gen data ----
+         matA = generate_dic(
+            setup.list_dic[i_dic],
+            setup.m,
+            setup.n,
+            setup.normalize
+         )
+
+         lip = np.linalg.norm(matA, ord=2)**2
+         vecy = np.random.randn(setup.m)
+         vecy /= np.linalg.norm(vecy, 2)
+
+
+         # ---- 2. Compute parameters ----
+         vec_gammas = gamma_sequence_generator(
+            setup.m, 
+            setup.n,
+            setup.list_sequence[i_seq][0],
+            setup.list_sequence[i_seq][1:]
+         )
+
+         lbd_max = get_lambda_max(vecy, matA, vec_gammas)
+
+
+         # ---- 3. XP ----
+         for i_ratio, ratio in enumerate(setup.list_ratio_lbd):
+            i_xp += 1
+
+            slopePb = SlopePb(matA, vecy, vec_gammas, ratio, lbdmax=lbd_max)
+            vecx_hat = mat_pvopt[i_dic, i_seq, i_ratio, rep, :]
+            vecu_hat = mat_dvopt[i_dic, i_seq, i_ratio, rep, :]
+
+            gap = slopePb.eval_gap(vecx_hat, vecu_hat)
+            print(f"xp {i_xp} / {nb_xp} --- (gap={gap})")
+
+            if (gap > xpparams.stopping_gap):
+               gap_old = gap
+
+               # ---- 3a. Find solution ----
+               params = SlopeParameters()
+               params.vecx_init = np.copy(vecx_hat)
+               params.lipchitz_constant = lip
+               params.lipchitz_update = EnumLipchitzOptions.EXACT
+               params.max_it = 1e7
+               params.gap_stopping = xpparams.stopping_gap
+               params.time_stopping = np.inf
+               params.screening1 = GapTestAll(vec_gammas)
+               params.eval_gap   = True
+               params.eval_gap_it = setup.eval_gap_it
+               params.accelerated = False
+               params.verbose = False
+               
+               out_slope = slope_gp(vecy, matA, ratio * lbd_max, vec_gammas, params)
+
+               dv  = vecy - matA @  out_slope["sol"]
+               dv = slopePb.make_dual_scaling(dv)
+               gap = slopePb.eval_gap(out_slope["sol"], dv)
+
+               if gap <= gap_old:
+                  mat_pvopt[i_dic, i_seq, i_ratio, rep, :] = out_slope["sol"]
+                  mat_dvopt[i_dic, i_seq, i_ratio, rep, :] = dv
+
+                  # Save
+                  np.savez(out_file_name,
+                     mat_seed=mat_seed,
+                     mat_pvopt=mat_pvopt,
+                     mat_dvopt=mat_dvopt,
+                     version = __version__,
+                     allow_pickle=True
+                  )
+
+               # if i_xp ==3:
+               #    exit()
\ No newline at end of file
diff --git a/experiments/SIAM/xp_0_balls/xp_b_screening.py b/experiments/SIAM/xp_0_balls/xp_b_screening.py
new file mode 100644
index 0000000000000000000000000000000000000000..d9b65cda52ceb40369447f0aeec104e71512f500
--- /dev/null
+++ b/experiments/SIAM/xp_0_balls/xp_b_screening.py
@@ -0,0 +1,122 @@
+# -*- coding: utf-8 -*-
+import argparse
+
+import numpy as np
+
+# Algorithm import
+from src import __version__
+from src.dictionaries import generate_dic
+from src.utils import get_lambda_max, gamma_sequence_generator
+
+# Screening
+from src.screening.gap_test_p_1 import GapTestPequalOne
+from src.screening.gap_test_p_q import GapTestPequalQ
+from src.screening.gap_test_all import GapTestAll
+
+# XP import
+from experiments.SIAM.slopepb import SlopePb
+from experiments.SIAM.setup import Setup
+from experiments.SIAM.xp_1_balls import xpparams
+
+
+parser=argparse.ArgumentParser()
+parser.add_argument('--id', help='setup id', type=str, default=1)
+args=parser.parse_args()
+
+# -------------------------
+#       Load  Setup
+# -------------------------
+
+setup = Setup(args.id)
+solutions_filename = f"results/xp_setup{args.id}.npz"
+
+try:
+   # Try to load existing results
+   load_results = np.load(solutions_filename, allow_pickle=True)
+   mat_seed  = load_results["mat_seed"]
+   mat_pvopt = load_results["mat_pvopt"]
+   mat_dvopt = load_results["mat_dvopt"]
+
+except FileNotFoundError:
+   print("No result found")
+   sys.exit(1)
+
+
+# -------------------------
+#        Screening Xp
+# -------------------------
+
+# percentage of detected zero
+mat_nb_zero_detected = np.full(
+    (3, xpparams.nb_point, setup.nb_dic, setup.nb_sequence, setup.nb_ratio_lbd, setup.n_rep),
+    np.nan,
+    dtype=int
+)
+screening_filename = f"results/xp_setup{args.id}_screening.npz"
+
+i_xp = 0
+nb_xp = setup.nb_dic * setup.n_rep * setup.nb_sequence * setup.nb_ratio_lbd
+for i_dic in range(setup.nb_dic):
+   for i_seq, seq in enumerate(setup.list_sequence):
+      for rep in range(setup.n_rep):
+
+         np.random.seed(mat_seed[i_dic, i_seq, rep])
+
+         # ---- 1. Gen data ----
+         matA = generate_dic(
+            setup.list_dic[i_dic],
+            setup.m,
+            setup.n,
+            setup.normalize
+         )
+
+         vecy = np.random.randn(setup.m)
+         vecy /= np.linalg.norm(vecy)
+
+
+         # ---- 2. Compute parameters ----
+         vec_gammas = gamma_sequence_generator(
+            setup.m, 
+            setup.n,
+            setup.list_sequence[i_seq][0],
+            setup.list_sequence[i_seq][1:]
+         )
+
+         lbd_max = get_lambda_max(vecy, matA, vec_gammas)
+
+         # ---- 3. XP ----
+         for i_ratio, ratio in enumerate(setup.list_ratio_lbd):
+            i_xp += 1
+            print(f"xp {i_xp} / {nb_xp}")
+
+            slopePb = SlopePb(matA, vecy, vec_gammas, ratio, lbdmax=lbd_max)
+            vecx_hat = mat_pvopt[i_dic, i_seq, i_ratio, rep, :]
+            vecu_hat = mat_dvopt[i_dic, i_seq, i_ratio, rep, :]
+            gap = slopePb.eval_gap(vecx_hat, vecu_hat)
+            Atu = matA.T @ vecu_hat
+
+            # ---- 3b. Build thin safe ball ----
+            rgap = np.sqrt(2 * gap)
+
+            # ---- 3c. Testing sphere ---- 
+            list_tests = [
+               GapTestPequalQ(),
+               GapTestAll(vec_gammas),
+               GapTestPequalOne(vec_gammas, np.arange(setup.n, dtype=np.double)),
+            ]
+
+
+            for i_offset, offset in enumerate(xpparams.vec_offsets):
+               for i_test, test in enumerate(list_tests):
+                  out = test.apply_test(np.abs(Atu), gap, ratio * lbd_max, vec_gammas, offset_radius=offset)
+                  mat_nb_zero_detected[i_test, i_offset, i_dic, i_seq, i_ratio, rep] = np.sum(out)
+
+
+         # Save
+         np.savez(
+            screening_filename,
+            mat_nb_zero_detected=mat_nb_zero_detected,
+            list_test = [test.get_name() for test in list_tests],
+            version = __version__,
+            allow_pickle=True
+         )
\ No newline at end of file
diff --git a/experiments/SIAM/xp_0_balls/xp_c_viz_paper.py b/experiments/SIAM/xp_0_balls/xp_c_viz_paper.py
new file mode 100644
index 0000000000000000000000000000000000000000..8ecde14a22a13d95982d5beb4790ef894fe2cd31
--- /dev/null
+++ b/experiments/SIAM/xp_0_balls/xp_c_viz_paper.py
@@ -0,0 +1,139 @@
+# -*- coding: utf-8 -*-
+from decimal import Decimal
+import json, argparse
+
+import numpy as np
+
+# XP import
+from xps.SIAM.slopepb import SlopePb
+from xps.SIAM.setup import Setup
+from xps.SIAM.xp_1_balls import xpparams
+from xps.SIAM.xp_1_balls.process_data import process
+
+
+parser=argparse.ArgumentParser()
+parser.add_argument('--noshow', help='do not display figures', action="store_true")
+parser.add_argument('--save', help='save figure', action="store_true")
+parser.add_argument('--id', help='setup id', type=str, default="SIAM")
+parser.add_argument('--i_seq', type=int, default=0)
+args=parser.parse_args()
+
+import matplotlib
+if args.noshow:
+   matplotlib.use('PS')
+else:
+   matplotlib.use("TkAgg")
+
+import matplotlib.pyplot as plt
+from matplotlib.legend_handler import HandlerBase
+import matplotlib.font_manager as font_manager
+
+# -------------------------
+#        Font stuff 
+# -------------------------
+
+fs = 20
+
+font_math = font_manager.FontProperties(
+   fname='../fonts/cmunrm.ttf',
+   # weight='normal',
+   # style='normal',
+   math_fontfamily="cm",
+   size=fs+2
+)
+
+font_text = font_manager.FontProperties(
+   fname='../fonts/cmunrm.ttf',
+   # weight='normal',
+   # style='normal',
+   math_fontfamily="cm",
+   size=fs+2
+)
+
+font_ttt = font_manager.FontProperties(
+   # fname='../fonts/ectt1000.ttf',
+   fname='../fonts/computer-modern/cmuntt.ttf',
+   weight='bold',
+   style='normal',
+   size=fs
+)
+
+
+# -------------------------
+#        Load Results
+# -------------------------
+
+
+setup_oscar = Setup(args.id)
+dic_process_oscar = process(setup_oscar)
+
+mat_pc_detected_oscar = dic_process_oscar["mat_pc_detected"]
+list_tests_oscar      = dic_process_oscar["list_tests"]
+
+
+
+# -------------------------
+#        Plot Results
+# -------------------------
+
+i_dic = 2
+i_lbd = 1
+
+fs=22
+fs_ylabels = 20
+list_colors  = ["tab:blue", "tab:orange", "tab:green"]
+list_legends = ["$p_q=q\\;\\forall q$", "all", "$p_q=1\\;\\forall q$"]
+# "best $r_q \\;\\forall q$ "
+
+print("printing xp_0_ball parameters with")
+print(f"- {setup_oscar.list_dic[i_dic]} dictionary")
+print(f"- lbd / lbd_max = {setup_oscar.list_ratio_lbd[i_lbd]}")
+
+
+f, ax = plt.subplots(1, 1, figsize=(.7*16, .6*9), sharex=True, sharey=True)
+
+# ax.set_title(f"OSCAR-{i_seq+1}", fontsize=fs+2)
+ax.set_xlabel(
+   r"$R$", 
+   fontsize=fs+2,
+   fontproperties=font_math,
+)
+ax.set_ylabel(
+   "% of zero entries detected",
+   fontsize=fs+2,
+   fontproperties=font_text
+)
+
+for tick in ax.xaxis.get_major_ticks():
+   tick.label.set_fontproperties(font_math)
+   tick.label.set_fontsize(20)
+
+for tick in ax.yaxis.get_major_ticks():
+   tick.label.set_fontproperties(font_math)
+   tick.label.set_fontsize(20)
+
+for i_test in [2, 0, 1]:
+   ax.plot(
+      xpparams.vec_offsets, 
+      100 * mat_pc_detected_oscar[i_test, :, i_dic, args.i_seq, i_lbd],
+      label = list_legends[i_test],
+      linewidth=4.,
+      alpha=.9,
+      color=list_colors[i_test]
+   )
+
+ax.legend(
+   fontsize=fs-2,
+   prop=font_math
+)
+
+ax.set_xscale("log")
+ax.set_xlim([1e-6, 1e0])
+ax.set_ylim([-2, 102])
+
+if args.save:
+   filename = f"figs/xp_illustration_screening{args.i_seq}.eps"
+   plt.savefig(filename, bbox_inches='tight')
+
+if not args.noshow:
+   plt.show()
\ No newline at end of file
diff --git a/src/screening/gap_test_all.c b/src/screening/gap_test_all.c
index 00a93b305924ba2297d346f92340f40b95d1d162..df0336502b4844e06da89b0964c46eba042dae9c 100644
--- a/src/screening/gap_test_all.c
+++ b/src/screening/gap_test_all.c
@@ -1181,6 +1181,7 @@ typedef npy_cdouble __pyx_t_5numpy_complex_t;
 struct __pyx_obj_3src_9screening_12gap_test_all_GapTestAll {
   PyObject_HEAD
   double *vec_cumsum_gammas;
+  char *name;
 };
 
 
@@ -1298,15 +1299,6 @@ static Py_ssize_t __Pyx_minusones[] = { -1, -1, -1, -1, -1, -1, -1, -1 };
 static Py_ssize_t __Pyx_zeros[] = { 0, 0, 0, 0, 0, 0, 0, 0 };
 
 #define __Pyx_BufPtrStrided1d(type, buf, i0, s0) (type)((char*)buf + i0 * s0)
-/* PyObjectSetAttrStr.proto */
-#if CYTHON_USE_TYPE_SLOTS
-#define __Pyx_PyObject_DelAttrStr(o,n) __Pyx_PyObject_SetAttrStr(o, n, NULL)
-static CYTHON_INLINE int __Pyx_PyObject_SetAttrStr(PyObject* obj, PyObject* attr_name, PyObject* value);
-#else
-#define __Pyx_PyObject_DelAttrStr(o,n)   PyObject_DelAttr(o,n)
-#define __Pyx_PyObject_SetAttrStr(o,n,v) PyObject_SetAttr(o,n,v)
-#endif
-
 /* PyThreadStateGet.proto */
 #if CYTHON_FAST_THREAD_STATE
 #define __Pyx_PyThreadState_declare  PyThreadState *__pyx_tstate;
@@ -1801,7 +1793,7 @@ static const char __pyx_k_np[] = "np";
 static const char __pyx_k_gap[] = "gap";
 static const char __pyx_k_lbd[] = "lbd";
 static const char __pyx_k_main[] = "__main__";
-static const char __pyx_k_name[] = "name";
+static const char __pyx_k_name[] = "__name__";
 static const char __pyx_k_sqrt[] = "sqrt";
 static const char __pyx_k_test[] = "__test__";
 static const char __pyx_k_dtype[] = "dtype";
@@ -1813,7 +1805,6 @@ static const char __pyx_k_Atcabs[] = "Atcabs";
 static const char __pyx_k_arange[] = "arange";
 static const char __pyx_k_cumsum[] = "cumsum";
 static const char __pyx_k_import[] = "__import__";
-static const char __pyx_k_name_2[] = "__name__";
 static const char __pyx_k_reduce[] = "__reduce__";
 static const char __pyx_k_argsort[] = "argsort";
 static const char __pyx_k_getstate[] = "__getstate__";
@@ -1828,7 +1819,6 @@ static const char __pyx_k_reduce_cython[] = "__reduce_cython__";
 static const char __pyx_k_setstate_cython[] = "__setstate_cython__";
 static const char __pyx_k_cline_in_traceback[] = "cline_in_traceback";
 static const char __pyx_k_coeff_dual_scaling[] = "coeff_dual_scaling";
-static const char __pyx_k_GAP_test_all_cython[] = "GAP test-all (cython)";
 static const char __pyx_k_AbstractGapScreening[] = "AbstractGapScreening";
 static const char __pyx_k_src_screening_ascreening[] = "src.screening.ascreening";
 static const char __pyx_k_numpy_core_multiarray_failed_to[] = "numpy.core.multiarray failed to import";
@@ -1836,7 +1826,6 @@ static const char __pyx_k_numpy_core_umath_failed_to_impor[] = "numpy.core.umath
 static const char __pyx_k_self_vec_cumsum_gammas_cannot_be[] = "self.vec_cumsum_gammas cannot be converted to a Python object for pickling";
 static PyObject *__pyx_n_s_AbstractGapScreening;
 static PyObject *__pyx_n_s_Atcabs;
-static PyObject *__pyx_kp_s_GAP_test_all_cython;
 static PyObject *__pyx_n_s_GapTestAll;
 static PyObject *__pyx_n_s_ImportError;
 static PyObject *__pyx_n_s_TypeError;
@@ -1853,7 +1842,6 @@ static PyObject *__pyx_n_s_index;
 static PyObject *__pyx_n_s_lbd;
 static PyObject *__pyx_n_s_main;
 static PyObject *__pyx_n_s_name;
-static PyObject *__pyx_n_s_name_2;
 static PyObject *__pyx_n_s_np;
 static PyObject *__pyx_n_s_numpy;
 static PyObject *__pyx_kp_s_numpy_core_multiarray_failed_to;
@@ -1887,7 +1875,7 @@ static PyObject *__pyx_tuple__4;
 static PyObject *__pyx_tuple__5;
 /* Late includes */
 
-/* "src/screening/gap_test_all.pyx":29
+/* "src/screening/gap_test_all.pyx":30
  *    @cython.boundscheck(False) # turn off bounds-checking for entire function
  *    @cython.wraparound(False)  # turn off negative index wrapping for entire function
  *    def __init__(             # <<<<<<<<<<<<<<
@@ -1924,7 +1912,7 @@ static int __pyx_pw_3src_9screening_12gap_test_all_10GapTestAll_1__init__(PyObje
         else goto __pyx_L5_argtuple_error;
       }
       if (unlikely(kw_args > 0)) {
-        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "__init__") < 0)) __PYX_ERR(0, 29, __pyx_L3_error)
+        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "__init__") < 0)) __PYX_ERR(0, 30, __pyx_L3_error)
       }
     } else if (PyTuple_GET_SIZE(__pyx_args) != 1) {
       goto __pyx_L5_argtuple_error;
@@ -1935,13 +1923,13 @@ static int __pyx_pw_3src_9screening_12gap_test_all_10GapTestAll_1__init__(PyObje
   }
   goto __pyx_L4_argument_unpacking_done;
   __pyx_L5_argtuple_error:;
-  __Pyx_RaiseArgtupleInvalid("__init__", 1, 1, 1, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 29, __pyx_L3_error)
+  __Pyx_RaiseArgtupleInvalid("__init__", 1, 1, 1, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 30, __pyx_L3_error)
   __pyx_L3_error:;
   __Pyx_AddTraceback("src.screening.gap_test_all.GapTestAll.__init__", __pyx_clineno, __pyx_lineno, __pyx_filename);
   __Pyx_RefNannyFinishContext();
   return -1;
   __pyx_L4_argument_unpacking_done:;
-  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_vec_gammas), __pyx_ptype_5numpy_ndarray, 1, "vec_gammas", 0))) __PYX_ERR(0, 30, __pyx_L1_error)
+  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_vec_gammas), __pyx_ptype_5numpy_ndarray, 1, "vec_gammas", 0))) __PYX_ERR(0, 31, __pyx_L1_error)
   __pyx_r = __pyx_pf_3src_9screening_12gap_test_all_10GapTestAll___init__(((struct __pyx_obj_3src_9screening_12gap_test_all_GapTestAll *)__pyx_v_self), __pyx_v_vec_gammas);
 
   /* function exit code */
@@ -1974,11 +1962,11 @@ static int __pyx_pf_3src_9screening_12gap_test_all_10GapTestAll___init__(struct
   __pyx_pybuffernd_vec_gammas.rcbuffer = &__pyx_pybuffer_vec_gammas;
   {
     __Pyx_BufFmt_StackElem __pyx_stack[1];
-    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_vec_gammas.rcbuffer->pybuffer, (PyObject*)__pyx_v_vec_gammas, &__Pyx_TypeInfo_nn_npy_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 29, __pyx_L1_error)
+    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_vec_gammas.rcbuffer->pybuffer, (PyObject*)__pyx_v_vec_gammas, &__Pyx_TypeInfo_nn_npy_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 30, __pyx_L1_error)
   }
   __pyx_pybuffernd_vec_gammas.diminfo[0].strides = __pyx_pybuffernd_vec_gammas.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_vec_gammas.diminfo[0].shape = __pyx_pybuffernd_vec_gammas.rcbuffer->pybuffer.shape[0];
 
-  /* "src/screening/gap_test_all.pyx":33
+  /* "src/screening/gap_test_all.pyx":34
  *    ):
  * 
  *       cdef int n = vec_gammas.shape[0]             # <<<<<<<<<<<<<<
@@ -1987,7 +1975,7 @@ static int __pyx_pf_3src_9screening_12gap_test_all_10GapTestAll___init__(struct
  */
   __pyx_v_n = (__pyx_v_vec_gammas->dimensions[0]);
 
-  /* "src/screening/gap_test_all.pyx":36
+  /* "src/screening/gap_test_all.pyx":37
  *       cdef int i
  * 
  *       self.vec_cumsum_gammas = <double*> PyMem_Malloc((n+1) * sizeof(double))             # <<<<<<<<<<<<<<
@@ -1996,7 +1984,7 @@ static int __pyx_pf_3src_9screening_12gap_test_all_10GapTestAll___init__(struct
  */
   __pyx_v_self->vec_cumsum_gammas = ((double *)PyMem_Malloc(((__pyx_v_n + 1) * (sizeof(double)))));
 
-  /* "src/screening/gap_test_all.pyx":37
+  /* "src/screening/gap_test_all.pyx":38
  * 
  *       self.vec_cumsum_gammas = <double*> PyMem_Malloc((n+1) * sizeof(double))
  *       self.vec_cumsum_gammas[0] = 0             # <<<<<<<<<<<<<<
@@ -2005,7 +1993,7 @@ static int __pyx_pf_3src_9screening_12gap_test_all_10GapTestAll___init__(struct
  */
   (__pyx_v_self->vec_cumsum_gammas[0]) = 0.0;
 
-  /* "src/screening/gap_test_all.pyx":39
+  /* "src/screening/gap_test_all.pyx":40
  *       self.vec_cumsum_gammas[0] = 0
  * 
  *       for i in range(n):             # <<<<<<<<<<<<<<
@@ -2017,7 +2005,7 @@ static int __pyx_pf_3src_9screening_12gap_test_all_10GapTestAll___init__(struct
   for (__pyx_t_3 = 0; __pyx_t_3 < __pyx_t_2; __pyx_t_3+=1) {
     __pyx_v_i = __pyx_t_3;
 
-    /* "src/screening/gap_test_all.pyx":40
+    /* "src/screening/gap_test_all.pyx":41
  * 
  *       for i in range(n):
  *          self.vec_cumsum_gammas[i+1] = self.vec_cumsum_gammas[i] + vec_gammas[i]             # <<<<<<<<<<<<<<
@@ -2028,16 +2016,16 @@ static int __pyx_pf_3src_9screening_12gap_test_all_10GapTestAll___init__(struct
     (__pyx_v_self->vec_cumsum_gammas[(__pyx_v_i + 1)]) = ((__pyx_v_self->vec_cumsum_gammas[__pyx_v_i]) + (*__Pyx_BufPtrStrided1d(npy_double *, __pyx_pybuffernd_vec_gammas.rcbuffer->pybuffer.buf, __pyx_t_4, __pyx_pybuffernd_vec_gammas.diminfo[0].strides)));
   }
 
-  /* "src/screening/gap_test_all.pyx":48
+  /* "src/screening/gap_test_all.pyx":49
  *       # self.vec_cumsum_gammas = &buf[0]
  * 
  *       self.name = "GAP test-all (cython)"             # <<<<<<<<<<<<<<
  * 
  * 
  */
-  if (__Pyx_PyObject_SetAttrStr(((PyObject *)__pyx_v_self), __pyx_n_s_name, __pyx_kp_s_GAP_test_all_cython) < 0) __PYX_ERR(0, 48, __pyx_L1_error)
+  __pyx_v_self->name = ((char *)"GAP test-all (cython)");
 
-  /* "src/screening/gap_test_all.pyx":29
+  /* "src/screening/gap_test_all.pyx":30
  *    @cython.boundscheck(False) # turn off bounds-checking for entire function
  *    @cython.wraparound(False)  # turn off negative index wrapping for entire function
  *    def __init__(             # <<<<<<<<<<<<<<
@@ -2065,7 +2053,7 @@ static int __pyx_pf_3src_9screening_12gap_test_all_10GapTestAll___init__(struct
   return __pyx_r;
 }
 
-/* "src/screening/gap_test_all.pyx":51
+/* "src/screening/gap_test_all.pyx":52
  * 
  * 
  *    def __dealloc__(self):             # <<<<<<<<<<<<<<
@@ -2088,7 +2076,7 @@ static void __pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_2__dealloc__(st
   __Pyx_RefNannyDeclarations
   __Pyx_RefNannySetupContext("__dealloc__", 0);
 
-  /* "src/screening/gap_test_all.pyx":52
+  /* "src/screening/gap_test_all.pyx":53
  * 
  *    def __dealloc__(self):
  *       PyMem_Free(self.vec_cumsum_gammas)  # no-op if self.data is NULL             # <<<<<<<<<<<<<<
@@ -2097,7 +2085,7 @@ static void __pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_2__dealloc__(st
  */
   PyMem_Free(__pyx_v_self->vec_cumsum_gammas);
 
-  /* "src/screening/gap_test_all.pyx":51
+  /* "src/screening/gap_test_all.pyx":52
  * 
  * 
  *    def __dealloc__(self):             # <<<<<<<<<<<<<<
@@ -2109,7 +2097,7 @@ static void __pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_2__dealloc__(st
   __Pyx_RefNannyFinishContext();
 }
 
-/* "src/screening/gap_test_all.pyx":57
+/* "src/screening/gap_test_all.pyx":58
  *    @cython.boundscheck(False) # turn off bounds-checking for entire function
  *    # @cython.wraparound(False)  # turn off negative index wrapping for entire function
  *    def apply_test(self,             # <<<<<<<<<<<<<<
@@ -2138,7 +2126,7 @@ static PyObject *__pyx_pw_3src_9screening_12gap_test_all_10GapTestAll_5apply_tes
     static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_Atcabs,&__pyx_n_s_gap,&__pyx_n_s_lbd,&__pyx_n_s_vec_gammas,&__pyx_n_s_coeff_dual_scaling,&__pyx_n_s_offset_radius,&__pyx_n_s_index,0};
     PyObject* values[7] = {0,0,0,0,0,0,0};
 
-    /* "src/screening/gap_test_all.pyx":64
+    /* "src/screening/gap_test_all.pyx":65
  *       double coeff_dual_scaling =1.,
  *       double offset_radius=0.,
  *       np.ndarray[long, ndim=1] index=None             # <<<<<<<<<<<<<<
@@ -2176,19 +2164,19 @@ static PyObject *__pyx_pw_3src_9screening_12gap_test_all_10GapTestAll_5apply_tes
         case  1:
         if (likely((values[1] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_gap)) != 0)) kw_args--;
         else {
-          __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, 1); __PYX_ERR(0, 57, __pyx_L3_error)
+          __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, 1); __PYX_ERR(0, 58, __pyx_L3_error)
         }
         CYTHON_FALLTHROUGH;
         case  2:
         if (likely((values[2] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_lbd)) != 0)) kw_args--;
         else {
-          __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, 2); __PYX_ERR(0, 57, __pyx_L3_error)
+          __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, 2); __PYX_ERR(0, 58, __pyx_L3_error)
         }
         CYTHON_FALLTHROUGH;
         case  3:
         if (likely((values[3] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_vec_gammas)) != 0)) kw_args--;
         else {
-          __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, 3); __PYX_ERR(0, 57, __pyx_L3_error)
+          __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, 3); __PYX_ERR(0, 58, __pyx_L3_error)
         }
         CYTHON_FALLTHROUGH;
         case  4:
@@ -2210,7 +2198,7 @@ static PyObject *__pyx_pw_3src_9screening_12gap_test_all_10GapTestAll_5apply_tes
         }
       }
       if (unlikely(kw_args > 0)) {
-        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "apply_test") < 0)) __PYX_ERR(0, 57, __pyx_L3_error)
+        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "apply_test") < 0)) __PYX_ERR(0, 58, __pyx_L3_error)
       }
     } else {
       switch (PyTuple_GET_SIZE(__pyx_args)) {
@@ -2229,16 +2217,16 @@ static PyObject *__pyx_pw_3src_9screening_12gap_test_all_10GapTestAll_5apply_tes
       }
     }
     __pyx_v_Atcabs = ((PyArrayObject *)values[0]);
-    __pyx_v_gap = __pyx_PyFloat_AsDouble(values[1]); if (unlikely((__pyx_v_gap == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 59, __pyx_L3_error)
-    __pyx_v_lbd = __pyx_PyFloat_AsDouble(values[2]); if (unlikely((__pyx_v_lbd == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 60, __pyx_L3_error)
+    __pyx_v_gap = __pyx_PyFloat_AsDouble(values[1]); if (unlikely((__pyx_v_gap == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 60, __pyx_L3_error)
+    __pyx_v_lbd = __pyx_PyFloat_AsDouble(values[2]); if (unlikely((__pyx_v_lbd == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 61, __pyx_L3_error)
     __pyx_v_vec_gammas = ((PyArrayObject *)values[3]);
     if (values[4]) {
-      __pyx_v_coeff_dual_scaling = __pyx_PyFloat_AsDouble(values[4]); if (unlikely((__pyx_v_coeff_dual_scaling == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 62, __pyx_L3_error)
+      __pyx_v_coeff_dual_scaling = __pyx_PyFloat_AsDouble(values[4]); if (unlikely((__pyx_v_coeff_dual_scaling == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 63, __pyx_L3_error)
     } else {
       __pyx_v_coeff_dual_scaling = ((double)1.);
     }
     if (values[5]) {
-      __pyx_v_offset_radius = __pyx_PyFloat_AsDouble(values[5]); if (unlikely((__pyx_v_offset_radius == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 63, __pyx_L3_error)
+      __pyx_v_offset_radius = __pyx_PyFloat_AsDouble(values[5]); if (unlikely((__pyx_v_offset_radius == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 64, __pyx_L3_error)
     } else {
       __pyx_v_offset_radius = ((double)0.);
     }
@@ -2246,18 +2234,18 @@ static PyObject *__pyx_pw_3src_9screening_12gap_test_all_10GapTestAll_5apply_tes
   }
   goto __pyx_L4_argument_unpacking_done;
   __pyx_L5_argtuple_error:;
-  __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 57, __pyx_L3_error)
+  __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 58, __pyx_L3_error)
   __pyx_L3_error:;
   __Pyx_AddTraceback("src.screening.gap_test_all.GapTestAll.apply_test", __pyx_clineno, __pyx_lineno, __pyx_filename);
   __Pyx_RefNannyFinishContext();
   return NULL;
   __pyx_L4_argument_unpacking_done:;
-  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_Atcabs), __pyx_ptype_5numpy_ndarray, 1, "Atcabs", 0))) __PYX_ERR(0, 58, __pyx_L1_error)
-  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_vec_gammas), __pyx_ptype_5numpy_ndarray, 1, "vec_gammas", 0))) __PYX_ERR(0, 61, __pyx_L1_error)
-  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_index), __pyx_ptype_5numpy_ndarray, 1, "index", 0))) __PYX_ERR(0, 64, __pyx_L1_error)
+  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_Atcabs), __pyx_ptype_5numpy_ndarray, 1, "Atcabs", 0))) __PYX_ERR(0, 59, __pyx_L1_error)
+  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_vec_gammas), __pyx_ptype_5numpy_ndarray, 1, "vec_gammas", 0))) __PYX_ERR(0, 62, __pyx_L1_error)
+  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_index), __pyx_ptype_5numpy_ndarray, 1, "index", 0))) __PYX_ERR(0, 65, __pyx_L1_error)
   __pyx_r = __pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_test(((struct __pyx_obj_3src_9screening_12gap_test_all_GapTestAll *)__pyx_v_self), __pyx_v_Atcabs, __pyx_v_gap, __pyx_v_lbd, __pyx_v_vec_gammas, __pyx_v_coeff_dual_scaling, __pyx_v_offset_radius, __pyx_v_index);
 
-  /* "src/screening/gap_test_all.pyx":57
+  /* "src/screening/gap_test_all.pyx":58
  *    @cython.boundscheck(False) # turn off bounds-checking for entire function
  *    # @cython.wraparound(False)  # turn off negative index wrapping for entire function
  *    def apply_test(self,             # <<<<<<<<<<<<<<
@@ -2372,35 +2360,35 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
   __pyx_pybuffernd_index.rcbuffer = &__pyx_pybuffer_index;
   {
     __Pyx_BufFmt_StackElem __pyx_stack[1];
-    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_Atcabs.rcbuffer->pybuffer, (PyObject*)__pyx_v_Atcabs, &__Pyx_TypeInfo_nn_npy_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 57, __pyx_L1_error)
+    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_Atcabs.rcbuffer->pybuffer, (PyObject*)__pyx_v_Atcabs, &__Pyx_TypeInfo_nn_npy_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 58, __pyx_L1_error)
   }
   __pyx_pybuffernd_Atcabs.diminfo[0].strides = __pyx_pybuffernd_Atcabs.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_Atcabs.diminfo[0].shape = __pyx_pybuffernd_Atcabs.rcbuffer->pybuffer.shape[0];
   {
     __Pyx_BufFmt_StackElem __pyx_stack[1];
-    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_vec_gammas.rcbuffer->pybuffer, (PyObject*)__pyx_v_vec_gammas, &__Pyx_TypeInfo_nn_npy_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 57, __pyx_L1_error)
+    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_vec_gammas.rcbuffer->pybuffer, (PyObject*)__pyx_v_vec_gammas, &__Pyx_TypeInfo_nn_npy_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 58, __pyx_L1_error)
   }
   __pyx_pybuffernd_vec_gammas.diminfo[0].strides = __pyx_pybuffernd_vec_gammas.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_vec_gammas.diminfo[0].shape = __pyx_pybuffernd_vec_gammas.rcbuffer->pybuffer.shape[0];
   {
     __Pyx_BufFmt_StackElem __pyx_stack[1];
-    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_index.rcbuffer->pybuffer, (PyObject*)__pyx_v_index, &__Pyx_TypeInfo_long, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 57, __pyx_L1_error)
+    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_index.rcbuffer->pybuffer, (PyObject*)__pyx_v_index, &__Pyx_TypeInfo_long, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 58, __pyx_L1_error)
   }
   __pyx_pybuffernd_index.diminfo[0].strides = __pyx_pybuffernd_index.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_index.diminfo[0].shape = __pyx_pybuffernd_index.rcbuffer->pybuffer.shape[0];
 
-  /* "src/screening/gap_test_all.pyx":106
+  /* "src/screening/gap_test_all.pyx":107
  *       # in order to match the paper indexation rules
  * 
  *       cdef double radius  = coeff_dual_scaling * np.sqrt(2 * gap) + offset_radius             # <<<<<<<<<<<<<<
  *       # cdef double lbd_aug = coeff_dual_scaling * lbd
  *       cdef np.ndarray[np.npy_double, ndim=1] coeff_lbd_gamma = coeff_dual_scaling * lbd * vec_gammas
  */
-  __pyx_t_1 = PyFloat_FromDouble(__pyx_v_coeff_dual_scaling); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 106, __pyx_L1_error)
+  __pyx_t_1 = PyFloat_FromDouble(__pyx_v_coeff_dual_scaling); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 107, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
-  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 106, __pyx_L1_error)
+  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 107, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_3);
-  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_sqrt); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 106, __pyx_L1_error)
+  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_sqrt); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 107, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_4);
   __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
-  __pyx_t_3 = PyFloat_FromDouble((2.0 * __pyx_v_gap)); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 106, __pyx_L1_error)
+  __pyx_t_3 = PyFloat_FromDouble((2.0 * __pyx_v_gap)); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 107, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_3);
   __pyx_t_5 = NULL;
   if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_4))) {
@@ -2415,42 +2403,42 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
   __pyx_t_2 = (__pyx_t_5) ? __Pyx_PyObject_Call2Args(__pyx_t_4, __pyx_t_5, __pyx_t_3) : __Pyx_PyObject_CallOneArg(__pyx_t_4, __pyx_t_3);
   __Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
   __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
-  if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 106, __pyx_L1_error)
+  if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 107, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_2);
   __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
-  __pyx_t_4 = PyNumber_Multiply(__pyx_t_1, __pyx_t_2); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 106, __pyx_L1_error)
+  __pyx_t_4 = PyNumber_Multiply(__pyx_t_1, __pyx_t_2); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 107, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_4);
   __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
   __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
-  __pyx_t_2 = PyFloat_FromDouble(__pyx_v_offset_radius); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 106, __pyx_L1_error)
+  __pyx_t_2 = PyFloat_FromDouble(__pyx_v_offset_radius); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 107, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_2);
-  __pyx_t_1 = PyNumber_Add(__pyx_t_4, __pyx_t_2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 106, __pyx_L1_error)
+  __pyx_t_1 = PyNumber_Add(__pyx_t_4, __pyx_t_2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 107, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
   __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
   __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
-  __pyx_t_6 = __pyx_PyFloat_AsDouble(__pyx_t_1); if (unlikely((__pyx_t_6 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 106, __pyx_L1_error)
+  __pyx_t_6 = __pyx_PyFloat_AsDouble(__pyx_t_1); if (unlikely((__pyx_t_6 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 107, __pyx_L1_error)
   __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
   __pyx_v_radius = __pyx_t_6;
 
-  /* "src/screening/gap_test_all.pyx":108
+  /* "src/screening/gap_test_all.pyx":109
  *       cdef double radius  = coeff_dual_scaling * np.sqrt(2 * gap) + offset_radius
  *       # cdef double lbd_aug = coeff_dual_scaling * lbd
  *       cdef np.ndarray[np.npy_double, ndim=1] coeff_lbd_gamma = coeff_dual_scaling * lbd * vec_gammas             # <<<<<<<<<<<<<<
  * 
  *       cdef int n = Atcabs.shape[0]
  */
-  __pyx_t_1 = PyFloat_FromDouble((__pyx_v_coeff_dual_scaling * __pyx_v_lbd)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 108, __pyx_L1_error)
+  __pyx_t_1 = PyFloat_FromDouble((__pyx_v_coeff_dual_scaling * __pyx_v_lbd)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 109, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
-  __pyx_t_2 = PyNumber_Multiply(__pyx_t_1, ((PyObject *)__pyx_v_vec_gammas)); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 108, __pyx_L1_error)
+  __pyx_t_2 = PyNumber_Multiply(__pyx_t_1, ((PyObject *)__pyx_v_vec_gammas)); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 109, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_2);
   __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
-  if (!(likely(((__pyx_t_2) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_2, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 108, __pyx_L1_error)
+  if (!(likely(((__pyx_t_2) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_2, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 109, __pyx_L1_error)
   __pyx_t_7 = ((PyArrayObject *)__pyx_t_2);
   {
     __Pyx_BufFmt_StackElem __pyx_stack[1];
     if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_coeff_lbd_gamma.rcbuffer->pybuffer, (PyObject*)__pyx_t_7, &__Pyx_TypeInfo_nn_npy_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
       __pyx_v_coeff_lbd_gamma = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_coeff_lbd_gamma.rcbuffer->pybuffer.buf = NULL;
-      __PYX_ERR(0, 108, __pyx_L1_error)
+      __PYX_ERR(0, 109, __pyx_L1_error)
     } else {__pyx_pybuffernd_coeff_lbd_gamma.diminfo[0].strides = __pyx_pybuffernd_coeff_lbd_gamma.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_coeff_lbd_gamma.diminfo[0].shape = __pyx_pybuffernd_coeff_lbd_gamma.rcbuffer->pybuffer.shape[0];
     }
   }
@@ -2458,7 +2446,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
   __pyx_v_coeff_lbd_gamma = ((PyArrayObject *)__pyx_t_2);
   __pyx_t_2 = 0;
 
-  /* "src/screening/gap_test_all.pyx":110
+  /* "src/screening/gap_test_all.pyx":111
  *       cdef np.ndarray[np.npy_double, ndim=1] coeff_lbd_gamma = coeff_dual_scaling * lbd * vec_gammas
  * 
  *       cdef int n = Atcabs.shape[0]             # <<<<<<<<<<<<<<
@@ -2467,7 +2455,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
  */
   __pyx_v_n = (__pyx_v_Atcabs->dimensions[0]);
 
-  /* "src/screening/gap_test_all.pyx":112
+  /* "src/screening/gap_test_all.pyx":113
  *       cdef int n = Atcabs.shape[0]
  *       cdef int k, q, r
  *       cdef double tau = 0             # <<<<<<<<<<<<<<
@@ -2476,40 +2464,40 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
  */
   __pyx_v_tau = 0.0;
 
-  /* "src/screening/gap_test_all.pyx":116
+  /* "src/screening/gap_test_all.pyx":117
  * 
  *       # arange because I want entries 0, 1 being 0 and 1!
  *       cdef np.ndarray[np.int_t, ndim=1] vec_p_star = np.arange(n+1, dtype=long)             # <<<<<<<<<<<<<<
  *       cdef np.ndarray[np.int_t, ndim=1] vec_q_star = np.arange(n+1, dtype=long)
  * 
  */
-  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 116, __pyx_L1_error)
+  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 117, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_2);
-  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_arange); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 116, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_arange); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 117, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
   __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
-  __pyx_t_2 = __Pyx_PyInt_From_long((__pyx_v_n + 1)); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 116, __pyx_L1_error)
+  __pyx_t_2 = __Pyx_PyInt_From_long((__pyx_v_n + 1)); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 117, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_2);
-  __pyx_t_4 = PyTuple_New(1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 116, __pyx_L1_error)
+  __pyx_t_4 = PyTuple_New(1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 117, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_4);
   __Pyx_GIVEREF(__pyx_t_2);
   PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_2);
   __pyx_t_2 = 0;
-  __pyx_t_2 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 116, __pyx_L1_error)
+  __pyx_t_2 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 117, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_2);
-  if (PyDict_SetItem(__pyx_t_2, __pyx_n_s_dtype, ((PyObject *)(&PyLong_Type))) < 0) __PYX_ERR(0, 116, __pyx_L1_error)
-  __pyx_t_3 = __Pyx_PyObject_Call(__pyx_t_1, __pyx_t_4, __pyx_t_2); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 116, __pyx_L1_error)
+  if (PyDict_SetItem(__pyx_t_2, __pyx_n_s_dtype, ((PyObject *)(&PyLong_Type))) < 0) __PYX_ERR(0, 117, __pyx_L1_error)
+  __pyx_t_3 = __Pyx_PyObject_Call(__pyx_t_1, __pyx_t_4, __pyx_t_2); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 117, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_3);
   __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
   __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
   __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
-  if (!(likely(((__pyx_t_3) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_3, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 116, __pyx_L1_error)
+  if (!(likely(((__pyx_t_3) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_3, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 117, __pyx_L1_error)
   __pyx_t_8 = ((PyArrayObject *)__pyx_t_3);
   {
     __Pyx_BufFmt_StackElem __pyx_stack[1];
     if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_vec_p_star.rcbuffer->pybuffer, (PyObject*)__pyx_t_8, &__Pyx_TypeInfo_nn___pyx_t_5numpy_int_t, PyBUF_FORMAT| PyBUF_STRIDES| PyBUF_WRITABLE, 1, 0, __pyx_stack) == -1)) {
       __pyx_v_vec_p_star = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_vec_p_star.rcbuffer->pybuffer.buf = NULL;
-      __PYX_ERR(0, 116, __pyx_L1_error)
+      __PYX_ERR(0, 117, __pyx_L1_error)
     } else {__pyx_pybuffernd_vec_p_star.diminfo[0].strides = __pyx_pybuffernd_vec_p_star.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_vec_p_star.diminfo[0].shape = __pyx_pybuffernd_vec_p_star.rcbuffer->pybuffer.shape[0];
     }
   }
@@ -2517,40 +2505,40 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
   __pyx_v_vec_p_star = ((PyArrayObject *)__pyx_t_3);
   __pyx_t_3 = 0;
 
-  /* "src/screening/gap_test_all.pyx":117
+  /* "src/screening/gap_test_all.pyx":118
  *       # arange because I want entries 0, 1 being 0 and 1!
  *       cdef np.ndarray[np.int_t, ndim=1] vec_p_star = np.arange(n+1, dtype=long)
  *       cdef np.ndarray[np.int_t, ndim=1] vec_q_star = np.arange(n+1, dtype=long)             # <<<<<<<<<<<<<<
  * 
  *       # output
  */
-  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 117, __pyx_L1_error)
+  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 118, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_3);
-  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_arange); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 117, __pyx_L1_error)
+  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_arange); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 118, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_2);
   __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
-  __pyx_t_3 = __Pyx_PyInt_From_long((__pyx_v_n + 1)); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 117, __pyx_L1_error)
+  __pyx_t_3 = __Pyx_PyInt_From_long((__pyx_v_n + 1)); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 118, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_3);
-  __pyx_t_4 = PyTuple_New(1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 117, __pyx_L1_error)
+  __pyx_t_4 = PyTuple_New(1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 118, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_4);
   __Pyx_GIVEREF(__pyx_t_3);
   PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_3);
   __pyx_t_3 = 0;
-  __pyx_t_3 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 117, __pyx_L1_error)
+  __pyx_t_3 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 118, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_3);
-  if (PyDict_SetItem(__pyx_t_3, __pyx_n_s_dtype, ((PyObject *)(&PyLong_Type))) < 0) __PYX_ERR(0, 117, __pyx_L1_error)
-  __pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_2, __pyx_t_4, __pyx_t_3); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 117, __pyx_L1_error)
+  if (PyDict_SetItem(__pyx_t_3, __pyx_n_s_dtype, ((PyObject *)(&PyLong_Type))) < 0) __PYX_ERR(0, 118, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_2, __pyx_t_4, __pyx_t_3); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 118, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
   __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
   __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
   __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
-  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 117, __pyx_L1_error)
+  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 118, __pyx_L1_error)
   __pyx_t_9 = ((PyArrayObject *)__pyx_t_1);
   {
     __Pyx_BufFmt_StackElem __pyx_stack[1];
     if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_vec_q_star.rcbuffer->pybuffer, (PyObject*)__pyx_t_9, &__Pyx_TypeInfo_nn___pyx_t_5numpy_int_t, PyBUF_FORMAT| PyBUF_STRIDES| PyBUF_WRITABLE, 1, 0, __pyx_stack) == -1)) {
       __pyx_v_vec_q_star = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_vec_q_star.rcbuffer->pybuffer.buf = NULL;
-      __PYX_ERR(0, 117, __pyx_L1_error)
+      __PYX_ERR(0, 118, __pyx_L1_error)
     } else {__pyx_pybuffernd_vec_q_star.diminfo[0].strides = __pyx_pybuffernd_vec_q_star.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_vec_q_star.diminfo[0].shape = __pyx_pybuffernd_vec_q_star.rcbuffer->pybuffer.shape[0];
     }
   }
@@ -2558,40 +2546,40 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
   __pyx_v_vec_q_star = ((PyArrayObject *)__pyx_t_1);
   __pyx_t_1 = 0;
 
-  /* "src/screening/gap_test_all.pyx":120
+  /* "src/screening/gap_test_all.pyx":121
  * 
  *       # output
  *       cdef np.ndarray[np.npy_bool, ndim=1] calI_screen = np.zeros(n, dtype=bool)             # <<<<<<<<<<<<<<
  * 
  *       # 1. Sort in descending order
  */
-  __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 120, __pyx_L1_error)
+  __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 121, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
-  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_zeros); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 120, __pyx_L1_error)
+  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_zeros); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 121, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_3);
   __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
-  __pyx_t_1 = __Pyx_PyInt_From_int(__pyx_v_n); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 120, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyInt_From_int(__pyx_v_n); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 121, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
-  __pyx_t_4 = PyTuple_New(1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 120, __pyx_L1_error)
+  __pyx_t_4 = PyTuple_New(1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 121, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_4);
   __Pyx_GIVEREF(__pyx_t_1);
   PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_1);
   __pyx_t_1 = 0;
-  __pyx_t_1 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 120, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 121, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
-  if (PyDict_SetItem(__pyx_t_1, __pyx_n_s_dtype, ((PyObject*)&PyBool_Type)) < 0) __PYX_ERR(0, 120, __pyx_L1_error)
-  __pyx_t_2 = __Pyx_PyObject_Call(__pyx_t_3, __pyx_t_4, __pyx_t_1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 120, __pyx_L1_error)
+  if (PyDict_SetItem(__pyx_t_1, __pyx_n_s_dtype, ((PyObject*)&PyBool_Type)) < 0) __PYX_ERR(0, 121, __pyx_L1_error)
+  __pyx_t_2 = __Pyx_PyObject_Call(__pyx_t_3, __pyx_t_4, __pyx_t_1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 121, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_2);
   __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
   __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
   __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
-  if (!(likely(((__pyx_t_2) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_2, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 120, __pyx_L1_error)
+  if (!(likely(((__pyx_t_2) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_2, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 121, __pyx_L1_error)
   __pyx_t_10 = ((PyArrayObject *)__pyx_t_2);
   {
     __Pyx_BufFmt_StackElem __pyx_stack[1];
     if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_calI_screen.rcbuffer->pybuffer, (PyObject*)__pyx_t_10, &__Pyx_TypeInfo_nn_npy_bool, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
       __pyx_v_calI_screen = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_calI_screen.rcbuffer->pybuffer.buf = NULL;
-      __PYX_ERR(0, 120, __pyx_L1_error)
+      __PYX_ERR(0, 121, __pyx_L1_error)
     } else {__pyx_pybuffernd_calI_screen.diminfo[0].strides = __pyx_pybuffernd_calI_screen.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_calI_screen.diminfo[0].shape = __pyx_pybuffernd_calI_screen.rcbuffer->pybuffer.shape[0];
     }
   }
@@ -2599,7 +2587,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
   __pyx_v_calI_screen = ((PyArrayObject *)__pyx_t_2);
   __pyx_t_2 = 0;
 
-  /* "src/screening/gap_test_all.pyx":123
+  /* "src/screening/gap_test_all.pyx":124
  * 
  *       # 1. Sort in descending order
  *       if index is None:             # <<<<<<<<<<<<<<
@@ -2610,16 +2598,16 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
   __pyx_t_12 = (__pyx_t_11 != 0);
   if (__pyx_t_12) {
 
-    /* "src/screening/gap_test_all.pyx":124
+    /* "src/screening/gap_test_all.pyx":125
  *       # 1. Sort in descending order
  *       if index is None:
  *          index = np.argsort(Atcabs)[::-1]             # <<<<<<<<<<<<<<
  * 
  *       # 2. Precomputing quantities
  */
-    __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 124, __pyx_L1_error)
+    __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 125, __pyx_L1_error)
     __Pyx_GOTREF(__pyx_t_1);
-    __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_argsort); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 124, __pyx_L1_error)
+    __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_argsort); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 125, __pyx_L1_error)
     __Pyx_GOTREF(__pyx_t_4);
     __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
     __pyx_t_1 = NULL;
@@ -2634,13 +2622,13 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
     }
     __pyx_t_2 = (__pyx_t_1) ? __Pyx_PyObject_Call2Args(__pyx_t_4, __pyx_t_1, ((PyObject *)__pyx_v_Atcabs)) : __Pyx_PyObject_CallOneArg(__pyx_t_4, ((PyObject *)__pyx_v_Atcabs));
     __Pyx_XDECREF(__pyx_t_1); __pyx_t_1 = 0;
-    if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 124, __pyx_L1_error)
+    if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 125, __pyx_L1_error)
     __Pyx_GOTREF(__pyx_t_2);
     __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
-    __pyx_t_4 = __Pyx_PyObject_GetItem(__pyx_t_2, __pyx_slice_); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 124, __pyx_L1_error)
+    __pyx_t_4 = __Pyx_PyObject_GetItem(__pyx_t_2, __pyx_slice_); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 125, __pyx_L1_error)
     __Pyx_GOTREF(__pyx_t_4);
     __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
-    if (!(likely(((__pyx_t_4) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_4, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 124, __pyx_L1_error)
+    if (!(likely(((__pyx_t_4) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_4, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 125, __pyx_L1_error)
     __pyx_t_13 = ((PyArrayObject *)__pyx_t_4);
     {
       __Pyx_BufFmt_StackElem __pyx_stack[1];
@@ -2657,13 +2645,13 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
         __pyx_t_15 = __pyx_t_16 = __pyx_t_17 = 0;
       }
       __pyx_pybuffernd_index.diminfo[0].strides = __pyx_pybuffernd_index.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_index.diminfo[0].shape = __pyx_pybuffernd_index.rcbuffer->pybuffer.shape[0];
-      if (unlikely(__pyx_t_14 < 0)) __PYX_ERR(0, 124, __pyx_L1_error)
+      if (unlikely(__pyx_t_14 < 0)) __PYX_ERR(0, 125, __pyx_L1_error)
     }
     __pyx_t_13 = 0;
     __Pyx_DECREF_SET(__pyx_v_index, ((PyArrayObject *)__pyx_t_4));
     __pyx_t_4 = 0;
 
-    /* "src/screening/gap_test_all.pyx":123
+    /* "src/screening/gap_test_all.pyx":124
  * 
  *       # 1. Sort in descending order
  *       if index is None:             # <<<<<<<<<<<<<<
@@ -2672,40 +2660,40 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
  */
   }
 
-  /* "src/screening/gap_test_all.pyx":127
+  /* "src/screening/gap_test_all.pyx":128
  * 
  *       # 2. Precomputing quantities
  *       cdef np.ndarray[np.npy_double, ndim=1] vec_f = np.cumsum(             # <<<<<<<<<<<<<<
  *          coeff_lbd_gamma[::-1] - Atcabs[index[::-1]] - radius
  *       )[::-1]
  */
-  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 127, __pyx_L1_error)
+  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 128, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_2);
-  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_cumsum); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 127, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_cumsum); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 128, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
   __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
 
-  /* "src/screening/gap_test_all.pyx":128
+  /* "src/screening/gap_test_all.pyx":129
  *       # 2. Precomputing quantities
  *       cdef np.ndarray[np.npy_double, ndim=1] vec_f = np.cumsum(
  *          coeff_lbd_gamma[::-1] - Atcabs[index[::-1]] - radius             # <<<<<<<<<<<<<<
  *       )[::-1]
  * 
  */
-  __pyx_t_2 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_coeff_lbd_gamma), __pyx_slice_); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 128, __pyx_L1_error)
+  __pyx_t_2 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_coeff_lbd_gamma), __pyx_slice_); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 129, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_2);
-  __pyx_t_3 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_index), __pyx_slice_); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 128, __pyx_L1_error)
+  __pyx_t_3 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_index), __pyx_slice_); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 129, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_3);
-  __pyx_t_5 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_Atcabs), __pyx_t_3); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 128, __pyx_L1_error)
+  __pyx_t_5 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_Atcabs), __pyx_t_3); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 129, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_5);
   __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
-  __pyx_t_3 = PyNumber_Subtract(__pyx_t_2, __pyx_t_5); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 128, __pyx_L1_error)
+  __pyx_t_3 = PyNumber_Subtract(__pyx_t_2, __pyx_t_5); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 129, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_3);
   __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
   __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
-  __pyx_t_5 = PyFloat_FromDouble(__pyx_v_radius); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 128, __pyx_L1_error)
+  __pyx_t_5 = PyFloat_FromDouble(__pyx_v_radius); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 129, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_5);
-  __pyx_t_2 = PyNumber_Subtract(__pyx_t_3, __pyx_t_5); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 128, __pyx_L1_error)
+  __pyx_t_2 = PyNumber_Subtract(__pyx_t_3, __pyx_t_5); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 129, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_2);
   __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
   __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
@@ -2722,27 +2710,27 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
   __pyx_t_4 = (__pyx_t_5) ? __Pyx_PyObject_Call2Args(__pyx_t_1, __pyx_t_5, __pyx_t_2) : __Pyx_PyObject_CallOneArg(__pyx_t_1, __pyx_t_2);
   __Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
   __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
-  if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 127, __pyx_L1_error)
+  if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 128, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_4);
   __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
 
-  /* "src/screening/gap_test_all.pyx":129
+  /* "src/screening/gap_test_all.pyx":130
  *       cdef np.ndarray[np.npy_double, ndim=1] vec_f = np.cumsum(
  *          coeff_lbd_gamma[::-1] - Atcabs[index[::-1]] - radius
  *       )[::-1]             # <<<<<<<<<<<<<<
  * 
  *       cdef double best_bound_q = vec_f[0] - coeff_lbd_gamma[0]
  */
-  __pyx_t_1 = __Pyx_PyObject_GetItem(__pyx_t_4, __pyx_slice_); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 129, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyObject_GetItem(__pyx_t_4, __pyx_slice_); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 130, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
   __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
-  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 129, __pyx_L1_error)
+  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 130, __pyx_L1_error)
   __pyx_t_18 = ((PyArrayObject *)__pyx_t_1);
   {
     __Pyx_BufFmt_StackElem __pyx_stack[1];
     if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_vec_f.rcbuffer->pybuffer, (PyObject*)__pyx_t_18, &__Pyx_TypeInfo_nn_npy_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
       __pyx_v_vec_f = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_vec_f.rcbuffer->pybuffer.buf = NULL;
-      __PYX_ERR(0, 127, __pyx_L1_error)
+      __PYX_ERR(0, 128, __pyx_L1_error)
     } else {__pyx_pybuffernd_vec_f.diminfo[0].strides = __pyx_pybuffernd_vec_f.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_vec_f.diminfo[0].shape = __pyx_pybuffernd_vec_f.rcbuffer->pybuffer.shape[0];
     }
   }
@@ -2750,7 +2738,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
   __pyx_v_vec_f = ((PyArrayObject *)__pyx_t_1);
   __pyx_t_1 = 0;
 
-  /* "src/screening/gap_test_all.pyx":131
+  /* "src/screening/gap_test_all.pyx":132
  *       )[::-1]
  * 
  *       cdef double best_bound_q = vec_f[0] - coeff_lbd_gamma[0]             # <<<<<<<<<<<<<<
@@ -2763,7 +2751,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
   if (__pyx_t_20 < 0) __pyx_t_20 += __pyx_pybuffernd_coeff_lbd_gamma.diminfo[0].shape;
   __pyx_v_best_bound_q = ((*__Pyx_BufPtrStrided1d(npy_double *, __pyx_pybuffernd_vec_f.rcbuffer->pybuffer.buf, __pyx_t_19, __pyx_pybuffernd_vec_f.diminfo[0].strides)) - (*__Pyx_BufPtrStrided1d(npy_double *, __pyx_pybuffernd_coeff_lbd_gamma.rcbuffer->pybuffer.buf, __pyx_t_20, __pyx_pybuffernd_coeff_lbd_gamma.diminfo[0].strides)));
 
-  /* "src/screening/gap_test_all.pyx":132
+  /* "src/screening/gap_test_all.pyx":133
  * 
  *       cdef double best_bound_q = vec_f[0] - coeff_lbd_gamma[0]
  *       cdef double curr_bound_q = 0.             # <<<<<<<<<<<<<<
@@ -2772,7 +2760,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
  */
   __pyx_v_curr_bound_q = 0.;
 
-  /* "src/screening/gap_test_all.pyx":133
+  /* "src/screening/gap_test_all.pyx":134
  *       cdef double best_bound_q = vec_f[0] - coeff_lbd_gamma[0]
  *       cdef double curr_bound_q = 0.
  *       cdef double best_bound_p = vec_f[0]             # <<<<<<<<<<<<<<
@@ -2783,7 +2771,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
   if (__pyx_t_20 < 0) __pyx_t_20 += __pyx_pybuffernd_vec_f.diminfo[0].shape;
   __pyx_v_best_bound_p = (*__Pyx_BufPtrStrided1d(npy_double *, __pyx_pybuffernd_vec_f.rcbuffer->pybuffer.buf, __pyx_t_20, __pyx_pybuffernd_vec_f.diminfo[0].strides));
 
-  /* "src/screening/gap_test_all.pyx":134
+  /* "src/screening/gap_test_all.pyx":135
  *       cdef double curr_bound_q = 0.
  *       cdef double best_bound_p = vec_f[0]
  *       cdef double curr_bound_p = 0.             # <<<<<<<<<<<<<<
@@ -2792,7 +2780,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
  */
   __pyx_v_curr_bound_p = 0.;
 
-  /* "src/screening/gap_test_all.pyx":135
+  /* "src/screening/gap_test_all.pyx":136
  *       cdef double best_bound_p = vec_f[0]
  *       cdef double curr_bound_p = 0.
  *       for k in range(2, n+1):             # <<<<<<<<<<<<<<
@@ -2804,7 +2792,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
   for (__pyx_t_14 = 2; __pyx_t_14 < __pyx_t_22; __pyx_t_14+=1) {
     __pyx_v_k = __pyx_t_14;
 
-    /* "src/screening/gap_test_all.pyx":139
+    /* "src/screening/gap_test_all.pyx":140
  * 
  *          # 1. Evaluate p*
  *          curr_bound_p = vec_f[k-1]             # <<<<<<<<<<<<<<
@@ -2815,7 +2803,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
     if (__pyx_t_20 < 0) __pyx_t_20 += __pyx_pybuffernd_vec_f.diminfo[0].shape;
     __pyx_v_curr_bound_p = (*__Pyx_BufPtrStrided1d(npy_double *, __pyx_pybuffernd_vec_f.rcbuffer->pybuffer.buf, __pyx_t_20, __pyx_pybuffernd_vec_f.diminfo[0].strides));
 
-    /* "src/screening/gap_test_all.pyx":140
+    /* "src/screening/gap_test_all.pyx":141
  *          # 1. Evaluate p*
  *          curr_bound_p = vec_f[k-1]
  *          if curr_bound_p > best_bound_p:             # <<<<<<<<<<<<<<
@@ -2825,7 +2813,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
     __pyx_t_12 = ((__pyx_v_curr_bound_p > __pyx_v_best_bound_p) != 0);
     if (__pyx_t_12) {
 
-      /* "src/screening/gap_test_all.pyx":141
+      /* "src/screening/gap_test_all.pyx":142
  *          curr_bound_p = vec_f[k-1]
  *          if curr_bound_p > best_bound_p:
  *             best_bound_p = curr_bound_p             # <<<<<<<<<<<<<<
@@ -2834,7 +2822,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
  */
       __pyx_v_best_bound_p = __pyx_v_curr_bound_p;
 
-      /* "src/screening/gap_test_all.pyx":142
+      /* "src/screening/gap_test_all.pyx":143
  *          if curr_bound_p > best_bound_p:
  *             best_bound_p = curr_bound_p
  *             vec_p_star[k] = k             # <<<<<<<<<<<<<<
@@ -2845,7 +2833,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
       if (__pyx_t_20 < 0) __pyx_t_20 += __pyx_pybuffernd_vec_p_star.diminfo[0].shape;
       *__Pyx_BufPtrStrided1d(__pyx_t_5numpy_int_t *, __pyx_pybuffernd_vec_p_star.rcbuffer->pybuffer.buf, __pyx_t_20, __pyx_pybuffernd_vec_p_star.diminfo[0].strides) = __pyx_v_k;
 
-      /* "src/screening/gap_test_all.pyx":140
+      /* "src/screening/gap_test_all.pyx":141
  *          # 1. Evaluate p*
  *          curr_bound_p = vec_f[k-1]
  *          if curr_bound_p > best_bound_p:             # <<<<<<<<<<<<<<
@@ -2855,7 +2843,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
       goto __pyx_L6;
     }
 
-    /* "src/screening/gap_test_all.pyx":144
+    /* "src/screening/gap_test_all.pyx":145
  *             vec_p_star[k] = k
  *          else:
  *             vec_p_star[k] = vec_p_star[k-1]             # <<<<<<<<<<<<<<
@@ -2871,7 +2859,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
     }
     __pyx_L6:;
 
-    /* "src/screening/gap_test_all.pyx":147
+    /* "src/screening/gap_test_all.pyx":148
  * 
  *          # 1. Evaluate q*
  *          curr_bound_q = vec_f[k-1] - coeff_lbd_gamma[k-1]             # <<<<<<<<<<<<<<
@@ -2884,7 +2872,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
     if (__pyx_t_19 < 0) __pyx_t_19 += __pyx_pybuffernd_coeff_lbd_gamma.diminfo[0].shape;
     __pyx_v_curr_bound_q = ((*__Pyx_BufPtrStrided1d(npy_double *, __pyx_pybuffernd_vec_f.rcbuffer->pybuffer.buf, __pyx_t_20, __pyx_pybuffernd_vec_f.diminfo[0].strides)) - (*__Pyx_BufPtrStrided1d(npy_double *, __pyx_pybuffernd_coeff_lbd_gamma.rcbuffer->pybuffer.buf, __pyx_t_19, __pyx_pybuffernd_coeff_lbd_gamma.diminfo[0].strides)));
 
-    /* "src/screening/gap_test_all.pyx":148
+    /* "src/screening/gap_test_all.pyx":149
  *          # 1. Evaluate q*
  *          curr_bound_q = vec_f[k-1] - coeff_lbd_gamma[k-1]
  *          if curr_bound_q > best_bound_q:             # <<<<<<<<<<<<<<
@@ -2894,7 +2882,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
     __pyx_t_12 = ((__pyx_v_curr_bound_q > __pyx_v_best_bound_q) != 0);
     if (__pyx_t_12) {
 
-      /* "src/screening/gap_test_all.pyx":149
+      /* "src/screening/gap_test_all.pyx":150
  *          curr_bound_q = vec_f[k-1] - coeff_lbd_gamma[k-1]
  *          if curr_bound_q > best_bound_q:
  *             best_bound_q = curr_bound_q             # <<<<<<<<<<<<<<
@@ -2903,7 +2891,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
  */
       __pyx_v_best_bound_q = __pyx_v_curr_bound_q;
 
-      /* "src/screening/gap_test_all.pyx":150
+      /* "src/screening/gap_test_all.pyx":151
  *          if curr_bound_q > best_bound_q:
  *             best_bound_q = curr_bound_q
  *             vec_q_star[k] = k             # <<<<<<<<<<<<<<
@@ -2914,7 +2902,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
       if (__pyx_t_19 < 0) __pyx_t_19 += __pyx_pybuffernd_vec_q_star.diminfo[0].shape;
       *__Pyx_BufPtrStrided1d(__pyx_t_5numpy_int_t *, __pyx_pybuffernd_vec_q_star.rcbuffer->pybuffer.buf, __pyx_t_19, __pyx_pybuffernd_vec_q_star.diminfo[0].strides) = __pyx_v_k;
 
-      /* "src/screening/gap_test_all.pyx":148
+      /* "src/screening/gap_test_all.pyx":149
  *          # 1. Evaluate q*
  *          curr_bound_q = vec_f[k-1] - coeff_lbd_gamma[k-1]
  *          if curr_bound_q > best_bound_q:             # <<<<<<<<<<<<<<
@@ -2924,7 +2912,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
       goto __pyx_L7;
     }
 
-    /* "src/screening/gap_test_all.pyx":152
+    /* "src/screening/gap_test_all.pyx":153
  *             vec_q_star[k] = k
  *          else:
  *             vec_q_star[k] = vec_q_star[k-1]             # <<<<<<<<<<<<<<
@@ -2941,16 +2929,16 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
     __pyx_L7:;
   }
 
-  /* "src/screening/gap_test_all.pyx":155
+  /* "src/screening/gap_test_all.pyx":156
  * 
  *       # 3. Tests
  *       for l in range(n, 0, -1):             # <<<<<<<<<<<<<<
  *          q = vec_q_star[l]
  *          while q >= 1:
  */
-  __pyx_t_1 = __Pyx_PyInt_From_int(__pyx_v_n); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 155, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyInt_From_int(__pyx_v_n); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 156, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
-  __pyx_t_4 = PyTuple_New(3); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 155, __pyx_L1_error)
+  __pyx_t_4 = PyTuple_New(3); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 156, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_4);
   __Pyx_GIVEREF(__pyx_t_1);
   PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_1);
@@ -2961,16 +2949,16 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
   __Pyx_GIVEREF(__pyx_int_neg_1);
   PyTuple_SET_ITEM(__pyx_t_4, 2, __pyx_int_neg_1);
   __pyx_t_1 = 0;
-  __pyx_t_1 = __Pyx_PyObject_Call(__pyx_builtin_range, __pyx_t_4, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 155, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyObject_Call(__pyx_builtin_range, __pyx_t_4, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 156, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
   __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
   if (likely(PyList_CheckExact(__pyx_t_1)) || PyTuple_CheckExact(__pyx_t_1)) {
     __pyx_t_4 = __pyx_t_1; __Pyx_INCREF(__pyx_t_4); __pyx_t_23 = 0;
     __pyx_t_24 = NULL;
   } else {
-    __pyx_t_23 = -1; __pyx_t_4 = PyObject_GetIter(__pyx_t_1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 155, __pyx_L1_error)
+    __pyx_t_23 = -1; __pyx_t_4 = PyObject_GetIter(__pyx_t_1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 156, __pyx_L1_error)
     __Pyx_GOTREF(__pyx_t_4);
-    __pyx_t_24 = Py_TYPE(__pyx_t_4)->tp_iternext; if (unlikely(!__pyx_t_24)) __PYX_ERR(0, 155, __pyx_L1_error)
+    __pyx_t_24 = Py_TYPE(__pyx_t_4)->tp_iternext; if (unlikely(!__pyx_t_24)) __PYX_ERR(0, 156, __pyx_L1_error)
   }
   __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
   for (;;) {
@@ -2978,17 +2966,17 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
       if (likely(PyList_CheckExact(__pyx_t_4))) {
         if (__pyx_t_23 >= PyList_GET_SIZE(__pyx_t_4)) break;
         #if CYTHON_ASSUME_SAFE_MACROS && !CYTHON_AVOID_BORROWED_REFS
-        __pyx_t_1 = PyList_GET_ITEM(__pyx_t_4, __pyx_t_23); __Pyx_INCREF(__pyx_t_1); __pyx_t_23++; if (unlikely(0 < 0)) __PYX_ERR(0, 155, __pyx_L1_error)
+        __pyx_t_1 = PyList_GET_ITEM(__pyx_t_4, __pyx_t_23); __Pyx_INCREF(__pyx_t_1); __pyx_t_23++; if (unlikely(0 < 0)) __PYX_ERR(0, 156, __pyx_L1_error)
         #else
-        __pyx_t_1 = PySequence_ITEM(__pyx_t_4, __pyx_t_23); __pyx_t_23++; if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 155, __pyx_L1_error)
+        __pyx_t_1 = PySequence_ITEM(__pyx_t_4, __pyx_t_23); __pyx_t_23++; if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 156, __pyx_L1_error)
         __Pyx_GOTREF(__pyx_t_1);
         #endif
       } else {
         if (__pyx_t_23 >= PyTuple_GET_SIZE(__pyx_t_4)) break;
         #if CYTHON_ASSUME_SAFE_MACROS && !CYTHON_AVOID_BORROWED_REFS
-        __pyx_t_1 = PyTuple_GET_ITEM(__pyx_t_4, __pyx_t_23); __Pyx_INCREF(__pyx_t_1); __pyx_t_23++; if (unlikely(0 < 0)) __PYX_ERR(0, 155, __pyx_L1_error)
+        __pyx_t_1 = PyTuple_GET_ITEM(__pyx_t_4, __pyx_t_23); __Pyx_INCREF(__pyx_t_1); __pyx_t_23++; if (unlikely(0 < 0)) __PYX_ERR(0, 156, __pyx_L1_error)
         #else
-        __pyx_t_1 = PySequence_ITEM(__pyx_t_4, __pyx_t_23); __pyx_t_23++; if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 155, __pyx_L1_error)
+        __pyx_t_1 = PySequence_ITEM(__pyx_t_4, __pyx_t_23); __pyx_t_23++; if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 156, __pyx_L1_error)
         __Pyx_GOTREF(__pyx_t_1);
         #endif
       }
@@ -2998,7 +2986,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
         PyObject* exc_type = PyErr_Occurred();
         if (exc_type) {
           if (likely(__Pyx_PyErr_GivenExceptionMatches(exc_type, PyExc_StopIteration))) PyErr_Clear();
-          else __PYX_ERR(0, 155, __pyx_L1_error)
+          else __PYX_ERR(0, 156, __pyx_L1_error)
         }
         break;
       }
@@ -3007,20 +2995,20 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
     __Pyx_XDECREF_SET(__pyx_v_l, __pyx_t_1);
     __pyx_t_1 = 0;
 
-    /* "src/screening/gap_test_all.pyx":156
+    /* "src/screening/gap_test_all.pyx":157
  *       # 3. Tests
  *       for l in range(n, 0, -1):
  *          q = vec_q_star[l]             # <<<<<<<<<<<<<<
  *          while q >= 1:
  *             # Implicit definition of the best value of p given q
  */
-    __pyx_t_1 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_vec_q_star), __pyx_v_l); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 156, __pyx_L1_error)
+    __pyx_t_1 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_vec_q_star), __pyx_v_l); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 157, __pyx_L1_error)
     __Pyx_GOTREF(__pyx_t_1);
-    __pyx_t_14 = __Pyx_PyInt_As_int(__pyx_t_1); if (unlikely((__pyx_t_14 == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 156, __pyx_L1_error)
+    __pyx_t_14 = __Pyx_PyInt_As_int(__pyx_t_1); if (unlikely((__pyx_t_14 == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 157, __pyx_L1_error)
     __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
     __pyx_v_q = __pyx_t_14;
 
-    /* "src/screening/gap_test_all.pyx":157
+    /* "src/screening/gap_test_all.pyx":158
  *       for l in range(n, 0, -1):
  *          q = vec_q_star[l]
  *          while q >= 1:             # <<<<<<<<<<<<<<
@@ -3031,7 +3019,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
       __pyx_t_12 = ((__pyx_v_q >= 1) != 0);
       if (!__pyx_t_12) break;
 
-      /* "src/screening/gap_test_all.pyx":159
+      /* "src/screening/gap_test_all.pyx":160
  *          while q >= 1:
  *             # Implicit definition of the best value of p given q
  *             p = vec_p_star[q]             # <<<<<<<<<<<<<<
@@ -3040,83 +3028,83 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
  */
       __pyx_t_19 = __pyx_v_q;
       if (__pyx_t_19 < 0) __pyx_t_19 += __pyx_pybuffernd_vec_p_star.diminfo[0].shape;
-      __pyx_t_1 = __Pyx_PyInt_From_npy_long((*__Pyx_BufPtrStrided1d(__pyx_t_5numpy_int_t *, __pyx_pybuffernd_vec_p_star.rcbuffer->pybuffer.buf, __pyx_t_19, __pyx_pybuffernd_vec_p_star.diminfo[0].strides))); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 159, __pyx_L1_error)
+      __pyx_t_1 = __Pyx_PyInt_From_npy_long((*__Pyx_BufPtrStrided1d(__pyx_t_5numpy_int_t *, __pyx_pybuffernd_vec_p_star.rcbuffer->pybuffer.buf, __pyx_t_19, __pyx_pybuffernd_vec_p_star.diminfo[0].strides))); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 160, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_1);
       __Pyx_XDECREF_SET(__pyx_v_p, __pyx_t_1);
       __pyx_t_1 = 0;
 
-      /* "src/screening/gap_test_all.pyx":162
+      /* "src/screening/gap_test_all.pyx":163
  * 
  *             # Evaluation of the threshold
  *             tau = vec_f[p-1] - vec_f[q-1] + (coeff_lbd_gamma[q-1] - radius)             # <<<<<<<<<<<<<<
  * 
  *             # Test
  */
-      __pyx_t_1 = __Pyx_PyInt_SubtractObjC(__pyx_v_p, __pyx_int_1, 1, 0, 0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 162, __pyx_L1_error)
+      __pyx_t_1 = __Pyx_PyInt_SubtractObjC(__pyx_v_p, __pyx_int_1, 1, 0, 0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 163, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_1);
-      __pyx_t_2 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_vec_f), __pyx_t_1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 162, __pyx_L1_error)
+      __pyx_t_2 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_vec_f), __pyx_t_1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 163, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_2);
       __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
       __pyx_t_19 = (__pyx_v_q - 1);
       if (__pyx_t_19 < 0) __pyx_t_19 += __pyx_pybuffernd_vec_f.diminfo[0].shape;
-      __pyx_t_1 = PyFloat_FromDouble((*__Pyx_BufPtrStrided1d(npy_double *, __pyx_pybuffernd_vec_f.rcbuffer->pybuffer.buf, __pyx_t_19, __pyx_pybuffernd_vec_f.diminfo[0].strides))); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 162, __pyx_L1_error)
+      __pyx_t_1 = PyFloat_FromDouble((*__Pyx_BufPtrStrided1d(npy_double *, __pyx_pybuffernd_vec_f.rcbuffer->pybuffer.buf, __pyx_t_19, __pyx_pybuffernd_vec_f.diminfo[0].strides))); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 163, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_1);
-      __pyx_t_5 = PyNumber_Subtract(__pyx_t_2, __pyx_t_1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 162, __pyx_L1_error)
+      __pyx_t_5 = PyNumber_Subtract(__pyx_t_2, __pyx_t_1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 163, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_5);
       __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
       __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
       __pyx_t_19 = (__pyx_v_q - 1);
       if (__pyx_t_19 < 0) __pyx_t_19 += __pyx_pybuffernd_coeff_lbd_gamma.diminfo[0].shape;
-      __pyx_t_1 = PyFloat_FromDouble(((*__Pyx_BufPtrStrided1d(npy_double *, __pyx_pybuffernd_coeff_lbd_gamma.rcbuffer->pybuffer.buf, __pyx_t_19, __pyx_pybuffernd_coeff_lbd_gamma.diminfo[0].strides)) - __pyx_v_radius)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 162, __pyx_L1_error)
+      __pyx_t_1 = PyFloat_FromDouble(((*__Pyx_BufPtrStrided1d(npy_double *, __pyx_pybuffernd_coeff_lbd_gamma.rcbuffer->pybuffer.buf, __pyx_t_19, __pyx_pybuffernd_coeff_lbd_gamma.diminfo[0].strides)) - __pyx_v_radius)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 163, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_1);
-      __pyx_t_2 = PyNumber_Add(__pyx_t_5, __pyx_t_1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 162, __pyx_L1_error)
+      __pyx_t_2 = PyNumber_Add(__pyx_t_5, __pyx_t_1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 163, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_2);
       __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
       __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
-      __pyx_t_6 = __pyx_PyFloat_AsDouble(__pyx_t_2); if (unlikely((__pyx_t_6 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 162, __pyx_L1_error)
+      __pyx_t_6 = __pyx_PyFloat_AsDouble(__pyx_t_2); if (unlikely((__pyx_t_6 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 163, __pyx_L1_error)
       __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
       __pyx_v_tau = __pyx_t_6;
 
-      /* "src/screening/gap_test_all.pyx":165
+      /* "src/screening/gap_test_all.pyx":166
  * 
  *             # Test
  *             if Atcabs[index[l-1]] >= tau:             # <<<<<<<<<<<<<<
  *                calI_screen[index[l:]] = True
  *                return calI_screen
  */
-      __pyx_t_2 = __Pyx_PyInt_SubtractObjC(__pyx_v_l, __pyx_int_1, 1, 0, 0); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 165, __pyx_L1_error)
+      __pyx_t_2 = __Pyx_PyInt_SubtractObjC(__pyx_v_l, __pyx_int_1, 1, 0, 0); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 166, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_2);
-      __pyx_t_1 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_index), __pyx_t_2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 165, __pyx_L1_error)
+      __pyx_t_1 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_index), __pyx_t_2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 166, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_1);
       __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
-      __pyx_t_2 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_Atcabs), __pyx_t_1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 165, __pyx_L1_error)
+      __pyx_t_2 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_Atcabs), __pyx_t_1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 166, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_2);
       __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
-      __pyx_t_1 = PyFloat_FromDouble(__pyx_v_tau); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 165, __pyx_L1_error)
+      __pyx_t_1 = PyFloat_FromDouble(__pyx_v_tau); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 166, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_1);
-      __pyx_t_5 = PyObject_RichCompare(__pyx_t_2, __pyx_t_1, Py_GE); __Pyx_XGOTREF(__pyx_t_5); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 165, __pyx_L1_error)
+      __pyx_t_5 = PyObject_RichCompare(__pyx_t_2, __pyx_t_1, Py_GE); __Pyx_XGOTREF(__pyx_t_5); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 166, __pyx_L1_error)
       __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
       __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
-      __pyx_t_12 = __Pyx_PyObject_IsTrue(__pyx_t_5); if (unlikely(__pyx_t_12 < 0)) __PYX_ERR(0, 165, __pyx_L1_error)
+      __pyx_t_12 = __Pyx_PyObject_IsTrue(__pyx_t_5); if (unlikely(__pyx_t_12 < 0)) __PYX_ERR(0, 166, __pyx_L1_error)
       __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
       if (__pyx_t_12) {
 
-        /* "src/screening/gap_test_all.pyx":166
+        /* "src/screening/gap_test_all.pyx":167
  *             # Test
  *             if Atcabs[index[l-1]] >= tau:
  *                calI_screen[index[l:]] = True             # <<<<<<<<<<<<<<
  *                return calI_screen
  * 
  */
-        __pyx_t_5 = PySlice_New(__pyx_v_l, Py_None, Py_None); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 166, __pyx_L1_error)
+        __pyx_t_5 = PySlice_New(__pyx_v_l, Py_None, Py_None); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 167, __pyx_L1_error)
         __Pyx_GOTREF(__pyx_t_5);
-        __pyx_t_1 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_index), __pyx_t_5); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 166, __pyx_L1_error)
+        __pyx_t_1 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_index), __pyx_t_5); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 167, __pyx_L1_error)
         __Pyx_GOTREF(__pyx_t_1);
         __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
-        if (unlikely(PyObject_SetItem(((PyObject *)__pyx_v_calI_screen), __pyx_t_1, Py_True) < 0)) __PYX_ERR(0, 166, __pyx_L1_error)
+        if (unlikely(PyObject_SetItem(((PyObject *)__pyx_v_calI_screen), __pyx_t_1, Py_True) < 0)) __PYX_ERR(0, 167, __pyx_L1_error)
         __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
 
-        /* "src/screening/gap_test_all.pyx":167
+        /* "src/screening/gap_test_all.pyx":168
  *             if Atcabs[index[l-1]] >= tau:
  *                calI_screen[index[l:]] = True
  *                return calI_screen             # <<<<<<<<<<<<<<
@@ -3129,7 +3117,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
         __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
         goto __pyx_L0;
 
-        /* "src/screening/gap_test_all.pyx":165
+        /* "src/screening/gap_test_all.pyx":166
  * 
  *             # Test
  *             if Atcabs[index[l-1]] >= tau:             # <<<<<<<<<<<<<<
@@ -3138,24 +3126,24 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
  */
       }
 
-      /* "src/screening/gap_test_all.pyx":170
+      /* "src/screening/gap_test_all.pyx":171
  * 
  *             # Next critical point
  *             q = vec_q_star[p-1]             # <<<<<<<<<<<<<<
  * 
  * 
  */
-      __pyx_t_1 = __Pyx_PyInt_SubtractObjC(__pyx_v_p, __pyx_int_1, 1, 0, 0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 170, __pyx_L1_error)
+      __pyx_t_1 = __Pyx_PyInt_SubtractObjC(__pyx_v_p, __pyx_int_1, 1, 0, 0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 171, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_1);
-      __pyx_t_5 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_vec_q_star), __pyx_t_1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 170, __pyx_L1_error)
+      __pyx_t_5 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_vec_q_star), __pyx_t_1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 171, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_5);
       __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
-      __pyx_t_14 = __Pyx_PyInt_As_int(__pyx_t_5); if (unlikely((__pyx_t_14 == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 170, __pyx_L1_error)
+      __pyx_t_14 = __Pyx_PyInt_As_int(__pyx_t_5); if (unlikely((__pyx_t_14 == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 171, __pyx_L1_error)
       __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
       __pyx_v_q = __pyx_t_14;
     }
 
-    /* "src/screening/gap_test_all.pyx":155
+    /* "src/screening/gap_test_all.pyx":156
  * 
  *       # 3. Tests
  *       for l in range(n, 0, -1):             # <<<<<<<<<<<<<<
@@ -3165,22 +3153,22 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
   }
   __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
 
-  /* "src/screening/gap_test_all.pyx":173
+  /* "src/screening/gap_test_all.pyx":174
  * 
  * 
  *       calI_screen[index[l:]] = True             # <<<<<<<<<<<<<<
  *       return calI_screen
  */
-  if (unlikely(!__pyx_v_l)) { __Pyx_RaiseUnboundLocalError("l"); __PYX_ERR(0, 173, __pyx_L1_error) }
-  __pyx_t_4 = PySlice_New(__pyx_v_l, Py_None, Py_None); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 173, __pyx_L1_error)
+  if (unlikely(!__pyx_v_l)) { __Pyx_RaiseUnboundLocalError("l"); __PYX_ERR(0, 174, __pyx_L1_error) }
+  __pyx_t_4 = PySlice_New(__pyx_v_l, Py_None, Py_None); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 174, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_4);
-  __pyx_t_5 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_index), __pyx_t_4); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 173, __pyx_L1_error)
+  __pyx_t_5 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_index), __pyx_t_4); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 174, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_5);
   __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
-  if (unlikely(PyObject_SetItem(((PyObject *)__pyx_v_calI_screen), __pyx_t_5, Py_True) < 0)) __PYX_ERR(0, 173, __pyx_L1_error)
+  if (unlikely(PyObject_SetItem(((PyObject *)__pyx_v_calI_screen), __pyx_t_5, Py_True) < 0)) __PYX_ERR(0, 174, __pyx_L1_error)
   __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
 
-  /* "src/screening/gap_test_all.pyx":174
+  /* "src/screening/gap_test_all.pyx":175
  * 
  *       calI_screen[index[l:]] = True
  *       return calI_screen             # <<<<<<<<<<<<<<
@@ -3190,7 +3178,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_all_10GapTestAll_4apply_tes
   __pyx_r = ((PyObject *)__pyx_v_calI_screen);
   goto __pyx_L0;
 
-  /* "src/screening/gap_test_all.pyx":57
+  /* "src/screening/gap_test_all.pyx":58
  *    @cython.boundscheck(False) # turn off bounds-checking for entire function
  *    # @cython.wraparound(False)  # turn off negative index wrapping for entire function
  *    def apply_test(self,             # <<<<<<<<<<<<<<
@@ -4526,7 +4514,6 @@ static struct PyModuleDef __pyx_moduledef = {
 static __Pyx_StringTabEntry __pyx_string_tab[] = {
   {&__pyx_n_s_AbstractGapScreening, __pyx_k_AbstractGapScreening, sizeof(__pyx_k_AbstractGapScreening), 0, 0, 1, 1},
   {&__pyx_n_s_Atcabs, __pyx_k_Atcabs, sizeof(__pyx_k_Atcabs), 0, 0, 1, 1},
-  {&__pyx_kp_s_GAP_test_all_cython, __pyx_k_GAP_test_all_cython, sizeof(__pyx_k_GAP_test_all_cython), 0, 0, 1, 0},
   {&__pyx_n_s_GapTestAll, __pyx_k_GapTestAll, sizeof(__pyx_k_GapTestAll), 0, 0, 1, 1},
   {&__pyx_n_s_ImportError, __pyx_k_ImportError, sizeof(__pyx_k_ImportError), 0, 0, 1, 1},
   {&__pyx_n_s_TypeError, __pyx_k_TypeError, sizeof(__pyx_k_TypeError), 0, 0, 1, 1},
@@ -4543,7 +4530,6 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = {
   {&__pyx_n_s_lbd, __pyx_k_lbd, sizeof(__pyx_k_lbd), 0, 0, 1, 1},
   {&__pyx_n_s_main, __pyx_k_main, sizeof(__pyx_k_main), 0, 0, 1, 1},
   {&__pyx_n_s_name, __pyx_k_name, sizeof(__pyx_k_name), 0, 0, 1, 1},
-  {&__pyx_n_s_name_2, __pyx_k_name_2, sizeof(__pyx_k_name_2), 0, 0, 1, 1},
   {&__pyx_n_s_np, __pyx_k_np, sizeof(__pyx_k_np), 0, 0, 1, 1},
   {&__pyx_n_s_numpy, __pyx_k_numpy, sizeof(__pyx_k_numpy), 0, 0, 1, 1},
   {&__pyx_kp_s_numpy_core_multiarray_failed_to, __pyx_k_numpy_core_multiarray_failed_to, sizeof(__pyx_k_numpy_core_multiarray_failed_to), 0, 0, 1, 0},
@@ -4564,7 +4550,7 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = {
   {0, 0, 0, 0, 0, 0, 0}
 };
 static CYTHON_SMALL_CODE int __Pyx_InitCachedBuiltins(void) {
-  __pyx_builtin_range = __Pyx_GetBuiltinName(__pyx_n_s_range); if (!__pyx_builtin_range) __PYX_ERR(0, 39, __pyx_L1_error)
+  __pyx_builtin_range = __Pyx_GetBuiltinName(__pyx_n_s_range); if (!__pyx_builtin_range) __PYX_ERR(0, 40, __pyx_L1_error)
   __pyx_builtin_TypeError = __Pyx_GetBuiltinName(__pyx_n_s_TypeError); if (!__pyx_builtin_TypeError) __PYX_ERR(1, 2, __pyx_L1_error)
   __pyx_builtin_ImportError = __Pyx_GetBuiltinName(__pyx_n_s_ImportError); if (!__pyx_builtin_ImportError) __PYX_ERR(2, 945, __pyx_L1_error)
   return 0;
@@ -4576,14 +4562,14 @@ static CYTHON_SMALL_CODE int __Pyx_InitCachedConstants(void) {
   __Pyx_RefNannyDeclarations
   __Pyx_RefNannySetupContext("__Pyx_InitCachedConstants", 0);
 
-  /* "src/screening/gap_test_all.pyx":124
+  /* "src/screening/gap_test_all.pyx":125
  *       # 1. Sort in descending order
  *       if index is None:
  *          index = np.argsort(Atcabs)[::-1]             # <<<<<<<<<<<<<<
  * 
  *       # 2. Precomputing quantities
  */
-  __pyx_slice_ = PySlice_New(Py_None, Py_None, __pyx_int_neg_1); if (unlikely(!__pyx_slice_)) __PYX_ERR(0, 124, __pyx_L1_error)
+  __pyx_slice_ = PySlice_New(Py_None, Py_None, __pyx_int_neg_1); if (unlikely(!__pyx_slice_)) __PYX_ERR(0, 125, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_slice_);
   __Pyx_GIVEREF(__pyx_slice_);
 
@@ -4950,7 +4936,7 @@ if (!__Pyx_RefNanny) {
   if (__Pyx_init_sys_getdefaultencoding_params() < 0) __PYX_ERR(0, 1, __pyx_L1_error)
   #endif
   if (__pyx_module_is_main_src__screening__gap_test_all) {
-    if (PyObject_SetAttr(__pyx_m, __pyx_n_s_name_2, __pyx_n_s_main) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
+    if (PyObject_SetAttr(__pyx_m, __pyx_n_s_name, __pyx_n_s_main) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
   }
   #if PY_MAJOR_VERSION >= 3
   {
@@ -5822,20 +5808,6 @@ fail:;
   return -1;
 }
 
-/* PyObjectSetAttrStr */
-  #if CYTHON_USE_TYPE_SLOTS
-static CYTHON_INLINE int __Pyx_PyObject_SetAttrStr(PyObject* obj, PyObject* attr_name, PyObject* value) {
-    PyTypeObject* tp = Py_TYPE(obj);
-    if (likely(tp->tp_setattro))
-        return tp->tp_setattro(obj, attr_name, value);
-#if PY_MAJOR_VERSION < 3
-    if (likely(tp->tp_setattr))
-        return tp->tp_setattr(obj, PyString_AS_STRING(attr_name), value);
-#endif
-    return PyObject_SetAttr(obj, attr_name, value);
-}
-#endif
-
 /* PyErrFetchRestore */
   #if CYTHON_FAST_THREAD_STATE
 static CYTHON_INLINE void __Pyx_ErrRestoreInState(PyThreadState *tstate, PyObject *type, PyObject *value, PyObject *tb) {
@@ -6826,7 +6798,7 @@ static CYTHON_INLINE PyObject* __Pyx_PyObject_GetAttrStrNoError(PyObject* obj, P
   static int __Pyx_setup_reduce_is_named(PyObject* meth, PyObject* name) {
   int ret;
   PyObject *name_attr;
-  name_attr = __Pyx_PyObject_GetAttrStr(meth, __pyx_n_s_name_2);
+  name_attr = __Pyx_PyObject_GetAttrStr(meth, __pyx_n_s_name);
   if (likely(name_attr)) {
       ret = PyObject_RichCompareBool(name_attr, name, Py_EQ);
   } else {
diff --git a/src/screening/gap_test_all.pyx b/src/screening/gap_test_all.pyx
index 4bd697c7d95a6381e1e662a89ab1d818c566277b..b172470ee2e0c9bb2de347a322e3d0494c92e3ac 100644
--- a/src/screening/gap_test_all.pyx
+++ b/src/screening/gap_test_all.pyx
@@ -23,6 +23,7 @@ cdef class GapTestAll:
    """
 
    cdef double *vec_cumsum_gammas
+   cdef char* name
 
    @cython.boundscheck(False) # turn off bounds-checking for entire function
    @cython.wraparound(False)  # turn off negative index wrapping for entire function
diff --git a/src/screening/gap_test_p_1.c b/src/screening/gap_test_p_1.c
index ade654c8ab195b25ba9c2e7816c72fbd84394692..595c7697e4685cd32e7c1986eff7d7d007a3e885 100644
--- a/src/screening/gap_test_p_1.c
+++ b/src/screening/gap_test_p_1.c
@@ -1182,6 +1182,7 @@ struct __pyx_obj_3src_9screening_12gap_test_p_1_GapTestPequalOne {
   PyObject_HEAD
   double *vec_cumsum_gammas;
   double *vec_kappa_q;
+  char *name;
 };
 
 
@@ -1394,15 +1395,6 @@ static CYTHON_INLINE PyObject* __Pyx_PyObject_CallMethO(PyObject *func, PyObject
 /* PyObjectCallOneArg.proto */
 static CYTHON_INLINE PyObject* __Pyx_PyObject_CallOneArg(PyObject *func, PyObject *arg);
 
-/* PyObjectSetAttrStr.proto */
-#if CYTHON_USE_TYPE_SLOTS
-#define __Pyx_PyObject_DelAttrStr(o,n) __Pyx_PyObject_SetAttrStr(o, n, NULL)
-static CYTHON_INLINE int __Pyx_PyObject_SetAttrStr(PyObject* obj, PyObject* attr_name, PyObject* value);
-#else
-#define __Pyx_PyObject_DelAttrStr(o,n)   PyObject_DelAttr(o,n)
-#define __Pyx_PyObject_SetAttrStr(o,n,v) PyObject_SetAttr(o,n,v)
-#endif
-
 /* PyThreadStateGet.proto */
 #if CYTHON_FAST_THREAD_STATE
 #define __Pyx_PyThreadState_declare  PyThreadState *__pyx_tstate;
@@ -1787,7 +1779,7 @@ static const char __pyx_k_np[] = "np";
 static const char __pyx_k_gap[] = "gap";
 static const char __pyx_k_lbd[] = "lbd";
 static const char __pyx_k_main[] = "__main__";
-static const char __pyx_k_name[] = "name";
+static const char __pyx_k_name[] = "__name__";
 static const char __pyx_k_sqrt[] = "sqrt";
 static const char __pyx_k_test[] = "__test__";
 static const char __pyx_k_dtype[] = "dtype";
@@ -1799,7 +1791,6 @@ static const char __pyx_k_Atcabs[] = "Atcabs";
 static const char __pyx_k_argmax[] = "argmax";
 static const char __pyx_k_cumsum[] = "cumsum";
 static const char __pyx_k_import[] = "__import__";
-static const char __pyx_k_name_2[] = "__name__";
 static const char __pyx_k_reduce[] = "__reduce__";
 static const char __pyx_k_argsort[] = "argsort";
 static const char __pyx_k_getstate[] = "__getstate__";
@@ -1814,7 +1805,6 @@ static const char __pyx_k_setstate_cython[] = "__setstate_cython__";
 static const char __pyx_k_GapTestPequalOne[] = "GapTestPequalOne";
 static const char __pyx_k_cline_in_traceback[] = "cline_in_traceback";
 static const char __pyx_k_coeff_dual_scaling[] = "coeff_dual_scaling";
-static const char __pyx_k_Gap_test_p_1_cython[] = "Gap test p=1 (cython)";
 static const char __pyx_k_AbstractGapScreening[] = "AbstractGapScreening";
 static const char __pyx_k_vec_coherence_function[] = "vec_coherence_function";
 static const char __pyx_k_src_screening_ascreening[] = "src.screening.ascreening";
@@ -1824,7 +1814,6 @@ static const char __pyx_k_numpy_core_umath_failed_to_impor[] = "numpy.core.umath
 static PyObject *__pyx_n_s_AbstractGapScreening;
 static PyObject *__pyx_n_s_Atcabs;
 static PyObject *__pyx_n_s_GapTestPequalOne;
-static PyObject *__pyx_kp_s_Gap_test_p_1_cython;
 static PyObject *__pyx_n_s_ImportError;
 static PyObject *__pyx_n_s_TypeError;
 static PyObject *__pyx_n_s_argmax;
@@ -1840,7 +1829,6 @@ static PyObject *__pyx_n_s_index;
 static PyObject *__pyx_n_s_lbd;
 static PyObject *__pyx_n_s_main;
 static PyObject *__pyx_n_s_name;
-static PyObject *__pyx_n_s_name_2;
 static PyObject *__pyx_n_s_np;
 static PyObject *__pyx_n_s_numpy;
 static PyObject *__pyx_kp_s_numpy_core_multiarray_failed_to;
@@ -1875,7 +1863,7 @@ static PyObject *__pyx_tuple__5;
 static PyObject *__pyx_tuple__6;
 /* Late includes */
 
-/* "src/screening/gap_test_p_1.pyx":38
+/* "src/screening/gap_test_p_1.pyx":39
  *    @cython.boundscheck(False) # turn off bounds-checking for entire function
  *    @cython.wraparound(False)  # turn off negative index wrapping for entire function
  *    def __init__(             # <<<<<<<<<<<<<<
@@ -1917,11 +1905,11 @@ static int __pyx_pw_3src_9screening_12gap_test_p_1_16GapTestPequalOne_1__init__(
         case  1:
         if (likely((values[1] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_vec_coherence_function)) != 0)) kw_args--;
         else {
-          __Pyx_RaiseArgtupleInvalid("__init__", 1, 2, 2, 1); __PYX_ERR(0, 38, __pyx_L3_error)
+          __Pyx_RaiseArgtupleInvalid("__init__", 1, 2, 2, 1); __PYX_ERR(0, 39, __pyx_L3_error)
         }
       }
       if (unlikely(kw_args > 0)) {
-        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "__init__") < 0)) __PYX_ERR(0, 38, __pyx_L3_error)
+        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "__init__") < 0)) __PYX_ERR(0, 39, __pyx_L3_error)
       }
     } else if (PyTuple_GET_SIZE(__pyx_args) != 2) {
       goto __pyx_L5_argtuple_error;
@@ -1934,14 +1922,14 @@ static int __pyx_pw_3src_9screening_12gap_test_p_1_16GapTestPequalOne_1__init__(
   }
   goto __pyx_L4_argument_unpacking_done;
   __pyx_L5_argtuple_error:;
-  __Pyx_RaiseArgtupleInvalid("__init__", 1, 2, 2, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 38, __pyx_L3_error)
+  __Pyx_RaiseArgtupleInvalid("__init__", 1, 2, 2, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 39, __pyx_L3_error)
   __pyx_L3_error:;
   __Pyx_AddTraceback("src.screening.gap_test_p_1.GapTestPequalOne.__init__", __pyx_clineno, __pyx_lineno, __pyx_filename);
   __Pyx_RefNannyFinishContext();
   return -1;
   __pyx_L4_argument_unpacking_done:;
-  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_vec_gammas), __pyx_ptype_5numpy_ndarray, 1, "vec_gammas", 0))) __PYX_ERR(0, 40, __pyx_L1_error)
-  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_vec_coherence_function), __pyx_ptype_5numpy_ndarray, 1, "vec_coherence_function", 0))) __PYX_ERR(0, 41, __pyx_L1_error)
+  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_vec_gammas), __pyx_ptype_5numpy_ndarray, 1, "vec_gammas", 0))) __PYX_ERR(0, 41, __pyx_L1_error)
+  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_vec_coherence_function), __pyx_ptype_5numpy_ndarray, 1, "vec_coherence_function", 0))) __PYX_ERR(0, 42, __pyx_L1_error)
   __pyx_r = __pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne___init__(((struct __pyx_obj_3src_9screening_12gap_test_p_1_GapTestPequalOne *)__pyx_v_self), __pyx_v_vec_gammas, __pyx_v_vec_coherence_function);
 
   /* function exit code */
@@ -1985,16 +1973,16 @@ static int __pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne___init__(s
   __pyx_pybuffernd_vec_coherence_function.rcbuffer = &__pyx_pybuffer_vec_coherence_function;
   {
     __Pyx_BufFmt_StackElem __pyx_stack[1];
-    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_vec_gammas.rcbuffer->pybuffer, (PyObject*)__pyx_v_vec_gammas, &__Pyx_TypeInfo_nn_npy_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 38, __pyx_L1_error)
+    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_vec_gammas.rcbuffer->pybuffer, (PyObject*)__pyx_v_vec_gammas, &__Pyx_TypeInfo_nn_npy_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 39, __pyx_L1_error)
   }
   __pyx_pybuffernd_vec_gammas.diminfo[0].strides = __pyx_pybuffernd_vec_gammas.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_vec_gammas.diminfo[0].shape = __pyx_pybuffernd_vec_gammas.rcbuffer->pybuffer.shape[0];
   {
     __Pyx_BufFmt_StackElem __pyx_stack[1];
-    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_vec_coherence_function.rcbuffer->pybuffer, (PyObject*)__pyx_v_vec_coherence_function, &__Pyx_TypeInfo_nn_npy_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 38, __pyx_L1_error)
+    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_vec_coherence_function.rcbuffer->pybuffer, (PyObject*)__pyx_v_vec_coherence_function, &__Pyx_TypeInfo_nn_npy_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 39, __pyx_L1_error)
   }
   __pyx_pybuffernd_vec_coherence_function.diminfo[0].strides = __pyx_pybuffernd_vec_coherence_function.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_vec_coherence_function.diminfo[0].shape = __pyx_pybuffernd_vec_coherence_function.rcbuffer->pybuffer.shape[0];
 
-  /* "src/screening/gap_test_p_1.pyx":44
+  /* "src/screening/gap_test_p_1.pyx":45
  *       ):
  * 
  *       cdef int n = vec_gammas.shape[0]             # <<<<<<<<<<<<<<
@@ -2003,7 +1991,7 @@ static int __pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne___init__(s
  */
   __pyx_v_n = (__pyx_v_vec_gammas->dimensions[0]);
 
-  /* "src/screening/gap_test_p_1.pyx":56
+  /* "src/screening/gap_test_p_1.pyx":57
  *       # self.vec_kappa_q = &buf2[0]
  * 
  *       self.vec_cumsum_gammas = <double*> PyMem_Malloc(n * sizeof(double))             # <<<<<<<<<<<<<<
@@ -2012,7 +2000,7 @@ static int __pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne___init__(s
  */
   __pyx_v_self->vec_cumsum_gammas = ((double *)PyMem_Malloc((__pyx_v_n * (sizeof(double)))));
 
-  /* "src/screening/gap_test_p_1.pyx":57
+  /* "src/screening/gap_test_p_1.pyx":58
  * 
  *       self.vec_cumsum_gammas = <double*> PyMem_Malloc(n * sizeof(double))
  *       self.vec_cumsum_gammas[0] = vec_gammas[0]             # <<<<<<<<<<<<<<
@@ -2022,7 +2010,7 @@ static int __pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne___init__(s
   __pyx_t_1 = 0;
   (__pyx_v_self->vec_cumsum_gammas[0]) = (*__Pyx_BufPtrStrided1d(npy_double *, __pyx_pybuffernd_vec_gammas.rcbuffer->pybuffer.buf, __pyx_t_1, __pyx_pybuffernd_vec_gammas.diminfo[0].strides));
 
-  /* "src/screening/gap_test_p_1.pyx":59
+  /* "src/screening/gap_test_p_1.pyx":60
  *       self.vec_cumsum_gammas[0] = vec_gammas[0]
  * 
  *       self.vec_kappa_q = <double*> PyMem_Malloc(n * sizeof(double))             # <<<<<<<<<<<<<<
@@ -2031,7 +2019,7 @@ static int __pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne___init__(s
  */
   __pyx_v_self->vec_kappa_q = ((double *)PyMem_Malloc((__pyx_v_n * (sizeof(double)))));
 
-  /* "src/screening/gap_test_p_1.pyx":61
+  /* "src/screening/gap_test_p_1.pyx":62
  *       self.vec_kappa_q = <double*> PyMem_Malloc(n * sizeof(double))
  * 
  *       cdef int i = 0             # <<<<<<<<<<<<<<
@@ -2040,7 +2028,7 @@ static int __pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne___init__(s
  */
   __pyx_v_i = 0;
 
-  /* "src/screening/gap_test_p_1.pyx":62
+  /* "src/screening/gap_test_p_1.pyx":63
  * 
  *       cdef int i = 0
  *       for i in range(1, n):             # <<<<<<<<<<<<<<
@@ -2052,7 +2040,7 @@ static int __pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne___init__(s
   for (__pyx_t_4 = 1; __pyx_t_4 < __pyx_t_3; __pyx_t_4+=1) {
     __pyx_v_i = __pyx_t_4;
 
-    /* "src/screening/gap_test_p_1.pyx":63
+    /* "src/screening/gap_test_p_1.pyx":64
  *       cdef int i = 0
  *       for i in range(1, n):
  *          self.vec_cumsum_gammas[i] = self.vec_cumsum_gammas[i-1] + vec_gammas[i]             # <<<<<<<<<<<<<<
@@ -2063,7 +2051,7 @@ static int __pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne___init__(s
     (__pyx_v_self->vec_cumsum_gammas[__pyx_v_i]) = ((__pyx_v_self->vec_cumsum_gammas[(__pyx_v_i - 1)]) + (*__Pyx_BufPtrStrided1d(npy_double *, __pyx_pybuffernd_vec_gammas.rcbuffer->pybuffer.buf, __pyx_t_1, __pyx_pybuffernd_vec_gammas.diminfo[0].strides)));
   }
 
-  /* "src/screening/gap_test_p_1.pyx":65
+  /* "src/screening/gap_test_p_1.pyx":66
  *          self.vec_cumsum_gammas[i] = self.vec_cumsum_gammas[i-1] + vec_gammas[i]
  * 
  *       for i in range(n):             # <<<<<<<<<<<<<<
@@ -2075,20 +2063,20 @@ static int __pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne___init__(s
   for (__pyx_t_4 = 0; __pyx_t_4 < __pyx_t_3; __pyx_t_4+=1) {
     __pyx_v_i = __pyx_t_4;
 
-    /* "src/screening/gap_test_p_1.pyx":66
+    /* "src/screening/gap_test_p_1.pyx":67
  * 
  *       for i in range(n):
  *          self.vec_kappa_q[i] = np.sqrt(             # <<<<<<<<<<<<<<
  *          (1. + vec_coherence_function[i]) * (i + 1.)
  *       )
  */
-    __Pyx_GetModuleGlobalName(__pyx_t_6, __pyx_n_s_np); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 66, __pyx_L1_error)
+    __Pyx_GetModuleGlobalName(__pyx_t_6, __pyx_n_s_np); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 67, __pyx_L1_error)
     __Pyx_GOTREF(__pyx_t_6);
-    __pyx_t_7 = __Pyx_PyObject_GetAttrStr(__pyx_t_6, __pyx_n_s_sqrt); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 66, __pyx_L1_error)
+    __pyx_t_7 = __Pyx_PyObject_GetAttrStr(__pyx_t_6, __pyx_n_s_sqrt); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 67, __pyx_L1_error)
     __Pyx_GOTREF(__pyx_t_7);
     __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
 
-    /* "src/screening/gap_test_p_1.pyx":67
+    /* "src/screening/gap_test_p_1.pyx":68
  *       for i in range(n):
  *          self.vec_kappa_q[i] = np.sqrt(
  *          (1. + vec_coherence_function[i]) * (i + 1.)             # <<<<<<<<<<<<<<
@@ -2096,7 +2084,7 @@ static int __pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne___init__(s
  * 
  */
     __pyx_t_1 = __pyx_v_i;
-    __pyx_t_6 = PyFloat_FromDouble(((1. + (*__Pyx_BufPtrStrided1d(npy_double *, __pyx_pybuffernd_vec_coherence_function.rcbuffer->pybuffer.buf, __pyx_t_1, __pyx_pybuffernd_vec_coherence_function.diminfo[0].strides))) * (__pyx_v_i + 1.))); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 67, __pyx_L1_error)
+    __pyx_t_6 = PyFloat_FromDouble(((1. + (*__Pyx_BufPtrStrided1d(npy_double *, __pyx_pybuffernd_vec_coherence_function.rcbuffer->pybuffer.buf, __pyx_t_1, __pyx_pybuffernd_vec_coherence_function.diminfo[0].strides))) * (__pyx_v_i + 1.))); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 68, __pyx_L1_error)
     __Pyx_GOTREF(__pyx_t_6);
     __pyx_t_8 = NULL;
     if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_7))) {
@@ -2111,32 +2099,32 @@ static int __pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne___init__(s
     __pyx_t_5 = (__pyx_t_8) ? __Pyx_PyObject_Call2Args(__pyx_t_7, __pyx_t_8, __pyx_t_6) : __Pyx_PyObject_CallOneArg(__pyx_t_7, __pyx_t_6);
     __Pyx_XDECREF(__pyx_t_8); __pyx_t_8 = 0;
     __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
-    if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 66, __pyx_L1_error)
+    if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 67, __pyx_L1_error)
     __Pyx_GOTREF(__pyx_t_5);
     __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
 
-    /* "src/screening/gap_test_p_1.pyx":66
+    /* "src/screening/gap_test_p_1.pyx":67
  * 
  *       for i in range(n):
  *          self.vec_kappa_q[i] = np.sqrt(             # <<<<<<<<<<<<<<
  *          (1. + vec_coherence_function[i]) * (i + 1.)
  *       )
  */
-    __pyx_t_9 = __pyx_PyFloat_AsDouble(__pyx_t_5); if (unlikely((__pyx_t_9 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 66, __pyx_L1_error)
+    __pyx_t_9 = __pyx_PyFloat_AsDouble(__pyx_t_5); if (unlikely((__pyx_t_9 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 67, __pyx_L1_error)
     __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
     (__pyx_v_self->vec_kappa_q[__pyx_v_i]) = __pyx_t_9;
   }
 
-  /* "src/screening/gap_test_p_1.pyx":75
+  /* "src/screening/gap_test_p_1.pyx":76
  *       # ).copy()
  * 
  *       self.name = "Gap test p=1 (cython)"             # <<<<<<<<<<<<<<
  * 
  * 
  */
-  if (__Pyx_PyObject_SetAttrStr(((PyObject *)__pyx_v_self), __pyx_n_s_name, __pyx_kp_s_Gap_test_p_1_cython) < 0) __PYX_ERR(0, 75, __pyx_L1_error)
+  __pyx_v_self->name = ((char *)"Gap test p=1 (cython)");
 
-  /* "src/screening/gap_test_p_1.pyx":38
+  /* "src/screening/gap_test_p_1.pyx":39
  *    @cython.boundscheck(False) # turn off bounds-checking for entire function
  *    @cython.wraparound(False)  # turn off negative index wrapping for entire function
  *    def __init__(             # <<<<<<<<<<<<<<
@@ -2170,7 +2158,7 @@ static int __pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne___init__(s
   return __pyx_r;
 }
 
-/* "src/screening/gap_test_p_1.pyx":78
+/* "src/screening/gap_test_p_1.pyx":79
  * 
  * 
  *    def __dealloc__(self):             # <<<<<<<<<<<<<<
@@ -2193,7 +2181,7 @@ static void __pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_2__deallo
   __Pyx_RefNannyDeclarations
   __Pyx_RefNannySetupContext("__dealloc__", 0);
 
-  /* "src/screening/gap_test_p_1.pyx":79
+  /* "src/screening/gap_test_p_1.pyx":80
  * 
  *    def __dealloc__(self):
  *       PyMem_Free(self.vec_cumsum_gammas)  # no-op if self.data is NULL             # <<<<<<<<<<<<<<
@@ -2202,7 +2190,7 @@ static void __pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_2__deallo
  */
   PyMem_Free(__pyx_v_self->vec_cumsum_gammas);
 
-  /* "src/screening/gap_test_p_1.pyx":80
+  /* "src/screening/gap_test_p_1.pyx":81
  *    def __dealloc__(self):
  *       PyMem_Free(self.vec_cumsum_gammas)  # no-op if self.data is NULL
  *       PyMem_Free(self.vec_kappa_q)  # no-op if self.data is NULL             # <<<<<<<<<<<<<<
@@ -2211,7 +2199,7 @@ static void __pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_2__deallo
  */
   PyMem_Free(__pyx_v_self->vec_kappa_q);
 
-  /* "src/screening/gap_test_p_1.pyx":78
+  /* "src/screening/gap_test_p_1.pyx":79
  * 
  * 
  *    def __dealloc__(self):             # <<<<<<<<<<<<<<
@@ -2223,7 +2211,7 @@ static void __pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_2__deallo
   __Pyx_RefNannyFinishContext();
 }
 
-/* "src/screening/gap_test_p_1.pyx":85
+/* "src/screening/gap_test_p_1.pyx":86
  *    @cython.boundscheck(False) # turn off bounds-checking for entire function
  *    @cython.wraparound(False)  # turn off negative index wrapping for entire function
  *    def apply_test(self,             # <<<<<<<<<<<<<<
@@ -2252,7 +2240,7 @@ static PyObject *__pyx_pw_3src_9screening_12gap_test_p_1_16GapTestPequalOne_5app
     static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_Atcabs,&__pyx_n_s_gap,&__pyx_n_s_lbd,&__pyx_n_s_vec_gammas,&__pyx_n_s_coeff_dual_scaling,&__pyx_n_s_offset_radius,&__pyx_n_s_index,0};
     PyObject* values[7] = {0,0,0,0,0,0,0};
 
-    /* "src/screening/gap_test_p_1.pyx":92
+    /* "src/screening/gap_test_p_1.pyx":93
  *       double coeff_dual_scaling = 1.,
  *       double offset_radius=0.,
  *       np.ndarray[long, ndim=1] index=None             # <<<<<<<<<<<<<<
@@ -2290,19 +2278,19 @@ static PyObject *__pyx_pw_3src_9screening_12gap_test_p_1_16GapTestPequalOne_5app
         case  1:
         if (likely((values[1] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_gap)) != 0)) kw_args--;
         else {
-          __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, 1); __PYX_ERR(0, 85, __pyx_L3_error)
+          __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, 1); __PYX_ERR(0, 86, __pyx_L3_error)
         }
         CYTHON_FALLTHROUGH;
         case  2:
         if (likely((values[2] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_lbd)) != 0)) kw_args--;
         else {
-          __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, 2); __PYX_ERR(0, 85, __pyx_L3_error)
+          __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, 2); __PYX_ERR(0, 86, __pyx_L3_error)
         }
         CYTHON_FALLTHROUGH;
         case  3:
         if (likely((values[3] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_vec_gammas)) != 0)) kw_args--;
         else {
-          __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, 3); __PYX_ERR(0, 85, __pyx_L3_error)
+          __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, 3); __PYX_ERR(0, 86, __pyx_L3_error)
         }
         CYTHON_FALLTHROUGH;
         case  4:
@@ -2324,7 +2312,7 @@ static PyObject *__pyx_pw_3src_9screening_12gap_test_p_1_16GapTestPequalOne_5app
         }
       }
       if (unlikely(kw_args > 0)) {
-        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "apply_test") < 0)) __PYX_ERR(0, 85, __pyx_L3_error)
+        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "apply_test") < 0)) __PYX_ERR(0, 86, __pyx_L3_error)
       }
     } else {
       switch (PyTuple_GET_SIZE(__pyx_args)) {
@@ -2343,16 +2331,16 @@ static PyObject *__pyx_pw_3src_9screening_12gap_test_p_1_16GapTestPequalOne_5app
       }
     }
     __pyx_v_Atcabs = ((PyArrayObject *)values[0]);
-    __pyx_v_gap = __pyx_PyFloat_AsDouble(values[1]); if (unlikely((__pyx_v_gap == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 87, __pyx_L3_error)
-    __pyx_v_lbd = __pyx_PyFloat_AsDouble(values[2]); if (unlikely((__pyx_v_lbd == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 88, __pyx_L3_error)
+    __pyx_v_gap = __pyx_PyFloat_AsDouble(values[1]); if (unlikely((__pyx_v_gap == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 88, __pyx_L3_error)
+    __pyx_v_lbd = __pyx_PyFloat_AsDouble(values[2]); if (unlikely((__pyx_v_lbd == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 89, __pyx_L3_error)
     __pyx_v_vec_gammas = ((PyArrayObject *)values[3]);
     if (values[4]) {
-      __pyx_v_coeff_dual_scaling = __pyx_PyFloat_AsDouble(values[4]); if (unlikely((__pyx_v_coeff_dual_scaling == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 90, __pyx_L3_error)
+      __pyx_v_coeff_dual_scaling = __pyx_PyFloat_AsDouble(values[4]); if (unlikely((__pyx_v_coeff_dual_scaling == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 91, __pyx_L3_error)
     } else {
       __pyx_v_coeff_dual_scaling = ((double)1.);
     }
     if (values[5]) {
-      __pyx_v_offset_radius = __pyx_PyFloat_AsDouble(values[5]); if (unlikely((__pyx_v_offset_radius == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 91, __pyx_L3_error)
+      __pyx_v_offset_radius = __pyx_PyFloat_AsDouble(values[5]); if (unlikely((__pyx_v_offset_radius == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 92, __pyx_L3_error)
     } else {
       __pyx_v_offset_radius = ((double)0.);
     }
@@ -2360,18 +2348,18 @@ static PyObject *__pyx_pw_3src_9screening_12gap_test_p_1_16GapTestPequalOne_5app
   }
   goto __pyx_L4_argument_unpacking_done;
   __pyx_L5_argtuple_error:;
-  __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 85, __pyx_L3_error)
+  __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 86, __pyx_L3_error)
   __pyx_L3_error:;
   __Pyx_AddTraceback("src.screening.gap_test_p_1.GapTestPequalOne.apply_test", __pyx_clineno, __pyx_lineno, __pyx_filename);
   __Pyx_RefNannyFinishContext();
   return NULL;
   __pyx_L4_argument_unpacking_done:;
-  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_Atcabs), __pyx_ptype_5numpy_ndarray, 1, "Atcabs", 0))) __PYX_ERR(0, 86, __pyx_L1_error)
-  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_vec_gammas), __pyx_ptype_5numpy_ndarray, 1, "vec_gammas", 0))) __PYX_ERR(0, 89, __pyx_L1_error)
-  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_index), __pyx_ptype_5numpy_ndarray, 1, "index", 0))) __PYX_ERR(0, 92, __pyx_L1_error)
+  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_Atcabs), __pyx_ptype_5numpy_ndarray, 1, "Atcabs", 0))) __PYX_ERR(0, 87, __pyx_L1_error)
+  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_vec_gammas), __pyx_ptype_5numpy_ndarray, 1, "vec_gammas", 0))) __PYX_ERR(0, 90, __pyx_L1_error)
+  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_index), __pyx_ptype_5numpy_ndarray, 1, "index", 0))) __PYX_ERR(0, 93, __pyx_L1_error)
   __pyx_r = __pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4apply_test(((struct __pyx_obj_3src_9screening_12gap_test_p_1_GapTestPequalOne *)__pyx_v_self), __pyx_v_Atcabs, __pyx_v_gap, __pyx_v_lbd, __pyx_v_vec_gammas, __pyx_v_coeff_dual_scaling, __pyx_v_offset_radius, __pyx_v_index);
 
-  /* "src/screening/gap_test_p_1.pyx":85
+  /* "src/screening/gap_test_p_1.pyx":86
  *    @cython.boundscheck(False) # turn off bounds-checking for entire function
  *    @cython.wraparound(False)  # turn off negative index wrapping for entire function
  *    def apply_test(self,             # <<<<<<<<<<<<<<
@@ -2458,21 +2446,21 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
   __pyx_pybuffernd_index.rcbuffer = &__pyx_pybuffer_index;
   {
     __Pyx_BufFmt_StackElem __pyx_stack[1];
-    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_Atcabs.rcbuffer->pybuffer, (PyObject*)__pyx_v_Atcabs, &__Pyx_TypeInfo_nn_npy_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 85, __pyx_L1_error)
+    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_Atcabs.rcbuffer->pybuffer, (PyObject*)__pyx_v_Atcabs, &__Pyx_TypeInfo_nn_npy_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 86, __pyx_L1_error)
   }
   __pyx_pybuffernd_Atcabs.diminfo[0].strides = __pyx_pybuffernd_Atcabs.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_Atcabs.diminfo[0].shape = __pyx_pybuffernd_Atcabs.rcbuffer->pybuffer.shape[0];
   {
     __Pyx_BufFmt_StackElem __pyx_stack[1];
-    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_vec_gammas.rcbuffer->pybuffer, (PyObject*)__pyx_v_vec_gammas, &__Pyx_TypeInfo_nn_npy_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 85, __pyx_L1_error)
+    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_vec_gammas.rcbuffer->pybuffer, (PyObject*)__pyx_v_vec_gammas, &__Pyx_TypeInfo_nn_npy_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 86, __pyx_L1_error)
   }
   __pyx_pybuffernd_vec_gammas.diminfo[0].strides = __pyx_pybuffernd_vec_gammas.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_vec_gammas.diminfo[0].shape = __pyx_pybuffernd_vec_gammas.rcbuffer->pybuffer.shape[0];
   {
     __Pyx_BufFmt_StackElem __pyx_stack[1];
-    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_index.rcbuffer->pybuffer, (PyObject*)__pyx_v_index, &__Pyx_TypeInfo_long, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 85, __pyx_L1_error)
+    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_index.rcbuffer->pybuffer, (PyObject*)__pyx_v_index, &__Pyx_TypeInfo_long, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 86, __pyx_L1_error)
   }
   __pyx_pybuffernd_index.diminfo[0].strides = __pyx_pybuffernd_index.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_index.diminfo[0].shape = __pyx_pybuffernd_index.rcbuffer->pybuffer.shape[0];
 
-  /* "src/screening/gap_test_p_1.pyx":133
+  /* "src/screening/gap_test_p_1.pyx":134
  *       """
  * 
  *       cdef int n = Atcabs.shape[0]             # <<<<<<<<<<<<<<
@@ -2481,7 +2469,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
  */
   __pyx_v_n = (__pyx_v_Atcabs->dimensions[0]);
 
-  /* "src/screening/gap_test_p_1.pyx":135
+  /* "src/screening/gap_test_p_1.pyx":136
  *       cdef int n = Atcabs.shape[0]
  *       cdef int q, l_q
  *       cdef int l_min = 0             # <<<<<<<<<<<<<<
@@ -2490,21 +2478,21 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
  */
   __pyx_v_l_min = 0;
 
-  /* "src/screening/gap_test_p_1.pyx":138
+  /* "src/screening/gap_test_p_1.pyx":139
  * 
  *       cdef double bound
  *       cdef double radius = coeff_dual_scaling * np.sqrt(2 * gap) + offset_radius             # <<<<<<<<<<<<<<
  *       cdef double coeff_lbd = coeff_dual_scaling * lbd
  * 
  */
-  __pyx_t_1 = PyFloat_FromDouble(__pyx_v_coeff_dual_scaling); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 138, __pyx_L1_error)
+  __pyx_t_1 = PyFloat_FromDouble(__pyx_v_coeff_dual_scaling); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 139, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
-  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 138, __pyx_L1_error)
+  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 139, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_3);
-  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_sqrt); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 138, __pyx_L1_error)
+  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_sqrt); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 139, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_4);
   __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
-  __pyx_t_3 = PyFloat_FromDouble((2.0 * __pyx_v_gap)); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 138, __pyx_L1_error)
+  __pyx_t_3 = PyFloat_FromDouble((2.0 * __pyx_v_gap)); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 139, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_3);
   __pyx_t_5 = NULL;
   if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_4))) {
@@ -2519,24 +2507,24 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
   __pyx_t_2 = (__pyx_t_5) ? __Pyx_PyObject_Call2Args(__pyx_t_4, __pyx_t_5, __pyx_t_3) : __Pyx_PyObject_CallOneArg(__pyx_t_4, __pyx_t_3);
   __Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
   __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
-  if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 138, __pyx_L1_error)
+  if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 139, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_2);
   __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
-  __pyx_t_4 = PyNumber_Multiply(__pyx_t_1, __pyx_t_2); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 138, __pyx_L1_error)
+  __pyx_t_4 = PyNumber_Multiply(__pyx_t_1, __pyx_t_2); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 139, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_4);
   __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
   __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
-  __pyx_t_2 = PyFloat_FromDouble(__pyx_v_offset_radius); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 138, __pyx_L1_error)
+  __pyx_t_2 = PyFloat_FromDouble(__pyx_v_offset_radius); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 139, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_2);
-  __pyx_t_1 = PyNumber_Add(__pyx_t_4, __pyx_t_2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 138, __pyx_L1_error)
+  __pyx_t_1 = PyNumber_Add(__pyx_t_4, __pyx_t_2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 139, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
   __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
   __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
-  __pyx_t_6 = __pyx_PyFloat_AsDouble(__pyx_t_1); if (unlikely((__pyx_t_6 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 138, __pyx_L1_error)
+  __pyx_t_6 = __pyx_PyFloat_AsDouble(__pyx_t_1); if (unlikely((__pyx_t_6 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 139, __pyx_L1_error)
   __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
   __pyx_v_radius = __pyx_t_6;
 
-  /* "src/screening/gap_test_p_1.pyx":139
+  /* "src/screening/gap_test_p_1.pyx":140
  *       cdef double bound
  *       cdef double radius = coeff_dual_scaling * np.sqrt(2 * gap) + offset_radius
  *       cdef double coeff_lbd = coeff_dual_scaling * lbd             # <<<<<<<<<<<<<<
@@ -2545,19 +2533,19 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
  */
   __pyx_v_coeff_lbd = (__pyx_v_coeff_dual_scaling * __pyx_v_lbd);
 
-  /* "src/screening/gap_test_p_1.pyx":141
+  /* "src/screening/gap_test_p_1.pyx":142
  *       cdef double coeff_lbd = coeff_dual_scaling * lbd
  * 
  *       cdef np.ndarray[np.npy_double, ndim=1] sigmas = np.zeros(n+1)             # <<<<<<<<<<<<<<
  * 
  *       # output
  */
-  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 141, __pyx_L1_error)
+  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 142, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_2);
-  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_zeros); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 141, __pyx_L1_error)
+  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_zeros); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 142, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_4);
   __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
-  __pyx_t_2 = __Pyx_PyInt_From_long((__pyx_v_n + 1)); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 141, __pyx_L1_error)
+  __pyx_t_2 = __Pyx_PyInt_From_long((__pyx_v_n + 1)); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 142, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_2);
   __pyx_t_3 = NULL;
   if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_4))) {
@@ -2572,16 +2560,16 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
   __pyx_t_1 = (__pyx_t_3) ? __Pyx_PyObject_Call2Args(__pyx_t_4, __pyx_t_3, __pyx_t_2) : __Pyx_PyObject_CallOneArg(__pyx_t_4, __pyx_t_2);
   __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
   __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
-  if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 141, __pyx_L1_error)
+  if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 142, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
   __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
-  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 141, __pyx_L1_error)
+  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 142, __pyx_L1_error)
   __pyx_t_7 = ((PyArrayObject *)__pyx_t_1);
   {
     __Pyx_BufFmt_StackElem __pyx_stack[1];
     if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_sigmas.rcbuffer->pybuffer, (PyObject*)__pyx_t_7, &__Pyx_TypeInfo_nn_npy_double, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
       __pyx_v_sigmas = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_sigmas.rcbuffer->pybuffer.buf = NULL;
-      __PYX_ERR(0, 141, __pyx_L1_error)
+      __PYX_ERR(0, 142, __pyx_L1_error)
     } else {__pyx_pybuffernd_sigmas.diminfo[0].strides = __pyx_pybuffernd_sigmas.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_sigmas.diminfo[0].shape = __pyx_pybuffernd_sigmas.rcbuffer->pybuffer.shape[0];
     }
   }
@@ -2589,40 +2577,40 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
   __pyx_v_sigmas = ((PyArrayObject *)__pyx_t_1);
   __pyx_t_1 = 0;
 
-  /* "src/screening/gap_test_p_1.pyx":144
+  /* "src/screening/gap_test_p_1.pyx":145
  * 
  *       # output
  *       cdef np.ndarray[np.npy_bool, ndim=1] calI_screen = np.zeros(n, dtype=bool)             # <<<<<<<<<<<<<<
  * 
  *       # 1. Sort in descending order
  */
-  __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 144, __pyx_L1_error)
+  __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 145, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
-  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_zeros); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 144, __pyx_L1_error)
+  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_zeros); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 145, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_4);
   __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
-  __pyx_t_1 = __Pyx_PyInt_From_int(__pyx_v_n); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 144, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyInt_From_int(__pyx_v_n); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 145, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
-  __pyx_t_2 = PyTuple_New(1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 144, __pyx_L1_error)
+  __pyx_t_2 = PyTuple_New(1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 145, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_2);
   __Pyx_GIVEREF(__pyx_t_1);
   PyTuple_SET_ITEM(__pyx_t_2, 0, __pyx_t_1);
   __pyx_t_1 = 0;
-  __pyx_t_1 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 144, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 145, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
-  if (PyDict_SetItem(__pyx_t_1, __pyx_n_s_dtype, ((PyObject*)&PyBool_Type)) < 0) __PYX_ERR(0, 144, __pyx_L1_error)
-  __pyx_t_3 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_t_2, __pyx_t_1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 144, __pyx_L1_error)
+  if (PyDict_SetItem(__pyx_t_1, __pyx_n_s_dtype, ((PyObject*)&PyBool_Type)) < 0) __PYX_ERR(0, 145, __pyx_L1_error)
+  __pyx_t_3 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_t_2, __pyx_t_1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 145, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_3);
   __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
   __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
   __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
-  if (!(likely(((__pyx_t_3) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_3, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 144, __pyx_L1_error)
+  if (!(likely(((__pyx_t_3) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_3, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 145, __pyx_L1_error)
   __pyx_t_8 = ((PyArrayObject *)__pyx_t_3);
   {
     __Pyx_BufFmt_StackElem __pyx_stack[1];
     if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_calI_screen.rcbuffer->pybuffer, (PyObject*)__pyx_t_8, &__Pyx_TypeInfo_nn_npy_bool, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
       __pyx_v_calI_screen = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_calI_screen.rcbuffer->pybuffer.buf = NULL;
-      __PYX_ERR(0, 144, __pyx_L1_error)
+      __PYX_ERR(0, 145, __pyx_L1_error)
     } else {__pyx_pybuffernd_calI_screen.diminfo[0].strides = __pyx_pybuffernd_calI_screen.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_calI_screen.diminfo[0].shape = __pyx_pybuffernd_calI_screen.rcbuffer->pybuffer.shape[0];
     }
   }
@@ -2630,7 +2618,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
   __pyx_v_calI_screen = ((PyArrayObject *)__pyx_t_3);
   __pyx_t_3 = 0;
 
-  /* "src/screening/gap_test_p_1.pyx":147
+  /* "src/screening/gap_test_p_1.pyx":148
  * 
  *       # 1. Sort in descending order
  *       if index is None:             # <<<<<<<<<<<<<<
@@ -2641,16 +2629,16 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
   __pyx_t_10 = (__pyx_t_9 != 0);
   if (__pyx_t_10) {
 
-    /* "src/screening/gap_test_p_1.pyx":148
+    /* "src/screening/gap_test_p_1.pyx":149
  *       # 1. Sort in descending order
  *       if index is None:
  *          index = np.argsort(Atcabs)[::-1]             # <<<<<<<<<<<<<<
  * 
  *          # Cdric' notations compliant
  */
-    __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 148, __pyx_L1_error)
+    __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 149, __pyx_L1_error)
     __Pyx_GOTREF(__pyx_t_1);
-    __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_argsort); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 148, __pyx_L1_error)
+    __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_argsort); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 149, __pyx_L1_error)
     __Pyx_GOTREF(__pyx_t_2);
     __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
     __pyx_t_1 = NULL;
@@ -2665,13 +2653,13 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
     }
     __pyx_t_3 = (__pyx_t_1) ? __Pyx_PyObject_Call2Args(__pyx_t_2, __pyx_t_1, ((PyObject *)__pyx_v_Atcabs)) : __Pyx_PyObject_CallOneArg(__pyx_t_2, ((PyObject *)__pyx_v_Atcabs));
     __Pyx_XDECREF(__pyx_t_1); __pyx_t_1 = 0;
-    if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 148, __pyx_L1_error)
+    if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 149, __pyx_L1_error)
     __Pyx_GOTREF(__pyx_t_3);
     __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
-    __pyx_t_2 = __Pyx_PyObject_GetItem(__pyx_t_3, __pyx_slice_); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 148, __pyx_L1_error)
+    __pyx_t_2 = __Pyx_PyObject_GetItem(__pyx_t_3, __pyx_slice_); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 149, __pyx_L1_error)
     __Pyx_GOTREF(__pyx_t_2);
     __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
-    if (!(likely(((__pyx_t_2) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_2, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 148, __pyx_L1_error)
+    if (!(likely(((__pyx_t_2) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_2, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 149, __pyx_L1_error)
     __pyx_t_11 = ((PyArrayObject *)__pyx_t_2);
     {
       __Pyx_BufFmt_StackElem __pyx_stack[1];
@@ -2688,13 +2676,13 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
         __pyx_t_13 = __pyx_t_14 = __pyx_t_15 = 0;
       }
       __pyx_pybuffernd_index.diminfo[0].strides = __pyx_pybuffernd_index.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_index.diminfo[0].shape = __pyx_pybuffernd_index.rcbuffer->pybuffer.shape[0];
-      if (unlikely(__pyx_t_12 < 0)) __PYX_ERR(0, 148, __pyx_L1_error)
+      if (unlikely(__pyx_t_12 < 0)) __PYX_ERR(0, 149, __pyx_L1_error)
     }
     __pyx_t_11 = 0;
     __Pyx_DECREF_SET(__pyx_v_index, ((PyArrayObject *)__pyx_t_2));
     __pyx_t_2 = 0;
 
-    /* "src/screening/gap_test_p_1.pyx":147
+    /* "src/screening/gap_test_p_1.pyx":148
  * 
  *       # 1. Sort in descending order
  *       if index is None:             # <<<<<<<<<<<<<<
@@ -2703,19 +2691,19 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
  */
   }
 
-  /* "src/screening/gap_test_p_1.pyx":151
+  /* "src/screening/gap_test_p_1.pyx":152
  * 
  *          # Cdric' notations compliant
  *       sigmas[1:] = np.cumsum(Atcabs[index])             # <<<<<<<<<<<<<<
  * 
  *       for q in range(n-1, -1, -1):
  */
-  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 151, __pyx_L1_error)
+  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 152, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_3);
-  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_cumsum); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 151, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_cumsum); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 152, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
   __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
-  __pyx_t_3 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_Atcabs), ((PyObject *)__pyx_v_index)); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 151, __pyx_L1_error)
+  __pyx_t_3 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_Atcabs), ((PyObject *)__pyx_v_index)); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 152, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_3);
   __pyx_t_4 = NULL;
   if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_1))) {
@@ -2730,13 +2718,13 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
   __pyx_t_2 = (__pyx_t_4) ? __Pyx_PyObject_Call2Args(__pyx_t_1, __pyx_t_4, __pyx_t_3) : __Pyx_PyObject_CallOneArg(__pyx_t_1, __pyx_t_3);
   __Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
   __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
-  if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 151, __pyx_L1_error)
+  if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 152, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_2);
   __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
-  if (unlikely(PyObject_SetItem(((PyObject *)__pyx_v_sigmas), __pyx_slice__2, __pyx_t_2) < 0)) __PYX_ERR(0, 151, __pyx_L1_error)
+  if (unlikely(PyObject_SetItem(((PyObject *)__pyx_v_sigmas), __pyx_slice__2, __pyx_t_2) < 0)) __PYX_ERR(0, 152, __pyx_L1_error)
   __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
 
-  /* "src/screening/gap_test_p_1.pyx":153
+  /* "src/screening/gap_test_p_1.pyx":154
  *       sigmas[1:] = np.cumsum(Atcabs[index])
  * 
  *       for q in range(n-1, -1, -1):             # <<<<<<<<<<<<<<
@@ -2746,7 +2734,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
   for (__pyx_t_12 = (__pyx_v_n - 1); __pyx_t_12 > -1; __pyx_t_12-=1) {
     __pyx_v_q = __pyx_t_12;
 
-    /* "src/screening/gap_test_p_1.pyx":154
+    /* "src/screening/gap_test_p_1.pyx":155
  * 
  *       for q in range(n-1, -1, -1):
  *          bound = coeff_lbd * self.vec_cumsum_gammas[q] - self.vec_kappa_q[q] * radius - sigmas[q]             # <<<<<<<<<<<<<<
@@ -2756,7 +2744,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
     __pyx_t_16 = __pyx_v_q;
     __pyx_v_bound = (((__pyx_v_coeff_lbd * (__pyx_v_self->vec_cumsum_gammas[__pyx_v_q])) - ((__pyx_v_self->vec_kappa_q[__pyx_v_q]) * __pyx_v_radius)) - (*__Pyx_BufPtrStrided1d(npy_double *, __pyx_pybuffernd_sigmas.rcbuffer->pybuffer.buf, __pyx_t_16, __pyx_pybuffernd_sigmas.diminfo[0].strides)));
 
-    /* "src/screening/gap_test_p_1.pyx":156
+    /* "src/screening/gap_test_p_1.pyx":157
  *          bound = coeff_lbd * self.vec_cumsum_gammas[q] - self.vec_kappa_q[q] * radius - sigmas[q]
  * 
  *          if Atcabs[index[q]] < bound:             # <<<<<<<<<<<<<<
@@ -2768,7 +2756,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
     __pyx_t_10 = (((*__Pyx_BufPtrStrided1d(npy_double *, __pyx_pybuffernd_Atcabs.rcbuffer->pybuffer.buf, __pyx_t_17, __pyx_pybuffernd_Atcabs.diminfo[0].strides)) < __pyx_v_bound) != 0);
     if (__pyx_t_10) {
 
-      /* "src/screening/gap_test_p_1.pyx":157
+      /* "src/screening/gap_test_p_1.pyx":158
  * 
  *          if Atcabs[index[q]] < bound:
  *             l_q = 0             # <<<<<<<<<<<<<<
@@ -2777,7 +2765,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
  */
       __pyx_v_l_q = 0;
 
-      /* "src/screening/gap_test_p_1.pyx":156
+      /* "src/screening/gap_test_p_1.pyx":157
  *          bound = coeff_lbd * self.vec_cumsum_gammas[q] - self.vec_kappa_q[q] * radius - sigmas[q]
  * 
  *          if Atcabs[index[q]] < bound:             # <<<<<<<<<<<<<<
@@ -2787,7 +2775,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
       goto __pyx_L6;
     }
 
-    /* "src/screening/gap_test_p_1.pyx":159
+    /* "src/screening/gap_test_p_1.pyx":160
  *             l_q = 0
  * 
  *          elif Atcabs[index[n-1]] >= bound:             # <<<<<<<<<<<<<<
@@ -2799,7 +2787,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
     __pyx_t_10 = (((*__Pyx_BufPtrStrided1d(npy_double *, __pyx_pybuffernd_Atcabs.rcbuffer->pybuffer.buf, __pyx_t_17, __pyx_pybuffernd_Atcabs.diminfo[0].strides)) >= __pyx_v_bound) != 0);
     if (__pyx_t_10) {
 
-      /* "src/screening/gap_test_p_1.pyx":160
+      /* "src/screening/gap_test_p_1.pyx":161
  * 
  *          elif Atcabs[index[n-1]] >= bound:
  *             return calI_screen             # <<<<<<<<<<<<<<
@@ -2811,7 +2799,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
       __pyx_r = ((PyObject *)__pyx_v_calI_screen);
       goto __pyx_L0;
 
-      /* "src/screening/gap_test_p_1.pyx":159
+      /* "src/screening/gap_test_p_1.pyx":160
  *             l_q = 0
  * 
  *          elif Atcabs[index[n-1]] >= bound:             # <<<<<<<<<<<<<<
@@ -2820,7 +2808,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
  */
     }
 
-    /* "src/screening/gap_test_p_1.pyx":163
+    /* "src/screening/gap_test_p_1.pyx":164
  * 
  *          else:
  *             l_q = q + np.argmax(Atcabs[index[q:]] < bound)             # <<<<<<<<<<<<<<
@@ -2828,27 +2816,27 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
  *          l_min = max(l_min, l_q)
  */
     /*else*/ {
-      __pyx_t_2 = __Pyx_PyInt_From_int(__pyx_v_q); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 163, __pyx_L1_error)
+      __pyx_t_2 = __Pyx_PyInt_From_int(__pyx_v_q); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 164, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_2);
-      __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 163, __pyx_L1_error)
+      __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 164, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_3);
-      __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_argmax); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 163, __pyx_L1_error)
+      __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_argmax); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 164, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_4);
       __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
-      __pyx_t_3 = __Pyx_PyInt_From_int(__pyx_v_q); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 163, __pyx_L1_error)
+      __pyx_t_3 = __Pyx_PyInt_From_int(__pyx_v_q); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 164, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_3);
-      __pyx_t_5 = PySlice_New(__pyx_t_3, Py_None, Py_None); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 163, __pyx_L1_error)
+      __pyx_t_5 = PySlice_New(__pyx_t_3, Py_None, Py_None); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 164, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_5);
       __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
-      __pyx_t_3 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_index), __pyx_t_5); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 163, __pyx_L1_error)
+      __pyx_t_3 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_index), __pyx_t_5); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 164, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_3);
       __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
-      __pyx_t_5 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_Atcabs), __pyx_t_3); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 163, __pyx_L1_error)
+      __pyx_t_5 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_Atcabs), __pyx_t_3); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 164, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_5);
       __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
-      __pyx_t_3 = PyFloat_FromDouble(__pyx_v_bound); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 163, __pyx_L1_error)
+      __pyx_t_3 = PyFloat_FromDouble(__pyx_v_bound); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 164, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_3);
-      __pyx_t_18 = PyObject_RichCompare(__pyx_t_5, __pyx_t_3, Py_LT); __Pyx_XGOTREF(__pyx_t_18); if (unlikely(!__pyx_t_18)) __PYX_ERR(0, 163, __pyx_L1_error)
+      __pyx_t_18 = PyObject_RichCompare(__pyx_t_5, __pyx_t_3, Py_LT); __Pyx_XGOTREF(__pyx_t_18); if (unlikely(!__pyx_t_18)) __PYX_ERR(0, 164, __pyx_L1_error)
       __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
       __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
       __pyx_t_3 = NULL;
@@ -2864,20 +2852,20 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
       __pyx_t_1 = (__pyx_t_3) ? __Pyx_PyObject_Call2Args(__pyx_t_4, __pyx_t_3, __pyx_t_18) : __Pyx_PyObject_CallOneArg(__pyx_t_4, __pyx_t_18);
       __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
       __Pyx_DECREF(__pyx_t_18); __pyx_t_18 = 0;
-      if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 163, __pyx_L1_error)
+      if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 164, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_1);
       __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
-      __pyx_t_4 = PyNumber_Add(__pyx_t_2, __pyx_t_1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 163, __pyx_L1_error)
+      __pyx_t_4 = PyNumber_Add(__pyx_t_2, __pyx_t_1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 164, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_4);
       __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
       __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
-      __pyx_t_19 = __Pyx_PyInt_As_int(__pyx_t_4); if (unlikely((__pyx_t_19 == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 163, __pyx_L1_error)
+      __pyx_t_19 = __Pyx_PyInt_As_int(__pyx_t_4); if (unlikely((__pyx_t_19 == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 164, __pyx_L1_error)
       __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
       __pyx_v_l_q = __pyx_t_19;
     }
     __pyx_L6:;
 
-    /* "src/screening/gap_test_p_1.pyx":165
+    /* "src/screening/gap_test_p_1.pyx":166
  *             l_q = q + np.argmax(Atcabs[index[q:]] < bound)
  * 
  *          l_min = max(l_min, l_q)             # <<<<<<<<<<<<<<
@@ -2894,24 +2882,24 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
     __pyx_v_l_min = __pyx_t_21;
   }
 
-  /* "src/screening/gap_test_p_1.pyx":167
+  /* "src/screening/gap_test_p_1.pyx":168
  *          l_min = max(l_min, l_q)
  * 
  *       calI_screen[index[l_min:]] = True             # <<<<<<<<<<<<<<
  *       return calI_screen
  */
-  __pyx_t_4 = __Pyx_PyInt_From_int(__pyx_v_l_min); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 167, __pyx_L1_error)
+  __pyx_t_4 = __Pyx_PyInt_From_int(__pyx_v_l_min); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 168, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_4);
-  __pyx_t_1 = PySlice_New(__pyx_t_4, Py_None, Py_None); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 167, __pyx_L1_error)
+  __pyx_t_1 = PySlice_New(__pyx_t_4, Py_None, Py_None); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 168, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
   __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
-  __pyx_t_4 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_index), __pyx_t_1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 167, __pyx_L1_error)
+  __pyx_t_4 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_index), __pyx_t_1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 168, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_4);
   __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
-  if (unlikely(PyObject_SetItem(((PyObject *)__pyx_v_calI_screen), __pyx_t_4, Py_True) < 0)) __PYX_ERR(0, 167, __pyx_L1_error)
+  if (unlikely(PyObject_SetItem(((PyObject *)__pyx_v_calI_screen), __pyx_t_4, Py_True) < 0)) __PYX_ERR(0, 168, __pyx_L1_error)
   __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
 
-  /* "src/screening/gap_test_p_1.pyx":168
+  /* "src/screening/gap_test_p_1.pyx":169
  * 
  *       calI_screen[index[l_min:]] = True
  *       return calI_screen             # <<<<<<<<<<<<<<
@@ -2921,7 +2909,7 @@ static PyObject *__pyx_pf_3src_9screening_12gap_test_p_1_16GapTestPequalOne_4app
   __pyx_r = ((PyObject *)__pyx_v_calI_screen);
   goto __pyx_L0;
 
-  /* "src/screening/gap_test_p_1.pyx":85
+  /* "src/screening/gap_test_p_1.pyx":86
  *    @cython.boundscheck(False) # turn off bounds-checking for entire function
  *    @cython.wraparound(False)  # turn off negative index wrapping for entire function
  *    def apply_test(self,             # <<<<<<<<<<<<<<
@@ -4248,7 +4236,6 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = {
   {&__pyx_n_s_AbstractGapScreening, __pyx_k_AbstractGapScreening, sizeof(__pyx_k_AbstractGapScreening), 0, 0, 1, 1},
   {&__pyx_n_s_Atcabs, __pyx_k_Atcabs, sizeof(__pyx_k_Atcabs), 0, 0, 1, 1},
   {&__pyx_n_s_GapTestPequalOne, __pyx_k_GapTestPequalOne, sizeof(__pyx_k_GapTestPequalOne), 0, 0, 1, 1},
-  {&__pyx_kp_s_Gap_test_p_1_cython, __pyx_k_Gap_test_p_1_cython, sizeof(__pyx_k_Gap_test_p_1_cython), 0, 0, 1, 0},
   {&__pyx_n_s_ImportError, __pyx_k_ImportError, sizeof(__pyx_k_ImportError), 0, 0, 1, 1},
   {&__pyx_n_s_TypeError, __pyx_k_TypeError, sizeof(__pyx_k_TypeError), 0, 0, 1, 1},
   {&__pyx_n_s_argmax, __pyx_k_argmax, sizeof(__pyx_k_argmax), 0, 0, 1, 1},
@@ -4264,7 +4251,6 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = {
   {&__pyx_n_s_lbd, __pyx_k_lbd, sizeof(__pyx_k_lbd), 0, 0, 1, 1},
   {&__pyx_n_s_main, __pyx_k_main, sizeof(__pyx_k_main), 0, 0, 1, 1},
   {&__pyx_n_s_name, __pyx_k_name, sizeof(__pyx_k_name), 0, 0, 1, 1},
-  {&__pyx_n_s_name_2, __pyx_k_name_2, sizeof(__pyx_k_name_2), 0, 0, 1, 1},
   {&__pyx_n_s_np, __pyx_k_np, sizeof(__pyx_k_np), 0, 0, 1, 1},
   {&__pyx_n_s_numpy, __pyx_k_numpy, sizeof(__pyx_k_numpy), 0, 0, 1, 1},
   {&__pyx_kp_s_numpy_core_multiarray_failed_to, __pyx_k_numpy_core_multiarray_failed_to, sizeof(__pyx_k_numpy_core_multiarray_failed_to), 0, 0, 1, 0},
@@ -4286,7 +4272,7 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = {
   {0, 0, 0, 0, 0, 0, 0}
 };
 static CYTHON_SMALL_CODE int __Pyx_InitCachedBuiltins(void) {
-  __pyx_builtin_range = __Pyx_GetBuiltinName(__pyx_n_s_range); if (!__pyx_builtin_range) __PYX_ERR(0, 62, __pyx_L1_error)
+  __pyx_builtin_range = __Pyx_GetBuiltinName(__pyx_n_s_range); if (!__pyx_builtin_range) __PYX_ERR(0, 63, __pyx_L1_error)
   __pyx_builtin_TypeError = __Pyx_GetBuiltinName(__pyx_n_s_TypeError); if (!__pyx_builtin_TypeError) __PYX_ERR(1, 2, __pyx_L1_error)
   __pyx_builtin_ImportError = __Pyx_GetBuiltinName(__pyx_n_s_ImportError); if (!__pyx_builtin_ImportError) __PYX_ERR(2, 945, __pyx_L1_error)
   return 0;
@@ -4298,25 +4284,25 @@ static CYTHON_SMALL_CODE int __Pyx_InitCachedConstants(void) {
   __Pyx_RefNannyDeclarations
   __Pyx_RefNannySetupContext("__Pyx_InitCachedConstants", 0);
 
-  /* "src/screening/gap_test_p_1.pyx":148
+  /* "src/screening/gap_test_p_1.pyx":149
  *       # 1. Sort in descending order
  *       if index is None:
  *          index = np.argsort(Atcabs)[::-1]             # <<<<<<<<<<<<<<
  * 
  *          # Cdric' notations compliant
  */
-  __pyx_slice_ = PySlice_New(Py_None, Py_None, __pyx_int_neg_1); if (unlikely(!__pyx_slice_)) __PYX_ERR(0, 148, __pyx_L1_error)
+  __pyx_slice_ = PySlice_New(Py_None, Py_None, __pyx_int_neg_1); if (unlikely(!__pyx_slice_)) __PYX_ERR(0, 149, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_slice_);
   __Pyx_GIVEREF(__pyx_slice_);
 
-  /* "src/screening/gap_test_p_1.pyx":151
+  /* "src/screening/gap_test_p_1.pyx":152
  * 
  *          # Cdric' notations compliant
  *       sigmas[1:] = np.cumsum(Atcabs[index])             # <<<<<<<<<<<<<<
  * 
  *       for q in range(n-1, -1, -1):
  */
-  __pyx_slice__2 = PySlice_New(__pyx_int_1, Py_None, Py_None); if (unlikely(!__pyx_slice__2)) __PYX_ERR(0, 151, __pyx_L1_error)
+  __pyx_slice__2 = PySlice_New(__pyx_int_1, Py_None, Py_None); if (unlikely(!__pyx_slice__2)) __PYX_ERR(0, 152, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_slice__2);
   __Pyx_GIVEREF(__pyx_slice__2);
 
@@ -4682,7 +4668,7 @@ if (!__Pyx_RefNanny) {
   if (__Pyx_init_sys_getdefaultencoding_params() < 0) __PYX_ERR(0, 1, __pyx_L1_error)
   #endif
   if (__pyx_module_is_main_src__screening__gap_test_p_1) {
-    if (PyObject_SetAttr(__pyx_m, __pyx_n_s_name_2, __pyx_n_s_main) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
+    if (PyObject_SetAttr(__pyx_m, __pyx_n_s_name, __pyx_n_s_main) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
   }
   #if PY_MAJOR_VERSION >= 3
   {
@@ -5866,20 +5852,6 @@ static CYTHON_INLINE PyObject* __Pyx_PyObject_CallOneArg(PyObject *func, PyObjec
 }
 #endif
 
-/* PyObjectSetAttrStr */
-  #if CYTHON_USE_TYPE_SLOTS
-static CYTHON_INLINE int __Pyx_PyObject_SetAttrStr(PyObject* obj, PyObject* attr_name, PyObject* value) {
-    PyTypeObject* tp = Py_TYPE(obj);
-    if (likely(tp->tp_setattro))
-        return tp->tp_setattro(obj, attr_name, value);
-#if PY_MAJOR_VERSION < 3
-    if (likely(tp->tp_setattr))
-        return tp->tp_setattr(obj, PyString_AS_STRING(attr_name), value);
-#endif
-    return PyObject_SetAttr(obj, attr_name, value);
-}
-#endif
-
 /* PyErrFetchRestore */
   #if CYTHON_FAST_THREAD_STATE
 static CYTHON_INLINE void __Pyx_ErrRestoreInState(PyThreadState *tstate, PyObject *type, PyObject *value, PyObject *tb) {
@@ -6429,7 +6401,7 @@ static CYTHON_INLINE PyObject* __Pyx_PyObject_GetAttrStrNoError(PyObject* obj, P
   static int __Pyx_setup_reduce_is_named(PyObject* meth, PyObject* name) {
   int ret;
   PyObject *name_attr;
-  name_attr = __Pyx_PyObject_GetAttrStr(meth, __pyx_n_s_name_2);
+  name_attr = __Pyx_PyObject_GetAttrStr(meth, __pyx_n_s_name);
   if (likely(name_attr)) {
       ret = PyObject_RichCompareBool(name_attr, name, Py_EQ);
   } else {
diff --git a/src/screening/gap_test_p_1.pyx b/src/screening/gap_test_p_1.pyx
index 369afdd87394a0bed898c3c988718783d08aaddd..9fc4aa1b89cebdf328d49436599e0dfc4d02977c 100644
--- a/src/screening/gap_test_p_1.pyx
+++ b/src/screening/gap_test_p_1.pyx
@@ -29,6 +29,7 @@ cdef class GapTestPequalOne:
 
    cdef double* vec_cumsum_gammas
    cdef double* vec_kappa_q
+   cdef char* name
 
    # cdef np.ndarray vec_cumsum_gammas
    # cdef np.ndarray vec_kappa_q