Add as_cvx_operator and example

examples/solvers/cvx_solvers.py

@@ -0,0 +1,61 @@
+# 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/>.
+
+"""Example of interfacing with cvx solvers.
+
+This example solves the problem
+
+    -Laplacian(x) = b
+
+where b is a gaussian peak at the origin.
+"""
+
+# 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 odl
+import proximal
+
+# Create discrete space, a square from [-1, 1] x [-1, 1] with (11 x 11) points
+space = odl.uniform_discr([-1, -1], [1, 1], [100, 100])
+
+# Create odl operator for negative laplacian
+laplacian = -odl.Laplacian(space)
+
+# Create right hand side
+rhs = laplacian(odl.phantom.shepp_logan(space))
+
+# Convert laplacian to cvx operator
+cvx_laplacian = odl.operator.oputils.as_cvx_operator(laplacian)
+
+# Convert to array
+rhs_arr = rhs.asarray()
+
+x = proximal.Variable(space.shape)
+funcs = [7.5 * proximal.sum_squares(cvx_laplacian(x) - rhs_arr),
+         proximal.norm1(proximal.grad(x))]
+prob = proximal.Problem(funcs)
+
+# Solve the problem
+with odl.util.Timer('solve problem'):
+    prob.solve(verbose=True, solver='pc', eps_abs=1e-3, eps_rel=1e-3)
+
+# Convert back to odl and display result
+result_odl = space.element(x.value)
+result_odl.show('result')

odl/operator/oputils.py

@@ -243,6 +243,52 @@ 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
+
+    Examples
+    --------
+    Wrap operator and solve simple problem (here toy problem ``Ix = b``)
+
+    >>> import odl
+    >>> op = odl.IdentityOperator(odl.rn(3))
+    >>> scipy_op = as_cvx_operator(op)
+
+    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