alternative solution

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

@@ -25,6 +25,19 @@ from .scaling_operator import ScalingOperator
 
 class SandwichOperator(EndomorphicOperator):
     """Operator which is equivalent to the expression `bun.adjoint*cheese*bun`.
+    """
+
+    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
         ----------
@@ -33,16 +46,16 @@ class SandwichOperator(EndomorphicOperator):
         cheese: EndomorphicOperator
             the cheese part
         """
-
-    def __init__(self, bun, cheese=None):
-        super(SandwichOperator, self).__init__()
-        self._bun = bun
         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):
@@ -56,10 +69,6 @@ class SandwichOperator(EndomorphicOperator):
         return self._op.apply(x, mode)
 
     def draw_sample(self, from_inverse=False, dtype=np.float64):
-        # Drawing samples from diagonal operators is easy (inverse is possible)
-        if isinstance(self._op, (ScalingOperator, DiagonalOperator)):
-            return self._op.draw_sample(from_inverse, dtype)
-
         # Inverse samples from general sandwiches is not possible
         if from_inverse:
             raise NotImplementedError(