Deal with interlacing more efficiently.

This moves the FFT out of the main particle loop for the interlaced
painter.

It fixes the previous wrong result. But I don't see why...

nbodykit/base/catalogmesh.py

@@ -330,12 +330,8 @@ def to_real_field(self, out=None, normalize=True):
         # since out may have non-zero elements, messing up our interlacing sum
         if self.interlaced:
 
-            # whether we can re-use "toret" workspace
-            if out is None:
-                real1 = toret
-            else:
-                real1 = RealField(pm)
-                real1[:] = 0
+            real1 = RealField(pm)
+            real1[:] = 0
 
             # the second, shifted mesh (always needed)
             real2 = RealField(pm)
@@ -407,21 +403,23 @@ def to_real_field(self, out=None, normalize=True):
                 # paint to two shifted meshes
                 pm.paint(p, mass=w * v, resampler=paintbrush, hold=True, out=real1)
                 pm.paint(p, mass=w * v, resampler=paintbrush, transform=shifted, hold=True, out=real2)
-                c1 = real1.r2c()
-                c2 = real2.r2c()
 
-                # and then combine
-                for k, s1, s2 in zip(c1.slabs.x, c1.slabs, c2.slabs):
-                    kH = sum(k[i] * H[i] for i in range(3))
-                    s1[...] = s1[...] * 0.5 + s2[...] * 0.5 * numpy.exp(0.5 * 1j * kH)
+        if self.interlaced:
+            # compose the two interlaced fields into the final result.
+            c1 = real1.r2c()
+            c2 = real2.r2c()
+
+            # and then combine
+            for k, s1, s2 in zip(c1.slabs.x, c1.slabs, c2.slabs):
+                kH = sum(k[i] * H[i] for i in range(3))
+                s1[...] = s1[...] * 0.5 + s2[...] * 0.5 * numpy.exp(0.5 * 1j * kH)
 
-                # FFT back to real-space
-                # NOTE: cannot use "toret" here in case user supplied "out"
-                c1.c2r(real1)
+            # FFT back to real-space
+            # NOTE: cannot use "toret" here in case user supplied "out"
+            c1.c2r(real1)
 
-                # need to add to the returned mesh if user supplied "out"
-                if real1 is not toret:
-                    toret[:] += real1[:]
+            # need to add to the returned mesh if user supplied "out"
+            toret[:] += real1[:]
 
         # unweighted number of objects
         N = pm.comm.allreduce(Nlocal)