Repository: SALib/SALib

Issue:
Migrate to numpy>=2.0.0
Numpy 2.0.0 was released at June 16. It seems that SALib is incompatible with the new numpy now.

E.g., the quick start example, if run this after `pip install SALib --force-reinstall`
```python
from SALib.sample import saltelli
from SALib.analyze import sobol
from SALib.test_functions import Ishigami
import numpy as np

problem = {
  'num_vars': 3,
  'names': ['x1', 'x2', 'x3'],
  'bounds': [[-np.pi, np.pi]]*3
}

# Generate samples
param_values = saltelli.sample(problem, 1024)

# Run model (example)
Y = Ishigami.evaluate(param_values)

# Perform analysis
Si = sobol.analyze(problem, Y, print_to_console=True)
# Returns a dictionary with keys 'S1', 'S1_conf', 'ST', and 'ST_conf'
# (first and total-order indices with bootstrap confidence intervals)
```

will yield the following error message,
```
Traceback (most recent call last):
  File "...", line 19, in <module>
    Si = sobol.analyze(problem, Y, print_to_console=True)
  File "python3.9/site-packages/SALib/analyze/sobol.py", line 151, in analyze
    S["S1"][j] = first_order(A, AB[:, j], B)
  File ".../python3.9/site-packages/SALib/analyze/sobol.py", line 215, in first_order
    if y.ptp() == 0:
AttributeError: `ptp` was removed from the ndarray class in NumPy 2.0. Use np.ptp(arr, ...) instead.
```

Any plan on migrating to numpy>=2.0.0?

Discussion:

Patch:
diff --git a/src/SALib/analyze/delta.py b/src/SALib/analyze/delta.py
index f80cad7..783b603 100644
--- a/src/SALib/analyze/delta.py
+++ b/src/SALib/analyze/delta.py
@@ -137,13 +137,13 @@ def calc_delta(Y, Ygrid, X, m):
 
         # if not np.all(np.equal(Y_ix, Y_ix[0])):
         Y_ix = Y[ix]
-        if Y_ix.ptp() != 0.0:
+        if np.ptp(Y_ix) != 0.0:
             fyc = gaussian_kde(Y_ix, bw_method="silverman")(Ygrid)
             fy_ = np.abs(fy - fyc)
         else:
             fy_ = np.abs(fy)
 
-        d_hat += (nm / (2 * N)) * np.trapz(fy_, Ygrid)
+        d_hat += (nm / (2 * N)) * np.trapezoid(fy_, Ygrid)
 
     return d_hat
 
@@ -168,7 +168,7 @@ def bias_reduced_delta(Y, Ygrid, X, m, num_resamples, conf_level, y_resamples):
 def sobol_first(Y, X, m):
     # pre-process to catch constant array
     # see: https://github.com/numpy/numpy/issues/9631
-    if Y.ptp() == 0.0:
+    if np.ptp(Y) == 0.0:
         # Catch constant results
         # If Y does not change then it is not sensitive to anything...
         return 0.0
diff --git a/src/SALib/analyze/morris.py b/src/SALib/analyze/morris.py
index 010ece4..c749a02 100644
--- a/src/SALib/analyze/morris.py
+++ b/src/SALib/analyze/morris.py
@@ -245,7 +245,7 @@ def _compute_grouped_sigma(
 
     sigma = np.zeros(groups.shape[1], dtype=float)
     np.copyto(sigma, sigma_agg, where=groups.sum(axis=0) == 1)
-    np.copyto(sigma, np.NAN, where=groups.sum(axis=0) != 1)
+    np.copyto(sigma, np.nan, where=groups.sum(axis=0) != 1)
 
     return sigma
 
diff --git a/src/SALib/analyze/sobol.py b/src/SALib/analyze/sobol.py
index 67c587e..d38101b 100644
--- a/src/SALib/analyze/sobol.py
+++ b/src/SALib/analyze/sobol.py
@@ -154,7 +154,7 @@ def analyze(
             if keep_resamples:
                 S["S1_conf_all"][:, j] = S1_conf_j
 
-            var_diff = np.r_[A[r], B[r]].ptp()
+            var_diff = np.ptp(np.r_[A[r], B[r]])
             if var_diff != 0.0:
                 S["S1_conf"][j] = Z * S1_conf_j.std(ddof=1)
             else:
@@ -212,7 +212,7 @@ def first_order(A, AB, B):
     sample variance
     """
     y = np.r_[A, B]
-    if y.ptp() == 0:
+    if np.ptp(y) == 0:
         warn(CONST_RESULT_MSG)
         return np.array([0.0])
 
@@ -225,7 +225,7 @@ def total_order(A, AB, B):
     sample variance
     """
     y = np.r_[A, B]
-    if y.ptp() == 0:
+    if np.ptp(y) == 0:
         warn(CONST_RESULT_MSG)
         return np.array([0.0])
 
@@ -235,7 +235,7 @@ def total_order(A, AB, B):
 def second_order(A, ABj, ABk, BAj, B):
     """Second order estimator following Saltelli 2002"""
     y = np.r_[A, B]
-    if y.ptp() == 0:
+    if np.ptp(y) == 0:
         warn(CONST_RESULT_MSG)
         return np.array([0.0])
 
diff --git a/src/SALib/util/__init__.py b/src/SALib/util/__init__.py
index 335b858..72fbe0f 100644
--- a/src/SALib/util/__init__.py
+++ b/src/SALib/util/__init__.py
@@ -299,7 +299,7 @@ def compute_groups_matrix(groups: List):
         groups and a list of unique group names
     """
     num_vars = len(groups)
-    unique_group_names = pd.unique(groups)
+    unique_group_names = pd.unique(np.array(groups))
     number_of_groups = len(unique_group_names)
 
     indices = dict([(x, i) for (i, x) in enumerate(unique_group_names)])


Test patch:
diff --git a/src/SALib/test_functions/Sobol_G.py b/src/SALib/test_functions/Sobol_G.py
index 1c5194f..2aedd46 100644
--- a/src/SALib/test_functions/Sobol_G.py
+++ b/src/SALib/test_functions/Sobol_G.py
@@ -81,7 +81,7 @@ def evaluate(values, a=None, delta=None, alpha=None):
         mod_x = shift_of_x - integral
         temp_y = np.abs(2 * mod_x - 1) ** alpha
         y_elements = ((1 + alpha) * temp_y + a) / (1 + a)
-        Y[i] = np.prod(y_elements)
+        Y[i] = y_elements.prod()
 
     return Y
 
@@ -103,7 +103,7 @@ def _total_variance(a=None, alpha=None):
         alpha = np.ones_like(a)
 
     a = np.array(a)
-    return np.add(-1, np.product(1 + _partial_first_order_variance(a, alpha), axis=0))
+    return np.add(-1, np.prod(1 + _partial_first_order_variance(a, alpha), axis=0))
 
 
 def sensitivity_index(a, alpha=None):
@@ -116,6 +116,6 @@ def total_sensitivity_index(a, alpha=None):
 
     pv = _partial_first_order_variance(a, alpha)
     tv = _total_variance(a, alpha)
-    product_pv = np.product(1 + pv, axis=0)
+    product_pv = np.prod(1 + pv, axis=0)
 
     return np.divide(pv * np.divide(product_pv, 1 + pv.T), tv)
diff --git a/tests/sample/test_latin.py b/tests/sample/test_latin.py
index f8e6e6b..a4b9e5f 100644
--- a/tests/sample/test_latin.py
+++ b/tests/sample/test_latin.py
@@ -137,8 +137,8 @@ class TestLatinSample:
         # Group samples should have the same values
         # Get (max - min) with the `ptp()` method, the result of which should be
         # an array of zeros
-        diff = samples[:, ::2].ptp(axis=1)
+        diff = np.ptp(samples[:, ::2], axis=1)
         assert np.all(diff == 0), "Grouped samples do not have the same values"
 
-        diff = samples[:, 1::2].ptp(axis=1)
+        diff = np.ptp(samples[:, 1::2], axis=1)
         assert np.all(diff == 0), "Grouped samples do not have the same values"
diff --git a/tests/test_analyze_morris.py b/tests/test_analyze_morris.py
index e5b7760..eb059ee 100644
--- a/tests/test_analyze_morris.py
+++ b/tests/test_analyze_morris.py
@@ -931,7 +931,7 @@ def test_compute_grouped_sigma():
     sigma = np.std(ee, axis=1, ddof=1)
 
     actual = _compute_grouped_sigma(sigma, group_matrix)
-    desired = np.array([1.79352911, np.NAN], dtype=float)
+    desired = np.array([1.79352911, np.nan], dtype=float)
     assert_allclose(actual, desired, rtol=1e-1)
 
 
diff --git a/tests/test_regression.py b/tests/test_regression.py
index a987013..61229bf 100644
--- a/tests/test_regression.py
+++ b/tests/test_regression.py
@@ -118,7 +118,7 @@ class TestMorris:
 
         assert_allclose(Si["mu"], [9.786986, -9.938717e-13], atol=0, rtol=1e-5)
 
-        assert_allclose(Si["sigma"], [6.453729, np.NaN], atol=0, rtol=1e-5)
+        assert_allclose(Si["sigma"], [6.453729, np.nan], atol=0, rtol=1e-5)
 
         assert_allclose(Si["mu_star"], [9.786986, 7.875], atol=0, rtol=1e-5)
 
diff --git a/tests/test_to_df.py b/tests/test_to_df.py
index b024cb8..e42755a 100644
--- a/tests/test_to_df.py
+++ b/tests/test_to_df.py
@@ -1,7 +1,7 @@
 import numpy as np
 import pandas as pd
 from SALib.sample import (
-    saltelli,
+    sobol as sobol_sample,
     morris as morris_sample,
     finite_diff,
     fast_sampler,
@@ -60,7 +60,7 @@ def test_sobol_to_df():
     params = ["x1", "x2", "x3"]
     problem = {"num_vars": 3, "names": params, "bounds": [[-np.pi, np.pi]] * 3}
 
-    X = saltelli.sample(problem, 512)
+    X = sobol_sample.sample(problem, 512)
     Y = Ishigami.evaluate(X)
     Si = sobol.analyze(problem, Y, print_to_console=False)
     total, first, second = Si.to_df()