Merge pull request #102 from KaveIO/scipy_upgrade_to_qmvn

[fix]: migrate from scipy mvn to qmvn functions

Author: mbaak

.github/workflows/test_matrix.yml

@@ -14,7 +14,7 @@ jobs:
       fail-fast: false
       matrix:
         platform: [windows-latest, macos-latest, ubuntu-latest]
-        python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"]
+        python-version: ["3.9", "3.10", "3.11", "3.12"]
 
     runs-on: ${{ matrix.platform }}
 

.github/workflows/tests.yml

@@ -10,7 +10,7 @@ jobs:
     - name: Set up Python
       uses: actions/setup-python@v5
       with:
-        python-version: 3.8
+        python-version: 3.9
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip

CHANGES.rst

@@ -2,6 +2,14 @@
 Release notes
 =============
 
+Version 0.12.5, Jul 2025
+------------------------
+
+- FIX: scipy 1.16.0 no longer supports mvn, code now migrated to qmvn.
+  https://github.com/KaveIO/PhiK/issues/101
+  https://github.com/KaveIO/PhiK/pull/102
+- Drop support for Python 3.8, has reached end of life.
+
 Version 0.12.4, Jan 2024
 ------------------------
 

README.rst

@@ -2,7 +2,7 @@
 Phi_K Correlation Constant
 ==========================
 
-* Version: 0.12.4. Released: Jan 2024
+* Version: 0.12.5. Released: Jul 2025
 * Release notes: https://github.com/KaveIO/PhiK/blob/master/CHANGES.rst
 * Repository: https://github.com/kaveio/phik
 * Documentation: https://phik.readthedocs.io

phik/bivariate.py

@@ -17,24 +17,42 @@
 import warnings
 
 import numpy as np
+import scipy
 from scipy import optimize
-from scipy.stats._mvn import mvnun
 
+_scipy_version = [int(v) for v in scipy.__version__.split('.')]
+USE_QMVN = True if _scipy_version[0] >= 1 and _scipy_version[1] >= 16 else False
+if USE_QMVN:
+    from scipy.stats._qmvnt import _qauto, _qmvn
+else:
+    from scipy.stats._mvn import mvnun
 
-def _mvn_un(rho: float, lower: tuple, upper: tuple) -> float:
+
+
+
+def _mvn_un(rho: float, lower: tuple, upper: tuple,
+            rng: np.random.Generator = np.random.default_rng(42)) -> float:
     """Perform integral of bivariate normal gauss with correlation
 
     Integral is performed using scipy's mvn library.
 
     :param float rho: tilt parameter
     :param tuple lower: tuple of lower corner of integral area
     :param tuple upper: tuple of upper corner of integral area
+    :param np.random.Generator rng: default_rng(42), optional
     :returns float: integral value
     """
     mu = np.array([0.0, 0.0])
     S = np.array([[1.0, rho], [rho, 1.0]])
-    p, i = mvnun(lower, upper, mu, S)
-    return p
+    return _calc_mvnun(lower=lower, upper=upper, mu=mu, S=S, rng=rng)
+
+
+def _calc_mvnun(lower, upper, mu, S, rng = np.random.default_rng(42)):
+    if USE_QMVN:
+        res = _qauto(_qmvn, S, lower, upper, rng)[0]
+    else:
+        res = mvnun(lower, upper, mu, S)[0]
+    return res
 
 
 def _mvn_array(rho: float, sx: np.ndarray, sy: np.ndarray) -> list:
@@ -55,7 +73,7 @@ def _mvn_array(rho: float, sx: np.ndarray, sy: np.ndarray) -> list:
     mu = np.array([0.0, 0.0])
     S = np.array([[1.0, rho], [rho, 1.0]])
 
-    # callling mvn.mvnun is expansive, so we only calculate half of the matrix, then symmetrize
+    # callling mvn.mvnun is expensive, so we only calculate half of the matrix, then symmetrize
     # add half block, which is symmetric in x
     odd_odd = False
     ranges = [
@@ -81,10 +99,6 @@ def _mvn_array(rho: float, sx: np.ndarray, sy: np.ndarray) -> list:
     return corr
 
 
-def _calc_mvnun(lower, upper, mu, S):
-    return mvnun(lower, upper, mu, S)[0]
-
-
 def bivariate_normal_theory(
     rho: float,
     nx: int = -1,

phik/simcore/__init__.py

@@ -1,6 +1,10 @@
 import importlib.util
 
-_ext_spec = importlib.util.find_spec("phik.lib._phik_simulation_core")
+try:
+    _ext_spec = importlib.util.find_spec("phik.lib._phik_simulation_core")
+except ModuleNotFoundError:
+    _ext_spec = None
+
 if _ext_spec is not None:
     from phik.lib._phik_simulation_core import _sim_2d_data_patefield
 

pyproject.toml

@@ -4,16 +4,15 @@ build-backend = "scikit_build_core.build"
 
 [project]
 name = "phik"
-version = "0.12.4"
+version = "0.12.5"
 description = "Phi_K correlation analyzer library"
 readme = "README.rst"
 authors = [{ name = "KPMG N.V. The Netherlands", email = "[email protected]" }]
-requires-python = ">=3.8"
+requires-python = ">=3.9"
 classifiers = [
   "Development Status :: 4 - Beta",
   "License :: OSI Approved :: MIT License",
   "Programming Language :: Python :: 3 :: Only",
-  "Programming Language :: Python :: 3.8",
   "Programming Language :: Python :: 3.9",
   "Programming Language :: Python :: 3.10",
   "Programming Language :: Python :: 3.11",

tests/test_phik.py

@@ -20,7 +20,7 @@
 import pandas as pd
 import numpy as np
 from phik import resources, bivariate
-from phik.simulation import sim_2d_data_patefield
+from phik.simulation import sim_2d_data_patefield, CPP_SUPPORT
 from phik.binning import auto_bin_data, bin_data
 from phik.phik import phik_observed_vs_expected_from_rebinned_df, phik_from_hist2d
 from phik.statistics import get_dependent_frequency_estimates
@@ -51,7 +51,7 @@ def test_phik_from_hist2d(self):
         observed = df[cols].hist2d(interval_cols=interval_cols)
 
         phik_value = phik_from_hist2d(observed)
-        self.assertAlmostEqual(phik_value, 0.7685888294891855)
+        self.assertAlmostEqual(phik_value, 0.7685888294891855, places=3)
 
     def test_phik_observed_vs_expected_from_hist2d(self):
         """Test the calculation of Phi_K value from hist2d"""
@@ -67,7 +67,7 @@ def test_phik_observed_vs_expected_from_hist2d(self):
         expected = get_dependent_frequency_estimates(observed)
 
         phik_value = phik_from_hist2d(observed=observed, expected=expected)
-        self.assertAlmostEqual(phik_value, 0.7685888294891855)
+        self.assertAlmostEqual(phik_value, 0.7685888294891855, places=3)
 
     def test_phik_matrix(self):
         """Test the calculation of Phi_K"""
@@ -79,29 +79,25 @@ def test_phik_matrix(self):
         interval_cols = ["driver_age", "mileage"]
         phik_corr = df.phik_matrix(interval_cols=interval_cols)
 
-        self.assertTrue(
-            np.isclose(
-                phik_corr.values[cols.index("car_color"), cols.index("area")],
-                0.5904561614620166,
-            )
+        self.assertAlmostEqual(
+            phik_corr.values[cols.index("car_color"), cols.index("area")],
+            0.5904561614620166,
+            places=3,
         )
-        self.assertTrue(
-            np.isclose(
-                phik_corr.values[cols.index("area"), cols.index("car_color")],
-                0.5904561614620166,
-            )
+        self.assertAlmostEqual(
+            phik_corr.values[cols.index("area"), cols.index("car_color")],
+            0.5904561614620166,
+            places=3,
         )
-        self.assertTrue(
-            np.isclose(
-                phik_corr.values[cols.index("mileage"), cols.index("car_size")],
-                0.768588987856336,
-            )
+        self.assertAlmostEqual(
+            phik_corr.values[cols.index("mileage"), cols.index("car_size")],
+            0.768588987856336,
+            places=3,
         )
-        self.assertTrue(
-            np.isclose(
-                phik_corr.values[cols.index("car_size"), cols.index("mileage")],
-                0.768588987856336,
-            )
+        self.assertAlmostEqual(
+            phik_corr.values[cols.index("car_size"), cols.index("mileage")],
+            0.768588987856336,
+            places=3,
         )
 
     def test_phik_matrix_observed_vs_expected(self):
@@ -154,9 +150,9 @@ def test_global_phik(self):
         car_size = (np.where(gk[1] == "car_size"))[0][0]
         mileage = (np.where(gk[1] == "mileage"))[0][0]
 
-        self.assertTrue(np.isclose(gk[0][area][0], 0.6057528003711345))
-        self.assertTrue(np.isclose(gk[0][car_size][0], 0.76858883))
-        self.assertTrue(np.isclose(gk[0][mileage][0], 0.768588987856336))
+        self.assertAlmostEqual(gk[0][area][0], 0.6057528003711345, places=3)
+        self.assertAlmostEqual(gk[0][car_size][0], 0.76858883, places=3)
+        self.assertAlmostEqual(gk[0][mileage][0], 0.768588987856336, places=3)
 
     def test_significance_matrix_asymptotic(self):
         """Test significance calculation"""
@@ -313,6 +309,7 @@ def test_outlier_significance_matrices(self):
 
         self.assertTrue(isinstance(om, dict))
 
+    @pytest.mark.skipif(not CPP_SUPPORT, reason="cpp not supported")
     def test_simulation_2d_patefield(self):
         """Test simulation code using patefield algorithm."""
         og_state = np.random.get_state()