From 3e609890c55f185be8e6812c2e2415fce8d64b35 Mon Sep 17 00:00:00 2001 From: Yu Feng Date: Mon, 23 Oct 2017 13:21:48 -0700 Subject: [PATCH 1/4] Add a test case that fails the cmean of the interlaced tomesh. It gives a bogus powerspectrum -- a few times larger. --- .../source/catalogmesh/tests/test_species.py | 39 +++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/nbodykit/source/catalogmesh/tests/test_species.py b/nbodykit/source/catalogmesh/tests/test_species.py index 14e47a4cd..8f10a56f5 100644 --- a/nbodykit/source/catalogmesh/tests/test_species.py +++ b/nbodykit/source/catalogmesh/tests/test_species.py @@ -77,3 +77,42 @@ def test_paint(comm): # must be the same assert_allclose(combined.value, (real1.value + real2.value)/norm, atol=1e-5) + +@MPITest([1, 4]) +def test_paint_interlaced(comm): + + CurrentMPIComm.set(comm) + + # the test case fails only if there is enough particles to trigger + # the second loop of the interlaced painter; these parameters will do it. + + # the catalog + source1 = UniformCatalog(nbar=1e-0, BoxSize=111, seed=111) + source2 = UniformCatalog(nbar=1e-0, BoxSize=111, seed=111) + source1['Weight'] = 1.0 + source2['Weight'] = 0.1 + cat = MultipleSpeciesCatalog(['data', 'randoms'], source1, source2, use_cache=False) + + # the meshes + mesh = cat.to_mesh(Nmesh=32, interlaced=True) + mesh1 = source1.to_mesh(Nmesh=32, interlaced=True) + mesh2 = source2.to_mesh(Nmesh=32, interlaced=True) + + # paint + real1 = mesh1.to_real_field() + real2 = mesh2.to_real_field() + assert_allclose(real1.cmean(), 1.0) + assert_allclose(real2.cmean(), 1.0) + + # un-normalize real1 and real2 + real1[:] *= real1.attrs['num_per_cell'] + real2[:] *= real2.attrs['num_per_cell'] + norm = real1.attrs['num_per_cell'] + real2.attrs['num_per_cell'] + + # the combined density field + #combined = mesh.to_real_field() + combined = mesh.paint() + + assert_allclose(combined.cmean(), 1.0) + # must be the same + assert_allclose(combined.value, (real1.value + real2.value)/norm, atol=1e-5) From a4bd776ade4a0d5ca1d0648807842c93817760d9 Mon Sep 17 00:00:00 2001 From: Yu Feng Date: Mon, 23 Oct 2017 13:22:55 -0700 Subject: [PATCH 2/4] 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 | 34 ++++++++++++++++------------------ 1 file changed, 16 insertions(+), 18 deletions(-) diff --git a/nbodykit/base/catalogmesh.py b/nbodykit/base/catalogmesh.py index 4b1d9e20c..6b646a1eb 100644 --- a/nbodykit/base/catalogmesh.py +++ b/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) From 387d4c0382e33b4ed23f52685f277b2c50914b49 Mon Sep 17 00:00:00 2001 From: Yu Feng Date: Mon, 23 Oct 2017 13:30:43 -0700 Subject: [PATCH 3/4] mark the end of the loop explicitly. --- nbodykit/base/catalogmesh.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/nbodykit/base/catalogmesh.py b/nbodykit/base/catalogmesh.py index 6b646a1eb..fca8e6e9d 100644 --- a/nbodykit/base/catalogmesh.py +++ b/nbodykit/base/catalogmesh.py @@ -404,7 +404,12 @@ def to_real_field(self, out=None, normalize=True): 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) - if self.interlaced: + # now the loop over particles is done + + if not self.interlaced: + # nothing to do, toret is already filled. + pass + else: # compose the two interlaced fields into the final result. c1 = real1.r2c() c2 = real2.r2c() From 13e7da94c8883a2820868a3b033db9cd31748467 Mon Sep 17 00:00:00 2001 From: Yu Feng Date: Mon, 23 Oct 2017 13:32:48 -0700 Subject: [PATCH 4/4] convert the error to a warning. Sometimes we do push zero particles through the pipe due to fluctuations or weird cuts. The code shall not crash due to this. --- nbodykit/base/catalogmesh.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/nbodykit/base/catalogmesh.py b/nbodykit/base/catalogmesh.py index fca8e6e9d..2767d363a 100644 --- a/nbodykit/base/catalogmesh.py +++ b/nbodykit/base/catalogmesh.py @@ -2,6 +2,7 @@ from nbodykit.base.catalog import CatalogSource, CatalogSourceBase import numpy import logging +import warnings # for converting from particle to mesh from pmesh import window @@ -437,8 +438,10 @@ def to_real_field(self, out=None, normalize=True): # make sure we painted something! if N == 0: - raise ValueError(("trying to paint particle source to mesh, " - "but no particles were found!")) + warnings.warn(("trying to paint particle source to mesh, " + "but no particles were found!"), + RuntimeWarning + ) # shot noise is volume / un-weighted number shotnoise = numpy.prod(pm.BoxSize) / N