Merge branch 'main' into sprint19_MetropolishHastings_renaming

README.md

@@ -28,14 +28,14 @@ import numpy as np
 import matplotlib.pyplot as plt
 from cuqi.testproblem import Deconvolution2D
 from cuqi.data import grains
-from cuqi.distribution import Laplace_diff, Gaussian
+from cuqi.distribution import LMRF, Gaussian
 from cuqi.problem import BayesianProblem
 
 # Step 1: Model and data, y = Ax
 A, y_data, info = Deconvolution2D.get_components(dim=128, phantom=grains())
 
-# Step 2: Prior, x ~ Laplace_diff(0, 0.01)
-x = Laplace_diff(location=np.zeros(A.domain_dim),
+# Step 2: Prior, x ~ LMRF(0, 0.01)
+x = LMRF(location=np.zeros(A.domain_dim),
                  scale=0.01,
                  bc_type='neumann',
                  physical_dim=2)

cuqi/distribution/__init__.py

@@ -1,14 +1,13 @@
 from ._distribution import Distribution
 from ._beta import Beta
 from ._cauchy import Cauchy
-from ._cauchy_diff import Cauchy_diff
+from ._cmrf import CMRF
 from ._gamma import Gamma
 from ._gaussian import Gaussian, JointGaussianSqrtPrec
 from ._gmrf import GMRF
 from ._inverse_gamma import InverseGamma
-from ._laplace_diff import Laplace_diff
-from ._laplace import Laplace
 from ._lmrf import LMRF
+from ._laplace import Laplace
 from ._lognormal import Lognormal
 from ._normal import Normal
 from ._posterior import Posterior

cuqi/distribution/_cauchy_diff.py → cuqi/distribution/_cmrf.py

@@ -5,7 +5,7 @@
 from cuqi.operator import FirstOrderFiniteDifference
 from cuqi.geometry import _DefaultGeometry, Image2D, _get_identity_geometries
 
-class Cauchy_diff(Distribution):
+class CMRF(Distribution):
     """Cauchy distribution on the difference between neighboring nodes.
 
     For 1D `(physical_dim=1)`, the Cauchy difference distribution assumes that
@@ -42,7 +42,7 @@ class Cauchy_diff(Distribution):
 
         import cuqi
         import numpy as np
-        prior = cuqi.distribution.Cauchy_diff(location=np.zeros(128), scale=0.1)
+        prior = cuqi.distribution.CMRF(location=np.zeros(128), scale=0.1)
  
     """
 
@@ -60,7 +60,7 @@ def __init__(self, location, scale, bc_type="zero", physical_dim=1, **kwargs):
             num_nodes = (N, N)
             if isinstance(self.geometry, _DefaultGeometry):
                 self.geometry = Image2D(num_nodes)
-            print("Warning: 2D Cauchy_diff is still experimental. Use at own risk.")
+            print("Warning: 2D CMRF is still experimental. Use at own risk.")
         elif physical_dim == 1:
             num_nodes = self.dim
         else:
@@ -89,7 +89,7 @@ def _gradient(self, val, **kwargs):
             warnings.warn('Gradient not implemented for {}'.format(type(self.location)))
 
     def _sample(self,N=1,rng=None):
-        raise NotImplementedError("'Cauchy_diff.sample' is not implemented. Sampling can be performed with the 'sampler' module.")
+        raise NotImplementedError("'CMRF.sample' is not implemented. Sampling can be performed with the 'sampler' module.")
 
     # def cdf(self, x):   # TODO
     #     return 1/np.pi * np.atan((x-self.loc)/self.scale)

cuqi/distribution/_lmrf.py

@@ -1,44 +1,88 @@
-from cuqi.distribution import Distribution
-from cuqi.operator import FirstOrderFiniteDifference
 import numpy as np
+from cuqi.geometry import _DefaultGeometry, Image2D
+from cuqi.operator import FirstOrderFiniteDifference
+from cuqi.distribution import Distribution
 
 class LMRF(Distribution):
-    """
-        Parameters
-        ----------
-        partition_size : int
-            The dimension of the distribution in one physical dimension. 
+    """Laplace distribution on the difference between neighboring nodes.
+
+    For 1D `(physical_dim=1)`, the Laplace difference distribution assumes that
+
+    .. math::
+
+        x_i-x_{i-1} \sim \mathrm{Laplace}(0, b),
+
+    where :math:`b` is the scale parameter.
+
+    For 2D `(physical_dim=2)` the differences are defined in both horizontal and vertical directions.
+
+    It is possible to define boundary conditions using the `bc_type` parameter.
+
+    The location parameter is a shift of the :math:`\mathbf{x}`.
 
-        physical_dim : int
-            The physical dimension of what the distribution represents (can take the values 1 or 2).
+    Parameters
+    ----------
+    location : scalar or ndarray
+        The location parameter of the distribution.
+
+    scale : scalar
+        The scale parameter of the distribution.
+
+    bc_type : string
+        The boundary conditions of the difference operator.
+
+    physical_dim : int
+        The physical dimension of what the distribution represents (can take the values 1 or 2).
+
+    Example
+    -------
+    .. code-block:: python
+
+        import cuqi
+        import numpy as np
+        prior = cuqi.distribution.LMRF(location=np.zeros(128), scale=0.1)
+ 
     """
-        
-    def __init__(self, mean, prec, partition_size, physical_dim, bc_type, **kwargs):
-        super().__init__(**kwargs)
-        self.mean = mean.reshape(len(mean), 1)
-        self.prec = prec
-        self._partition_size = partition_size          # partition size
-        self._bc_type = bc_type      # boundary conditions
+    def __init__(self, location, scale, bc_type="zero", physical_dim=1, **kwargs):
+        # Init from abstract distribution class
+        super().__init__(**kwargs) 
+
+        self.location = location
+        self.scale = scale
+        self._bc_type = bc_type
         self._physical_dim = physical_dim
-        if physical_dim == 1: 
-            num_nodes = (partition_size,) 
+
+        if physical_dim == 2:
+            N = int(np.sqrt(self.dim))
+            num_nodes = (N, N)
+            if isinstance(self.geometry, _DefaultGeometry):
+                self.geometry = Image2D(num_nodes)
+
+        elif physical_dim == 1:
+            num_nodes = self.dim
         else:
-            num_nodes = (partition_size,partition_size)
+            raise ValueError("Only physical dimension 1 or 2 supported.")
+
+        self._diff_op = FirstOrderFiniteDifference(num_nodes=num_nodes, bc_type=bc_type)
 
-        self._diff_op = FirstOrderFiniteDifference( num_nodes, bc_type= bc_type) 
+        # Check if location parameter is non-zero vector (not supported)
+        if callable(location) or np.linalg.norm(location) > 0:
+            raise ValueError("Non-zero location parameter not supported.")
 
-    @property
-    def dim(self):
-        return self._diff_op.dim
+    def pdf(self, x):
+        Dx = self._diff_op @ (x-self.location)  # np.diff(X)
+        return (1/(2*self.scale))**(len(Dx)) * np.exp(-np.linalg.norm(Dx, ord=1, axis=0)/self.scale)
 
     def logpdf(self, x):
+        Dx = self._diff_op @ (x-self.location)
+        return len(Dx)*(-(np.log(2)+np.log(self.scale))) - np.linalg.norm(Dx, ord=1, axis=0)/self.scale
 
-        if self._physical_dim == 1 or self._physical_dim == 2:
-            const = self.dim *(np.log(self.prec)-np.log(2)) 
-            y = const -  self.prec*(np.linalg.norm(self._diff_op@x, ord=1))
-        else:
-            raise NotImplementedError
-        return y
+    def _sample(self,N=1,rng=None):
+        raise NotImplementedError("'LMRF.sample' is not implemented. Sampling can be performed with the 'sampler' module.")
+
+    # def cdf(self, x):   # TODO
+    #     return 1/2 + 1/2*np.sign(x-self.loc)*(1-np.exp(-np.linalg.norm(x, ord=1, axis=0)/self.scale))
 
-    def _sample(self, N):
-        raise NotImplementedError
+    # def sample(self):   # TODO
+    #     p = np.random.rand(self.dim)
+    #     return self.loc - self.scale*np.sign(p-1/2)*np.log(1-2*abs(p-1/2))

cuqi/problem/_problem.py

@@ -5,7 +5,7 @@
 
 import cuqi
 from cuqi import config
-from cuqi.distribution import Distribution, Gaussian, InverseGamma, Laplace_diff, GMRF, Lognormal, Posterior, LMRF, Beta, JointDistribution, Gamma, Cauchy_diff
+from cuqi.distribution import Distribution, Gaussian, InverseGamma, LMRF, GMRF, Lognormal, Posterior, Beta, JointDistribution, Gamma, CMRF
 from cuqi.density import Density
 from cuqi.model import LinearModel, Model
 from cuqi.likelihood import Likelihood
@@ -330,8 +330,8 @@ def sample_posterior(self, Ns, callback=None) -> cuqi.samples.Samples:
         elif hasattr(self.prior,"sqrtprecTimesMean") and hasattr(self.likelihood.distribution,"sqrtprec") and isinstance(self.model,LinearModel):
             return self._sampleLinearRTO(Ns, callback)
 
-        # For Laplace_diff we use our awesome unadjusted Laplace approximation!
-        elif self._check_posterior(self, Laplace_diff, Gaussian):
+        # For LMRF we use our awesome unadjusted Laplace approximation!
+        elif self._check_posterior(self, LMRF, Gaussian):
             return self._sampleUGLA(Ns, callback)
 
         # If we have gradients, use NUTS!
@@ -343,10 +343,6 @@ def sample_posterior(self, Ns, callback=None) -> cuqi.samples.Samples:
         elif self._check_posterior(self, (Gaussian, GMRF), Gaussian):
             return self._samplepCN(Ns, callback)
 
-        # For the remainder of valid cases we use CWMH
-        elif self._check_posterior(self, LMRF):
-            return self._sampleCWMH(Ns, callback)
-
         else:
             raise NotImplementedError(f"Automatic sampler choice is not implemented for model: {type(self.model)}, likelihood: {type(self.likelihood.distribution)} and prior: {type(self.prior)} and dim {self.prior.dim}. Manual sampler choice can be done via the 'sampler' module. Posterior distribution can be extracted via '.posterior' of any testproblem (BayesianProblem).")
 
@@ -581,7 +577,7 @@ def gradfunc(x): return -density.gradient(x)
             if disp: print("Optimizing with approximate gradients.") 
 
         # Compute point estimate
-        if self._check_posterior(self, Cauchy_diff, must_have_gradient=True): # Use L-BFGS-B for Cauchy_diff prior as it has better performance for this multi-modal posterior
+        if self._check_posterior(self, CMRF, must_have_gradient=True): # Use L-BFGS-B for CMRF prior as it has better performance for this multi-modal posterior
             if disp: print(f"Using scipy.optimize.L_BFGS_B on negative log of {density.__class__.__name__}")
             if disp: print("x0: ones vector")
             solver = cuqi.solver.L_BFGS_B(func, x0, gradfunc=gradfunc)
@@ -717,16 +713,16 @@ def _determine_sampling_strategy(self):
             if self._check_posterior(cond_target, Gamma, (Gaussian, GMRF)): 
                 sampling_strategy[par_name] = cuqi.sampler.Conjugate
 
-            # Gamma prior, Laplace_diff likelihood -> ConjugateApprox
-            elif self._check_posterior(cond_target, Gamma, Laplace_diff):
+            # Gamma prior, LMRF likelihood -> ConjugateApprox
+            elif self._check_posterior(cond_target, Gamma, LMRF):
                 sampling_strategy[par_name] = cuqi.sampler.ConjugateApprox
 
             # Gaussian prior, Gaussian likelihood, Linear model -> Linear_RTO
             elif self._check_posterior(cond_target, (Gaussian, GMRF), Gaussian, LinearModel):
                 sampling_strategy[par_name] = cuqi.sampler.Linear_RTO
 
-            # Laplace_diff prior, Gaussian likelihood, Linear model -> UGLA
-            elif self._check_posterior(cond_target, Laplace_diff, Gaussian, LinearModel):
+            # LMRF prior, Gaussian likelihood, Linear model -> UGLA
+            elif self._check_posterior(cond_target, LMRF, Gaussian, LinearModel):
                 sampling_strategy[par_name] = cuqi.sampler.UGLA
 
             else:

cuqi/sampler/_conjugate_approx.py

@@ -1,4 +1,4 @@
-from cuqi.distribution import Posterior, Laplace_diff, Gamma
+from cuqi.distribution import Posterior, LMRF, Gamma
 import numpy as np
 import scipy as sp
 
@@ -9,15 +9,15 @@ class ConjugateApprox: # TODO: Subclass from Sampler once updated
     by a conjugate pair.
 
     Currently supported pairs are:
-    - (Laplace_diff, Gamma): Approximated by (Gaussian, Gamma)
+    - (LMRF, Gamma): Approximated by (Gaussian, Gamma)
 
     For more information on conjugate pairs, see https://en.wikipedia.org/wiki/Conjugate_prior.
 
     """
 
 
     def __init__(self, target: Posterior):
-        if not isinstance(target.likelihood.distribution, Laplace_diff):
+        if not isinstance(target.likelihood.distribution, LMRF):
             raise ValueError("Conjugate sampler only works with Laplace diff likelihood function")
         if not isinstance(target.prior, Gamma):
             raise ValueError("Conjugate sampler only works with Gamma prior")
@@ -31,7 +31,7 @@ def step(self, x=None):
         D = self.target.likelihood.distribution._diff_op
         n = D.shape[0]
 
-        # Gaussian approximation of Laplace_diff prior as function of x_k
+        # Gaussian approximation of LMRF prior as function of x_k
         # See Uribe et al. (2022) for details
         # Current has a zero mean assumption on likelihood! TODO
         beta=1e-5