Merge branch 'more_samples_for_sandwiches' into 'NIFTy_4'

Allow inverse samples from sandwiches which are diagonal

See merge request ift/NIFTy!258

demos/krylov_sampling.py

@@ -31,7 +31,7 @@ d = R(s_x) + n
 
 R_p = R * FFT * A
 j = R_p.adjoint(N.inverse(d))
-D_inv = ift.SandwichOperator(R_p, N.inverse) + S.inverse
+D_inv = ift.SandwichOperator.make(R_p, N.inverse) + S.inverse
 
 
 N_samps = 200

demos/poisson_demo.py

@@ -80,12 +80,12 @@ if __name__ == "__main__":
     IC = ift.GradientNormController(name="inverter", iteration_limit=500,
                                     tol_abs_gradnorm=1e-3)
     inverter = ift.ConjugateGradient(controller=IC)
-    D = (ift.SandwichOperator(R, N.inverse) + Phi_h.inverse).inverse
+    D = (ift.SandwichOperator.make(R, N.inverse) + Phi_h.inverse).inverse
     D = ift.InversionEnabler(D, inverter, approximation=Phi_h)
     m = HT(D(j))
 
     # Uncertainty
-    D = ift.SandwichOperator(aHT, D)  # real space propagator
+    D = ift.SandwichOperator.make(aHT, D)  # real space propagator
     Dhat = ift.probe_with_posterior_samples(D.inverse, None,
                                             nprobes=nprobes)[1]
     sig = ift.sqrt(Dhat)

demos/wiener_filter_easy.py

@@ -51,7 +51,7 @@ if __name__ == "__main__":
     IC = ift.GradientNormController(name="inverter", iteration_limit=500,
                                     tol_abs_gradnorm=0.1)
     inverter = ift.ConjugateGradient(controller=IC)
-    D = (ift.SandwichOperator(R, N.inverse) + Sh.inverse).inverse
+    D = (ift.SandwichOperator.make(R, N.inverse) + Sh.inverse).inverse
     D = ift.InversionEnabler(D, inverter, approximation=Sh)
     m = D(j)
 

nifty4/library/poisson_energy.py

@@ -46,7 +46,7 @@ class PoissonEnergy(Energy):
 
         R1 = Instrument*Rho*ht
         self._grad = (phipos + R1.adjoint_times((lam-d)/(lam+eps))).lock()
-        self._curv = Phi_h.inverse + SandwichOperator(R1, W)
+        self._curv = Phi_h.inverse + SandwichOperator.make(R1, W)
 
     def at(self, position):
         return self.__class__(position, self._d, self._Instrument,

nifty4/library/wiener_filter_curvature.py

@@ -39,5 +39,5 @@ def WienerFilterCurvature(R, N, S, inverter):
     inverter : Minimizer
         The minimizer to use during numerical inversion
     """
-    op = SandwichOperator(R, N.inverse) + S.inverse
+    op = SandwichOperator.make(R, N.inverse) + S.inverse
     return InversionEnabler(op, inverter, S.inverse)

nifty4/operators/sandwich_operator.py

@@ -16,31 +16,46 @@
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
 # and financially supported by the Studienstiftung des deutschen Volkes.
 
+import numpy as np
+
+from .diagonal_operator import DiagonalOperator
 from .endomorphic_operator import EndomorphicOperator
 from .scaling_operator import ScalingOperator
-import numpy as np
 
 
 class SandwichOperator(EndomorphicOperator):
     """Operator which is equivalent to the expression `bun.adjoint*cheese*bun`.
-
-    Parameters
-    ----------
-    bun: LinearOperator
-        the bun part
-    cheese: EndomorphicOperator
-        the cheese part
     """
 
-    def __init__(self, bun, cheese=None):
+    def __init__(self, bun, cheese, op, _callingfrommake=False):
+        if not _callingfrommake:
+            raise NotImplementedError
         super(SandwichOperator, self).__init__()
         self._bun = bun
+        self._cheese = cheese
+        self._op = op
+
+    @staticmethod
+    def make(bun, cheese=None):
+        """Build a SandwichOperator (or something simpler if possible)
+
+        Parameters
+        ----------
+        bun: LinearOperator
+            the bun part
+        cheese: EndomorphicOperator
+            the cheese part
+        """
         if cheese is None:
-            self._cheese = ScalingOperator(1., bun.target)
-            self._op = bun.adjoint*bun
+            cheese = ScalingOperator(1., bun.target)
+            op = bun.adjoint*bun
         else:
-            self._cheese = cheese
-            self._op = bun.adjoint*cheese*bun
+            op = bun.adjoint*cheese*bun
+
+        # if our sandwich is diagonal, we can return immediately
+        if isinstance(op, (ScalingOperator, DiagonalOperator)):
+            return op
+        return SandwichOperator(bun, cheese, op, _callingfrommake=True)
 
     @property
     def domain(self):
@@ -54,8 +69,11 @@ class SandwichOperator(EndomorphicOperator):
         return self._op.apply(x, mode)
 
     def draw_sample(self, from_inverse=False, dtype=np.float64):
+        # Inverse samples from general sandwiches is not possible
         if from_inverse:
             raise NotImplementedError(
                 "cannot draw from inverse of this operator")
+
+        # Samples from general sandwiches
         return self._bun.adjoint_times(
             self._cheese.draw_sample(from_inverse, dtype))