Repository: CUQI-DTU/CUQIpy

Issue:
Bump up scipy version requirement
**General description of task**

1. The recently introduced `isotonic_regression` (in Regularized LinearRTO) requires scipy>=1.12, while the current requirements.txt specifies scipy<1.13. As a result, the only version(s) of scipy that is compatible with CUQIpy now is 1.12.*. This is rather restrictive, particularly when working with other packages.
2. The latest arviz (0.21.0) doesn't work under the environment of Python 3.9 and scipy 1.13.
3. Functions, like `asfarray`, `(sparse matrix).A` and `.NaN`, have been removed from numpy 2.0.
4. Unit test `test_TP_BayesianProblem_sample` in `tests/test_bayesian_inversion.py` compares sample statistics with those created under scipy 1.12. However, the default optimization behaviour has changed from 1.13, so `RegularizedLinearRTO` will behave differently under scipy 1.12 and >=1.13. In the end, the unit test will fail under scipy>=1.13.

**Definition of Done (Feature/change)**
- [x] For Python<=3.9, we stick with scipy 1.12 and for Python 3.10, 1.11 and 3.12, we remove scipy's version requirement. (to solve 1 and 2)
- [x] Raise error if "increasing" or "decreasing" type of constraints is specified when `isotonic_regression` doesn't exist in the installed scipy (to solve 1)
- [x] Update codes involving removed functions from numpy and scipy. (to solve 1 and 3)
- [x] Skip regularized priors in `test_TP_BayesianProblem_sample` if installed scipy>=1.13. (to solve 4) 

**Definition of Done (Mandatory)**
- [x] Documentation added (docstrings on all public methods/classes)
- [x] Unit tests added/updated (and passing!)
- [x] Code reviewed and approved by 2 CUQI-maintainers
- [x] Online documentation is rendered correctly (esp. math) 

**Defintion of Done (Optional)**
- [x] HowTo/tutorial added

Discussion:
@chaozg: Hi @jakobsj , could you help review the dod here?@jakobsj: Thanks @chaozg these will be great fixes to include. For DoD can you please: 
1. For the isotonic, I agree with your proposed solution to omit access to this very niche feature from CUQIpy (replace by an error message or similar) to allow installing CUQIpy with more versions of scipy. Please add this to the description and DoD.
2. Please make sure that the DoD lists solutions to all the points/problems mentioned under "general description of task" for example the asfarray point is not mentioned in DoD.
@chaozg: Hi @jakobsj thanks for your suggestion. Could you have a look at the DoD again?@jakobsj: Thank you @chaozg, DoD approved!

Solution patch:
diff --git a/cuqi/distribution/__init__.py b/cuqi/distribution/__init__.py
index 3a51f5472..4842e2cb4 100644
--- a/cuqi/distribution/__init__.py
+++ b/cuqi/distribution/__init__.py
@@ -14,6 +14,6 @@
 from ._normal import Normal
 from ._truncated_normal import TruncatedNormal
 from ._posterior import Posterior
-from ._uniform import Uniform
+from ._uniform import Uniform, UnboundedUniform
 from ._custom import UserDefinedDistribution, DistributionGallery
 from ._joint_distribution import JointDistribution, _StackedJointDistribution, MultipleLikelihoodPosterior
diff --git a/cuqi/distribution/_truncated_normal.py b/cuqi/distribution/_truncated_normal.py
index 3f2ade3c5..51e97eb4d 100644
--- a/cuqi/distribution/_truncated_normal.py
+++ b/cuqi/distribution/_truncated_normal.py
@@ -107,7 +107,7 @@ def _gradient(self, x, *args, **kwargs):
         """
         # check if x falls in the range between np.array a and b
         if np.any(x < self.low) or np.any(x > self.high):
-            return np.NaN*np.ones_like(x)
+            return np.nan*np.ones_like(x)
         else:
             return self._normal.gradient(x, *args, **kwargs)
 
diff --git a/cuqi/distribution/_uniform.py b/cuqi/distribution/_uniform.py
index a61ae71e6..87d3054b1 100644
--- a/cuqi/distribution/_uniform.py
+++ b/cuqi/distribution/_uniform.py
@@ -1,5 +1,6 @@
 import numpy as np
 from cuqi.distribution import Distribution
+from cuqi.geometry import Geometry
 
 class Uniform(Distribution):
 
@@ -46,7 +47,7 @@ def gradient(self, x):
         Computes the gradient of logpdf at the given values of x.
         """
         if np.any(x < self.low) or np.any(x > self.high):
-            return np.NaN*np.ones_like(x)
+            return np.nan*np.ones_like(x)
         else:
             return np.zeros_like(x)
 
@@ -57,4 +58,37 @@ def _sample(self,N=1, rng=None):
         else:
             s = np.random.uniform(self.low, self.high, (N,self.dim)).T
 
-        return s
\ No newline at end of file
+        return s
+
+class UnboundedUniform(Distribution):
+    """
+    Unbounded uniform distribution. This is a special case of the
+    Uniform distribution, where the lower and upper bounds are set to
+    -inf and inf, respectively. This distribution is not normalizable,
+    and therefore cannot be sampled from. It is mainly used for
+    initializing non-informative priors.
+    Parameters
+    ----------
+    geometry : int or Geometry
+        The geometry of the distribution. If an integer is given, it is
+        interpreted as the dimension of the distribution. If a
+        Geometry object is given, its par_dim attribute is used.
+    """
+    def __init__(self, geometry, is_symmetric=True, **kwargs):
+        super().__init__(geometry=geometry, is_symmetric=is_symmetric, **kwargs) 
+
+    def logpdf(self, x):
+        """
+        Evaluate the logarithm of the unnormalized PDF at the given values of x.
+        """
+        # Always return 1.0 (the unnormalized log PDF)
+        return 1.0
+
+    def gradient(self, x):
+        """
+        Computes the gradient of logpdf at the given values of x.
+        """
+        return np.zeros_like(x)
+
+    def _sample(self, N=1, rng=None):
+        raise NotImplementedError("Cannot sample from UnboundedUniform distribution")
\ No newline at end of file
diff --git a/cuqi/implicitprior/_regularizedGaussian.py b/cuqi/implicitprior/_regularizedGaussian.py
index a768ccf57..38a493aee 100644
--- a/cuqi/implicitprior/_regularizedGaussian.py
+++ b/cuqi/implicitprior/_regularizedGaussian.py
@@ -189,12 +189,18 @@ def _parse_preset_constraint_input(self, constraint, optional_regularization_par
             elif c_lower == "increasing":
                 if not isinstance(self.geometry, Continuous1D):
                     raise ValueError("Geometry not supported for " + c_lower)
-                self._constraint_prox = lambda z, _: spoptimize.isotonic_regression(z, increasing=True).x
+                if hasattr(spoptimize, 'isotonic_regression'):
+                    self._constraint_prox = lambda z, _: spoptimize.isotonic_regression(z, increasing=True).x
+                else:
+                    raise AttributeError(f"The function 'isotonic_regression' does not exist in scipy.optimize. Installed scipy version: {spoptimize.__version__}. You need to install a scipy >= 1.12.0")
                 self._preset["constraint"] = "increasing"
             elif c_lower == "decreasing":
                 if not isinstance(self.geometry, Continuous1D):
                     raise ValueError("Geometry not supported for " + c_lower)
-                self._constraint_prox = lambda z, _: spoptimize.isotonic_regression(z, increasing=False).x
+                if hasattr(spoptimize, 'isotonic_regression'):
+                    self._constraint_prox = lambda z, _: spoptimize.isotonic_regression(z, increasing=False).x
+                else:
+                    raise AttributeError(f"The function 'isotonic_regression' does not exist in scipy.optimize. Installed scipy version: {spoptimize.__version__}. You need to install a scipy >= 1.12.0")
                 self._preset["constraint"] = "decreasing"
             elif c_lower == "convex":
                 if not isinstance(self.geometry, Continuous1D):
diff --git a/cuqi/utilities/_utilities.py b/cuqi/utilities/_utilities.py
index d43639dca..a246d2865 100644
--- a/cuqi/utilities/_utilities.py
+++ b/cuqi/utilities/_utilities.py
@@ -188,7 +188,7 @@ def approx_derivative(func, wrt, direction=None, epsilon=np.sqrt(np.finfo(float)
     # We compute the Jacobian matrix of func using forward differences.
     # If the function is scalar-valued, we compute the gradient instead.
     # If the direction is provided, we compute the direction-Jacobian product.
-    wrt = np.asfarray(wrt)
+    wrt = np.asarray(wrt)
     f0 = func(wrt)
     Matr = np.zeros([infer_len(wrt), infer_len(f0)])
     dx = np.zeros(len(wrt))
diff --git a/requirements.txt b/requirements.txt
index 3ae2538d6..540fe9b91 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,5 +1,6 @@
 matplotlib
 numpy>=1.17.0
-scipy<1.13
+scipy==1.12.0; python_version <= '3.9'
+scipy; python_version > '3.9'
 arviz
 tqdm


Test patch:
diff --git a/tests/test_bayesian_inversion.py b/tests/test_bayesian_inversion.py
index 97d27e854..4ea132186 100644
--- a/tests/test_bayesian_inversion.py
+++ b/tests/test_bayesian_inversion.py
@@ -2,7 +2,9 @@
 
 import pytest
 import numpy as np
+import scipy
 import sys
+from packaging import version
 
 from cuqi.testproblem import Deconvolution1D
 from cuqi.distribution import Gaussian, GMRF, CMRF, LMRF, Gamma
@@ -25,6 +27,14 @@ def test_TP_BayesianProblem_sample(copy_reference, TP_type, phantom, prior, Ns,
     # SKIP NUTS test if not windows (for now)
     if isinstance(prior, CMRF) and not sys.platform.startswith('win'):
         pytest.skip("NUTS(CMRF) regression test is not implemented for this platform")
+    elif (
+        isinstance(prior, RegularizedGaussian) 
+        or isinstance(prior, RegularizedGMRF)
+    ) and version.parse(scipy.__version__) >= version.parse('1.13.0'):
+        pytest.skip(
+            "Saved reference data is not consistent with scipy>=1.13.0 "
+            "(due to different default behaviour of scipy.optimize)"
+        )
 
     np.random.seed(19937)
 
diff --git a/tests/test_distribution.py b/tests/test_distribution.py
index 358d4a265..961afbad3 100644
--- a/tests/test_distribution.py
+++ b/tests/test_distribution.py
@@ -203,7 +203,6 @@ def test_Uniform_logpdf(low, high, toeval, expected):
     UD = cuqi.distribution.Uniform(low, high)
     assert np.allclose(UD.logpdf(toeval), expected) 
 
-
 @pytest.mark.parametrize("low,high,expected",[(np.array([1, .5]), 
                                                np.array([2, 2.5]),
                                                np.array([[1.5507979 , 1.29090474, 1.89294695],
@@ -217,6 +216,19 @@ def test_Uniform_sample(low, high, expected):
     print(cuqi_samples)
     assert np.allclose(cuqi_samples.samples, expected) 
 
+def test_UnboundedUniform_logpdf_and_gradient():
+    x1 = cuqi.distribution.UnboundedUniform(geometry=1)
+    assert np.allclose(x1.logpdf(np.array([1])), 1.0)
+    assert np.allclose(x1.gradient(np.array([0])), np.array([0]))
+    x2 = cuqi.distribution.UnboundedUniform(geometry=cuqi.geometry.Discrete(2))
+    assert np.allclose(x2.logpdf(np.array([1, 2])), 1.0)
+    assert np.allclose(x2.gradient(np.array([0, 0])), np.array([0, 0]))
+
+def test_UnboundedUniform_sample():
+    x = cuqi.distribution.UnboundedUniform(geometry=2)
+    with pytest.raises(NotImplementedError, match="Cannot sample from UnboundedUniform distribution"):
+        x.sample()
+
 @pytest.mark.parametrize("distribution, kwargs",
                          [(cuqi.distribution.Uniform, 
                           {'low':np.array([2, 2.5, 3,5]),
diff --git a/tests/test_distributions_shape.py b/tests/test_distributions_shape.py
index 7d0f69d1c..7a5b13780 100644
--- a/tests/test_distributions_shape.py
+++ b/tests/test_distributions_shape.py
@@ -13,6 +13,7 @@
     "Posterior",
     "Distribution",
     "JointGaussianSqrtPrec",
+    "UnboundedUniform",
 ]
 
 # Define cases to skip (these are TODO)
diff --git a/tests/test_implicit_priors.py b/tests/test_implicit_priors.py
index ed07f5896..2a8b4efa6 100644
--- a/tests/test_implicit_priors.py
+++ b/tests/test_implicit_priors.py
@@ -1,6 +1,8 @@
 import cuqi
 import numpy as np
+import scipy
 import pytest
+from packaging import version
 
 def test_RegularizedGaussian_default_init():
     """ Test that the implicit regularized Gaussian requires at least 1 regularization argument """
@@ -186,6 +188,7 @@ def test_RegularizedGaussian_double_preset():
     assert len(x.proximal[0]) == 2
     assert len(x.proximal[1]) == 2
 
+@pytest.mark.skipif(version.parse(scipy.__version__) < version.parse('1.12.0'), reason="isotonic_regression in not available in scipy < 1.12.0")
 def test_regression_increasing():
     """ Regression test for the increasing constraints in RegularizedGaussian"""
 
diff --git a/tests/test_model.py b/tests/test_model.py
index 0bb44beb4..77fd6cd43 100644
--- a/tests/test_model.py
+++ b/tests/test_model.py
@@ -19,7 +19,7 @@ def test_LinearModel_getMatrix(seed):
     mat1 = model1.get_matrix() #Normal matrix
     mat2 = model2.get_matrix() #Sparse matrix (generated from functions)
 
-    assert np.allclose(mat1,mat2.A)
+    assert np.allclose(mat1,mat2.toarray())
 
 def test_initialize_model_dim():
     model1 = cuqi.model.Model(lambda x:x, range_geometry=4, domain_geometry=4)