Merge branch 'cleanup_my_ops' into 'NIFTy_6'

Cleanup my ops for NIFTy6 release

See merge request !469

ChangeLog

@@ -16,7 +16,6 @@ In addition to the below changes, the following operators were introduced:
 * PartialConjugate: Conjugates parts of a multi-field
 * SliceOperator: Geometry preserving mask operator
 * SplitOperator: Splits a single field into a multi-field
-* SwitchSpacesOperator: Permutes the domain entries of fields
 
 FFT convention adjusted
 =======================

nifty6/__init__.py

@@ -43,9 +43,9 @@ from .operators.selection_operators import SliceOperator, SplitOperator
 from .operators.block_diagonal_operator import BlockDiagonalOperator
 from .operators.outer_product_operator import OuterProduct
 from .operators.simple_linear_operators import (
-    VdotOperator, ConjugationOperator, Realizer,
-    FieldAdapter, ducktape, GeometryRemover, NullOperator,
-    MatrixProductOperator, PartialExtractor, SwitchSpacesOperator)
+    VdotOperator, ConjugationOperator, Realizer, FieldAdapter, ducktape,
+    GeometryRemover, NullOperator, PartialExtractor)
+from .operators.matrix_product_operator import MatrixProductOperator
 from .operators.value_inserter import ValueInserter
 from .operators.energy_operators import (
     EnergyOperator, GaussianEnergy, PoissonianEnergy, InverseGammaLikelihood,

nifty6/operators/matrix_product_operator.py

+# This program 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.
+#
+# This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Copyright(C) 2013-2019 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+from ..domain_tuple import DomainTuple
+from ..field import Field
+from .endomorphic_operator import EndomorphicOperator
+from .. import utilities
+import numpy as np
+
+
+class MatrixProductOperator(EndomorphicOperator):
+    """Endomorphic matrix multiplication with input field.
+
+    This operator supports scipy.sparse matrices and numpy arrays
+    as the matrix to be applied.
+
+    For numpy array matrices, can apply the matrix over a subspace
+    of the input.
+
+    If the input arrays have more than one dimension, for
+    scipy.sparse matrices the `flatten` keyword argument must be
+    set to true. This means that the input field will be flattened
+    before applying the matrix and reshaped to its original shape
+    afterwards.
+
+    Matrices are tested regarding their compatibility with the
+    called for application method.
+
+    Flattening and subspace application are mutually exclusive.
+
+    Parameters
+    ----------
+    domain: :class:`Domain` or :class:`DomainTuple`
+        Domain of the operator.
+        If :class:`DomainTuple` it is assumed to have only one entry.
+    matrix: scipy.sparse matrix or numpy array
+        Quadratic matrix of shape `(domain.shape, domain.shape)`
+        (if `not flatten`) that supports `matrix.transpose()`.
+        If it is not a numpy array, needs to be applicable to the val
+        array of input fields by `matrix.dot()`.
+    spaces: int or tuple of int, optional
+        The subdomain(s) of "domain" which the operator acts on.
+        If None, it acts on all elements.
+        Only possible for numpy array matrices.
+        If `len(domain) > 1` and `flatten=False`, this parameter is
+        mandatory.
+    flatten: boolean, optional
+        Whether the input value array should be flattened before
+        applying the matrix and reshaped to its original shape
+        afterwards.
+        Needed for scipy.sparse matrices if `len(domain) > 1`.
+    """
+    def __init__(self, domain, matrix, spaces=None, flatten=False):
+        self._capability = self.TIMES | self.ADJOINT_TIMES
+        self._domain = DomainTuple.make(domain)
+
+        mat_dim = len(matrix.shape)
+
+        if mat_dim % 2 != 0 or \
+           matrix.shape != (matrix.shape[:mat_dim//2] + matrix.shape[:mat_dim//2]):
+            raise ValueError("Matrix must be quadratic.")
+        appl_dim = mat_dim // 2  # matrix application space dimension
+
+        # take shortcut for trivial case
+        if spaces is not None:
+            if len(self._domain.shape) == 1 and spaces == (0, ):
+                spaces = None
+
+        if spaces is None:
+            self._spaces = None
+            self._active_axes = utilities.my_sum(self._domain.axes)
+            appl_space_shape = self._domain.shape
+            if flatten:
+                appl_space_shape = (utilities.my_product(appl_space_shape), )
+        else:
+            if flatten:
+                raise ValueError(
+                    "Cannot flatten input AND apply to a subspace")
+            if not isinstance(matrix, np.ndarray):
+                raise ValueError(
+                    "Application to subspaces only supported for numpy array matrices."
+                )
+            self._spaces = utilities.parse_spaces(spaces, len(self._domain))
+            appl_space_shape = []
+            active_axes = []
+            for space_idx in spaces:
+                appl_space_shape += self._domain[space_idx].shape
+                active_axes += self._domain.axes[space_idx]
+            appl_space_shape = tuple(appl_space_shape)
+            self._active_axes = tuple(active_axes)
+
+            self._mat_last_n = tuple([-appl_dim + i for i in range(appl_dim)])
+            self._mat_first_n = np.arange(appl_dim)
+
+        # Test if the matrix and the array it will be applied to fit
+        if matrix.shape[:appl_dim] != appl_space_shape:
+            raise ValueError(
+                "Matrix and domain shapes are incompatible under the requested "
+                + "application scheme.\n" +
+                f"Matrix appl shape: {matrix.shape[:appl_dim]}, " +
+                f"appl_space_shape: {appl_space_shape}.")
+
+        self._mat = matrix
+        self._mat_tr = matrix.transpose().conjugate()
+        self._flatten = flatten
+
+    def apply(self, x, mode):
+        self._check_input(x, mode)
+        times = (mode == self.TIMES)
+        m = self._mat if times else self._mat_tr
+
+        if self._spaces is None:
+            if not self._flatten:
+                res = m.dot(x.val)
+            else:
+                res = m.dot(x.val.flatten()).reshape(self._domain.shape)
+            return Field(self._domain, res)
+
+        mat_axes = self._mat_last_n if times else np.flip(self._mat_last_n)
+        move_axes = self._mat_first_n if times else np.flip(self._mat_first_n)
+        res = np.tensordot(m, x.val, axes=(mat_axes, self._active_axes))
+        res = np.moveaxis(res, move_axes, self._active_axes)
+        return Field(self._domain, res)

nifty6/operators/simple_linear_operators.py

@@ -349,156 +349,3 @@ class PartialExtractor(LinearOperator):
         res0 = MultiField.from_dict({key: x[key] for key in x.domain.keys()})
         res1 = MultiField.full(self._compldomain, 0.)
         return res0.unite(res1)
-
-
-class MatrixProductOperator(EndomorphicOperator):
-    """Endomorphic matrix multiplication with input field.
-
-    This operator supports scipy.sparse matrices and numpy arrays
-    as the matrix to be applied.
-
-    For numpy array matrices, can apply the matrix over a subspace
-    of the input.
-
-    If the input arrays have more than one dimension, for
-    scipy.sparse matrices the `flatten` keyword argument must be
-    set to true. This means that the input field will be flattened
-    before applying the matrix and reshaped to its original shape
-    afterwards.
-
-    Matrices are tested regarding their compatibility with the
-    called for application method.
-
-    Flattening and subspace application are mutually exclusive.
-
-    Parameters
-    ----------
-    domain: :class:`Domain` or :class:`DomainTuple`
-        Domain of the operator.
-        If :class:`DomainTuple` it is assumed to have only one entry.
-    matrix: scipy.sparse matrix or numpy array
-        Quadratic matrix of shape `(domain.shape, domain.shape)`
-        (if `not flatten`) that supports `matrix.transpose()`.
-        If it is not a numpy array, needs to be applicable to the val
-        array of input fields by `matrix.dot()`.
-    spaces: int or tuple of int, optional
-        The subdomain(s) of "domain" which the operator acts on.
-        If None, it acts on all elements.
-        Only possible for numpy array matrices.
-        If `len(domain) > 1` and `flatten=False`, this parameter is
-        mandatory.
-    flatten: boolean, optional
-        Whether the input value array should be flattened before
-        applying the matrix and reshaped to its original shape
-        afterwards.
-        Needed for scipy.sparse matrices if `len(domain) > 1`.
-    """
-    def __init__(self, domain, matrix, spaces=None, flatten=False):
-        self._capability = self.TIMES | self.ADJOINT_TIMES
-        self._domain = DomainTuple.make(domain)
-
-        mat_dim = len(matrix.shape)
-
-        if mat_dim % 2 != 0 or \
-           matrix.shape != (matrix.shape[:mat_dim//2] + matrix.shape[:mat_dim//2]):
-            raise ValueError("Matrix must be quadratic.")
-        appl_dim = mat_dim // 2  # matrix application space dimension
-
-        # take shortcut for trivial case
-        if spaces is not None:
-            if len(self._domain.shape) == 1 and spaces == (0, ):
-                spaces = None
-
-        if spaces is None:
-            self._spaces = None
-            self._active_axes = utilities.my_sum(self._domain.axes)
-            appl_space_shape = self._domain.shape
-            if flatten:
-                appl_space_shape = (utilities.my_product(appl_space_shape), )
-        else:
-            if flatten:
-                raise ValueError(
-                    "Cannot flatten input AND apply to a subspace")
-            if not isinstance(matrix, np.ndarray):
-                raise ValueError(
-                    "Application to subspaces only supported for numpy array matrices."
-                )
-            self._spaces = utilities.parse_spaces(spaces, len(self._domain))
-            appl_space_shape = []
-            active_axes = []
-            for space_idx in spaces:
-                appl_space_shape += self._domain[space_idx].shape
-                active_axes += self._domain.axes[space_idx]
-            appl_space_shape = tuple(appl_space_shape)
-            self._active_axes = tuple(active_axes)
-
-            self._mat_last_n = tuple([-appl_dim + i for i in range(appl_dim)])
-            self._mat_first_n = np.arange(appl_dim)
-
-        # Test if the matrix and the array it will be applied to fit
-        if matrix.shape[:appl_dim] != appl_space_shape:
-            raise ValueError(
-                "Matrix and domain shapes are incompatible under the requested "
-                + "application scheme.\n" +
-                f"Matrix appl shape: {matrix.shape[:appl_dim]}, " +
-                f"appl_space_shape: {appl_space_shape}.")
-
-        self._mat = matrix
-        self._mat_tr = matrix.transpose().conjugate()
-        self._flatten = flatten
-
-    def apply(self, x, mode):
-        self._check_input(x, mode)
-        times = (mode == self.TIMES)
-        m = self._mat if times else self._mat_tr
-
-        if self._spaces is None:
-            if not self._flatten:
-                res = m.dot(x.val)
-            else:
-                res = m.dot(x.val.flatten()).reshape(self._domain.shape)
-            return Field(self._domain, res)
-
-        mat_axes = self._mat_last_n if times else np.flip(self._mat_last_n)
-        move_axes = self._mat_first_n if times else np.flip(self._mat_first_n)
-        res = np.tensordot(m, x.val, axes=(mat_axes, self._active_axes))
-        res = np.moveaxis(res, move_axes, self._active_axes)
-        return Field(self._domain, res)
-
-
-class SwitchSpacesOperator(LinearOperator):
-    """Operator to permutate the domain entries of fields.
-
-    Exchanges the entries `space1` and `space2` of the input's domain.
-    """
-    def __init__(self, domain, space1, space2=0):
-        self._capability = self.TIMES | self.ADJOINT_TIMES
-        self._domain = DomainTuple.make(domain)
-
-        n_spaces = len(self._domain)
-        if space1 >= n_spaces or space1 < 0 \
-          or space2 >= n_spaces or space2 < 0:
-            raise ValueError("invalid space value")
-
-        tgt = list(self._domain)
-        tgt[space2] = self._domain[space1]
-        tgt[space1] = self._domain[space2]
-        self._target = DomainTuple.make(tgt)
-
-        dom_axes = self._domain.axes
-        tgt_axes = self._target.axes
-
-        self._axes_dom = dom_axes[space1] + dom_axes[space2]
-        self._axes_tgt = tgt_axes[space2] + tgt_axes[space1]
-
-    def apply(self, x, mode):
-        self._check_input(x, mode)
-
-        if mode == self.TIMES:
-            val = np.moveaxis(x.val, self._axes_dom, self._axes_tgt)
-            dom = self._target
-        else:
-            val = np.moveaxis(x.val, self._axes_tgt, self._axes_dom)
-            dom = self._domain
-
-        return Field(dom, val)

test/test_operators/test_adjoint.py

@@ -326,20 +326,3 @@ def testSlowFieldAdapter(seed):
     dom = {'a': ift.RGSpace(1), 'b': ift.RGSpace(2)}
     op = ift.operators.simple_linear_operators._SlowFieldAdapter(dom, 'a')
     ift.extra.consistency_check(op)
-
-@pmp('sp1', [0, 2])
-@pmp('sp2', [1])
-@pmp('seed', [12, 3])
-def testSwitchSpacesOperator(sp1, sp2, seed):
-    with ift.random.Context(seed):
-        dom1 = ift.RGSpace(1)
-        dom2 = ift.RGSpace((2, 2))
-        dom3 = ift.RGSpace(3)
-        dom = ift.DomainTuple.make([dom1, dom2, dom3])
-        op = ift.SwitchSpacesOperator(dom, sp1, sp2)
-
-        tgt = list(dom)
-        tgt[sp1] = dom[sp2]
-        tgt[sp2] = dom[sp1]
-        assert op.target == ift.DomainTuple.make(tgt)
-        ift.extra.consistency_check(op)

test/test_operators/test_einsum.py

@@ -92,6 +92,44 @@ def test_linear_einsum_contraction(space1, space2, dtype, n_invocations=10):
         assert_allclose(le.adjoint(r_adj).val, le_ift.adjoint(r_adj).val)
 
 
+class _SwitchSpacesOperator(ift.LinearOperator):
+    """Operator to permutate the domain entries of fields.
+
+    Exchanges the entries `space1` and `space2` of the input's domain.
+    """
+    def __init__(self, domain, space1, space2=0):
+        self._capability = self.TIMES | self.ADJOINT_TIMES
+        self._domain = ift.DomainTuple.make(domain)
+
+        n_spaces = len(self._domain)
+        if space1 >= n_spaces or space1 < 0 \
+          or space2 >= n_spaces or space2 < 0:
+            raise ValueError("invalid space value")
+
+        tgt = list(self._domain)
+        tgt[space2] = self._domain[space1]
+        tgt[space1] = self._domain[space2]
+        self._target = ift.DomainTuple.make(tgt)
+
+        dom_axes = self._domain.axes
+        tgt_axes = self._target.axes
+
+        self._axes_dom = dom_axes[space1] + dom_axes[space2]
+        self._axes_tgt = tgt_axes[space2] + tgt_axes[space1]
+
+    def apply(self, x, mode):
+        self._check_input(x, mode)
+
+        if mode == self.TIMES:
+            val = np.moveaxis(x.val, self._axes_dom, self._axes_tgt)
+            dom = self._target
+        else:
+            val = np.moveaxis(x.val, self._axes_tgt, self._axes_dom)
+            dom = self._domain
+
+        return ift.Field(dom, val)
+
+
 def test_multi_linear_einsum_outer(
     space1, space2, dtype, n_invocations=10, ntries=100
 ):
@@ -116,7 +154,7 @@ def test_multi_linear_einsum_outer(
         ift.full(mf_dom["dom01"], 1.), ift.DomainTuple.make(mf_dom["dom02"][1])
     )
     # SwitchSpacesOperator is equivalent to LinearEinsum with "ij->ji"
-    mle_ift = ift.SwitchSpacesOperator(
+    mle_ift = _SwitchSpacesOperator(
         outer_i.target, 1
     ) @ outer_i @ ift.FieldAdapter(mf_dom["dom01"], "dom01") * ift.FieldAdapter(
         mf_dom["dom02"], "dom02"