From 041b976b41c282f1635c39026039be8febc395c4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Cl=C3=A9ment=20Elvira?= <clement.elvira@centralesupelec.fr>
Date: Thu, 10 Nov 2022 16:16:38 +0100
Subject: [PATCH] Robustify SIAM xp 1

---
 experiments/SIAM/xp_1_balls/process_data.py   |   4 +-
 experiments/SIAM/xp_1_balls/xp_b_screening.py |   3 +-
 experiments/SIAM/xp_1_balls/xp_c_viz.py       | 117 +++++
 ...xp_c_viz_fig1.py => xp_c_viz_SIAM_fig1.py} |  15 +-
 ...xp_c_viz_fig2.py => xp_c_viz_SIMA_fig2.py} |  10 +-
 slopescreening/screening/ascreening.py        |   5 +-
 slopescreening/screening/gap_test_all.c       | 414 ++++++++++--------
 slopescreening/screening/gap_test_all.pyx     |   4 +
 slopescreening/screening/gap_test_p_1.c       | 290 +++++++-----
 slopescreening/screening/gap_test_p_1.pyx     |   4 +
 slopescreening/screening/gap_test_p_q.py      |   1 +
 11 files changed, 563 insertions(+), 304 deletions(-)
 create mode 100644 experiments/SIAM/xp_1_balls/xp_c_viz.py
 rename experiments/SIAM/xp_1_balls/{xp_c_viz_fig1.py => xp_c_viz_SIAM_fig1.py} (87%)
 rename experiments/SIAM/xp_1_balls/{xp_c_viz_fig2.py => xp_c_viz_SIMA_fig2.py} (93%)

diff --git a/experiments/SIAM/xp_1_balls/process_data.py b/experiments/SIAM/xp_1_balls/process_data.py
index 60b3a9a..3cf5430 100644
--- a/experiments/SIAM/xp_1_balls/process_data.py
+++ b/experiments/SIAM/xp_1_balls/process_data.py
@@ -19,7 +19,8 @@ def process(setup, log=True):
 
    mat_nb_zero_detected = screenings['mat_nb_zero_detected']
    nb_test = mat_nb_zero_detected.shape[0]
-   list_tests = screenings["list_test"]
+   list_tests   = screenings["list_test"]
+   list_legends = screenings["list_legends"]
 
    # ---- log ----
    if log:
@@ -48,6 +49,7 @@ def process(setup, log=True):
    return {
       "mat_pc_detected": np.mean(mat_pc_detected, axis=5),
       "list_tests": list_tests,
+      "list_legends": list_legends,
    }
 
    # {
diff --git a/experiments/SIAM/xp_1_balls/xp_b_screening.py b/experiments/SIAM/xp_1_balls/xp_b_screening.py
index 56691e1..f651127 100644
--- a/experiments/SIAM/xp_1_balls/xp_b_screening.py
+++ b/experiments/SIAM/xp_1_balls/xp_b_screening.py
@@ -120,7 +120,8 @@ for i_dic in range(setup.nb_dic):
          np.savez(
             screening_filename,
             mat_nb_zero_detected=mat_nb_zero_detected,
-            list_test = [test.get_name() for test in list_tests],
+            list_test    = [test.get_name() for test in list_tests],
+            list_legends = [test.get_legend_name() for test in list_tests],
             version = __version__,
             allow_pickle=True
          )
\ No newline at end of file
diff --git a/experiments/SIAM/xp_1_balls/xp_c_viz.py b/experiments/SIAM/xp_1_balls/xp_c_viz.py
new file mode 100644
index 0000000..d1ebf7a
--- /dev/null
+++ b/experiments/SIAM/xp_1_balls/xp_c_viz.py
@@ -0,0 +1,117 @@
+# -*- coding: utf-8 -*-
+import argparse
+from pathlib import Path
+
+import numpy as np
+
+# XP import
+from experiments.SIAM.slopepb import SlopePb
+from experiments.SIAM.setup import Setup
+
+import xpparams
+from process_data import process
+
+
+parser=argparse.ArgumentParser()
+parser.add_argument('--noshow', help='do not display figures', action="store_true")
+parser.add_argument('--noverbose', help='disable printing if true', 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
+
+
+# -------------------------
+#        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 = 0
+i_lbd = 0
+
+fs=22
+fs_ylabels = 20
+list_colors  = ["tab:blue", "tab:orange", "tab:green"]
+list_legends = dic_process_oscar["list_legends"]
+# list_legends = ["$p_q=q\\;\\forall q$", "all", "$p_q=1\\;\\forall q$"]
+# "best $r_q \\;\\forall q$ "
+
+if not args.noverbose:
+   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,
+)
+ax.set_ylabel(
+   "% of zero entries detected",
+   fontsize=fs+2,
+)
+
+for tick in ax.xaxis.get_major_ticks():
+   tick.label.set_fontsize(20)
+
+for tick in ax.yaxis.get_major_ticks():
+   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,
+)
+
+ax.set_xscale("log")
+ax.set_xlim([1e-6, 1e0])
+ax.set_ylim([-2, 102])
+
+if args.save:
+   folderfig = "figs"
+   Path(folderfig).mkdir(parents=True, exist_ok=True)
+
+   filename = folderfig + f"/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/experiments/SIAM/xp_1_balls/xp_c_viz_fig1.py b/experiments/SIAM/xp_1_balls/xp_c_viz_SIAM_fig1.py
similarity index 87%
rename from experiments/SIAM/xp_1_balls/xp_c_viz_fig1.py
rename to experiments/SIAM/xp_1_balls/xp_c_viz_SIAM_fig1.py
index db26ba1..9fff40f 100644
--- a/experiments/SIAM/xp_1_balls/xp_c_viz_fig1.py
+++ b/experiments/SIAM/xp_1_balls/xp_c_viz_SIAM_fig1.py
@@ -20,6 +20,11 @@ parser.add_argument('--id', help='setup id', type=str, default="SIAM")
 parser.add_argument('--i_seq', type=int, default=0)
 args=parser.parse_args()
 
+
+if args.id != "SIAM":
+   raise Exception("this file has been created to reproduce the SIAM figures. Use the SIAM id instead")
+
+
 import matplotlib
 if args.noshow:
    matplotlib.use('PS')
@@ -84,12 +89,14 @@ 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$"]
+list_legends = dic_process_oscar["list_legends"]
+# 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]}")
+if not args.noverbose:
+   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)
diff --git a/experiments/SIAM/xp_1_balls/xp_c_viz_fig2.py b/experiments/SIAM/xp_1_balls/xp_c_viz_SIMA_fig2.py
similarity index 93%
rename from experiments/SIAM/xp_1_balls/xp_c_viz_fig2.py
rename to experiments/SIAM/xp_1_balls/xp_c_viz_SIMA_fig2.py
index ac1599d..ac1879f 100644
--- a/experiments/SIAM/xp_1_balls/xp_c_viz_fig2.py
+++ b/experiments/SIAM/xp_1_balls/xp_c_viz_SIMA_fig2.py
@@ -19,6 +19,9 @@ parser.add_argument('--save', help='save figure', action="store_true")
 parser.add_argument('--id', help='setup id', type=str, default="SIAM")
 args=parser.parse_args()
 
+if args.id != "SIAM":
+   raise Exception("this file has been created to reproduce the SIAM figures. Use the SIAM id instead")
+
 import matplotlib
 if args.noshow:
    matplotlib.use("ps")
@@ -52,6 +55,7 @@ list_tests_oscar      = dic_process_oscar["list_tests"]
 
 i_dic = 0
 i_lbd = 1
+n_lbd = len(setup_oscar.list_ratio_lbd)
 
 fs=22
 fs_ylabels = 20
@@ -104,7 +108,7 @@ font_ttt = font_manager.FontProperties(
    size=fs+6
 )
 
-f, ax = plt.subplots(1, 3, figsize=(18, 6), sharex=True, sharey=True)
+f, ax = plt.subplots(1, n_lbd, figsize=(18, 6), sharex=True, sharey=True)
 
 
 for i_seq, _ in enumerate(setup_oscar.list_sequence):
@@ -120,7 +124,7 @@ for i_seq, _ in enumerate(setup_oscar.list_sequence):
          tick.label.set_fontsize(18)
 
       
-      for i_lbd in range(3):
+      for i_lbd in range(n_lbd):
          ax[i_seq].plot(
             xpparams.vec_offsets, 
             100 * mat_pc_detected_oscar[i_test, :, i_dic, i_seq, i_lbd],
@@ -144,7 +148,7 @@ ax[0].set_ylabel(
    fontsize=fs_ylabels+6,
 )
 
-ax[0].legend([object, object, object], ['$\\lambda / \\lambda_{\\max}=' + str(x) + '$' for x in [.3, .5, .8]],
+ax[0].legend([object, object, object], ['$\\lambda / \\lambda_{\\max}=' + str(x) + '$' for x in setup_oscar.list_ratio_lbd],
            handler_map={object: AnyObjectHandler()},
            prop=font_math
 )
diff --git a/slopescreening/screening/ascreening.py b/slopescreening/screening/ascreening.py
index 77cd284..503f99c 100644
--- a/slopescreening/screening/ascreening.py
+++ b/slopescreening/screening/ascreening.py
@@ -9,11 +9,14 @@ class AbstractGapScreening(ABC):
    def __init__(self):
       super(AbstractGapScreening, self).__init__()
 
-      self.name = "NONAME"
+      self.name = "void"
+      self.legend_name = "void"
 
    def get_name(self):
       return self.name
    
+   def get_legend_name(self):
+      return self.legend_name
 
    @abstractmethod
    def apply_test(self, Atu_abs, gap, lbd, vec_gammas, coeff_dual_scaling=1., offset_radius=0, index=None) -> np.ndarray:
diff --git a/slopescreening/screening/gap_test_all.c b/slopescreening/screening/gap_test_all.c
index 64673e1..62048d6 100644
--- a/slopescreening/screening/gap_test_all.c
+++ b/slopescreening/screening/gap_test_all.c
@@ -1925,6 +1925,7 @@ static PyObject *__pyx_builtin_range;
 static PyObject *__pyx_builtin_TypeError;
 static PyObject *__pyx_builtin_ImportError;
 static const char __pyx_k_np[] = "np";
+static const char __pyx_k_all[] = "all";
 static const char __pyx_k_gap[] = "gap";
 static const char __pyx_k_lbd[] = "lbd";
 static const char __pyx_k_main[] = "__main__";
@@ -1964,6 +1965,7 @@ static PyObject *__pyx_n_s_Atcabs;
 static PyObject *__pyx_n_s_GapTestAll;
 static PyObject *__pyx_n_s_ImportError;
 static PyObject *__pyx_n_s_TypeError;
+static PyObject *__pyx_n_s_all;
 static PyObject *__pyx_n_s_arange;
 static PyObject *__pyx_n_s_argsort;
 static PyObject *__pyx_n_s_cline_in_traceback;
@@ -1997,9 +1999,10 @@ static PyObject *__pyx_n_s_zeros;
 static int __pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll___init__(struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *__pyx_v_self, PyArrayObject *__pyx_v_vec_gammas); /* proto */
 static void __pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_2__dealloc__(struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *__pyx_v_self); /* proto */
 static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_4get_name(struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *__pyx_v_self); /* proto */
-static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_6apply_test(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *__pyx_v_self, PyArrayObject *__pyx_v_Atcabs, double __pyx_v_gap, double __pyx_v_lbd, PyArrayObject *__pyx_v_vec_gammas, double __pyx_v_coeff_dual_scaling, double __pyx_v_offset_radius, PyArrayObject *__pyx_v_index); /* proto */
-static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_8__reduce_cython__(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *__pyx_v_self); /* proto */
-static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_10__setstate_cython__(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *__pyx_v_self, CYTHON_UNUSED PyObject *__pyx_v___pyx_state); /* proto */
+static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_6get_legend_name(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *__pyx_v_self); /* proto */
+static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_8apply_test(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *__pyx_v_self, PyArrayObject *__pyx_v_Atcabs, double __pyx_v_gap, double __pyx_v_lbd, PyArrayObject *__pyx_v_vec_gammas, double __pyx_v_coeff_dual_scaling, double __pyx_v_offset_radius, PyArrayObject *__pyx_v_index); /* proto */
+static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_10__reduce_cython__(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *__pyx_v_self); /* proto */
+static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_12__setstate_cython__(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *__pyx_v_self, CYTHON_UNUSED PyObject *__pyx_v___pyx_state); /* proto */
 static PyObject *__pyx_tp_new_14slopescreening_9screening_12gap_test_all_GapTestAll(PyTypeObject *t, PyObject *a, PyObject *k); /*proto*/
 static PyObject *__pyx_int_0;
 static PyObject *__pyx_int_1;
@@ -2268,7 +2271,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
  *    def get_name(self):
  *       return self.name             # <<<<<<<<<<<<<<
  * 
- * 
+ *    @cython.boundscheck(False) # turn off bounds-checking for entire function
  */
   __Pyx_XDECREF(__pyx_r);
   __pyx_t_1 = __Pyx_PyBytes_FromString(__pyx_v_self->name); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 59, __pyx_L1_error)
@@ -2296,7 +2299,60 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
   return __pyx_r;
 }
 
-/* "slopescreening/screening/gap_test_all.pyx":64
+/* "slopescreening/screening/gap_test_all.pyx":63
+ *    @cython.boundscheck(False) # turn off bounds-checking for entire function
+ *    @cython.wraparound(False)  # turn off negative index wrapping for entire function
+ *    def get_legend_name(self):             # <<<<<<<<<<<<<<
+ *       return "all"
+ * 
+ */
+
+/* Python wrapper */
+static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_7get_legend_name(PyObject *__pyx_v_self, CYTHON_UNUSED PyObject *unused); /*proto*/
+static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_7get_legend_name(PyObject *__pyx_v_self, CYTHON_UNUSED PyObject *unused) {
+  PyObject *__pyx_r = 0;
+  __Pyx_RefNannyDeclarations
+  __Pyx_RefNannySetupContext("get_legend_name (wrapper)", 0);
+  __pyx_r = __pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_6get_legend_name(((struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *)__pyx_v_self));
+
+  /* function exit code */
+  __Pyx_RefNannyFinishContext();
+  return __pyx_r;
+}
+
+static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_6get_legend_name(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *__pyx_v_self) {
+  PyObject *__pyx_r = NULL;
+  __Pyx_RefNannyDeclarations
+  __Pyx_RefNannySetupContext("get_legend_name", 0);
+
+  /* "slopescreening/screening/gap_test_all.pyx":64
+ *    @cython.wraparound(False)  # turn off negative index wrapping for entire function
+ *    def get_legend_name(self):
+ *       return "all"             # <<<<<<<<<<<<<<
+ * 
+ *    @cython.boundscheck(False) # turn off bounds-checking for entire function
+ */
+  __Pyx_XDECREF(__pyx_r);
+  __Pyx_INCREF(__pyx_n_s_all);
+  __pyx_r = __pyx_n_s_all;
+  goto __pyx_L0;
+
+  /* "slopescreening/screening/gap_test_all.pyx":63
+ *    @cython.boundscheck(False) # turn off bounds-checking for entire function
+ *    @cython.wraparound(False)  # turn off negative index wrapping for entire function
+ *    def get_legend_name(self):             # <<<<<<<<<<<<<<
+ *       return "all"
+ * 
+ */
+
+  /* function exit code */
+  __pyx_L0:;
+  __Pyx_XGIVEREF(__pyx_r);
+  __Pyx_RefNannyFinishContext();
+  return __pyx_r;
+}
+
+/* "slopescreening/screening/gap_test_all.pyx":68
  *    @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,             # <<<<<<<<<<<<<<
@@ -2305,9 +2361,9 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
  */
 
 /* Python wrapper */
-static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_7apply_test(PyObject *__pyx_v_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
-static char __pyx_doc_14slopescreening_9screening_12gap_test_all_10GapTestAll_6apply_test[] = " Apply the Gap safe sphere screening test\n\n      Implementation of the safe screening test\n      \n      Parameters\n      ----------\n      Atcabs : np.ndarray\n         vector matA.T @ vecc where vec is dual admissible in absolute value\n         size [n]\n      gap : positive float\n         duality gap\n      lbd : positive float\n         regularization parameters\n      vec_gammas : np.ndarray\n         slope parameters\n         size [n,]\n      coeff_dual_scaling : positif float\n         If coeff_dual_scaling is not feasible, dual scaling factor\n         such taht vecu / coeff_dual_scaling os dual feasible\n         Here for code optimization purposes\n         Default value is 1. (vecu is feasible)\n      offset_radius : float\n         additive term added to the redius\n         default is 0\n      index : np.ndarray\n         Array of indices that sort Atu in absolute value\n         default is None\n      \n      Returns\n      -------\n      calI_screen : np.ndarray\n         vector of boolean\n         True if screening test passes and False otherwise\n         size [n,]\n      ";
-static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_7apply_test(PyObject *__pyx_v_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
+static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_9apply_test(PyObject *__pyx_v_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
+static char __pyx_doc_14slopescreening_9screening_12gap_test_all_10GapTestAll_8apply_test[] = " Apply the Gap safe sphere screening test\n\n      Implementation of the safe screening test\n      \n      Parameters\n      ----------\n      Atcabs : np.ndarray\n         vector matA.T @ vecc where vec is dual admissible in absolute value\n         size [n]\n      gap : positive float\n         duality gap\n      lbd : positive float\n         regularization parameters\n      vec_gammas : np.ndarray\n         slope parameters\n         size [n,]\n      coeff_dual_scaling : positif float\n         If coeff_dual_scaling is not feasible, dual scaling factor\n         such taht vecu / coeff_dual_scaling os dual feasible\n         Here for code optimization purposes\n         Default value is 1. (vecu is feasible)\n      offset_radius : float\n         additive term added to the redius\n         default is 0\n      index : np.ndarray\n         Array of indices that sort Atu in absolute value\n         default is None\n      \n      Returns\n      -------\n      calI_screen : np.ndarray\n         vector of boolean\n         True if screening test passes and False otherwise\n         size [n,]\n      ";
+static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_9apply_test(PyObject *__pyx_v_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
   PyArrayObject *__pyx_v_Atcabs = 0;
   double __pyx_v_gap;
   double __pyx_v_lbd;
@@ -2325,7 +2381,7 @@ static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAl
     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};
 
-    /* "slopescreening/screening/gap_test_all.pyx":71
+    /* "slopescreening/screening/gap_test_all.pyx":75
  *       double coeff_dual_scaling =1.,
  *       double offset_radius=0.,
  *       np.ndarray[long, ndim=1] index=None             # <<<<<<<<<<<<<<
@@ -2363,19 +2419,19 @@ static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAl
         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, 64, __pyx_L3_error)
+          __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, 1); __PYX_ERR(0, 68, __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, 64, __pyx_L3_error)
+          __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, 2); __PYX_ERR(0, 68, __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, 64, __pyx_L3_error)
+          __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, 3); __PYX_ERR(0, 68, __pyx_L3_error)
         }
         CYTHON_FALLTHROUGH;
         case  4:
@@ -2397,7 +2453,7 @@ static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAl
         }
       }
       if (unlikely(kw_args > 0)) {
-        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "apply_test") < 0)) __PYX_ERR(0, 64, __pyx_L3_error)
+        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "apply_test") < 0)) __PYX_ERR(0, 68, __pyx_L3_error)
       }
     } else {
       switch (PyTuple_GET_SIZE(__pyx_args)) {
@@ -2416,16 +2472,16 @@ static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAl
       }
     }
     __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, 66, __pyx_L3_error)
-    __pyx_v_lbd = __pyx_PyFloat_AsDouble(values[2]); if (unlikely((__pyx_v_lbd == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 67, __pyx_L3_error)
+    __pyx_v_gap = __pyx_PyFloat_AsDouble(values[1]); if (unlikely((__pyx_v_gap == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 70, __pyx_L3_error)
+    __pyx_v_lbd = __pyx_PyFloat_AsDouble(values[2]); if (unlikely((__pyx_v_lbd == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 71, __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, 69, __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, 73, __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, 70, __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, 74, __pyx_L3_error)
     } else {
       __pyx_v_offset_radius = ((double)0.);
     }
@@ -2433,18 +2489,18 @@ static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAl
   }
   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, 64, __pyx_L3_error)
+  __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 68, __pyx_L3_error)
   __pyx_L3_error:;
   __Pyx_AddTraceback("slopescreening.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, 65, __pyx_L1_error)
-  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_vec_gammas), __pyx_ptype_5numpy_ndarray, 1, "vec_gammas", 0))) __PYX_ERR(0, 68, __pyx_L1_error)
-  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_index), __pyx_ptype_5numpy_ndarray, 1, "index", 0))) __PYX_ERR(0, 71, __pyx_L1_error)
-  __pyx_r = __pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_6apply_test(((struct __pyx_obj_14slopescreening_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);
+  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_Atcabs), __pyx_ptype_5numpy_ndarray, 1, "Atcabs", 0))) __PYX_ERR(0, 69, __pyx_L1_error)
+  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_vec_gammas), __pyx_ptype_5numpy_ndarray, 1, "vec_gammas", 0))) __PYX_ERR(0, 72, __pyx_L1_error)
+  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_index), __pyx_ptype_5numpy_ndarray, 1, "index", 0))) __PYX_ERR(0, 75, __pyx_L1_error)
+  __pyx_r = __pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_8apply_test(((struct __pyx_obj_14slopescreening_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);
 
-  /* "slopescreening/screening/gap_test_all.pyx":64
+  /* "slopescreening/screening/gap_test_all.pyx":68
  *    @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,             # <<<<<<<<<<<<<<
@@ -2461,7 +2517,7 @@ static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAl
   return __pyx_r;
 }
 
-static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_6apply_test(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *__pyx_v_self, PyArrayObject *__pyx_v_Atcabs, double __pyx_v_gap, double __pyx_v_lbd, PyArrayObject *__pyx_v_vec_gammas, double __pyx_v_coeff_dual_scaling, double __pyx_v_offset_radius, PyArrayObject *__pyx_v_index) {
+static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_8apply_test(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *__pyx_v_self, PyArrayObject *__pyx_v_Atcabs, double __pyx_v_gap, double __pyx_v_lbd, PyArrayObject *__pyx_v_vec_gammas, double __pyx_v_coeff_dual_scaling, double __pyx_v_offset_radius, PyArrayObject *__pyx_v_index) {
   double __pyx_v_radius;
   PyArrayObject *__pyx_v_coeff_lbd_gamma = 0;
   int __pyx_v_n;
@@ -2559,35 +2615,35 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
   __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, 64, __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, 68, __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, 64, __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, 68, __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, 64, __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, 68, __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];
 
-  /* "slopescreening/screening/gap_test_all.pyx":113
+  /* "slopescreening/screening/gap_test_all.pyx":117
  *       # 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, 113, __pyx_L1_error)
+  __pyx_t_1 = PyFloat_FromDouble(__pyx_v_coeff_dual_scaling); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 117, __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, 113, __pyx_L1_error)
+  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 117, __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, 113, __pyx_L1_error)
+  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_sqrt); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 117, __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, 113, __pyx_L1_error)
+  __pyx_t_3 = PyFloat_FromDouble((2.0 * __pyx_v_gap)); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 117, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_3);
   __pyx_t_5 = NULL;
   if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_4))) {
@@ -2602,42 +2658,42 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
   __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, 113, __pyx_L1_error)
+  if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 117, __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, 113, __pyx_L1_error)
+  __pyx_t_4 = PyNumber_Multiply(__pyx_t_1, __pyx_t_2); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 117, __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, 113, __pyx_L1_error)
+  __pyx_t_2 = PyFloat_FromDouble(__pyx_v_offset_radius); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 117, __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, 113, __pyx_L1_error)
+  __pyx_t_1 = PyNumber_Add(__pyx_t_4, __pyx_t_2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 117, __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, 113, __pyx_L1_error)
+  __pyx_t_6 = __pyx_PyFloat_AsDouble(__pyx_t_1); if (unlikely((__pyx_t_6 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 117, __pyx_L1_error)
   __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
   __pyx_v_radius = __pyx_t_6;
 
-  /* "slopescreening/screening/gap_test_all.pyx":115
+  /* "slopescreening/screening/gap_test_all.pyx":119
  *       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, 115, __pyx_L1_error)
+  __pyx_t_1 = PyFloat_FromDouble((__pyx_v_coeff_dual_scaling * __pyx_v_lbd)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 119, __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, 115, __pyx_L1_error)
+  __pyx_t_2 = PyNumber_Multiply(__pyx_t_1, ((PyObject *)__pyx_v_vec_gammas)); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 119, __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, 115, __pyx_L1_error)
+  if (!(likely(((__pyx_t_2) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_2, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 119, __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, 115, __pyx_L1_error)
+      __PYX_ERR(0, 119, __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];
     }
   }
@@ -2645,7 +2701,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
   __pyx_v_coeff_lbd_gamma = ((PyArrayObject *)__pyx_t_2);
   __pyx_t_2 = 0;
 
-  /* "slopescreening/screening/gap_test_all.pyx":117
+  /* "slopescreening/screening/gap_test_all.pyx":121
  *       cdef np.ndarray[np.npy_double, ndim=1] coeff_lbd_gamma = coeff_dual_scaling * lbd * vec_gammas
  * 
  *       cdef int n = Atcabs.shape[0]             # <<<<<<<<<<<<<<
@@ -2654,7 +2710,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
  */
   __pyx_v_n = (__pyx_v_Atcabs->dimensions[0]);
 
-  /* "slopescreening/screening/gap_test_all.pyx":119
+  /* "slopescreening/screening/gap_test_all.pyx":123
  *       cdef int n = Atcabs.shape[0]
  *       cdef int k, q, r
  *       cdef double tau = 0             # <<<<<<<<<<<<<<
@@ -2663,40 +2719,40 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
  */
   __pyx_v_tau = 0.0;
 
-  /* "slopescreening/screening/gap_test_all.pyx":123
+  /* "slopescreening/screening/gap_test_all.pyx":127
  * 
  *       # 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, 123, __pyx_L1_error)
+  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 127, __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, 123, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_arange); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 127, __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, 123, __pyx_L1_error)
+  __pyx_t_2 = __Pyx_PyInt_From_long((__pyx_v_n + 1)); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 127, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_2);
-  __pyx_t_4 = PyTuple_New(1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 123, __pyx_L1_error)
+  __pyx_t_4 = PyTuple_New(1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 127, __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, 123, __pyx_L1_error)
+  __pyx_t_2 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 127, __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, 123, __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, 123, __pyx_L1_error)
+  if (PyDict_SetItem(__pyx_t_2, __pyx_n_s_dtype, ((PyObject *)(&PyLong_Type))) < 0) __PYX_ERR(0, 127, __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, 127, __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, 123, __pyx_L1_error)
+  if (!(likely(((__pyx_t_3) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_3, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 127, __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, 123, __pyx_L1_error)
+      __PYX_ERR(0, 127, __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];
     }
   }
@@ -2704,40 +2760,40 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
   __pyx_v_vec_p_star = ((PyArrayObject *)__pyx_t_3);
   __pyx_t_3 = 0;
 
-  /* "slopescreening/screening/gap_test_all.pyx":124
+  /* "slopescreening/screening/gap_test_all.pyx":128
  *       # 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, 124, __pyx_L1_error)
+  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 128, __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, 124, __pyx_L1_error)
+  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_arange); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 128, __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, 124, __pyx_L1_error)
+  __pyx_t_3 = __Pyx_PyInt_From_long((__pyx_v_n + 1)); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 128, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_3);
-  __pyx_t_4 = PyTuple_New(1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 124, __pyx_L1_error)
+  __pyx_t_4 = PyTuple_New(1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 128, __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, 124, __pyx_L1_error)
+  __pyx_t_3 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 128, __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, 124, __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, 124, __pyx_L1_error)
+  if (PyDict_SetItem(__pyx_t_3, __pyx_n_s_dtype, ((PyObject *)(&PyLong_Type))) < 0) __PYX_ERR(0, 128, __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, 128, __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, 124, __pyx_L1_error)
+  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 128, __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, 124, __pyx_L1_error)
+      __PYX_ERR(0, 128, __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];
     }
   }
@@ -2745,40 +2801,40 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
   __pyx_v_vec_q_star = ((PyArrayObject *)__pyx_t_1);
   __pyx_t_1 = 0;
 
-  /* "slopescreening/screening/gap_test_all.pyx":127
+  /* "slopescreening/screening/gap_test_all.pyx":131
  * 
  *       # 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, 127, __pyx_L1_error)
+  __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 131, __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, 127, __pyx_L1_error)
+  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_zeros); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 131, __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, 127, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyInt_From_int(__pyx_v_n); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 131, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
-  __pyx_t_4 = PyTuple_New(1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 127, __pyx_L1_error)
+  __pyx_t_4 = PyTuple_New(1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 131, __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, 127, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 131, __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, 127, __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, 127, __pyx_L1_error)
+  if (PyDict_SetItem(__pyx_t_1, __pyx_n_s_dtype, ((PyObject*)&PyBool_Type)) < 0) __PYX_ERR(0, 131, __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, 131, __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, 127, __pyx_L1_error)
+  if (!(likely(((__pyx_t_2) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_2, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 131, __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, 127, __pyx_L1_error)
+      __PYX_ERR(0, 131, __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];
     }
   }
@@ -2786,7 +2842,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
   __pyx_v_calI_screen = ((PyArrayObject *)__pyx_t_2);
   __pyx_t_2 = 0;
 
-  /* "slopescreening/screening/gap_test_all.pyx":130
+  /* "slopescreening/screening/gap_test_all.pyx":134
  * 
  *       # 1. Sort in descending order
  *       if index is None:             # <<<<<<<<<<<<<<
@@ -2797,16 +2853,16 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
   __pyx_t_12 = (__pyx_t_11 != 0);
   if (__pyx_t_12) {
 
-    /* "slopescreening/screening/gap_test_all.pyx":131
+    /* "slopescreening/screening/gap_test_all.pyx":135
  *       # 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, 131, __pyx_L1_error)
+    __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 135, __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, 131, __pyx_L1_error)
+    __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_argsort); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 135, __pyx_L1_error)
     __Pyx_GOTREF(__pyx_t_4);
     __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
     __pyx_t_1 = NULL;
@@ -2821,13 +2877,13 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
     }
     __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, 131, __pyx_L1_error)
+    if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 135, __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, 131, __pyx_L1_error)
+    __pyx_t_4 = __Pyx_PyObject_GetItem(__pyx_t_2, __pyx_slice_); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 135, __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, 131, __pyx_L1_error)
+    if (!(likely(((__pyx_t_4) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_4, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 135, __pyx_L1_error)
     __pyx_t_13 = ((PyArrayObject *)__pyx_t_4);
     {
       __Pyx_BufFmt_StackElem __pyx_stack[1];
@@ -2844,13 +2900,13 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
         __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, 131, __pyx_L1_error)
+      if (unlikely(__pyx_t_14 < 0)) __PYX_ERR(0, 135, __pyx_L1_error)
     }
     __pyx_t_13 = 0;
     __Pyx_DECREF_SET(__pyx_v_index, ((PyArrayObject *)__pyx_t_4));
     __pyx_t_4 = 0;
 
-    /* "slopescreening/screening/gap_test_all.pyx":130
+    /* "slopescreening/screening/gap_test_all.pyx":134
  * 
  *       # 1. Sort in descending order
  *       if index is None:             # <<<<<<<<<<<<<<
@@ -2859,40 +2915,40 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
  */
   }
 
-  /* "slopescreening/screening/gap_test_all.pyx":134
+  /* "slopescreening/screening/gap_test_all.pyx":138
  * 
  *       # 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, 134, __pyx_L1_error)
+  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 138, __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, 134, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_cumsum); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 138, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
   __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
 
-  /* "slopescreening/screening/gap_test_all.pyx":135
+  /* "slopescreening/screening/gap_test_all.pyx":139
  *       # 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, 135, __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, 139, __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, 135, __pyx_L1_error)
+  __pyx_t_3 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_index), __pyx_slice_); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 139, __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, 135, __pyx_L1_error)
+  __pyx_t_5 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_Atcabs), __pyx_t_3); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 139, __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, 135, __pyx_L1_error)
+  __pyx_t_3 = PyNumber_Subtract(__pyx_t_2, __pyx_t_5); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 139, __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, 135, __pyx_L1_error)
+  __pyx_t_5 = PyFloat_FromDouble(__pyx_v_radius); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 139, __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, 135, __pyx_L1_error)
+  __pyx_t_2 = PyNumber_Subtract(__pyx_t_3, __pyx_t_5); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 139, __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;
@@ -2909,27 +2965,27 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
   __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, 134, __pyx_L1_error)
+  if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 138, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_4);
   __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
 
-  /* "slopescreening/screening/gap_test_all.pyx":136
+  /* "slopescreening/screening/gap_test_all.pyx":140
  *       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, 136, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyObject_GetItem(__pyx_t_4, __pyx_slice_); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 140, __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, 136, __pyx_L1_error)
+  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 140, __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, 134, __pyx_L1_error)
+      __PYX_ERR(0, 138, __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];
     }
   }
@@ -2937,7 +2993,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
   __pyx_v_vec_f = ((PyArrayObject *)__pyx_t_1);
   __pyx_t_1 = 0;
 
-  /* "slopescreening/screening/gap_test_all.pyx":138
+  /* "slopescreening/screening/gap_test_all.pyx":142
  *       )[::-1]
  * 
  *       cdef double best_bound_q = vec_f[0] - coeff_lbd_gamma[0]             # <<<<<<<<<<<<<<
@@ -2950,7 +3006,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
   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)));
 
-  /* "slopescreening/screening/gap_test_all.pyx":139
+  /* "slopescreening/screening/gap_test_all.pyx":143
  * 
  *       cdef double best_bound_q = vec_f[0] - coeff_lbd_gamma[0]
  *       cdef double curr_bound_q = 0.             # <<<<<<<<<<<<<<
@@ -2959,7 +3015,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
  */
   __pyx_v_curr_bound_q = 0.;
 
-  /* "slopescreening/screening/gap_test_all.pyx":140
+  /* "slopescreening/screening/gap_test_all.pyx":144
  *       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]             # <<<<<<<<<<<<<<
@@ -2970,7 +3026,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
   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));
 
-  /* "slopescreening/screening/gap_test_all.pyx":141
+  /* "slopescreening/screening/gap_test_all.pyx":145
  *       cdef double curr_bound_q = 0.
  *       cdef double best_bound_p = vec_f[0]
  *       cdef double curr_bound_p = 0.             # <<<<<<<<<<<<<<
@@ -2979,7 +3035,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
  */
   __pyx_v_curr_bound_p = 0.;
 
-  /* "slopescreening/screening/gap_test_all.pyx":142
+  /* "slopescreening/screening/gap_test_all.pyx":146
  *       cdef double best_bound_p = vec_f[0]
  *       cdef double curr_bound_p = 0.
  *       for k in range(2, n+1):             # <<<<<<<<<<<<<<
@@ -2991,7 +3047,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
   for (__pyx_t_14 = 2; __pyx_t_14 < __pyx_t_22; __pyx_t_14+=1) {
     __pyx_v_k = __pyx_t_14;
 
-    /* "slopescreening/screening/gap_test_all.pyx":146
+    /* "slopescreening/screening/gap_test_all.pyx":150
  * 
  *          # 1. Evaluate p*
  *          curr_bound_p = vec_f[k-1]             # <<<<<<<<<<<<<<
@@ -3002,7 +3058,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
     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));
 
-    /* "slopescreening/screening/gap_test_all.pyx":147
+    /* "slopescreening/screening/gap_test_all.pyx":151
  *          # 1. Evaluate p*
  *          curr_bound_p = vec_f[k-1]
  *          if curr_bound_p > best_bound_p:             # <<<<<<<<<<<<<<
@@ -3012,7 +3068,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
     __pyx_t_12 = ((__pyx_v_curr_bound_p > __pyx_v_best_bound_p) != 0);
     if (__pyx_t_12) {
 
-      /* "slopescreening/screening/gap_test_all.pyx":148
+      /* "slopescreening/screening/gap_test_all.pyx":152
  *          curr_bound_p = vec_f[k-1]
  *          if curr_bound_p > best_bound_p:
  *             best_bound_p = curr_bound_p             # <<<<<<<<<<<<<<
@@ -3021,7 +3077,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
  */
       __pyx_v_best_bound_p = __pyx_v_curr_bound_p;
 
-      /* "slopescreening/screening/gap_test_all.pyx":149
+      /* "slopescreening/screening/gap_test_all.pyx":153
  *          if curr_bound_p > best_bound_p:
  *             best_bound_p = curr_bound_p
  *             vec_p_star[k] = k             # <<<<<<<<<<<<<<
@@ -3032,7 +3088,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
       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;
 
-      /* "slopescreening/screening/gap_test_all.pyx":147
+      /* "slopescreening/screening/gap_test_all.pyx":151
  *          # 1. Evaluate p*
  *          curr_bound_p = vec_f[k-1]
  *          if curr_bound_p > best_bound_p:             # <<<<<<<<<<<<<<
@@ -3042,7 +3098,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
       goto __pyx_L6;
     }
 
-    /* "slopescreening/screening/gap_test_all.pyx":151
+    /* "slopescreening/screening/gap_test_all.pyx":155
  *             vec_p_star[k] = k
  *          else:
  *             vec_p_star[k] = vec_p_star[k-1]             # <<<<<<<<<<<<<<
@@ -3058,7 +3114,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
     }
     __pyx_L6:;
 
-    /* "slopescreening/screening/gap_test_all.pyx":154
+    /* "slopescreening/screening/gap_test_all.pyx":158
  * 
  *          # 1. Evaluate q*
  *          curr_bound_q = vec_f[k-1] - coeff_lbd_gamma[k-1]             # <<<<<<<<<<<<<<
@@ -3071,7 +3127,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
     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)));
 
-    /* "slopescreening/screening/gap_test_all.pyx":155
+    /* "slopescreening/screening/gap_test_all.pyx":159
  *          # 1. Evaluate q*
  *          curr_bound_q = vec_f[k-1] - coeff_lbd_gamma[k-1]
  *          if curr_bound_q > best_bound_q:             # <<<<<<<<<<<<<<
@@ -3081,7 +3137,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
     __pyx_t_12 = ((__pyx_v_curr_bound_q > __pyx_v_best_bound_q) != 0);
     if (__pyx_t_12) {
 
-      /* "slopescreening/screening/gap_test_all.pyx":156
+      /* "slopescreening/screening/gap_test_all.pyx":160
  *          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             # <<<<<<<<<<<<<<
@@ -3090,7 +3146,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
  */
       __pyx_v_best_bound_q = __pyx_v_curr_bound_q;
 
-      /* "slopescreening/screening/gap_test_all.pyx":157
+      /* "slopescreening/screening/gap_test_all.pyx":161
  *          if curr_bound_q > best_bound_q:
  *             best_bound_q = curr_bound_q
  *             vec_q_star[k] = k             # <<<<<<<<<<<<<<
@@ -3101,7 +3157,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
       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;
 
-      /* "slopescreening/screening/gap_test_all.pyx":155
+      /* "slopescreening/screening/gap_test_all.pyx":159
  *          # 1. Evaluate q*
  *          curr_bound_q = vec_f[k-1] - coeff_lbd_gamma[k-1]
  *          if curr_bound_q > best_bound_q:             # <<<<<<<<<<<<<<
@@ -3111,7 +3167,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
       goto __pyx_L7;
     }
 
-    /* "slopescreening/screening/gap_test_all.pyx":159
+    /* "slopescreening/screening/gap_test_all.pyx":163
  *             vec_q_star[k] = k
  *          else:
  *             vec_q_star[k] = vec_q_star[k-1]             # <<<<<<<<<<<<<<
@@ -3128,16 +3184,16 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
     __pyx_L7:;
   }
 
-  /* "slopescreening/screening/gap_test_all.pyx":162
+  /* "slopescreening/screening/gap_test_all.pyx":166
  * 
  *       # 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, 162, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyInt_From_int(__pyx_v_n); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 166, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
-  __pyx_t_4 = PyTuple_New(3); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 162, __pyx_L1_error)
+  __pyx_t_4 = PyTuple_New(3); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 166, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_4);
   __Pyx_GIVEREF(__pyx_t_1);
   PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_1);
@@ -3148,16 +3204,16 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
   __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, 162, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyObject_Call(__pyx_builtin_range, __pyx_t_4, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 166, __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, 162, __pyx_L1_error)
+    __pyx_t_23 = -1; __pyx_t_4 = PyObject_GetIter(__pyx_t_1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 166, __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, 162, __pyx_L1_error)
+    __pyx_t_24 = Py_TYPE(__pyx_t_4)->tp_iternext; if (unlikely(!__pyx_t_24)) __PYX_ERR(0, 166, __pyx_L1_error)
   }
   __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
   for (;;) {
@@ -3165,17 +3221,17 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
       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, 162, __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, 166, __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, 162, __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, 166, __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, 162, __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, 166, __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, 162, __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, 166, __pyx_L1_error)
         __Pyx_GOTREF(__pyx_t_1);
         #endif
       }
@@ -3185,7 +3241,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
         PyObject* exc_type = PyErr_Occurred();
         if (exc_type) {
           if (likely(__Pyx_PyErr_GivenExceptionMatches(exc_type, PyExc_StopIteration))) PyErr_Clear();
-          else __PYX_ERR(0, 162, __pyx_L1_error)
+          else __PYX_ERR(0, 166, __pyx_L1_error)
         }
         break;
       }
@@ -3194,20 +3250,20 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
     __Pyx_XDECREF_SET(__pyx_v_l, __pyx_t_1);
     __pyx_t_1 = 0;
 
-    /* "slopescreening/screening/gap_test_all.pyx":163
+    /* "slopescreening/screening/gap_test_all.pyx":167
  *       # 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, 163, __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, 167, __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, 163, __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, 167, __pyx_L1_error)
     __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
     __pyx_v_q = __pyx_t_14;
 
-    /* "slopescreening/screening/gap_test_all.pyx":164
+    /* "slopescreening/screening/gap_test_all.pyx":168
  *       for l in range(n, 0, -1):
  *          q = vec_q_star[l]
  *          while q >= 1:             # <<<<<<<<<<<<<<
@@ -3218,7 +3274,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
       __pyx_t_12 = ((__pyx_v_q >= 1) != 0);
       if (!__pyx_t_12) break;
 
-      /* "slopescreening/screening/gap_test_all.pyx":166
+      /* "slopescreening/screening/gap_test_all.pyx":170
  *          while q >= 1:
  *             # Implicit definition of the best value of p given q
  *             p = vec_p_star[q]             # <<<<<<<<<<<<<<
@@ -3227,83 +3283,83 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
  */
       __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, 166, __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, 170, __pyx_L1_error)
       __Pyx_GOTREF(__pyx_t_1);
       __Pyx_XDECREF_SET(__pyx_v_p, __pyx_t_1);
       __pyx_t_1 = 0;
 
-      /* "slopescreening/screening/gap_test_all.pyx":169
+      /* "slopescreening/screening/gap_test_all.pyx":173
  * 
  *             # 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, 169, __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, 173, __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, 169, __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, 173, __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, 169, __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, 173, __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, 169, __pyx_L1_error)
+      __pyx_t_5 = PyNumber_Subtract(__pyx_t_2, __pyx_t_1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 173, __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, 169, __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, 173, __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, 169, __pyx_L1_error)
+      __pyx_t_2 = PyNumber_Add(__pyx_t_5, __pyx_t_1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 173, __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, 169, __pyx_L1_error)
+      __pyx_t_6 = __pyx_PyFloat_AsDouble(__pyx_t_2); if (unlikely((__pyx_t_6 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 173, __pyx_L1_error)
       __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
       __pyx_v_tau = __pyx_t_6;
 
-      /* "slopescreening/screening/gap_test_all.pyx":172
+      /* "slopescreening/screening/gap_test_all.pyx":176
  * 
  *             # 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, 172, __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, 176, __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, 172, __pyx_L1_error)
+      __pyx_t_1 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_index), __pyx_t_2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 176, __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, 172, __pyx_L1_error)
+      __pyx_t_2 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_Atcabs), __pyx_t_1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 176, __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, 172, __pyx_L1_error)
+      __pyx_t_1 = PyFloat_FromDouble(__pyx_v_tau); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 176, __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, 172, __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, 176, __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, 172, __pyx_L1_error)
+      __pyx_t_12 = __Pyx_PyObject_IsTrue(__pyx_t_5); if (unlikely(__pyx_t_12 < 0)) __PYX_ERR(0, 176, __pyx_L1_error)
       __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
       if (__pyx_t_12) {
 
-        /* "slopescreening/screening/gap_test_all.pyx":173
+        /* "slopescreening/screening/gap_test_all.pyx":177
  *             # 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, 173, __pyx_L1_error)
+        __pyx_t_5 = PySlice_New(__pyx_v_l, Py_None, Py_None); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 177, __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, 173, __pyx_L1_error)
+        __pyx_t_1 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_index), __pyx_t_5); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 177, __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, 173, __pyx_L1_error)
+        if (unlikely(PyObject_SetItem(((PyObject *)__pyx_v_calI_screen), __pyx_t_1, Py_True) < 0)) __PYX_ERR(0, 177, __pyx_L1_error)
         __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
 
-        /* "slopescreening/screening/gap_test_all.pyx":174
+        /* "slopescreening/screening/gap_test_all.pyx":178
  *             if Atcabs[index[l-1]] >= tau:
  *                calI_screen[index[l:]] = True
  *                return calI_screen             # <<<<<<<<<<<<<<
@@ -3316,7 +3372,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
         __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
         goto __pyx_L0;
 
-        /* "slopescreening/screening/gap_test_all.pyx":172
+        /* "slopescreening/screening/gap_test_all.pyx":176
  * 
  *             # Test
  *             if Atcabs[index[l-1]] >= tau:             # <<<<<<<<<<<<<<
@@ -3325,24 +3381,24 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
  */
       }
 
-      /* "slopescreening/screening/gap_test_all.pyx":177
+      /* "slopescreening/screening/gap_test_all.pyx":181
  * 
  *             # 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, 177, __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, 181, __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, 177, __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, 181, __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, 177, __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, 181, __pyx_L1_error)
       __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
       __pyx_v_q = __pyx_t_14;
     }
 
-    /* "slopescreening/screening/gap_test_all.pyx":162
+    /* "slopescreening/screening/gap_test_all.pyx":166
  * 
  *       # 3. Tests
  *       for l in range(n, 0, -1):             # <<<<<<<<<<<<<<
@@ -3352,22 +3408,22 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
   }
   __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
 
-  /* "slopescreening/screening/gap_test_all.pyx":180
+  /* "slopescreening/screening/gap_test_all.pyx":184
  * 
  * 
  *       calI_screen[index[l:]] = True             # <<<<<<<<<<<<<<
  *       return calI_screen
  */
-  if (unlikely(!__pyx_v_l)) { __Pyx_RaiseUnboundLocalError("l"); __PYX_ERR(0, 180, __pyx_L1_error) }
-  __pyx_t_4 = PySlice_New(__pyx_v_l, Py_None, Py_None); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 180, __pyx_L1_error)
+  if (unlikely(!__pyx_v_l)) { __Pyx_RaiseUnboundLocalError("l"); __PYX_ERR(0, 184, __pyx_L1_error) }
+  __pyx_t_4 = PySlice_New(__pyx_v_l, Py_None, Py_None); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 184, __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, 180, __pyx_L1_error)
+  __pyx_t_5 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_index), __pyx_t_4); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 184, __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, 180, __pyx_L1_error)
+  if (unlikely(PyObject_SetItem(((PyObject *)__pyx_v_calI_screen), __pyx_t_5, Py_True) < 0)) __PYX_ERR(0, 184, __pyx_L1_error)
   __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
 
-  /* "slopescreening/screening/gap_test_all.pyx":181
+  /* "slopescreening/screening/gap_test_all.pyx":185
  * 
  *       calI_screen[index[l:]] = True
  *       return calI_screen             # <<<<<<<<<<<<<<
@@ -3377,7 +3433,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
   __pyx_r = ((PyObject *)__pyx_v_calI_screen);
   goto __pyx_L0;
 
-  /* "slopescreening/screening/gap_test_all.pyx":64
+  /* "slopescreening/screening/gap_test_all.pyx":68
  *    @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,             # <<<<<<<<<<<<<<
@@ -3438,19 +3494,19 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
  */
 
 /* Python wrapper */
-static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_9__reduce_cython__(PyObject *__pyx_v_self, CYTHON_UNUSED PyObject *unused); /*proto*/
-static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_9__reduce_cython__(PyObject *__pyx_v_self, CYTHON_UNUSED PyObject *unused) {
+static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_11__reduce_cython__(PyObject *__pyx_v_self, CYTHON_UNUSED PyObject *unused); /*proto*/
+static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_11__reduce_cython__(PyObject *__pyx_v_self, CYTHON_UNUSED PyObject *unused) {
   PyObject *__pyx_r = 0;
   __Pyx_RefNannyDeclarations
   __Pyx_RefNannySetupContext("__reduce_cython__ (wrapper)", 0);
-  __pyx_r = __pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_8__reduce_cython__(((struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *)__pyx_v_self));
+  __pyx_r = __pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_10__reduce_cython__(((struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *)__pyx_v_self));
 
   /* function exit code */
   __Pyx_RefNannyFinishContext();
   return __pyx_r;
 }
 
-static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_8__reduce_cython__(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *__pyx_v_self) {
+static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_10__reduce_cython__(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *__pyx_v_self) {
   PyObject *__pyx_r = NULL;
   __Pyx_RefNannyDeclarations
   PyObject *__pyx_t_1 = NULL;
@@ -3495,19 +3551,19 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAl
  */
 
 /* Python wrapper */
-static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_11__setstate_cython__(PyObject *__pyx_v_self, PyObject *__pyx_v___pyx_state); /*proto*/
-static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_11__setstate_cython__(PyObject *__pyx_v_self, PyObject *__pyx_v___pyx_state) {
+static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_13__setstate_cython__(PyObject *__pyx_v_self, PyObject *__pyx_v___pyx_state); /*proto*/
+static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_13__setstate_cython__(PyObject *__pyx_v_self, PyObject *__pyx_v___pyx_state) {
   PyObject *__pyx_r = 0;
   __Pyx_RefNannyDeclarations
   __Pyx_RefNannySetupContext("__setstate_cython__ (wrapper)", 0);
-  __pyx_r = __pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_10__setstate_cython__(((struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *)__pyx_v_self), ((PyObject *)__pyx_v___pyx_state));
+  __pyx_r = __pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_12__setstate_cython__(((struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *)__pyx_v_self), ((PyObject *)__pyx_v___pyx_state));
 
   /* function exit code */
   __Pyx_RefNannyFinishContext();
   return __pyx_r;
 }
 
-static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_10__setstate_cython__(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *__pyx_v_self, CYTHON_UNUSED PyObject *__pyx_v___pyx_state) {
+static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_all_10GapTestAll_12__setstate_cython__(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_all_GapTestAll *__pyx_v_self, CYTHON_UNUSED PyObject *__pyx_v___pyx_state) {
   PyObject *__pyx_r = NULL;
   __Pyx_RefNannyDeclarations
   PyObject *__pyx_t_1 = NULL;
@@ -4591,9 +4647,10 @@ static void __pyx_tp_dealloc_14slopescreening_9screening_12gap_test_all_GapTestA
 
 static PyMethodDef __pyx_methods_14slopescreening_9screening_12gap_test_all_GapTestAll[] = {
   {"get_name", (PyCFunction)__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_5get_name, METH_NOARGS, 0},
-  {"apply_test", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_7apply_test, METH_VARARGS|METH_KEYWORDS, __pyx_doc_14slopescreening_9screening_12gap_test_all_10GapTestAll_6apply_test},
-  {"__reduce_cython__", (PyCFunction)__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_9__reduce_cython__, METH_NOARGS, 0},
-  {"__setstate_cython__", (PyCFunction)__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_11__setstate_cython__, METH_O, 0},
+  {"get_legend_name", (PyCFunction)__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_7get_legend_name, METH_NOARGS, 0},
+  {"apply_test", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_9apply_test, METH_VARARGS|METH_KEYWORDS, __pyx_doc_14slopescreening_9screening_12gap_test_all_10GapTestAll_8apply_test},
+  {"__reduce_cython__", (PyCFunction)__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_11__reduce_cython__, METH_NOARGS, 0},
+  {"__setstate_cython__", (PyCFunction)__pyx_pw_14slopescreening_9screening_12gap_test_all_10GapTestAll_13__setstate_cython__, METH_O, 0},
   {0, 0, 0, 0}
 };
 
@@ -4720,6 +4777,7 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = {
   {&__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},
+  {&__pyx_n_s_all, __pyx_k_all, sizeof(__pyx_k_all), 0, 0, 1, 1},
   {&__pyx_n_s_arange, __pyx_k_arange, sizeof(__pyx_k_arange), 0, 0, 1, 1},
   {&__pyx_n_s_argsort, __pyx_k_argsort, sizeof(__pyx_k_argsort), 0, 0, 1, 1},
   {&__pyx_n_s_cline_in_traceback, __pyx_k_cline_in_traceback, sizeof(__pyx_k_cline_in_traceback), 0, 0, 1, 1},
@@ -4765,14 +4823,14 @@ static CYTHON_SMALL_CODE int __Pyx_InitCachedConstants(void) {
   __Pyx_RefNannyDeclarations
   __Pyx_RefNannySetupContext("__Pyx_InitCachedConstants", 0);
 
-  /* "slopescreening/screening/gap_test_all.pyx":131
+  /* "slopescreening/screening/gap_test_all.pyx":135
  *       # 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, 131, __pyx_L1_error)
+  __pyx_slice_ = PySlice_New(Py_None, Py_None, __pyx_int_neg_1); if (unlikely(!__pyx_slice_)) __PYX_ERR(0, 135, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_slice_);
   __Pyx_GIVEREF(__pyx_slice_);
 
diff --git a/slopescreening/screening/gap_test_all.pyx b/slopescreening/screening/gap_test_all.pyx
index b9bdf97..457b339 100644
--- a/slopescreening/screening/gap_test_all.pyx
+++ b/slopescreening/screening/gap_test_all.pyx
@@ -58,6 +58,10 @@ cdef class GapTestAll:
    def get_name(self):
       return self.name
 
+   @cython.boundscheck(False) # turn off bounds-checking for entire function
+   @cython.wraparound(False)  # turn off negative index wrapping for entire function
+   def get_legend_name(self):
+      return "all"
 
    @cython.boundscheck(False) # turn off bounds-checking for entire function
    # @cython.wraparound(False)  # turn off negative index wrapping for entire function
diff --git a/slopescreening/screening/gap_test_p_1.c b/slopescreening/screening/gap_test_p_1.c
index e757c50..78887c1 100644
--- a/slopescreening/screening/gap_test_p_1.c
+++ b/slopescreening/screening/gap_test_p_1.c
@@ -1936,6 +1936,7 @@ static const char __pyx_k_vec_gammas[] = "vec_gammas";
 static const char __pyx_k_ImportError[] = "ImportError";
 static const char __pyx_k_offset_radius[] = "offset_radius";
 static const char __pyx_k_reduce_cython[] = "__reduce_cython__";
+static const char __pyx_k_p_q_1_forall_q[] = "$p_q=1\\;\\forall q$";
 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";
@@ -1969,6 +1970,7 @@ static PyObject *__pyx_n_s_numpy;
 static PyObject *__pyx_kp_s_numpy_core_multiarray_failed_to;
 static PyObject *__pyx_kp_s_numpy_core_umath_failed_to_impor;
 static PyObject *__pyx_n_s_offset_radius;
+static PyObject *__pyx_kp_s_p_q_1_forall_q;
 static PyObject *__pyx_n_s_range;
 static PyObject *__pyx_n_s_reduce;
 static PyObject *__pyx_n_s_reduce_cython;
@@ -1985,9 +1987,10 @@ static PyObject *__pyx_n_s_zeros;
 static int __pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne___init__(struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *__pyx_v_self, PyArrayObject *__pyx_v_vec_gammas, PyArrayObject *__pyx_v_vec_coherence_function); /* proto */
 static void __pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_2__dealloc__(struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *__pyx_v_self); /* proto */
 static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_4get_name(struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *__pyx_v_self); /* proto */
-static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_6apply_test(struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *__pyx_v_self, PyArrayObject *__pyx_v_Atcabs, double __pyx_v_gap, double __pyx_v_lbd, CYTHON_UNUSED PyArrayObject *__pyx_v_vec_gammas, double __pyx_v_coeff_dual_scaling, double __pyx_v_offset_radius, PyArrayObject *__pyx_v_index); /* proto */
-static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_8__reduce_cython__(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *__pyx_v_self); /* proto */
-static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_10__setstate_cython__(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *__pyx_v_self, CYTHON_UNUSED PyObject *__pyx_v___pyx_state); /* proto */
+static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_6get_legend_name(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *__pyx_v_self); /* proto */
+static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_8apply_test(struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *__pyx_v_self, PyArrayObject *__pyx_v_Atcabs, double __pyx_v_gap, double __pyx_v_lbd, CYTHON_UNUSED PyArrayObject *__pyx_v_vec_gammas, double __pyx_v_coeff_dual_scaling, double __pyx_v_offset_radius, PyArrayObject *__pyx_v_index); /* proto */
+static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_10__reduce_cython__(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *__pyx_v_self); /* proto */
+static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_12__setstate_cython__(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *__pyx_v_self, CYTHON_UNUSED PyObject *__pyx_v___pyx_state); /* proto */
 static PyObject *__pyx_tp_new_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne(PyTypeObject *t, PyObject *a, PyObject *k); /*proto*/
 static PyObject *__pyx_int_1;
 static PyObject *__pyx_int_neg_1;
@@ -2382,7 +2385,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
  *    def get_name(self):
  *       return self.name             # <<<<<<<<<<<<<<
  * 
- * 
+ *    @cython.boundscheck(False) # turn off bounds-checking for entire function
  */
   __Pyx_XDECREF(__pyx_r);
   __pyx_t_1 = __Pyx_PyBytes_FromString(__pyx_v_self->name); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 87, __pyx_L1_error)
@@ -2410,7 +2413,60 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
   return __pyx_r;
 }
 
-/* "slopescreening/screening/gap_test_p_1.pyx":92
+/* "slopescreening/screening/gap_test_p_1.pyx":91
+ *    @cython.boundscheck(False) # turn off bounds-checking for entire function
+ *    @cython.wraparound(False)  # turn off negative index wrapping for entire function
+ *    def get_legend_name(self):             # <<<<<<<<<<<<<<
+ *       return "$p_q=1\\;\\forall q$"
+ * 
+ */
+
+/* Python wrapper */
+static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_7get_legend_name(PyObject *__pyx_v_self, CYTHON_UNUSED PyObject *unused); /*proto*/
+static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_7get_legend_name(PyObject *__pyx_v_self, CYTHON_UNUSED PyObject *unused) {
+  PyObject *__pyx_r = 0;
+  __Pyx_RefNannyDeclarations
+  __Pyx_RefNannySetupContext("get_legend_name (wrapper)", 0);
+  __pyx_r = __pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_6get_legend_name(((struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *)__pyx_v_self));
+
+  /* function exit code */
+  __Pyx_RefNannyFinishContext();
+  return __pyx_r;
+}
+
+static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_6get_legend_name(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *__pyx_v_self) {
+  PyObject *__pyx_r = NULL;
+  __Pyx_RefNannyDeclarations
+  __Pyx_RefNannySetupContext("get_legend_name", 0);
+
+  /* "slopescreening/screening/gap_test_p_1.pyx":92
+ *    @cython.wraparound(False)  # turn off negative index wrapping for entire function
+ *    def get_legend_name(self):
+ *       return "$p_q=1\\;\\forall q$"             # <<<<<<<<<<<<<<
+ * 
+ *    @cython.boundscheck(False) # turn off bounds-checking for entire function
+ */
+  __Pyx_XDECREF(__pyx_r);
+  __Pyx_INCREF(__pyx_kp_s_p_q_1_forall_q);
+  __pyx_r = __pyx_kp_s_p_q_1_forall_q;
+  goto __pyx_L0;
+
+  /* "slopescreening/screening/gap_test_p_1.pyx":91
+ *    @cython.boundscheck(False) # turn off bounds-checking for entire function
+ *    @cython.wraparound(False)  # turn off negative index wrapping for entire function
+ *    def get_legend_name(self):             # <<<<<<<<<<<<<<
+ *       return "$p_q=1\\;\\forall q$"
+ * 
+ */
+
+  /* function exit code */
+  __pyx_L0:;
+  __Pyx_XGIVEREF(__pyx_r);
+  __Pyx_RefNannyFinishContext();
+  return __pyx_r;
+}
+
+/* "slopescreening/screening/gap_test_p_1.pyx":96
  *    @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,             # <<<<<<<<<<<<<<
@@ -2419,9 +2475,9 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
  */
 
 /* Python wrapper */
-static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_7apply_test(PyObject *__pyx_v_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
-static char __pyx_doc_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_6apply_test[] = " Apply the Gap safe sphere screening test\n\n      Implementation of the safe screening test\n      \n      Parameters\n      ----------\n      Atcabs : np.ndarray\n         vector matA.T @ vecc where vec is dual admissible in absolute value\n         size [n]\n      gap : positive float\n         duality gap\n      lbd : positive float\n         regularization parameters\n      vec_gammas : np.ndarray\n         slope parameters\n         size [n,]\n      vec_gammas : np.ndarray\n         slope parameters\n         size [n,]\n      coeff_dual_scaling : positif float\n         If coeff_dual_scaling is not feasible, dual scaling factor\n         such taht vecu / coeff_dual_scaling os dual feasible\n         Here for code optimization purposes\n         Default value is 1. (vecu is feasible)\n      offset_radius : float\n         additive term added to the redius\n         default is 0\n      index : np.ndarray\n         Array of indices that sort Atu in absolute value\n         default is None\n      \n      Returns\n      -------\n      calI_screen : np.ndarray\n         vector of boolean\n         True if screening test passes and False otherwise\n         size [n,]\n      ";
-static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_7apply_test(PyObject *__pyx_v_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
+static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_9apply_test(PyObject *__pyx_v_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
+static char __pyx_doc_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_8apply_test[] = " Apply the Gap safe sphere screening test\n\n      Implementation of the safe screening test\n      \n      Parameters\n      ----------\n      Atcabs : np.ndarray\n         vector matA.T @ vecc where vec is dual admissible in absolute value\n         size [n]\n      gap : positive float\n         duality gap\n      lbd : positive float\n         regularization parameters\n      vec_gammas : np.ndarray\n         slope parameters\n         size [n,]\n      vec_gammas : np.ndarray\n         slope parameters\n         size [n,]\n      coeff_dual_scaling : positif float\n         If coeff_dual_scaling is not feasible, dual scaling factor\n         such taht vecu / coeff_dual_scaling os dual feasible\n         Here for code optimization purposes\n         Default value is 1. (vecu is feasible)\n      offset_radius : float\n         additive term added to the redius\n         default is 0\n      index : np.ndarray\n         Array of indices that sort Atu in absolute value\n         default is None\n      \n      Returns\n      -------\n      calI_screen : np.ndarray\n         vector of boolean\n         True if screening test passes and False otherwise\n         size [n,]\n      ";
+static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_9apply_test(PyObject *__pyx_v_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
   PyArrayObject *__pyx_v_Atcabs = 0;
   double __pyx_v_gap;
   double __pyx_v_lbd;
@@ -2439,7 +2495,7 @@ static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
     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};
 
-    /* "slopescreening/screening/gap_test_p_1.pyx":99
+    /* "slopescreening/screening/gap_test_p_1.pyx":103
  *       double coeff_dual_scaling = 1.,
  *       double offset_radius=0.,
  *       np.ndarray[long, ndim=1] index=None             # <<<<<<<<<<<<<<
@@ -2477,19 +2533,19 @@ static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
         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, 92, __pyx_L3_error)
+          __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, 1); __PYX_ERR(0, 96, __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, 92, __pyx_L3_error)
+          __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, 2); __PYX_ERR(0, 96, __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, 92, __pyx_L3_error)
+          __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, 3); __PYX_ERR(0, 96, __pyx_L3_error)
         }
         CYTHON_FALLTHROUGH;
         case  4:
@@ -2511,7 +2567,7 @@ static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
         }
       }
       if (unlikely(kw_args > 0)) {
-        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "apply_test") < 0)) __PYX_ERR(0, 92, __pyx_L3_error)
+        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "apply_test") < 0)) __PYX_ERR(0, 96, __pyx_L3_error)
       }
     } else {
       switch (PyTuple_GET_SIZE(__pyx_args)) {
@@ -2530,16 +2586,16 @@ static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
       }
     }
     __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, 94, __pyx_L3_error)
-    __pyx_v_lbd = __pyx_PyFloat_AsDouble(values[2]); if (unlikely((__pyx_v_lbd == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 95, __pyx_L3_error)
+    __pyx_v_gap = __pyx_PyFloat_AsDouble(values[1]); if (unlikely((__pyx_v_gap == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 98, __pyx_L3_error)
+    __pyx_v_lbd = __pyx_PyFloat_AsDouble(values[2]); if (unlikely((__pyx_v_lbd == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 99, __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, 97, __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, 101, __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, 98, __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, 102, __pyx_L3_error)
     } else {
       __pyx_v_offset_radius = ((double)0.);
     }
@@ -2547,18 +2603,18 @@ static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
   }
   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, 92, __pyx_L3_error)
+  __Pyx_RaiseArgtupleInvalid("apply_test", 0, 4, 7, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 96, __pyx_L3_error)
   __pyx_L3_error:;
   __Pyx_AddTraceback("slopescreening.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, 93, __pyx_L1_error)
-  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_vec_gammas), __pyx_ptype_5numpy_ndarray, 1, "vec_gammas", 0))) __PYX_ERR(0, 96, __pyx_L1_error)
-  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_index), __pyx_ptype_5numpy_ndarray, 1, "index", 0))) __PYX_ERR(0, 99, __pyx_L1_error)
-  __pyx_r = __pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_6apply_test(((struct __pyx_obj_14slopescreening_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);
+  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_Atcabs), __pyx_ptype_5numpy_ndarray, 1, "Atcabs", 0))) __PYX_ERR(0, 97, __pyx_L1_error)
+  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_vec_gammas), __pyx_ptype_5numpy_ndarray, 1, "vec_gammas", 0))) __PYX_ERR(0, 100, __pyx_L1_error)
+  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_index), __pyx_ptype_5numpy_ndarray, 1, "index", 0))) __PYX_ERR(0, 103, __pyx_L1_error)
+  __pyx_r = __pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_8apply_test(((struct __pyx_obj_14slopescreening_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);
 
-  /* "slopescreening/screening/gap_test_p_1.pyx":92
+  /* "slopescreening/screening/gap_test_p_1.pyx":96
  *    @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,             # <<<<<<<<<<<<<<
@@ -2575,7 +2631,7 @@ static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
   return __pyx_r;
 }
 
-static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_6apply_test(struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *__pyx_v_self, PyArrayObject *__pyx_v_Atcabs, double __pyx_v_gap, double __pyx_v_lbd, CYTHON_UNUSED PyArrayObject *__pyx_v_vec_gammas, double __pyx_v_coeff_dual_scaling, double __pyx_v_offset_radius, PyArrayObject *__pyx_v_index) {
+static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_8apply_test(struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *__pyx_v_self, PyArrayObject *__pyx_v_Atcabs, double __pyx_v_gap, double __pyx_v_lbd, CYTHON_UNUSED PyArrayObject *__pyx_v_vec_gammas, double __pyx_v_coeff_dual_scaling, double __pyx_v_offset_radius, PyArrayObject *__pyx_v_index) {
   int __pyx_v_n;
   int __pyx_v_q;
   int __pyx_v_l_q;
@@ -2645,21 +2701,21 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
   __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, 92, __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, 96, __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, 92, __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, 96, __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, 92, __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, 96, __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];
 
-  /* "slopescreening/screening/gap_test_p_1.pyx":140
+  /* "slopescreening/screening/gap_test_p_1.pyx":144
  *       """
  * 
  *       cdef int n = Atcabs.shape[0]             # <<<<<<<<<<<<<<
@@ -2668,7 +2724,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
  */
   __pyx_v_n = (__pyx_v_Atcabs->dimensions[0]);
 
-  /* "slopescreening/screening/gap_test_p_1.pyx":142
+  /* "slopescreening/screening/gap_test_p_1.pyx":146
  *       cdef int n = Atcabs.shape[0]
  *       cdef int q, l_q
  *       cdef int l_min = 0             # <<<<<<<<<<<<<<
@@ -2677,21 +2733,21 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
  */
   __pyx_v_l_min = 0;
 
-  /* "slopescreening/screening/gap_test_p_1.pyx":145
+  /* "slopescreening/screening/gap_test_p_1.pyx":149
  * 
  *       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, 145, __pyx_L1_error)
+  __pyx_t_1 = PyFloat_FromDouble(__pyx_v_coeff_dual_scaling); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 149, __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, 145, __pyx_L1_error)
+  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 149, __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, 145, __pyx_L1_error)
+  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_sqrt); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 149, __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, 145, __pyx_L1_error)
+  __pyx_t_3 = PyFloat_FromDouble((2.0 * __pyx_v_gap)); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 149, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_3);
   __pyx_t_5 = NULL;
   if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_4))) {
@@ -2706,24 +2762,24 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
   __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, 145, __pyx_L1_error)
+  if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 149, __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, 145, __pyx_L1_error)
+  __pyx_t_4 = PyNumber_Multiply(__pyx_t_1, __pyx_t_2); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 149, __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, 145, __pyx_L1_error)
+  __pyx_t_2 = PyFloat_FromDouble(__pyx_v_offset_radius); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 149, __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, 145, __pyx_L1_error)
+  __pyx_t_1 = PyNumber_Add(__pyx_t_4, __pyx_t_2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 149, __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, 145, __pyx_L1_error)
+  __pyx_t_6 = __pyx_PyFloat_AsDouble(__pyx_t_1); if (unlikely((__pyx_t_6 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 149, __pyx_L1_error)
   __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
   __pyx_v_radius = __pyx_t_6;
 
-  /* "slopescreening/screening/gap_test_p_1.pyx":146
+  /* "slopescreening/screening/gap_test_p_1.pyx":150
  *       cdef double bound
  *       cdef double radius = coeff_dual_scaling * np.sqrt(2 * gap) + offset_radius
  *       cdef double coeff_lbd = coeff_dual_scaling * lbd             # <<<<<<<<<<<<<<
@@ -2732,19 +2788,19 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
  */
   __pyx_v_coeff_lbd = (__pyx_v_coeff_dual_scaling * __pyx_v_lbd);
 
-  /* "slopescreening/screening/gap_test_p_1.pyx":148
+  /* "slopescreening/screening/gap_test_p_1.pyx":152
  *       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, 148, __pyx_L1_error)
+  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 152, __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, 148, __pyx_L1_error)
+  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_zeros); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 152, __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, 148, __pyx_L1_error)
+  __pyx_t_2 = __Pyx_PyInt_From_long((__pyx_v_n + 1)); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 152, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_2);
   __pyx_t_3 = NULL;
   if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_4))) {
@@ -2759,16 +2815,16 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
   __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, 148, __pyx_L1_error)
+  if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 152, __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, 148, __pyx_L1_error)
+  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 152, __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, 148, __pyx_L1_error)
+      __PYX_ERR(0, 152, __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];
     }
   }
@@ -2776,40 +2832,40 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
   __pyx_v_sigmas = ((PyArrayObject *)__pyx_t_1);
   __pyx_t_1 = 0;
 
-  /* "slopescreening/screening/gap_test_p_1.pyx":151
+  /* "slopescreening/screening/gap_test_p_1.pyx":155
  * 
  *       # 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, 151, __pyx_L1_error)
+  __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 155, __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, 151, __pyx_L1_error)
+  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_zeros); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 155, __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, 151, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyInt_From_int(__pyx_v_n); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 155, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_1);
-  __pyx_t_2 = PyTuple_New(1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 151, __pyx_L1_error)
+  __pyx_t_2 = PyTuple_New(1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 155, __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, 151, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 155, __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, 151, __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, 151, __pyx_L1_error)
+  if (PyDict_SetItem(__pyx_t_1, __pyx_n_s_dtype, ((PyObject*)&PyBool_Type)) < 0) __PYX_ERR(0, 155, __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, 155, __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, 151, __pyx_L1_error)
+  if (!(likely(((__pyx_t_3) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_3, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 155, __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, 151, __pyx_L1_error)
+      __PYX_ERR(0, 155, __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];
     }
   }
@@ -2817,7 +2873,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
   __pyx_v_calI_screen = ((PyArrayObject *)__pyx_t_3);
   __pyx_t_3 = 0;
 
-  /* "slopescreening/screening/gap_test_p_1.pyx":154
+  /* "slopescreening/screening/gap_test_p_1.pyx":158
  * 
  *       # 1. Sort in descending order
  *       if index is None:             # <<<<<<<<<<<<<<
@@ -2828,16 +2884,16 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
   __pyx_t_10 = (__pyx_t_9 != 0);
   if (__pyx_t_10) {
 
-    /* "slopescreening/screening/gap_test_p_1.pyx":155
+    /* "slopescreening/screening/gap_test_p_1.pyx":159
  *       # 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, 155, __pyx_L1_error)
+    __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 159, __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, 155, __pyx_L1_error)
+    __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_argsort); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 159, __pyx_L1_error)
     __Pyx_GOTREF(__pyx_t_2);
     __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
     __pyx_t_1 = NULL;
@@ -2852,13 +2908,13 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
     }
     __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, 155, __pyx_L1_error)
+    if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 159, __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, 155, __pyx_L1_error)
+    __pyx_t_2 = __Pyx_PyObject_GetItem(__pyx_t_3, __pyx_slice_); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 159, __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, 155, __pyx_L1_error)
+    if (!(likely(((__pyx_t_2) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_2, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 159, __pyx_L1_error)
     __pyx_t_11 = ((PyArrayObject *)__pyx_t_2);
     {
       __Pyx_BufFmt_StackElem __pyx_stack[1];
@@ -2875,13 +2931,13 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
         __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, 155, __pyx_L1_error)
+      if (unlikely(__pyx_t_12 < 0)) __PYX_ERR(0, 159, __pyx_L1_error)
     }
     __pyx_t_11 = 0;
     __Pyx_DECREF_SET(__pyx_v_index, ((PyArrayObject *)__pyx_t_2));
     __pyx_t_2 = 0;
 
-    /* "slopescreening/screening/gap_test_p_1.pyx":154
+    /* "slopescreening/screening/gap_test_p_1.pyx":158
  * 
  *       # 1. Sort in descending order
  *       if index is None:             # <<<<<<<<<<<<<<
@@ -2890,19 +2946,19 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
  */
   }
 
-  /* "slopescreening/screening/gap_test_p_1.pyx":158
+  /* "slopescreening/screening/gap_test_p_1.pyx":162
  * 
  *          # 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, 158, __pyx_L1_error)
+  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 162, __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, 158, __pyx_L1_error)
+  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_cumsum); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 162, __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, 158, __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, 162, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_t_3);
   __pyx_t_4 = NULL;
   if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_1))) {
@@ -2917,13 +2973,13 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
   __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, 158, __pyx_L1_error)
+  if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 162, __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, 158, __pyx_L1_error)
+  if (unlikely(PyObject_SetItem(((PyObject *)__pyx_v_sigmas), __pyx_slice__2, __pyx_t_2) < 0)) __PYX_ERR(0, 162, __pyx_L1_error)
   __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
 
-  /* "slopescreening/screening/gap_test_p_1.pyx":160
+  /* "slopescreening/screening/gap_test_p_1.pyx":164
  *       sigmas[1:] = np.cumsum(Atcabs[index])
  * 
  *       for q in range(n-1, -1, -1):             # <<<<<<<<<<<<<<
@@ -2933,7 +2989,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
   for (__pyx_t_12 = (__pyx_v_n - 1); __pyx_t_12 > -1; __pyx_t_12-=1) {
     __pyx_v_q = __pyx_t_12;
 
-    /* "slopescreening/screening/gap_test_p_1.pyx":161
+    /* "slopescreening/screening/gap_test_p_1.pyx":165
  * 
  *       for q in range(n-1, -1, -1):
  *          bound = coeff_lbd * self.vec_cumsum_gammas[q] - self.vec_kappa_q[q] * radius - sigmas[q]             # <<<<<<<<<<<<<<
@@ -2943,7 +2999,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
     __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)));
 
-    /* "slopescreening/screening/gap_test_p_1.pyx":163
+    /* "slopescreening/screening/gap_test_p_1.pyx":167
  *          bound = coeff_lbd * self.vec_cumsum_gammas[q] - self.vec_kappa_q[q] * radius - sigmas[q]
  * 
  *          if Atcabs[index[q]] < bound:             # <<<<<<<<<<<<<<
@@ -2955,7 +3011,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
     __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) {
 
-      /* "slopescreening/screening/gap_test_p_1.pyx":164
+      /* "slopescreening/screening/gap_test_p_1.pyx":168
  * 
  *          if Atcabs[index[q]] < bound:
  *             l_q = 0             # <<<<<<<<<<<<<<
@@ -2964,7 +3020,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
  */
       __pyx_v_l_q = 0;
 
-      /* "slopescreening/screening/gap_test_p_1.pyx":163
+      /* "slopescreening/screening/gap_test_p_1.pyx":167
  *          bound = coeff_lbd * self.vec_cumsum_gammas[q] - self.vec_kappa_q[q] * radius - sigmas[q]
  * 
  *          if Atcabs[index[q]] < bound:             # <<<<<<<<<<<<<<
@@ -2974,7 +3030,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
       goto __pyx_L6;
     }
 
-    /* "slopescreening/screening/gap_test_p_1.pyx":166
+    /* "slopescreening/screening/gap_test_p_1.pyx":170
  *             l_q = 0
  * 
  *          elif Atcabs[index[n-1]] >= bound:             # <<<<<<<<<<<<<<
@@ -2986,7 +3042,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
     __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) {
 
-      /* "slopescreening/screening/gap_test_p_1.pyx":167
+      /* "slopescreening/screening/gap_test_p_1.pyx":171
  * 
  *          elif Atcabs[index[n-1]] >= bound:
  *             return calI_screen             # <<<<<<<<<<<<<<
@@ -2998,7 +3054,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
       __pyx_r = ((PyObject *)__pyx_v_calI_screen);
       goto __pyx_L0;
 
-      /* "slopescreening/screening/gap_test_p_1.pyx":166
+      /* "slopescreening/screening/gap_test_p_1.pyx":170
  *             l_q = 0
  * 
  *          elif Atcabs[index[n-1]] >= bound:             # <<<<<<<<<<<<<<
@@ -3007,7 +3063,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
  */
     }
 
-    /* "slopescreening/screening/gap_test_p_1.pyx":170
+    /* "slopescreening/screening/gap_test_p_1.pyx":174
  * 
  *          else:
  *             l_q = q + np.argmax(Atcabs[index[q:]] < bound)             # <<<<<<<<<<<<<<
@@ -3015,27 +3071,27 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
  *          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, 170, __pyx_L1_error)
+      __pyx_t_2 = __Pyx_PyInt_From_int(__pyx_v_q); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 174, __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, 170, __pyx_L1_error)
+      __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 174, __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, 170, __pyx_L1_error)
+      __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_argmax); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 174, __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, 170, __pyx_L1_error)
+      __pyx_t_3 = __Pyx_PyInt_From_int(__pyx_v_q); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 174, __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, 170, __pyx_L1_error)
+      __pyx_t_5 = PySlice_New(__pyx_t_3, Py_None, Py_None); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 174, __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, 170, __pyx_L1_error)
+      __pyx_t_3 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_index), __pyx_t_5); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 174, __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, 170, __pyx_L1_error)
+      __pyx_t_5 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_Atcabs), __pyx_t_3); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 174, __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, 170, __pyx_L1_error)
+      __pyx_t_3 = PyFloat_FromDouble(__pyx_v_bound); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 174, __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, 170, __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, 174, __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;
@@ -3051,20 +3107,20 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
       __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, 170, __pyx_L1_error)
+      if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 174, __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, 170, __pyx_L1_error)
+      __pyx_t_4 = PyNumber_Add(__pyx_t_2, __pyx_t_1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 174, __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, 170, __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, 174, __pyx_L1_error)
       __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
       __pyx_v_l_q = __pyx_t_19;
     }
     __pyx_L6:;
 
-    /* "slopescreening/screening/gap_test_p_1.pyx":172
+    /* "slopescreening/screening/gap_test_p_1.pyx":176
  *             l_q = q + np.argmax(Atcabs[index[q:]] < bound)
  * 
  *          l_min = max(l_min, l_q)             # <<<<<<<<<<<<<<
@@ -3081,24 +3137,24 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
     __pyx_v_l_min = __pyx_t_21;
   }
 
-  /* "slopescreening/screening/gap_test_p_1.pyx":174
+  /* "slopescreening/screening/gap_test_p_1.pyx":178
  *          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, 174, __pyx_L1_error)
+  __pyx_t_4 = __Pyx_PyInt_From_int(__pyx_v_l_min); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 178, __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, 174, __pyx_L1_error)
+  __pyx_t_1 = PySlice_New(__pyx_t_4, Py_None, Py_None); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 178, __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, 174, __pyx_L1_error)
+  __pyx_t_4 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_index), __pyx_t_1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 178, __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, 174, __pyx_L1_error)
+  if (unlikely(PyObject_SetItem(((PyObject *)__pyx_v_calI_screen), __pyx_t_4, Py_True) < 0)) __PYX_ERR(0, 178, __pyx_L1_error)
   __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
 
-  /* "slopescreening/screening/gap_test_p_1.pyx":175
+  /* "slopescreening/screening/gap_test_p_1.pyx":179
  * 
  *       calI_screen[index[l_min:]] = True
  *       return calI_screen             # <<<<<<<<<<<<<<
@@ -3108,7 +3164,7 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
   __pyx_r = ((PyObject *)__pyx_v_calI_screen);
   goto __pyx_L0;
 
-  /* "slopescreening/screening/gap_test_p_1.pyx":92
+  /* "slopescreening/screening/gap_test_p_1.pyx":96
  *    @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,             # <<<<<<<<<<<<<<
@@ -3159,19 +3215,19 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
  */
 
 /* Python wrapper */
-static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_9__reduce_cython__(PyObject *__pyx_v_self, CYTHON_UNUSED PyObject *unused); /*proto*/
-static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_9__reduce_cython__(PyObject *__pyx_v_self, CYTHON_UNUSED PyObject *unused) {
+static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_11__reduce_cython__(PyObject *__pyx_v_self, CYTHON_UNUSED PyObject *unused); /*proto*/
+static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_11__reduce_cython__(PyObject *__pyx_v_self, CYTHON_UNUSED PyObject *unused) {
   PyObject *__pyx_r = 0;
   __Pyx_RefNannyDeclarations
   __Pyx_RefNannySetupContext("__reduce_cython__ (wrapper)", 0);
-  __pyx_r = __pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_8__reduce_cython__(((struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *)__pyx_v_self));
+  __pyx_r = __pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_10__reduce_cython__(((struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *)__pyx_v_self));
 
   /* function exit code */
   __Pyx_RefNannyFinishContext();
   return __pyx_r;
 }
 
-static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_8__reduce_cython__(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *__pyx_v_self) {
+static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_10__reduce_cython__(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *__pyx_v_self) {
   PyObject *__pyx_r = NULL;
   __Pyx_RefNannyDeclarations
   PyObject *__pyx_t_1 = NULL;
@@ -3216,19 +3272,19 @@ static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPe
  */
 
 /* Python wrapper */
-static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_11__setstate_cython__(PyObject *__pyx_v_self, PyObject *__pyx_v___pyx_state); /*proto*/
-static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_11__setstate_cython__(PyObject *__pyx_v_self, PyObject *__pyx_v___pyx_state) {
+static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_13__setstate_cython__(PyObject *__pyx_v_self, PyObject *__pyx_v___pyx_state); /*proto*/
+static PyObject *__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_13__setstate_cython__(PyObject *__pyx_v_self, PyObject *__pyx_v___pyx_state) {
   PyObject *__pyx_r = 0;
   __Pyx_RefNannyDeclarations
   __Pyx_RefNannySetupContext("__setstate_cython__ (wrapper)", 0);
-  __pyx_r = __pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_10__setstate_cython__(((struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *)__pyx_v_self), ((PyObject *)__pyx_v___pyx_state));
+  __pyx_r = __pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_12__setstate_cython__(((struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *)__pyx_v_self), ((PyObject *)__pyx_v___pyx_state));
 
   /* function exit code */
   __Pyx_RefNannyFinishContext();
   return __pyx_r;
 }
 
-static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_10__setstate_cython__(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *__pyx_v_self, CYTHON_UNUSED PyObject *__pyx_v___pyx_state) {
+static PyObject *__pyx_pf_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_12__setstate_cython__(CYTHON_UNUSED struct __pyx_obj_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne *__pyx_v_self, CYTHON_UNUSED PyObject *__pyx_v___pyx_state) {
   PyObject *__pyx_r = NULL;
   __Pyx_RefNannyDeclarations
   PyObject *__pyx_t_1 = NULL;
@@ -4312,9 +4368,10 @@ static void __pyx_tp_dealloc_14slopescreening_9screening_12gap_test_p_1_GapTestP
 
 static PyMethodDef __pyx_methods_14slopescreening_9screening_12gap_test_p_1_GapTestPequalOne[] = {
   {"get_name", (PyCFunction)__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_5get_name, METH_NOARGS, 0},
-  {"apply_test", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_7apply_test, METH_VARARGS|METH_KEYWORDS, __pyx_doc_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_6apply_test},
-  {"__reduce_cython__", (PyCFunction)__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_9__reduce_cython__, METH_NOARGS, 0},
-  {"__setstate_cython__", (PyCFunction)__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_11__setstate_cython__, METH_O, 0},
+  {"get_legend_name", (PyCFunction)__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_7get_legend_name, METH_NOARGS, 0},
+  {"apply_test", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_9apply_test, METH_VARARGS|METH_KEYWORDS, __pyx_doc_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_8apply_test},
+  {"__reduce_cython__", (PyCFunction)__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_11__reduce_cython__, METH_NOARGS, 0},
+  {"__setstate_cython__", (PyCFunction)__pyx_pw_14slopescreening_9screening_12gap_test_p_1_16GapTestPequalOne_13__setstate_cython__, METH_O, 0},
   {0, 0, 0, 0}
 };
 
@@ -4459,6 +4516,7 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = {
   {&__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},
   {&__pyx_kp_s_numpy_core_umath_failed_to_impor, __pyx_k_numpy_core_umath_failed_to_impor, sizeof(__pyx_k_numpy_core_umath_failed_to_impor), 0, 0, 1, 0},
   {&__pyx_n_s_offset_radius, __pyx_k_offset_radius, sizeof(__pyx_k_offset_radius), 0, 0, 1, 1},
+  {&__pyx_kp_s_p_q_1_forall_q, __pyx_k_p_q_1_forall_q, sizeof(__pyx_k_p_q_1_forall_q), 0, 0, 1, 0},
   {&__pyx_n_s_range, __pyx_k_range, sizeof(__pyx_k_range), 0, 0, 1, 1},
   {&__pyx_n_s_reduce, __pyx_k_reduce, sizeof(__pyx_k_reduce), 0, 0, 1, 1},
   {&__pyx_n_s_reduce_cython, __pyx_k_reduce_cython, sizeof(__pyx_k_reduce_cython), 0, 0, 1, 1},
@@ -4487,25 +4545,25 @@ static CYTHON_SMALL_CODE int __Pyx_InitCachedConstants(void) {
   __Pyx_RefNannyDeclarations
   __Pyx_RefNannySetupContext("__Pyx_InitCachedConstants", 0);
 
-  /* "slopescreening/screening/gap_test_p_1.pyx":155
+  /* "slopescreening/screening/gap_test_p_1.pyx":159
  *       # 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, 155, __pyx_L1_error)
+  __pyx_slice_ = PySlice_New(Py_None, Py_None, __pyx_int_neg_1); if (unlikely(!__pyx_slice_)) __PYX_ERR(0, 159, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_slice_);
   __Pyx_GIVEREF(__pyx_slice_);
 
-  /* "slopescreening/screening/gap_test_p_1.pyx":158
+  /* "slopescreening/screening/gap_test_p_1.pyx":162
  * 
  *          # 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, 158, __pyx_L1_error)
+  __pyx_slice__2 = PySlice_New(__pyx_int_1, Py_None, Py_None); if (unlikely(!__pyx_slice__2)) __PYX_ERR(0, 162, __pyx_L1_error)
   __Pyx_GOTREF(__pyx_slice__2);
   __Pyx_GIVEREF(__pyx_slice__2);
 
diff --git a/slopescreening/screening/gap_test_p_1.pyx b/slopescreening/screening/gap_test_p_1.pyx
index 5b798fb..17cdf58 100644
--- a/slopescreening/screening/gap_test_p_1.pyx
+++ b/slopescreening/screening/gap_test_p_1.pyx
@@ -86,6 +86,10 @@ cdef class GapTestPequalOne:
    def get_name(self):
       return self.name
 
+   @cython.boundscheck(False) # turn off bounds-checking for entire function
+   @cython.wraparound(False)  # turn off negative index wrapping for entire function
+   def get_legend_name(self):
+      return "$p_q=1\\;\\forall q$"
 
    @cython.boundscheck(False) # turn off bounds-checking for entire function
    @cython.wraparound(False)  # turn off negative index wrapping for entire function
diff --git a/slopescreening/screening/gap_test_p_q.py b/slopescreening/screening/gap_test_p_q.py
index 1923359..83e88c0 100644
--- a/slopescreening/screening/gap_test_p_q.py
+++ b/slopescreening/screening/gap_test_p_q.py
@@ -11,6 +11,7 @@ class GapTestPequalQ(AbstractGapScreening):
       super(GapTestPequalQ, self).__init__()
 
       self.name = "Gap test p=q"
+      self.legend_name = "$p_q=q\\;\\forall q$"
 
    def apply_test(self, Atu_abs, gap, lbd, vec_gammas, coeff_dual_scaling=1., offset_radius=0, index=None) -> np.ndarray:
       """ GAP Sphere screening test detailed in Cédric's node
-- 
GitLab