ENH: add decorator to fix adjoint weightings

examples/operator/auto_adjoint_weighting.py

@@ -0,0 +1,75 @@
+"""Example demonstrating the usage of ``auto_adjoint_weighting``.
+
+This is an advanced example showing how the `auto_adjoint_weighting`
+can be used to automatically perform the correct weighting of the
+adjoint of an `Operator`, depending on the weightings used in the
+operator domain and range.
+
+The general idea is that users implement the adjoint for the operator
+variant that maps between ``R^n`` or ``C^n`` type spaces that have no
+weightings associated with them. The `auto_adjoint_weighting` decorator then
+implements the adjoint weighting for operator variants between spaces
+that are "similar" to ``R^n`` or ``C^n``, but are weighted, for instance
+discretized ``L^p`` function spaces.
+
+See the `auto_adjoint_weighting` documentation for more details.
+"""
+
+import odl
+from odl.operator.oputils import auto_adjoint_weighting
+
+
+class ScalingOp(odl.Operator):
+
+    """Operator that scales its input by a constant."""
+
+    def __init__(self, dom, ran, c):
+        super(ScalingOp, self).__init__(dom, ran, linear=True)
+        self.c = c
+
+    def _call(self, x):
+        return self.c * x
+
+    @property
+    @auto_adjoint_weighting
+    def adjoint(self):
+        """Adjoint of the scaling operator.
+
+        Note that we return the adjoint for the scaling operator
+        ``S: R^n -> R^n``, and the decorator implements the other cases.
+        """
+        return ScalingOp(self.range, self.domain, self.c)
+
+
+rn = odl.rn(2)  # Constant weight 1
+discr = odl.uniform_discr(0, 4, 2)  # Constant weight 2
+print('Checking scaling operators X -> X, Y -> Y and X -> Y with ')
+print('X = {!r}, Y ={!r}'.format(rn, discr))
+print('')
+
+op1 = ScalingOp(rn, rn, 2)  # Same weightings, no scaling in adjoint
+op2 = ScalingOp(discr, discr, 2)  # Same weightings, no scaling in adjoint
+op3 = ScalingOp(rn, discr, 2)  # Different weightings, adjoint scales
+
+# Look at output of adjoint
+print('X -> X adjoint at one:', op1.adjoint(op1.range.one()))
+print('Y -> Y adjoint at one:', op2.adjoint(op2.range.one()))
+print('X -> Y adjoint at one:', op3.adjoint(op3.range.one()))
+print('')
+
+# Check adjointness
+print('Adjointness check:')
+inner1_dom = op1.domain.one().inner(op1.adjoint(op1.range.one()))
+inner1_ran = op1(op1.domain.one()).inner(op1.range.one())
+print('X -> X:    <Sx, y> = {},  <x, S^*y> = {}'
+      ''.format(inner1_ran, inner1_dom))
+
+inner2_dom = op2.domain.one().inner(op2.adjoint(op2.range.one()))
+inner2_ran = op2(op2.domain.one()).inner(op2.range.one())
+print('Y -> Y:    <Sx, y> = {},  <x, S^*y> = {}'
+      ''.format(inner2_ran, inner2_dom))
+
+inner3_dom = op3.domain.one().inner(op3.adjoint(op3.range.one()))
+inner3_ran = op3(op3.domain.one()).inner(op3.range.one())
+print('X -> Y:    <Sx, y> = {},  <x, S^*y> = {}'
+      ''.format(inner3_ran, inner3_dom))

odl/discr/diff_ops.py

@@ -6,12 +6,13 @@
 # v. 2.0. If a copy of the MPL was not distributed with this file, You can
 # obtain one at https://mozilla.org/MPL/2.0/.
 
-"""Operators defined for tensor fields."""
+"""Differential operators."""
 
 from __future__ import print_function, division, absolute_import
 import numpy as np
 
 from odl.discr.lp_discr import DiscreteLp
+from odl.operator.oputils import auto_adjoint_weighting
 from odl.operator.tensor_ops import PointwiseTensorFieldOperator
 from odl.space import ProductSpace
 from odl.util import writable_array, signature_string, indent
@@ -164,6 +165,7 @@ def derivative(self, point=None):
             return self
 
     @property
+    @auto_adjoint_weighting
     def adjoint(self):
         """Return the adjoint operator."""
         if not self.is_linear:
@@ -380,6 +382,7 @@ def derivative(self, point=None):
             return self
 
     @property
+    @auto_adjoint_weighting
     def adjoint(self):
         """Adjoint of this operator.
 
@@ -594,6 +597,7 @@ def derivative(self, point=None):
             return self
 
     @property
+    @auto_adjoint_weighting
     def adjoint(self):
         """Adjoint of this operator.
 
@@ -644,6 +648,9 @@ def __init__(self, domain, range=None, pad_mode='constant', pad_const=0):
         ----------
         domain : `DiscreteLp`
             Space of elements which the operator is acting on.
+        range : `DiscreteLp`, optional
+            Space of elements to which the operator maps.
+            Default: ``domain``
         pad_mode : string, optional
             The padding mode to use outside the domain.
 
@@ -692,9 +699,6 @@ def __init__(self, domain, range=None, pad_mode='constant', pad_const=0):
         if range is None:
             range = domain
 
-        super(Laplacian, self).__init__(
-            domain, range, base_space=domain, linear=True)
-
         self.pad_mode, pad_mode_in = str(pad_mode).lower(), pad_mode
         if pad_mode not in _SUPPORTED_PAD_MODES:
             raise ValueError('`pad_mode` {} not understood'
@@ -705,6 +709,10 @@ def __init__(self, domain, range=None, pad_mode='constant', pad_const=0):
             raise ValueError('`pad_mode` {} not implemented for Laplacian.'
                              ''.format(pad_mode_in))
 
+        linear = not (self.pad_mode == 'constant' and pad_const != 0)
+        super(Laplacian, self).__init__(
+            domain, range, base_space=domain, linear=linear)
+
         self.pad_const = self.domain.field.element(pad_const)
 
     def _call(self, x, out=None):
@@ -760,11 +768,16 @@ def derivative(self, point=None):
             return self
 
     @property
+    @auto_adjoint_weighting
     def adjoint(self):
         """Return the adjoint operator.
 
         The laplacian is self-adjoint, so this returns ``self``.
         """
+        if not self.is_linear:
+            raise ValueError('operator with nonzero pad_const ({}) is not'
+                             ' linear and has no adjoint'
+                             ''.format(self.pad_const))
         return Laplacian(self.range, self.domain,
                          pad_mode=self.pad_mode, pad_const=0)
 

odl/discr/discr_ops.py

@@ -13,6 +13,7 @@
 
 from odl.discr import DiscreteLp, uniform_partition
 from odl.operator import Operator
+from odl.operator.oputils import auto_adjoint_weighting
 from odl.set import IntervalProd
 from odl.space import FunctionSpace, tensor_space
 from odl.util import (
@@ -98,24 +99,6 @@ def inverse(self):
         The returned operator is resampling defined in the opposite
         direction.
 
-        See Also
-        --------
-        adjoint : resampling is unitary, so the adjoint is the inverse.
-        """
-        return Resampling(self.range, self.domain)
-
-    @property
-    def adjoint(self):
-        """Return an (approximate) adjoint.
-
-        The result is only exact if the interpolation and sampling
-        operators of the underlying spaces match exactly.
-
-        Returns
-        -------
-        adjoint : Resampling
-            Resampling operator defined in the opposite direction.
-
         Examples
         --------
         Create resampling operator and inverse:
@@ -137,30 +120,61 @@ def adjoint(self):
         >>> y = [0.0, 0.0, 0.0, 1.0, 0.0, 0.0]
         >>> print(resampling(resampling_inv(y)))
         [ 0.,  0.,  0.,  0.,  0.,  0.]
+
+        See Also
+        --------
+        adjoint : resampling is unitary, so the adjoint is the inverse.
+        """
+        return Resampling(self.range, self.domain)
+
+    @property
+    @auto_adjoint_weighting
+    def adjoint(self):
+        """Return an (approximate) adjoint.
+
+        The result is only exact if the interpolation and sampling
+        operators of the underlying spaces match exactly.
+
+        Returns
+        -------
+        adjoint : Resampling
+            Resampling operator defined in the opposite direction.
         """
         return self.inverse
 
 
-class ResizingOperatorBase(Operator):
 
-    """Base class for `ResizingOperator` and its adjoint.
+class ResizingOperator(Operator):
 
-    This is an abstract class used to share code between the forward and
-    adjoint variants of the resizing operator.
+    """Operator mapping a discretized function to a new domain.
+
+    This operator is a mapping between uniformly discretized
+    `DiscreteLp` spaces with the same `DiscreteLp.cell_sides`,
+    but different `DiscreteLp.shape`. The underlying operation is array
+    resizing, i.e. no resampling is performed.
+    In axes where the domain is enlarged, the new entries are filled
+    ("padded") according to a provided parameter ``pad_mode``.
+
+    All resizing operator variants are linear, except constant padding
+    with constant != 0.
+
+    See `the online documentation
+    <https://odlgroup.github.io/odl/math/resizing_ops.html>`_
+    on resizing operators for mathematical details.
     """
 
     def __init__(self, domain, range=None, ran_shp=None, **kwargs):
         """Initialize a new instance.
 
         Parameters
         ----------
-        domain : uniform `DiscreteLp`
-            Uniformly discretized space, the operator can be applied
-            to its elements.
-        range : uniform `DiscreteLp`, optional
-            Uniformly discretized space in which the result of the
-            application of this operator lies.
-            For the default ``None``, a space with the same attributes
+        domain : `DiscreteLp`
+            Space of discretized functions to which the operator can be
+            applied. It must be uniformly discretized in axes where
+            resizing is applied.
+        range : `DiscreteLp`, optional
+            Space in which the result of the application of this operator
+            lies. For the default ``None``, a space with the same attributes
             as ``domain`` is used, except for its shape, which is set
             to ``ran_shp``.
         ran_shp : sequence of ints, optional
@@ -313,8 +327,7 @@ def __init__(self, domain, range=None, ran_shp=None, **kwargs):
         # padding mode 'constant' with `pad_const != 0` is not linear
         linear = (self.pad_mode != 'constant' or self.pad_const == 0.0)
 
-        super(ResizingOperatorBase, self).__init__(
-            domain, ran, linear=linear)
+        super(ResizingOperator, self).__init__(domain, ran, linear=linear)
 
     @property
     def offset(self):
@@ -337,26 +350,6 @@ def axes(self):
         return tuple(i for i in range(self.domain.ndim)
                      if self.domain.shape[i] != self.range.shape[i])
 
-
-class ResizingOperator(ResizingOperatorBase):
-
-    """Operator mapping a discretized function to a new domain.
-
-    This operator is a mapping between uniformly discretized
-    `DiscreteLp` spaces with the same `DiscreteLp.cell_sides`,
-    but different `DiscreteLp.shape`. The underlying operation is array
-    resizing, i.e. no resampling is performed.
-    In axes where the domain is enlarged, the new entries are filled
-    ("padded") according to a provided parameter ``pad_mode``.
-
-    All resizing operator variants are linear, except constant padding
-    with constant != 0.
-
-    See `the online documentation
-    <https://odlgroup.github.io/odl/math/resizing_ops.html>`_
-    on resizing operators for mathematical details.
-    """
-
     def _call(self, x, out):
         """Implement ``self(x, out)``."""
         with writable_array(out) as out_arr:
@@ -380,15 +373,16 @@ def derivative(self, point):
             return self
 
     @property
+    @auto_adjoint_weighting
     def adjoint(self):
         """Adjoint of this operator."""
         if not self.is_linear:
             raise NotImplementedError('this operator is not linear and '
                                       'thus has no adjoint')
 
-        forward_op = self
+        op = self
 
-        class ResizingOperatorAdjoint(ResizingOperatorBase):
+        class ResizingOperatorAdjoint(Operator):
 
             """Adjoint of `ResizingOperator`.
 
@@ -397,18 +391,23 @@ class ResizingOperatorAdjoint(ResizingOperatorBase):
             on resizing operators for mathematical details.
             """
 
+            def __init__(self):
+                """Initialize a new instance."""
+                super(ResizingOperatorAdjoint, self).__init__(
+                    op.range, op.domain, linear=True)
+
             def _call(self, x, out):
                 """Implement ``self(x, out)``."""
                 with writable_array(out) as out_arr:
-                    resize_array(x.asarray(), self.range.shape,
-                                 offset=self.offset, pad_mode=self.pad_mode,
+                    resize_array(x.asarray(), op.domain.shape,
+                                 offset=op.offset, pad_mode=op.pad_mode,
                                  pad_const=0, direction='adjoint',
                                  out=out_arr)
 
             @property
             def adjoint(self):
                 """Adjoint of the adjoint, i.e. the original operator."""
-                return forward_op
+                return op
 
             @property
             def inverse(self):
@@ -422,8 +421,7 @@ def inverse(self):
                     domain=self.range, range=self.domain,
                     pad_mode=self.pad_mode)
 
-        return ResizingOperatorAdjoint(domain=self.range, range=self.domain,
-                                       pad_mode=self.pad_mode)
+        return ResizingOperatorAdjoint()
 
     @property
     def inverse(self):