Updating figures and tests

Author: Carole Sudre

MetricsReloaded/metrics/calibration_measures.py

@@ -31,7 +31,7 @@
 import numpy as np
 import math
 from scipy.special import gamma
-
+import warnings
 # from metrics.pairwise_measures import CacheFunctionOutput
 from MetricsReloaded.utility.utils import (
     CacheFunctionOutput,
@@ -150,6 +150,7 @@ def expectation_calibration_error(self):
         if "bins_ece" in self.dict_args:
             nbins = self.dict_args["bins_ece"]
         else:
+            warnings.warn("Bins ECE not specified in optional arguments dictionary - default set to 10")
             nbins = 10
         step = 1.0 / nbins
         range_values = np.arange(0, 1.00001, step)
@@ -176,7 +177,55 @@ def expectation_calibration_error(self):
             else:
                 list_values.append(nsamples * np.abs(prop - np.mean(pred_sel)))
             numb_samples += nsamples
-        return np.sum(np.asarray(list_values)) / numb_samples
+        ece = np.sum(np.asarray(list_values)) / numb_samples
+        return ece
+    
+
+    def maximum_calibration_error(self):
+        """
+        Derives the maximum calibration error in the case of binary task
+        bins_mce is the key in the dictionary for the number of bins to consider
+        Default is 10
+
+        .. math::
+
+            MCE = max(|\dfrac{1}{|B_m|}\sum_{i \in B_m}1(pred_ik==ref_ik)-\dfrac{1}{|B_m|}\sum_{i \in B_m}pred_i|)
+
+        :return: mce
+
+        """
+        if "bins_mce" in self.dict_args:
+            nbins = self.dict_args["bins_mce"]
+        else:
+            warnings.warn("Bins MCE not specified in optional arguments dictionary - default set to 10")
+            nbins = 10
+        step = 1.0 / nbins
+        range_values = np.arange(0, 1.00001, step)
+        list_values = []
+        numb_samples = 0
+        pred_prob = self.pred[:,1]
+        for (l, u) in zip(range_values[:-1], range_values[1:]):
+            ref_tmp = np.where(
+                np.logical_and(pred_prob > l, pred_prob <= u),
+                self.ref,
+                np.ones_like(self.ref) * -1,
+            )
+            ref_sel = ref_tmp[ref_tmp > -1]
+            nsamples = np.size(ref_sel)
+            prop = np.sum(ref_sel) / nsamples
+            pred_tmp = np.where(
+                np.logical_and(pred_prob > l, pred_prob <= u),
+                pred_prob,
+                np.ones_like(pred_prob) * -1,
+            )
+            pred_sel = pred_tmp[pred_tmp > -1]
+            if nsamples == 0:
+                list_values.append(0)
+            else:
+                list_values.append(np.abs(prop - np.mean(pred_sel)))
+        mce = np.max(np.asarray(list_values))
+        return mce
+
 
     def brier_score(self):
         """

MetricsReloaded/metrics/pairwise_measures.py

@@ -265,6 +265,11 @@ def __init__(
             "fbeta": (self.fbeta, "FBeta"),
             "dsc":(self.dsc, "DSC"),
             "youden_ind": (self.youden_index, "YoudenInd"),
+            "ppv":(self.positive_predictive_value,'PPV'),
+            "npv":(self.negative_predictive_value,'NPV'),
+            "ior":(self.intersection_over_reference,"IoR"),
+            "sensitivity":(self.sensitivity,"Sens"),
+            "specificity":(self.specificity,"Spec"),
             "mcc": (self.matthews_correlation_coefficient, "MCC"),
             "cldice": (self.centreline_dsc, "CentreLineDSC"),
             "assd": (self.measured_average_distance, "ASSD"),
@@ -693,7 +698,7 @@ def pred_in_ref(self):
         else:
             return 0
 
-    def positive_predictive_values(self):
+    def positive_predictive_value(self):
         """
         Calculates the positive predictive value
 
@@ -785,10 +790,10 @@ def fbeta(self):
             warnings.warn("beta value not specified in option - default set to 1")
             beta = 1
         numerator = (
-            (1 + np.square(beta)) * self.positive_predictive_values() * self.recall()
+            (1 + np.square(beta)) * self.positive_predictive_value() * self.recall()
         )
         denominator = (
-            np.square(beta) * self.positive_predictive_values() + self.recall()
+            np.square(beta) * self.positive_predictive_value() + self.recall()
         )
         if np.isnan(denominator):
             if self.fp() + self.fn() > 0:
@@ -830,7 +835,7 @@ def net_benefit_treated(self):
         net_benefit = tp / n - fp / n * er
         return net_benefit
 
-    def negative_predictive_values(self):
+    def negative_predictive_value(self):
         """
         This function calculates the negative predictive value ratio between
         the number of true negatives and the total number of negative elements

test/test_metrics/test_calibration_metrics.py

@@ -4,9 +4,7 @@
 from scipy.special import gamma
 from MetricsReloaded.utility.utils import median_heuristic
 
-
-def test_expected_calibration_error():
-    f40_pred = [[1-0.22, 0.22 ],
+pred_224 = [[1-0.22, 0.22 ],
                 [1-0.48, 0.48],
                 [0.51,0.49],
                 [0.04, 0.96],
@@ -17,15 +15,38 @@ def test_expected_calibration_error():
                 [0.66, 0.34],
                 [0.13, 0.87]]
     #f40_pred = [0.22, 0.48, 0.49, 0.96, 0.55, 0.64, 0.78, 0.82, 0.34, 0.87]
-    f40_ref = [0, 1, 0, 0, 1, 1, 1, 1, 1, 0]
-    ppm = CalibrationMeasures(f40_pred, f40_ref)
-    ppm1 = CalibrationMeasures(f40_pred, f40_ref, dict_args={"bins_ece": 2})
-    value_test2 = ppm.expectation_calibration_error()
+ref_224 = [0, 1, 0, 0, 1, 1, 1, 1, 1, 0]
+
+def test_expected_calibration_error():
+    """
+    Using as reference SN 2.24 p67
+    """
+    ppm1 = CalibrationMeasures(pred_224, ref_224, dict_args={"bins_ece": 2})
+    ppm2 = CalibrationMeasures(pred_224, ref_224, dict_args={'bins_ece':5})
+    ppm3 = CalibrationMeasures(pred_224, ref_224)
     value_test1 = ppm1.expectation_calibration_error()
+    value_test2 = ppm2.expectation_calibration_error()
+    value_test3 = ppm3.expectation_calibration_error()
     expected_ece1 = 0.11
-    expected_ece2 = 0.36
+    expected_ece2 = 0.32
+    expected_ece3 = 0.36
+    assert_allclose(value_test1, expected_ece1, atol=0.01)
+    assert_allclose(value_test2, expected_ece2, atol=0.01)
+    assert_allclose(value_test3, expected_ece3, atol=0.01)
+
+def test_maximum_calibration_error():
+    ppm1 = CalibrationMeasures(pred_224, ref_224, dict_args={"bins_mce": 2})
+    ppm2 = CalibrationMeasures(pred_224, ref_224, dict_args={'bins_mce':5})
+    ppm3 = CalibrationMeasures(pred_224, ref_224)
+    value_test1 = ppm1.maximum_calibration_error()
+    value_test2 = ppm2.maximum_calibration_error()
+    value_test3 = ppm3.maximum_calibration_error()
+    expected_ece1 = 0.12
+    expected_ece2 = 0.55
+    expected_ece3 = 0.96
     assert_allclose(value_test1, expected_ece1, atol=0.01)
     assert_allclose(value_test2, expected_ece2, atol=0.01)
+    assert_allclose(value_test3, expected_ece3, atol=0.01)
 
 
 def test_logarithmic_score():

test/test_metrics/test_pairwise_measures.py

@@ -46,6 +46,14 @@
 ppm212_1 = PM(pred212, ref212)
 ppm212_2 = PM(pred212,ref212,dict_args={'boundary_dist':2})
 
+#Data for figure 5c (Hausdoff with annotation error p14 Pitfalls)
+ref5c = np.zeros([14, 14])
+ref5c[1, 1] = 1
+ref5c[9:12, 9:12] = 1
+pred5c = np.zeros([14, 14])
+pred5c [9:12, 9:12] = 1
+bpm5c = PM(pred5c, ref5c, dict_args={'hd_perc':95})
+
 ### Small size of structures relative to pixel/voxel size (DSC)
 ## Larger structure
 p_large_ref = np.zeros((11, 11))
@@ -77,7 +85,7 @@
 f27_ref2 = f27_pred1
 f27_pred2 = f27_ref1
 
-# Figure ClDice p 53 S2.14
+# Figure ClDice p 53 S2.14 pitfalls paper
 ref214 = np.zeros([24,24])
 ref214[1:10,7:12]=1
 ref214[10:12,3:19]=1
@@ -116,26 +124,26 @@
 
 
 
-# panoptic quality
-pq_pred1 = np.zeros([21, 21])
-pq_pred1[5:7, 2:5] = 1
-pq_pred2 = np.zeros([21, 21])
-pq_pred2[14:18, 4:6] = 1
-pq_pred2[16, 3] = 1
-pq_pred3 = np.zeros([21, 21])
-pq_pred3[14:18, 7:12] = 1
-pq_pred4 = np.zeros([21, 21])
-pq_pred4[2:8, 13:16] = 1
-pq_pred4[2:4, 12] = 1
-
-pq_ref1 = np.zeros([21, 21])
-pq_ref1[8:11, 3] = 1
-pq_ref1[9, 2:5] = 1
-pq_ref2 = np.zeros([21, 21])
-pq_ref2[14:19, 7:13] = 1
-pq_ref3 = np.zeros([21, 21])
-pq_ref3[2:7, 14:17] = 1
-pq_ref3[2:4, 12:14] = 1
+# panoptic quality Figure 3.51 p96
+pq_pred1 = np.zeros([18, 18])
+pq_pred1[ 3:7,1:3] = 1
+pq_pred1[3:6,3:7]=1
+pq_pred2 = np.zeros([18, 18])
+pq_pred2[13:16,4:6] = 1
+pq_pred3 = np.zeros([18, 18])
+pq_pred3[7:12,13:17] = 1
+pq_pred4 = np.zeros([18, 18])
+pq_pred4[13:15,13:17] = 1
+pq_pred4[15,15] = 1
+
+pq_ref1 = np.zeros([18, 18])
+pq_ref1[2:7, 1:3] = 1
+pq_ref1[2:5,3:6] = 1
+pq_ref2 = np.zeros([18, 18])
+pq_ref2[6:12,12:17] = 1
+pq_ref3 = np.zeros([18, 18])
+pq_ref3[14:15:,7:10] = 1
+pq_ref3[13:16,8:9] = 1
 
 f27_pred = np.concatenate([np.ones([81]), np.zeros([9]), np.ones([2]), np.zeros([8])])
 f27_ref = np.concatenate([np.ones([90]), np.zeros([10])])
@@ -324,16 +332,6 @@ def test_fn_map():
     fn2 = ppm210_2.fn()
     expected_fn1 = 12
     expected_fn2 = 0
-    # fn_map_1 = ppm210_1.__fn_map()
-    # expected_fn_map1 = np.zeros([14,14])
-    # expected_fn_map1[5:6,5:9] = 1
-    # expected_fn_map1[8:9,5:9] = 1 
-    # expected_fn_map1[5:9,5:6] = 1
-    # expected_fn_map1[5:9,8:9] = 1
-    # fn_map_2 = ppm210_2.__fn_map()
-    # expected_fn_map2 = np.zeros([14,14])
-    # assert_array_equal(fn_map_1, expected_fn_map1)
-    # assert_array_equal(fn_map_2, expected_fn_map2)
     assert fn1 == 12
     assert fn2 == 0 
 
@@ -553,8 +551,8 @@ def test_negative_predictive_value():
     """
     Taking figure SN 2.9 as inspiration p49 Pitfalls
     """
-    value_test1 = ppm29_1.negative_predictive_values()
-    value_test2 = ppm29_2.negative_predictive_values()
+    value_test1 = ppm29_1.negative_predictive_value()
+    value_test2 = ppm29_2.negative_predictive_value()
     expected_npv1 = 0.889
     expected_npv2 = 0.47
     assert_allclose(value_test1, expected_npv1, atol=0.001)
@@ -699,23 +697,27 @@ def test_nsd2():
     assert_allclose(value_test, expected_nsd2, atol=0.01)
 
 
-def test_iou():
+def test_intersection_over_union():
     bpm = PM(p_pred, p_ref)
     value_test = bpm.intersection_over_union()
     print("IoU ", value_test)
     expected_iou = 0.76
     assert_allclose(value_test, expected_iou, atol=0.01)
 
 
-def test_fbeta():
-    pm = PM(p_large_pred1, p_large_ref)
-    pm2 = PM(p_large_pred1, p_large_ref, dict_args={"beta": 1})
-    value_test = pm.fbeta()
-    value_test2 = pm2.fbeta()
-    print(value_test)
-    expected_fbeta = 0.986
-    assert_allclose(value_test, expected_fbeta, atol=0.001)
-    assert_allclose(value_test2, expected_fbeta, atol=0.001)
+def test_fbeta_beta_value():
+    """
+    Taking inspiration from SN 2.9 - p49 Pitfalls
+    """
+    expected_f11 = 0.86
+    expected_f12 = 0.94
+    ppm29_1.dict_args={'beta':1}
+    ppm29_2.dict_args={'beta':1}
+    value_test1 = ppm29_1.fbeta()
+    value_test2 = ppm29_2.fbeta()
+    assert_allclose(value_test1, expected_f11, atol=0.01)
+    assert_allclose(value_test2, expected_f12, atol=0.01)
+
 
 def test_sensitivity():
     """
@@ -749,13 +751,13 @@ def test_sens():
     assert_allclose(value_test, expected_sens, atol=0.01)
 
 
-def test_ppv():
+def test_positive_predictive_value():
     """
     Taking as inspiration figure SN2.9 p49 Pitfalls
     """
 
-    value_test1 = ppm29_1.positive_predictive_values()
-    value_test2 = ppm29_2.positive_predictive_values()
+    value_test1 = ppm29_1.positive_predictive_value()
+    value_test2 = ppm29_2.positive_predictive_value()
     expected_ppv1 = 0.82
     expected_ppv2 = 0.98
     assert_allclose(value_test1, expected_ppv1, atol=0.01)
@@ -817,15 +819,12 @@ def test_nsd_s210():
     assert_allclose(nsd_1,expected_nsd1,atol=0.01)
     assert_allclose(nsd_2,expected_nsd2,atol=0.01)
 
-def test_hd():
-    f20_ref = np.zeros([14, 14])
-    f20_ref[1, 1] = 1
-    f20_ref[9:12, 9:12] = 1
-    f20_pred = np.zeros([14, 14])
-    f20_pred[9:12, 9:12] = 1
-    bpm = PM(f20_pred, f20_ref, dict_args={"hd_perc": 95})
-    hausdorff_distance = bpm.measured_hausdorff_distance()
-    hausdorff_distance_perc = bpm.measured_hausdorff_distance_perc()
+def test_hausdorff_distance_5c():
+    """
+    Using figure 5c p14 as illustration for calculation of HD and HD95
+    """
+    hausdorff_distance = bpm5c.measured_hausdorff_distance()
+    hausdorff_distance_perc = bpm5c.measured_hausdorff_distance_perc()
     print(hausdorff_distance_perc)
     expected_hausdorff_distance = 11.31
     expected_hausdorff_distance_perc = 6.79

test/test_metrics/test_prob_pairwise_measures.py

@@ -14,7 +14,10 @@
 from MetricsReloaded.metrics.prob_pairwise_measures import ProbabilityPairwiseMeasures
 
 
-def test_auc():
+def test_auroc():
+    """
+    Based on SN2.18 p60 of Pitfalls paper
+    """
     ref = np.asarray([0, 0, 0, 1, 1, 1])
     pred_proba = np.asarray([0.21, 0.35, 0.63, 0.92, 0.32, 0.79])
     ppm = ProbabilityPairwiseMeasures(pred_proba, ref)
@@ -24,7 +27,10 @@ def test_auc():
     assert_allclose(value_test, expected_auc, atol=0.01)
 
 
-def test_ap():
+def test_average_precision():
+    """
+    Based on SN2.18 p60 of pitfalls paper
+    """
     ref = np.asarray([0, 0, 0, 1, 1, 1])
     pred_proba = np.asarray([0.21, 0.35, 0.63, 0.92, 0.32, 0.79])
     ppm = ProbabilityPairwiseMeasures(pred_proba, ref)