Add as_cvx_operator and example

examples/solvers/cvx_solvers_poisson.py

@@ -0,0 +1,64 @@
+# 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/>.
+
+"""Poissons 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.
+"""
+
+# Imports for common Python 2/3 codebase
+from __future__ import print_function, division, absolute_import
+from future import standard_library
+standard_library.install_aliases()
+
+import numpy as np
+import odl
+import proximal
+
+# Create space, 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 laplacian
+laplacian = odl.Laplacian(space)
+
+# Create right hand side
+rhs = laplacian(odl.phantom.shepp_logan(space, modified=True))
+rhs += odl.phantom.white_noise(space) * np.std(rhs) * 0.1
+
+# Convert laplacian to cvx operator
+cvx_laplacian = odl.operator.oputils.as_cvx_operator(laplacian)
+
+# Convert to array
+rhs_arr = rhs.asarray()
+
+# Set up optimization problem
+x = proximal.Variable(space.shape)
+funcs = [10 * proximal.sum_squares(cvx_laplacian(x) - rhs_arr),
+         proximal.norm1(proximal.grad(x))]
+prob = proximal.Problem(funcs)
+
+# Solve the problem
+prob.solve(verbose=True, solver='pc', eps_abs=1e-5, eps_rel=1e-5)
+
+# Convert back to odl and display result
+result_odl = space.element(x.value)
+result_odl.show('result')

examples/solvers/cvx_solvers_tomography.py

@@ -0,0 +1,90 @@
+# 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/>.
+
+"""Total variation tomography using the proximal solver.
+
+Solves the optimization problem
+
+    min_x  1/2 ||A(x) - g||_2^2 + lam || |grad(x)| ||_1
+
+Where ``A`` is a parallel beam forward projector, ``grad`` the spatial
+gradient and ``g`` is given noisy data.
+"""
+
+# Imports for common Python 2/3 codebase
+from __future__ import print_function, division, absolute_import
+from future import standard_library
+standard_library.install_aliases()
+
+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'   Require astra tomography to be installed.
+#                             Astra is much faster than scikit. Webpage:
+#                             https://github.com/astra-toolbox/astra-toolbox
+impl = 'astra_cuda'
+
+# Ray transform aka forward projection.
+ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl=impl)
+cvx_ray_trafo = odl.operator.oputils.as_cvx_operator(ray_trafo)
+
+# --- Generate artificial data --- #
+
+
+# Create phantom
+discr_phantom = odl.phantom.shepp_logan(reco_space, modified=True)
+
+# Create sinogram of forward projected phantom with noise
+data = ray_trafo(discr_phantom)
+data += odl.phantom.white_noise(ray_trafo.range) * np.mean(data) * 0.1
+rhs_arr = data.asarray()
+
+# --- Set up the inverse problem --- #
+
+# Note that proximal is not aware of the problem scaling
+x = proximal.Variable(reco_space.shape)
+funcs = [proximal.sum_squares(cvx_ray_trafo(x) - rhs_arr),
+         0.5 * proximal.norm1(proximal.grad(x))]
+
+# Solve the problem
+prob = proximal.Problem(funcs)
+prob.solve(verbose=True, solver='pc', eps_abs=1e-3, eps_rel=1e-3)
+
+# Convert back to odl and display result
+result_odl = reco_space.element(x.value)
+result_odl.show('result')

odl/operator/oputils.py

@@ -243,6 +243,44 @@ def rmatvec(v):
                                               rmatvec=rmatvec,
                                               dtype=dtype)
 
+
+def as_cvx_operator(op, norm_bound=None):
+    """Wrap ``op`` as a ``proximal.BlackBox``.
+
+    This is intended to be used with the cvx solvers.
+
+    Parameters
+    ----------
+    op : `Operator`
+        A linear operator that should be wrapped.
+    norm_bound : float, optional
+        An upper bound on the spectral norm of the operator.
+
+    Returns
+    -------
+    ``proximal.BlackBox`` : cvx_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 other nonlocal type, the overhead is significant.
+    """
+
+    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