Merge pull request #618 from jeverink/RegularizedGaussian-Additional-Constraints

Regularized Gaussian - additional constraints

Author: chaozg

cuqi/experimental/mcmc/_conjugate.py

@@ -4,7 +4,7 @@
 from cuqi.experimental.mcmc import Sampler
 from cuqi.distribution import Posterior, Gaussian, Gamma, GMRF, ModifiedHalfNormal
 from cuqi.implicitprior import RegularizedGaussian, RegularizedGMRF, RegularizedUnboundedUniform
-from cuqi.utilities import get_non_default_args, count_nonzero, count_constant_components_1D, count_constant_components_2D
+from cuqi.utilities import get_non_default_args, count_nonzero, count_within_bounds, count_constant_components_1D, count_constant_components_2D, piecewise_linear_1D_DoF
 from cuqi.geometry import Continuous1D, Continuous2D, Image2D
 
 class Conjugate(Sampler):
@@ -17,8 +17,8 @@ class Conjugate(Sampler):
     - (GMRF, Gamma) where Gamma is defined on the precision parameter of the GMRF
     - (RegularizedGaussian, Gamma) with preset constraints only and Gamma is defined on the precision parameter of the RegularizedGaussian
     - (RegularizedGMRF, Gamma) with preset constraints only and Gamma is defined on the precision parameter of the RegularizedGMRF
-    - (RegularizedGaussian, ModifiedHalfNormal) with preset constraints and regularization only
-    - (RegularizedGMRF, ModifiedHalfNormal) with preset constraints and regularization only
+    - (RegularizedGaussian, ModifiedHalfNormal) with most of the preset constraints and regularization
+    - (RegularizedGMRF, ModifiedHalfNormal) with most of the preset constraints and regularization
 
     Currently the Gamma and ModifiedHalfNormal distribution must be univariate.
 
@@ -147,8 +147,8 @@ def validate_target(self):
         if self.target.prior.dim != 1:
             raise ValueError("RegularizedGaussian-Gamma conjugacy only works with univariate ModifiedHalfNormal prior")
 
-        if self.target.likelihood.distribution.preset["constraint"] not in ["nonnegativity"]:
-            raise ValueError("RegularizedGaussian-Gamma conjugacy only works with implicit regularized Gaussian likelihood with nonnegativity constraints")
+        # Raises error if preset is not supported
+        _compute_sparsity_level(self.target)
 
         key_value_pairs = _get_conjugate_parameter(self.target)
         if len(key_value_pairs) != 1:
@@ -166,7 +166,7 @@ def validate_target(self):
     def conjugate_distribution(self):
         # Extract variables
         b = self.target.likelihood.data                                 # mu
-        m = np.count_nonzero(b)                                         # n
+        m = _compute_sparsity_level(self.target)
         Ax = self.target.likelihood.distribution.mean                   # x_i
         L = self.target.likelihood.distribution(np.array([1])).sqrtprec # L
         alpha = self.target.prior.shape                                 # alpha
@@ -183,9 +183,9 @@ def validate_target(self):
         if self.target.prior.dim != 1:
             raise ValueError("RegularizedUnboundedUniform-Gamma conjugacy only works with univariate Gamma prior")
 
-        if self.target.likelihood.distribution.preset["regularization"] not in ["l1", "tv"]:
-            raise ValueError("RegularizedUnboundedUniform-Gamma conjugacy only works with implicit regularized Gaussian likelihood with l1 or tv regularization")
-        
+        # Raises error if preset is not supported
+        _compute_sparsity_level(self.target)
+
         key_value_pairs = _get_conjugate_parameter(self.target)
         if len(key_value_pairs) != 1:
             raise ValueError(f"Multiple references to conjugate parameter {self.target.prior.name} found in likelihood. Only one occurance is supported.")
@@ -219,8 +219,8 @@ def validate_target(self):
         if self.target.prior.dim != 1:
             raise ValueError("RegularizedGaussian-ModifiedHalfNormal conjugacy only works with univariate ModifiedHalfNormal prior")
 
-        if self.target.likelihood.distribution.preset["regularization"] not in ["l1", "tv"]:
-            raise ValueError("RegularizedGaussian-ModifiedHalfNormal conjugacy only works with implicit regularized Gaussian likelihood with l1 or tv regularization")
+        # Raises error if preset is not supported
+        _compute_sparsity_level(self.target)
 
         key_value_pairs = _get_conjugate_parameter(self.target)
         if len(key_value_pairs) != 2:
@@ -266,23 +266,74 @@ def conjugate_distribution(self):
 
 
 def _compute_sparsity_level(target):
-    """Computes the sparsity level in accordance with Section 4 from [2],"""
+    """Computes the sparsity level in accordance with Section 4 from [2],
+    this can be interpreted as the number of degrees of freedom, that is,
+    the number of components n minus the dimension the of the subdifferential of the regularized.
+    """
     x = target.likelihood.data
-    if target.likelihood.distribution.preset["constraint"] == "nonnegativity":
-        if target.likelihood.distribution.preset["regularization"] == "l1":
-            m = count_nonzero(x)
-        elif target.likelihood.distribution.preset["regularization"] == "tv" and isinstance(target.likelihood.distribution.geometry, Continuous1D):
-            m = count_constant_components_1D(x, lower = 0.0)
-        elif target.likelihood.distribution.preset["regularization"] == "tv" and isinstance(target.likelihood.distribution.geometry, (Continuous2D, Image2D)):
-            m = count_constant_components_2D(target.likelihood.distribution.geometry.par2fun(x), lower = 0.0)
-    else: # No constraints, only regularization
-        if target.likelihood.distribution.preset["regularization"] == "l1":
-            m = count_nonzero(x)
-        elif target.likelihood.distribution.preset["regularization"] == "tv" and isinstance(target.likelihood.distribution.geometry, Continuous1D):
-            m = count_constant_components_1D(x)
-        elif target.likelihood.distribution.preset["regularization"] == "tv" and isinstance(target.likelihood.distribution.geometry, (Continuous2D, Image2D)):
-            m = count_constant_components_2D(target.likelihood.distribution.geometry.par2fun(x))
-    return m
+
+    constraint = target.likelihood.distribution.preset["constraint"]
+    regularization = target.likelihood.distribution.preset["regularization"]
+
+    # There is no reference for some of these conjugacy rules
+    if constraint == "nonnegativity":
+        if regularization in [None, "l1"]:
+            # Number of non-zero components in x
+            return count_nonzero(x)
+        elif regularization == "tv" and isinstance(target.likelihood.distribution.geometry, Continuous1D):
+            # Number of non-zero constant components in x
+            return count_constant_components_1D(x, lower = 0.0)
+        elif regularization == "tv" and isinstance(target.likelihood.distribution.geometry, (Continuous2D, Image2D)):
+            # Number of non-zero constant components in x
+            return count_constant_components_2D(target.likelihood.distribution.geometry.par2fun(x), lower = 0.0)
+    elif constraint == "box":
+        bounds = target.likelihood.distribution._box_bounds
+        if regularization is None:
+            # Number of components in x that are strictly between the lower and upper bound
+            return count_within_bounds(x, bounds[0], bounds[1])
+        elif regularization == "l1":
+            # Number of components in x that are strictly between the lower and upper bound and are not zero
+            return count_within_bounds(x, bounds[0], bounds[1], exception = 0.0)
+        elif regularization == "tv" and isinstance(target.likelihood.distribution.geometry, Continuous1D):
+            # Number of constant components in x between are strictly between the lower and upper bound
+            return count_constant_components_1D(x, lower = bounds[0], upper = bounds[1])
+        elif regularization == "tv" and isinstance(target.likelihood.distribution.geometry, (Continuous2D, Image2D)):
+            # Number of constant components in x between are strictly between the lower and upper bound
+            return count_constant_components_2D(target.likelihood.distribution.geometry.par2fun(x), lower = bounds[0], upper = bounds[1])
+    elif constraint in ["increasing", "decreasing"]:
+        if regularization is None:
+            # Number of constant components in x
+            return count_constant_components_1D(x)
+        elif regularization == "l1":
+            # Number of constant components in x that are not zero
+            return count_constant_components_1D(x, exception = 0.0)
+        elif regularization == "tv" and isinstance(target.likelihood.distribution.geometry, Continuous1D):
+            # Number of constant components in x
+            return count_constant_components_1D(x)
+        # Increasing and decreasing cannot be done in 2D
+    elif constraint in ["convex", "concave"]:
+        if regularization is None:
+            # Number of piecewise linear components in x
+            return piecewise_linear_1D_DoF(x)
+        elif regularization == "l1":
+            # Number of piecewise linear components in x that are not zero
+            return piecewise_linear_1D_DoF(x, exception_zero = True)
+        elif regularization == "tv" and isinstance(target.likelihood.distribution.geometry, Continuous1D):
+            # Number of piecewise linear components in x that are not flat
+            return piecewise_linear_1D_DoF(x, exception_flat = True)
+        # convex and concave has only been implemented in 1D
+    elif constraint == None: 
+        if regularization == "l1":
+            # Number of non-zero components in x
+            return count_nonzero(x)
+        elif regularization == "tv" and isinstance(target.likelihood.distribution.geometry, Continuous1D):
+            # Number of non-zero constant components in x
+            return count_constant_components_1D(x)
+        elif regularization == "tv" and isinstance(target.likelihood.distribution.geometry, (Continuous2D, Image2D)):
+            # Number of non-zero constant components in x
+            return count_constant_components_2D(target.likelihood.distribution.geometry.par2fun(x))
+
+    raise ValueError("RegularizedGaussian preset constraint and regularization choice is currently not supported with conjugacy.")
 
 
 def _get_conjugate_parameter(target):

cuqi/implicitprior/_regularizedGaussian.py

@@ -2,10 +2,12 @@
 from cuqi.distribution import Distribution, Gaussian
 from cuqi.solver import ProjectNonnegative, ProjectBox, ProximalL1
 from cuqi.geometry import Continuous1D, Continuous2D, Image2D
-from cuqi.operator import FirstOrderFiniteDifference, Operator
+from cuqi.operator import FirstOrderFiniteDifference, SecondOrderFiniteDifference, Operator
 
 import numpy as np
 import scipy.sparse as sparse
+import scipy.optimize as spoptimize
+
 from copy import copy
 
 
@@ -60,12 +62,18 @@ class RegularizedGaussian(Distribution):
         min_(z in C) 0.5||x-z||_2^2.
 
     constraint : string or None
-        Preset constraints that generate the corresponding proximal parameter. Can be set to "nonnegativity" and "box". Required for use in Gibbs.
-        For "box", the following additional parameters can be passed:
+        Preset constraints that generate the corresponding proximal parameter. Required for use in Gibbs. For any geometry the following can be chosen:
+        - "nonnegativity"
+        - "box", the following additional parameters can be passed:
             lower_bound : array_like or None
                 Lower bound of box, defaults to zero
             upper_bound : array_like
                 Upper bound of box, defaults to one
+        Additionally, for Continuous1D geometry the following options can be chosen:
+        - "increasing"
+        - "decreasing"
+        - "convex"
+        - "concave"
 
     regularization : string or None
         Preset regularization that generate the corresponding proximal parameter. Can be set to "l1" or 'tv'. Required for use in Gibbs in future update.
@@ -178,6 +186,28 @@ def _parse_preset_constraint_input(self, constraint, optional_regularization_par
                 self._proximal = lambda z, _: ProjectBox(z, _box_lower, _box_upper)
                 self._box_bounds = (np.ones(self.dim)*_box_lower, np.ones(self.dim)*_box_upper)
                 self._preset["constraint"] = "box"
+            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
+                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
+                self._preset["constraint"] = "decreasing"
+            elif c_lower == "convex":
+                if not isinstance(self.geometry, Continuous1D):
+                    raise ValueError("Geometry not supported for " + c_lower)
+                self._constraint_prox = lambda z, _: -ProjectNonnegative(-z)
+                self._constraint_oper = SecondOrderFiniteDifference(self.geometry.fun_shape, bc_type='neumann')
+                self._preset["constraint"] = "convex"
+            elif c_lower == "concave":
+                if not isinstance(self.geometry, Continuous1D):
+                    raise ValueError("Geometry not supported for " + c_lower)
+                self._constraint_prox = lambda z, _: ProjectNonnegative(z)
+                self._constraint_oper = SecondOrderFiniteDifference(self.geometry.fun_shape, bc_type='neumann')
+                self._preset["constraint"] = "concave"
             else:
                 raise ValueError("Constraint not supported.")
 
@@ -289,7 +319,7 @@ def _sample(self, N, rng=None):
 
     @staticmethod
     def constraint_options():
-        return ["nonnegativity", "box"]
+        return ["nonnegativity", "box", "increasing", "decreasing", "convex", "concave"]
 
     @staticmethod
     def regularization_options():

cuqi/utilities/__init__.py

@@ -14,8 +14,10 @@
     plot_1D_density,
     plot_2D_density,
     count_nonzero,
+    count_within_bounds,
     count_constant_components_1D,
-    count_constant_components_2D
+    count_constant_components_2D,
+    piecewise_linear_1D_DoF
 )
 
 from ._get_python_variable_name import _get_python_variable_name

cuqi/utilities/_utilities.py

@@ -360,8 +360,29 @@ def count_nonzero(x, threshold = 1e-6):
             Theshold for considering a value as nonzero.
         """
         return np.sum([np.abs(v) >= threshold for v in x])
+
+        
+def count_within_bounds(x, lower_bounds, upper_bounds, threshold = 1e-6, exception = np.nan):
+        """ Returns the number of values in an array whose value lies between the provided lower and upper bounds.
+
+        Parameters
+        ----------
+        x : `np.ndarray` 
+            Array to count elements of.
+
+        lower_bounds : `np.ndarray` 
+            Lower bound on values to disregard when counting.
+
+        upper_bounds : `np.ndarray` 
+            Upper bound on values to disregard.
+
+        threshold : float
+            Theshold for considering a value as nonzero.
+        """
+        return np.sum([l + threshold <= v and v <= u - threshold and not (exception - threshold <= v and v <= exception + threshold) for v, l, u in zip(x, lower_bounds, upper_bounds)])
+        
 
-def count_constant_components_1D(x, threshold = 1e-2, lower = -np.inf, upper = np.inf):
+def count_constant_components_1D(x, threshold = 1e-2, lower = -np.inf, upper = np.inf, exception = np.nan):
         """ Returns the number of piecewise constant components in a one-dimensional array
 
         Parameters
@@ -372,29 +393,40 @@ def count_constant_components_1D(x, threshold = 1e-2, lower = -np.inf, upper = n
         threshold : float
             Strict theshold on when the difference of neighbouring values is considered zero.
 
-        lower : float
+        lower : float or numpy.ndarray
             Piecewise constant components below this value are not counted.
 
-        upper : float
+        upper : float or numpy.ndarray
             Piecewise constant components above this value are not counted.
         """
 
+        if not isinstance(lower, np.ndarray):
+            lower = lower*np.ones_like(x)
+
+        if not isinstance(upper, np.ndarray):
+            upper = upper*np.ones_like(x)
+
         counter = 0
-        if x[0] > lower and x[0] < upper:
+        index = 0
+        if (x[index] > lower[index] and
+            x[index] < upper[index] and
+            not (exception - threshold <= x[index] and x[index] <= exception + threshold)):
             counter += 1
 
         x_previous = x[0]
 
         for x_current in x[1:]:
+            index += 1
             if (abs(x_previous - x_current) >= threshold and
-                x_current > lower and
-                x_current < upper):
+                x_current > lower[index] and
+                x_current < upper[index] and
+                not (exception - threshold <= x_current and x_current <= exception + threshold)):
                     counter += 1
 
             x_previous = x_current
 
         return counter
-        
+
 def count_constant_components_2D(x, threshold = 1e-2, lower = -np.inf, upper = np.inf):
         """ Returns the number of piecewise constant components in a two-dimensional array
 
@@ -406,12 +438,19 @@ def count_constant_components_2D(x, threshold = 1e-2, lower = -np.inf, upper = n
         threshold : float
             Strict theshold on when the difference of neighbouring values is considered zero.
 
-        lower : float
+        lower : float or numpy.ndarray
             Piecewise constant components below this value are not counted.
 
-        upper : float
+        upper : float or numpy.ndarray
             Piecewise constant components above this value are not counted.
         """
+
+        if not isinstance(lower, np.ndarray):
+            lower = lower*np.ones_like(x)
+
+        if not isinstance(upper, np.ndarray):
+            upper = upper*np.ones_like(x)
+
         filled = np.zeros_like(x, dtype = int)
         counter = 0
 
@@ -438,8 +477,51 @@ def process(i, j):
         for i in range(x.shape[0]):
             for j in range(x.shape[1]):
                 if filled[i,j] == 0:
-                    if x[i,j] > lower and x[i,j] < upper:
+                    if x[i,j] > lower[i,j] and x[i,j] < upper[i,j]:
                         counter += 1
                     process(i, j)
         return counter
-                    
+                    
+
+
+def piecewise_linear_1D_DoF(x, threshold = 1e-5, exception_zero = False, exception_flat = False):
+        """ Returns the degrees of freedom of a piecewise linear signal.
+        Assuming linear interpolation, this corresponds to the number of non-differentiable points, including end-points.
+
+        Parameters
+        ----------
+        x : `np.ndarray` 
+            1D Array to compute degrees of freedom of.
+
+        threshold : float
+            Strict theshold on when values are considered zero.
+
+        exception_zero : Boolean
+            Whether a zero piecewise linear components should be considered.
+
+        exception_flat : Boolean
+            Whether a flat piecewise linear components should be considered.
+        """
+
+        differences = x[1:] - x[:-1]
+        double_differences = differences[1:] - differences[:-1]
+
+        joints = [True] + [np.abs(d) >= threshold for d in double_differences] + [True]
+        
+        if not exception_zero and not exception_flat:
+            return np.sum(joints)
+        elif exception_zero:
+            return np.sum(joints and [np.abs(v) < threshold for v in x])
+        elif exception_flat:
+            prev_joint = None
+            counter = 0
+            for i in range(len(joints)):
+                if joints[i]:
+                    counter += 1
+                    if prev_joint is None:
+                        prev_joint = i
+                    else:
+                        if np.abs(x[i] - x[prev_joint]) < threshold:
+                            counter -= 1
+                        prev_joint = i
+            return counter