Fix correctness bugs across the core modules

glmdisc/__init__.py

@@ -215,8 +215,14 @@ def __init__(self, algorithm="SEM", test=True, validation=True, criterion="bic",
             raise ValueError(msg)
 
         # m_start doit être pas déconnant
-        if not 2 <= m_start <= 50:
-            msg = 'Please set 2 <= m_start <= 50'
+        if not isinstance(m_start, (int, np.integer)) or not 2 <= m_start <= 50:
+            msg = 'Please set an integer 2 <= m_start <= 50'
+            logger.error(msg)
+            raise ValueError(msg)
+
+        # burn_in must be a non-negative integer
+        if not isinstance(burn_in, (int, np.integer)) or isinstance(burn_in, bool) or burn_in < 0:
+            msg = 'Please set an integer burn_in >= 0'
             logger.error(msg)
             raise ValueError(msg)
 
@@ -266,6 +272,8 @@ def _check_is_fitted(self):
         """
         if self.algorithm == "SEM":
             try:
+                if self.best_reglog is None:
+                    raise sk.exceptions.NotFittedError("This Glmdisc instance is not fitted yet.")
                 sk.utils.validation.check_is_fitted(self.best_reglog)
                 for link in self.best_link:
                     if isinstance(link, sk.linear_model.LogisticRegression):

glmdisc/_bestFormula.py

@@ -16,9 +16,14 @@ def _best_formula_cont_sem(self, emap_best):
     for j in range(self.d_cont):
         best_disc.append([])
         for k in np.unique(emap_best[:, j]):
-            best_disc[j].append(
-                np.nanmin(self.predictors_cont[emap_best[:, j] == k, j]))
-        del best_disc[j][np.where(best_disc[j] == np.nanmin(best_disc[j]))[0][0]]
+            lower_bound = self.predictors_cont[emap_best[:, j] == k, j]
+            lower_bound = lower_bound[~np.isnan(lower_bound)]
+            # the level reserved for missing values has no lower bound
+            if lower_bound.size > 0:
+                best_disc[j].append(lower_bound.min())
+        # the smallest lower bound is not a cut-point
+        if best_disc[j]:
+            del best_disc[j][np.argmin(best_disc[j])]
         if len(best_disc[j]) < 1:
             logger.info("No cut-points found for continuous variable " + str(j) + ".")
         else:

glmdisc/_discretize.py

@@ -36,7 +36,7 @@ def _check_args_discretize_nn(self, predictors_cont, predictors_qual):
         raise ValueError(msg)
     if d_2 != self.d_qual:
         msg = ('Shape of ' + str(d_2) +
-               ' for predictors_cont does not match training set of size ' + str(self.d_qual) + '.')
+               ' for predictors_qual does not match training set of size ' + str(self.d_qual) + '.')
         logger.error(msg)
         raise ValueError(msg)
 
@@ -70,7 +70,7 @@ def _check_args_discretize_sem(self, predictors_cont, predictors_qual):
         raise ValueError(msg)
     if d_2 != sum(d_2bis) or d_2 != self.d_qual:
         msg = ('Shape of ' + str(d_2) +
-               ' for predictors_cont does not match provided link function '
+               ' for predictors_qual does not match provided link function '
                'of size ' + str(sum(d_2bis)) + ' and/or training set of size ' + str(self.d_qual) + '.')
         logger.error(msg)
         raise ValueError(msg)

glmdisc/_fit.py

@@ -5,7 +5,6 @@
 import numpy as np
 import sklearn as sk
 from loguru import logger
-from scipy import stats
 
 from glmdisc._fitNN import _fit_nn
 from glmdisc._fitSEM import _fit_sem
@@ -96,7 +95,7 @@ def _init_disc(self, continu_complete_case):
     initial quantization and transformation of categorical features to integers
     """
     self.affectations = [None] * (self.d_cont + self.d_qual)
-    edisc = np.random.choice(list(range(self.m_start)), size=(self.n, self.d_cont + self.d_qual))
+    edisc = np.random.randint(0, self.m_start, size=(self.n, self.d_cont + self.d_qual))
 
     for j in range(self.d_cont):
         edisc[np.invert(continu_complete_case[:, j]), j] = self.m_start
@@ -105,16 +104,9 @@ def _init_disc(self, continu_complete_case):
 
     for j in range(self.d_qual):
         self.affectations[j + self.d_cont] = sk.preprocessing.LabelEncoder().fit(self.predictors_qual[:, j])
-        if (self.m_start > stats.describe(sk.preprocessing.LabelEncoder().fit(self.predictors_qual[:, j]).transform(
-                self.predictors_qual[:, j])).minmax[1] + 1):
-            edisc[:, j + self.d_cont] = np.random.choice(list(range(
-                stats.describe(sk.preprocessing.LabelEncoder().fit(
-                    self.predictors_qual[:, j]).transform(
-                        self.predictors_qual[:, j])).minmax[1])),
-                size=self.n)
-        else:
-            edisc[:, j + self.d_cont] = np.random.choice(list(range(self.m_start)),
-                                                         size=self.n)
+        # a feature cannot start with more groups than it has levels
+        m_init = min(self.m_start, len(self.affectations[j + self.d_cont].classes_))
+        edisc[:, j + self.d_cont] = np.random.randint(0, m_init, size=self.n)
 
         predictors_trans[:, j] = (self.affectations[j + self.d_cont].transform(
             self.predictors_qual[:, j])).astype(int)

glmdisc/_fitNN.py

@@ -12,7 +12,6 @@
 import tensorflow as tf
 import tensorflow.keras.optimizers
 from loguru import logger
-from scipy import stats
 from tensorflow.keras import Model
 from tensorflow.keras.callbacks import Callback, ReduceLROnPlateau
 from tensorflow.keras.layers import Dense
@@ -152,8 +151,7 @@ def _initialize_neural_net(self):
     m_qual = [self.m_start] * (self.d_cont + self.d_qual)
 
     for j in range(self.d_qual):
-        num_levels = stats.describe(sk.preprocessing.LabelEncoder().fit(self.predictors_qual[:, j]).transform(
-            self.predictors_qual[:, j])).minmax[1]
+        num_levels = len(self.affectations[j + self.d_cont].classes_) - 1
         if self.m_start > num_levels:
             m_qual[j] = num_levels
 

glmdisc/_generateData.py

@@ -36,19 +36,20 @@ def generate_data(n, d, theta=None, plot=False):
     for i in range(d):
         xd[:, i] = pd.cut(x[:, i], bins=cuts, labels=[0, 1, 2])
 
-    log_odd = np.array([0] * n)
+    log_odd = np.zeros(n)
     for i in range(n):
         for j in range(d):
             log_odd[i] += theta[int(xd[i, j]), j]
 
     p = 1 / (1 + np.exp(- log_odd))
-    y = np.random.binomial(1, p)
+    y = np.random.binomial(1, p, size=n)
 
     if plot and d > 1:
-        for unique_pair_of_feature_as_tuple in combinations(range(d), 2):
-            plt.scatter(x[:, unique_pair_of_feature_as_tuple[0]],
-                        x[:, unique_pair_of_feature_as_tuple[1]],
-                        c=y,
-                        s=[1 for _ in range(len(x))])
+        for j1, j2 in combinations(range(d), 2):
+            plt.figure()
+            plt.scatter(x[:, j1], x[:, j2], c=y, s=1)
+            plt.xlabel("x" + str(j1))
+            plt.ylabel("x" + str(j2))
+            plt.show(block=False)
 
     return [x, y, theta]

glmdisc/_gini_utils.py

@@ -24,15 +24,15 @@ def compute_midrank(x):
     J = np.argsort(x)
     Z = x[J]
     N = len(x)
-    T = np.zeros(N, dtype=np.float)
+    T = np.zeros(N, dtype=float)
     i = 0
     while i < N:
         j = i
         while j < N and Z[j] == Z[i]:
             j += 1
         T[i:j] = 0.5 * (i + j - 1)
         i = j
-    T2 = np.empty(N, dtype=np.float)
+    T2 = np.empty(N, dtype=float)
     # Note(kazeevn) +1 is due to Python using 0-based indexing
     # instead of 1-based in the AUC formula in the paper
     T2[J] = T + 1
@@ -50,15 +50,15 @@ def compute_midrank_weight(x, sample_weight):
     Z = x[J]
     cumulative_weight = np.cumsum(sample_weight[J])
     N = len(x)
-    T = np.zeros(N, dtype=np.float)
+    T = np.zeros(N, dtype=float)
     i = 0
     while i < N:
         j = i
         while j < N and Z[j] == Z[i]:
             j += 1
         T[i:j] = cumulative_weight[i:j].mean()
         i = j
-    T2 = np.empty(N, dtype=np.float)
+    T2 = np.empty(N, dtype=float)
     T2[J] = T
     return T2
 
@@ -99,9 +99,9 @@ def fast_de_long_weights(predictions_sorted_transposed, label_1_count, sample_we
     negative_examples = predictions_sorted_transposed[:, m:]
     k = predictions_sorted_transposed.shape[0]
 
-    tx = np.empty([k, m], dtype=np.float)
-    ty = np.empty([k, n], dtype=np.float)
-    tz = np.empty([k, m + n], dtype=np.float)
+    tx = np.empty([k, m], dtype=float)
+    ty = np.empty([k, n], dtype=float)
+    tz = np.empty([k, m + n], dtype=float)
     for r in range(k):
         tx[r, :] = compute_midrank_weight(positive_examples[r, :], sample_weight[:m])
         ty[r, :] = compute_midrank_weight(negative_examples[r, :], sample_weight[m:])
@@ -149,9 +149,9 @@ def fast_de_long_no_weights(predictions_sorted_transposed, label_1_count):
     negative_examples = predictions_sorted_transposed[:, m:]
     k = predictions_sorted_transposed.shape[0]
 
-    tx = np.empty([k, m], dtype=np.float)
-    ty = np.empty([k, n], dtype=np.float)
-    tz = np.empty([k, m + n], dtype=np.float)
+    tx = np.empty([k, m], dtype=float)
+    ty = np.empty([k, n], dtype=float)
+    tz = np.empty([k, m + n], dtype=float)
     for r in range(k):
         tx[r, :] = compute_midrank(positive_examples[r, :])
         ty[r, :] = compute_midrank(negative_examples[r, :])
@@ -201,7 +201,7 @@ def gini(ground_truth, predictions, sample_weight=None, level=.95):
        sample_weight: np.array of floats of the sample weights
        level: float (between 0 and 1) describing the level required for the confidence interval
     """
-    if (not isinstance(level, float)) or (not level < 1 and level > 0):
+    if not isinstance(level, float) or not 0 < level < 1:
         raise ValueError("level should be a single float between 0 and 1")
     auc, auc_cov = delong_roc_variance(ground_truth, predictions, sample_weight)
     auc_std = np.sqrt(auc_cov)

glmdisc/_plot.py

@@ -50,24 +50,31 @@ def _plot_disc(self, predictors_cont_number, predictors_qual_number):
 
 
 def _plot_logodd(self):
-    # Gérer les manquants dans le GAM
-    lignes_completes = np.invert(
-        np.isnan(self.predictors_cont).sum(axis=1).astype(bool))
+    # Handle missing values in the GAM: keep the complete rows of the continuous features
+    if self.d_cont > 0:
+        lignes_completes = np.invert(np.isnan(self.predictors_cont).any(axis=1))
+    else:
+        lignes_completes = np.ones(self.n, dtype=bool)
+
+    features = []
+    if self.d_cont > 0:
+        features.append(pd.DataFrame(self.predictors_cont[lignes_completes, :]).apply(
+            lambda x: x.astype('float')))
+    if self.d_qual > 0:
+        features.append(pd.DataFrame(self.predictors_qual[lignes_completes, :]).apply(
+            lambda x: x.astype('category')))
 
-    # Fit du GAM sur tout le monde
+    # Fit the GAM on everyone
     gam = LogisticGAM(dtype=['numerical' for _ in range(self.d_cont)] + ['categorical' for _ in range(
-        self.d_qual)]).fit(
-        pd.concat([pd.DataFrame(self.predictors_cont[lignes_completes, :]).apply(
-            lambda x: x.astype('float')),
-            pd.DataFrame(self.predictors_qual[lignes_completes, :]).apply(
-                lambda x: x.astype('category'))], axis=1), self.labels[lignes_completes])
-
-    # Quelles que soient les valeurs de predictors_cont_number et
-    # predictors_qual_number, on plot tout pour l'instant
+        self.d_qual)]).fit(pd.concat(features, axis=1), self.labels[lignes_completes])
+
+    # Whatever the values of predictors_cont_number and
+    # predictors_qual_number, everything is plotted for now
     plt.figure()
-    fig, axs = plt.subplots(1, self.d_cont + self.d_qual)
+    # squeeze=False keeps axs a 2D array even with a single feature
+    fig, axs = plt.subplots(1, self.d_cont + self.d_qual, squeeze=False)
     plt.rcParams['figure.figsize'] = (28, 8)
-    for i, ax in enumerate(axs):
+    for i, ax in enumerate(axs[0]):
         try:
             XX = gam.generate_X_grid(term=i)
             ax.plot(XX[:, i], gam.partial_dependence(term=i, X=XX))

tests/test_gini.py

@@ -0,0 +1,50 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+import numpy as np
+import pytest
+from sklearn.metrics import roc_auc_score
+
+import glmdisc
+
+
+def _scores(n=200, seed=0):
+    rng = np.random.RandomState(seed)
+    y = rng.binomial(1, 0.4, size=n)
+    p = np.clip(0.3 * y + rng.uniform(size=n) * 0.7, 0, 1)
+    return y, p
+
+
+def test_gini_matches_auc():
+    y, p = _scores()
+    value, ci = glmdisc.gini(y, p)
+    assert value == pytest.approx(2 * roc_auc_score(y, p) - 1)
+    assert ci[0] < value < ci[1] <= 1
+
+
+def test_gini_with_ties():
+    y, p = _scores()
+    p = np.round(p, 1)
+    value, _ = glmdisc.gini(y, p)
+    assert value == pytest.approx(2 * roc_auc_score(y, p) - 1)
+
+
+def test_gini_with_sample_weight():
+    y, p = _scores()
+    w = np.random.RandomState(1).uniform(0.5, 2, size=len(y))
+    value, _ = glmdisc.gini(y, p, sample_weight=w)
+    assert value == pytest.approx(2 * roc_auc_score(y, p, sample_weight=w) - 1)
+
+
+@pytest.mark.parametrize("level", [0., -.5, 1., 1.5, 1, "0.95", None])
+def test_gini_invalid_level(level):
+    y, p = _scores()
+    with pytest.raises(ValueError):
+        glmdisc.gini(y, p, level=level)
+
+
+def test_gini_level_widens_interval():
+    y, p = _scores()
+    _, ci_90 = glmdisc.gini(y, p, level=.9)
+    _, ci_99 = glmdisc.gini(y, p, level=.99)
+    assert ci_99[0] < ci_90[0]
+    assert ci_99[1] >= ci_90[1]

tests/test_init.py

@@ -64,6 +64,17 @@ def test_m_start():
     with pytest.raises(ValueError):
         glmdisc.Glmdisc(m_start=1)
 
+    with pytest.raises(ValueError):
+        glmdisc.Glmdisc(m_start=2.5)
+
+
+def test_burn_in():
+    assert glmdisc.Glmdisc(burn_in=0).burn_in == 0
+
+    for burn_in in [-1, 2.5, "5", None, True]:
+        with pytest.raises(ValueError):
+            glmdisc.Glmdisc(burn_in=burn_in)
+
 
 def test_gini_penalized(caplog):
     glmdisc.Glmdisc(validation=False,
@@ -83,3 +94,11 @@ def test_validation_criterion(caplog):
 
     assert caplog.records[0].message == ('No need to penalize the log-likelihood when a validation set is used. '
                                          'Using log-likelihood instead.')
+
+
+def test_not_fitted_sem():
+    model = glmdisc.Glmdisc()
+    with pytest.raises(glmdisc.NotFittedError):
+        model._check_is_fitted()
+    with pytest.raises(glmdisc.NotFittedError):
+        model.predict(np.zeros((3, 1)), None)

tests/tests_NN/test_NN_fit.py

@@ -5,6 +5,8 @@
 import pandas as pd
 import random
 import glmdisc
+import sklearn as sk
+from glmdisc._fitNN import _initialize_neural_net
 import tensorflow
 from tensorflow.keras.optimizers import Adagrad
 from tensorflow.keras.callbacks import EarlyStopping
@@ -209,3 +211,20 @@ def test_fit():
     model.fit(predictors_cont=x, predictors_qual=None, labels=y, iter=20)
     model = glmdisc.Glmdisc(algorithm="NN", criterion="gini", test=False, validation=False)
     model.fit(predictors_cont=x, predictors_qual=None, labels=y, iter=20)
+
+
+def test_initialize_neural_net_reuses_label_encoders(monkeypatch):
+    np.random.seed(0)
+    x, y, theta = glmdisc.Glmdisc.generate_data(100, 2)
+    xd = np.digitize(x, [0.333, 0.666]).astype(float)
+    model = glmdisc.Glmdisc(algorithm="NN", validation=False, test=False)
+    model.predictors_cont = x
+    model.predictors_qual = xd
+    model.labels = y
+    model._init_disc(model._calculate_shape())
+
+    def fail_fit(self, y):
+        raise AssertionError("LabelEncoder fitted again")
+
+    monkeypatch.setattr(sk.preprocessing.LabelEncoder, "fit", fail_fit)
+    _initialize_neural_net(model)