Merge pull request #494 from odlgroup/cvx_interface

Add as_cvx_operator and example

README.md

@@ -56,6 +56,7 @@ where FOO, BAR, etc are the dependencies. Others require more complicated instal
 | [FFTW](https://github.com/pyFFTW/pyFFTW)   | Accelerated [FourierTransform](http://odl.readthedocs.io/generated/odl.trafos.fourier.FourierTransform.html) | fftw |
 | [PyWavelets](https://github.com/PyWavelets/pywt)   | Computation of the  [WaveletTransform](http://odl.readthedocs.io/generated/odl.trafos.wavelet.WaveletTransform.html) | pywavelets |
 | [matplotlib](http://matplotlib.org/)   | Visualization through the [show](http://odl.readthedocs.io/generated/odl.discr.lp_discr.DiscreteLpVector.show.html) command | show |
+| [proximal](http://github.com/comp-imaging/ProxImaL)   | Solution of some convex optimization problems | proximal |
 | [pytest](http://pytest.org/latest/)   | Unit tests | testing |
 
 

doc/source/guide/glossary.rst

@@ -76,17 +76,20 @@ Glossary
 
     proximal
         Given a closed proper convex function :math:`f`, the proximal operator is defined by
-        
+
             .. math::
-                
+
                 \text{prox}_f(v) = \arg\min_x f(x) + (1/2)||x - v||_2^2
-        
+
+        proximal is also occationally used instead of ProxImaL, then refering to the proximal
+        modelling language for the solution of convex optimization problems.
+
     proximal factory
-        A proximal factory associated with a function :math:`f` is a `callable`, 
+        A proximal factory associated with a function :math:`f` is a `callable`,
         which returns the proximal of the scaled function :math:`\sigma f` when called with a
-        scalar :math:`\sigma`. This is used due to the fact that optimization methods often use 
+        scalar :math:`\sigma`. This is used due to the fact that optimization methods often use
         :math:`\text{prox}_{\sigma f}` for varying :math:`\sigma`.
-        
+
     range
         Set of elements to which an operator maps, i.e. in which the result of
         an operator evaluation lies.

doc/source/guide/in_depth/index.rst

@@ -11,4 +11,5 @@ This is a more in depth guide to the different parts of ODL.
     linearspace_guide
     vectorization_guide
     numpy_guide
+    proximal_lang_guide
     chambolle_pock_guide

doc/source/guide/in_depth/proximal_lang_guide.rst

@@ -0,0 +1,52 @@
+.. _proximal_lang_in_depth:
+
+#######################
+Using ODL with ProxImaL
+#######################
+
+`Proximal
+<http://www.proximal-lang.org/en/latest/>`_ is a Python-embedded modeling language for image optimization problems and can be used with ODL to solve typical inverse problems phrased as optimization problems. The package is especially suited for non-differentiable problems such as total variance denoising.
+
+Here is a minimal example of solving Poisson's equation equation on an interval with a TV type regularizer (:math:`\min_x \ 10||-\Delta x - rhs||_2^2 + ||\nabla x||_1`)::
+
+   >>> space = odl.uniform_discr(0, 1, 5)
+   >>> op = -odl.Laplacian(space)
+   >>> proximal_lang_op = odl.as_proximal_lang_operator(op)
+   >>> rhs = space.element(lambda x: (x>0.4) & (x<0.6))  # indicator function on [0.4, 0.6]
+   >>> x = proximal.Variable(space.shape)
+   >>> prob = proximal.Problem([10 * proximal.sum_squares(x - rhs.asarray()),
+   >>>                          proximal.norm1(proximal.grad(x))])
+   >>> prob.solve()
+   >>> x.value
+   array([ 0.02352054,  0.02647946,  0.9       ,  0.02647946,  0.02352054])
+
+Notable differences between ODL and ProxImaL
+============================================
+
+It may be tempting to try to convert an arbitrary problem from ODL into ProxImaL, but some differences exist.
+
+Norms
+-----
+Norms in ODL are scaled according to the underlying function space. Hence a sequence of statements converging discretizations give rise to a converging norm::
+
+   >>> for n in range(2, 10000):
+   ...     X = odl.uniform_discr(0, 1, n)
+   ...     print(X.element(lambda x: x).norm())
+   0.559016994375
+   0.576628129734
+   0.577343052266
+   0.577350268468
+   >>> 1 / np.sqrt(3)  # exact result
+   0.577350269189
+
+this is not the case in ProxImaL, where the norm depends on the number of discretization points. Hence a scaling that is correct for a problem in ODL needs not be correct in proximal. This also changes the definition of things like the operator norm.
+
+This also has the added effect of changing the definition of derived features, like the spectral norm of operators.
+
+Spaces
+------
+ODL can represent some complicated spaces, like :math:`\mathbb{R}^3 \times \mathbb{C}^2` through the `ProductSpace` class::
+
+   >>> space = odl.ProductSpace(odl.rn(3), odl.cn(2))
+
+This can then be used in solvers and other structures. ProxImaL currently lacks an equivalent structure.

doc/source/prs.inc

@@ -58,6 +58,7 @@
 .. _PR 491: https://github.com/odlgroup/odl/pull/491
 .. _PR 492: https://github.com/odlgroup/odl/pull/492
 .. _PR 493: https://github.com/odlgroup/odl/pull/493
+.. _PR 494: https://github.com/odlgroup/odl/pull/494
 .. _PR 495: https://github.com/odlgroup/odl/pull/495
 .. _PR 500: https://github.com/odlgroup/odl/pull/500
 .. _PR 502: https://github.com/odlgroup/odl/pull/502

doc/source/refs.rst

@@ -42,6 +42,9 @@ References
 .. [GNS2009] Griva, I, Nash, S G, and Sofer, A. *Linear and nonlinear
    optimization*. Siam, 2009.
 
+.. [Hei+2016] Heide, F et al. *ProxImaL: Efficient Image Optimization using
+   Proximal Algorithms*. ACM Transactions on Graphics (TOG), 2016.
+
 .. [Kva1991] Kvaalen, E. *A faster Broyden method*. BIT Numerical
    Mathematics 31 (1991), pp 369--372.
 

doc/source/release_notes.rst

@@ -11,8 +11,7 @@ Next release
 
 New features
 ------------
-- First implementation of a deformation model based on linearized deformations using displacement
-  vector fields. (`PR 488`)
+- Add option to interface with ProxImaL solvers using ODL operators. (`PR 494`)
 
 ODL 0.3.1 Release Notes (2016-08-15)
 ====================================
@@ -23,7 +22,7 @@ splitting.
 
 New features
 ------------
-- New solvers based on the Douglas-Rachford and forward-backward splitting schemes. (`PR 478`_, 
+- New solvers based on the Douglas-Rachford and forward-backward splitting schemes. (`PR 478`_,
   `PR 480`_)
 - ``NormOperator`` and ``DistOperator`` added. (`PR 487`_)
 - Single-element ``NtuplesBase`` vectors can now be converted to ``float``, ``complex`` etc.
@@ -38,7 +37,7 @@ Improvements
   (`PR 489`_, `PR 482`_, `PR 491`_)
 - Clearer separation between attributes that are intended as part of the subclassing API and those
   that are not. (`PR 471`_)
-- Chambolle-Pock solver accepts also non-linear operators and has better documentation now. 
+- Chambolle-Pock solver accepts also non-linear operators and has better documentation now.
   (`PR 490`_)
 - Clean-up of imports. (`PR 492`_)
 - All solvers now check that the given start value ``x`` is in ``op.domain``. (`PR 502`_)

examples/solvers/proximal_lang_poisson.py

@@ -0,0 +1,63 @@
+# Copyright 2014-2016 The ODL development group
+#
+# This file is part of ODL.
+#
+# ODL is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ODL is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with ODL.  If not, see <http://www.gnu.org/licenses/>.
+
+"""Poisson's problem using the ProxImaL solver.
+
+Solves the optimization problem
+
+    min_x  10 ||laplacian(x) - g||_2^2 + || |grad(x)| ||_1
+
+Where ``laplacian`` is the spatial Laplacian, ``grad`` the spatial
+gradient and ``g`` is given noisy data.
+"""
+
+import numpy as np
+import odl
+import proximal
+
+# Create space defined on a square from [0, 0] to [100, 100] with (100 x 100)
+# points
+space = odl.uniform_discr([0, 0], [100, 100], [100, 100])
+
+# Create ODL operator for the Laplacian
+laplacian = odl.Laplacian(space)
+
+# Create right hand side
+phantom = odl.phantom.shepp_logan(space, modified=True)
+phantom.show('original image')
+rhs = laplacian(phantom)
+rhs += odl.phantom.white_noise(space) * np.std(rhs) * 0.1
+rhs.show('rhs')
+
+# Convert laplacian to ProxImaL operator
+proximal_lang_laplacian = odl.as_proximal_lang_operator(laplacian)
+
+# Convert to array
+rhs_arr = rhs.asarray()
+
+# Set up optimization problem
+x = proximal.Variable(space.shape)
+funcs = [10 * proximal.sum_squares(proximal_lang_laplacian(x) - rhs_arr),
+         proximal.norm1(proximal.grad(x))]
+
+# Solve the problem using ProxImaL
+prob = proximal.Problem(funcs)
+prob.solve(verbose=True)
+
+# Convert back to odl and display result
+result_odl = space.element(x.value)
+result_odl.show('result from ProxImaL')

examples/solvers/proximal_lang_tomography.py

@@ -0,0 +1,89 @@
+# Copyright 2014-2016 The ODL development group
+#
+# This file is part of ODL.
+#
+# ODL is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ODL is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with ODL.  If not, see <http://www.gnu.org/licenses/>.
+
+"""Tomography with TV regularization using the ProxImaL solver.
+
+Solves the optimization problem
+
+    min_{0 <= x <= 1}  ||A(x) - g||_2^2 + 0.2 || |grad(x)| ||_1
+
+Where ``A`` is a parallel beam forward projector, ``grad`` the spatial
+gradient and ``g`` is given noisy data.
+"""
+
+import numpy as np
+import odl
+import proximal
+
+
+# --- Set up the forward operator (ray transform) --- #
+
+
+# Discrete reconstruction space: discretized functions on the rectangle
+# [-20, 20]^2 with 300 samples per dimension.
+reco_space = odl.uniform_discr(
+    min_corner=[-20, -20], max_corner=[20, 20], nsamples=[300, 300],
+    dtype='float32')
+
+# Make a parallel beam geometry with flat detector
+# Angles: uniformly spaced, n = 360, min = 0, max = 2 * pi
+angle_partition = odl.uniform_partition(0, 2 * np.pi, 360)
+# Detector: uniformly sampled, n = 558, min = -30, max = 30
+detector_partition = odl.uniform_partition(-30, 30, 558)
+geometry = odl.tomo.Parallel2dGeometry(angle_partition, detector_partition)
+
+# The implementation of the ray transform to use, options:
+# 'scikit'                    Requires scikit-image (can be installed by
+#                             running ``pip install scikit-image``).
+# 'astra_cpu', 'astra_cuda'   Requires astra tomography to be installed.
+#                             Astra is much faster than scikit. Webpage:
+#                             https://github.com/astra-toolbox/astra-toolbox
+impl = 'astra_cuda'
+
+# Initialize the ray transform (forward projection).
+ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl=impl)
+
+# Convert ray transform to proximal language operator
+proximal_lang_ray_trafo = odl.as_proximal_lang_operator(ray_trafo)
+
+# Create sinogram of forward projected phantom with noise
+phantom = odl.phantom.shepp_logan(reco_space, modified=True)
+phantom.show('phantom')
+data = ray_trafo(phantom)
+data += odl.phantom.white_noise(ray_trafo.range) * np.mean(data) * 0.1
+data.show('noisy data')
+
+# Convert to array for ProxImaL
+rhs_arr = data.asarray()
+
+# Set up optimization problem
+# Note that proximal is not aware of the underlying space and only works with
+# matrices. Hence the norm in proximal does not match the norm in the ODL space
+# exactly.
+x = proximal.Variable(reco_space.shape)
+funcs = [proximal.sum_squares(proximal_lang_ray_trafo(x) - rhs_arr),
+         0.2 * proximal.norm1(proximal.grad(x)),
+         proximal.nonneg(x),
+         proximal.nonneg(1 - x)]
+
+# Solve the problem using ProxImaL
+prob = proximal.Problem(funcs)
+prob.solve(verbose=True)
+
+# Convert back to odl and display result
+result_odl = reco_space.element(x.value)
+result_odl.show('ProxImaL result')

odl/operator/operator.py

@@ -538,7 +538,7 @@ def _call(self, x, out=None, **kwargs):
         pattern.
 
         See the `documentation
-        <https://odl.readthedocs.org/guide/in_depth/operator_guide.html>`_
+        <https://odl.readthedocs.io/guide/in_depth/operator_guide.html>`_
         for more info on in-place vs. out-of-place evaluation.
 
         Parameters

odl/operator/oputils.py

@@ -28,7 +28,8 @@
 from odl.space.base_ntuples import FnBase
 from odl.space import ProductSpace
 
-__all__ = ('matrix_representation', 'power_method_opnorm', 'as_scipy_operator')
+__all__ = ('matrix_representation', 'power_method_opnorm', 'as_scipy_operator',
+           'as_proximal_lang_operator')
 
 
 def matrix_representation(op):
@@ -238,7 +239,7 @@ def as_scipy_operator(op):
     -----
     If the data representation of ``op``'s domain and range is of type
     `NumpyFn` this incurs no significant overhead. If the data type is `CudaFn`
-    or other nonlocal type, the overhead is significant.
+    or some other nonlocal type, the overhead is significant.
     """
     if not op.is_linear:
         raise ValueError('`op` needs to be linear')
@@ -268,6 +269,53 @@ def rmatvec(v):
                                               rmatvec=rmatvec,
                                               dtype=dtype)
 
+
+def as_proximal_lang_operator(op, norm_bound=None):
+    """Wrap ``op`` as a ``proximal.BlackBox``.
+
+    This is intended to be used with the `ProxImaL language solvers.
+    <https://github.com/comp-imaging/proximal>`_
+
+    Parameters
+    ----------
+    op : `Operator`
+        Linear operator to be wrapped. Its domain and range must implement
+        ``shape``, and elements in these need to implement ``asarray``.
+    norm_bound : float, optional
+        An upper bound on the spectral norm of the operator. Note that this is
+        the norm as defined by ProxImaL, and hence use the unweighted spaces.
+
+    Returns
+    -------
+    ``proximal.BlackBox`` : proximal_lang_operator
+        The wrapped operator.
+
+    Notes
+    -----
+    If the data representation of ``op``'s domain and range is of type
+    `NumpyFn` this incurs no significant overhead. If the data type is `CudaFn`
+    or some other nonlocal type, the overhead is significant.
+
+    References
+    ----------
+    For documentation on the proximal language (ProxImaL) see [Hei+2016]_.
+    """
+
+    # TODO: use out parameter once "as editable array" is added
+
+    def forward(inp, out):
+        out[:] = op(inp).asarray()
+
+    def adjoint(inp, out):
+        out[:] = op.adjoint(inp).asarray()
+
+    import proximal
+    return proximal.LinOpFactory(input_shape=op.domain.shape,
+                                 output_shape=op.range.shape,
+                                 forward=forward,
+                                 adjoint=adjoint,
+                                 norm_bound=norm_bound)
+
 if __name__ == '__main__':
     # pylint: disable=wrong-import-position
     from odl.util.testutils import run_doctests

setup.py

@@ -138,7 +138,8 @@ def run_tests(self):
         'fftw': 'pyfftw',
         'pywavelets': 'Pywavelets',
         'scikit' : 'scikit-image',
-        'all': test_requires + ['matplotlib', 'pyfftw', 'Pywavelets', 'scikit-image']
+        'proximal' : 'proximal',
+        'all': test_requires + ['matplotlib', 'pyfftw', 'Pywavelets', 'scikit-image', 'proximal']
     },
 
     cmdclass={'test': PyTest},