Generating images as tangos properties

[Martin-Rey]

I have been trying to generate dark matter projected images around my halos of interest (with pynbody handler). Following the baryonic example given in https://github.com/pynbody/tangos/blob/811444efd5b06c080470e0672da59784fffc87e9/tangos/properties/pynbody/images.py#L4

I derived the following class:

class DarkMatterCloseUpImages(PynbodyPropertyCalculation):
    names = "dm_closeup_image_x", "dm_closeup_image_x"

    @classmethod
    def plot_extent(cls):
        return 1000.0

    @classmethod
    def plot_xlabel(cls):
        return "x/a kpc"

    @classmethod
    def plot_ylabel(cls):
        return "y/a kpc"

    def plot_clabel(self):
         return r"M$_{\odot}$ kpc$^{-2}$"

    def requires_property(self):
        return ["shrink_center", "max_radius"]

    def region_specification(self, db_data):
        import pynbody
        side = pynbody.array.SimArray(4 * db_data['max_radius'], "kpc")
        xcenter = db_data['shrink_center'][0]
        ycenter = db_data['shrink_center'][1]
        zcenter = db_data['shrink_center'][2]

        x1 = xcenter - side/2
        y1 = ycenter - side/2
        z1 = zcenter - side/2
        x2 = xcenter + side/2
        y2 = ycenter + side/2
        z2 = zcenter + side/2
        return pynbody.filt.Cuboid(x1, y1=y1, z1=z1, x2=x2, y2=y2, z2=z2)

    def _render_projected(self, f, size):
        import pynbody.plot
        im = pynbody.plot.sph.image(f, qty='rho', width=size, units="Msol kpc^-2", noplot=True)
        return im

    @centred_calculation
    def calculate(self, particle_data, existing_properties):
        plot_size = 4 * existing_properties['max_radius']

        # This is the key. Pynbody tries to calculate the smooth array on the full snapshot for some reason ....
        particle_data.d['smooth']

        im_z = self._render_projected(particle_data.dm, plot_size)
        with particle_data.rotate_y(90):
            im_x = self._render_projected(particle_data.dm, plot_size)
        return im_z, im_x

My first question is on plot_extent(). Let's say I want to make a plot of the surroundings of the halo up to 4 radii. How can this be incorporated in the plot_extent that is a class method ?

The second and main problem is that this calculation fails every time throwing:

Traceback (most recent call last):
  File "/Users/mrey/Documents/tangos/tangos/tools/property_writer.py", line 382, in _get_property_value
    result = property_calculator.calculate(snapshot_data, db_data)
  File "/Users/mrey/Documents/tangos/tangos/properties/pynbody/centring.py", line 60, in new_fn
    return fn(self, halo, existing_properties)
  File "/Users/mrey/Documents/tangos/tangos/properties/pynbody/images.py", line 103, in calculate
    particle_data.d['smooth']
  File "/Users/mrey/anaconda3/lib/python3.6/site-packages/pynbody-0.43-py3.6-macosx-10.7-x86_64.egg/pynbody/snapshot/__init__.py", line 264, in __getitem__
    return self._get_array_with_lazy_actions(i)
  File "/Users/mrey/anaconda3/lib/python3.6/site-packages/pynbody-0.43-py3.6-macosx-10.7-x86_64.egg/pynbody/snapshot/__init__.py", line 357, in _get_array_with_lazy_actions
    self.__derive_if_required(name)
  File "/Users/mrey/anaconda3/lib/python3.6/site-packages/pynbody-0.43-py3.6-macosx-10.7-x86_64.egg/pynbody/snapshot/__init__.py", line 372, in __derive_if_required
    self._derive_array(name)
  File "/Users/mrey/anaconda3/lib/python3.6/site-packages/pynbody-0.43-py3.6-macosx-10.7-x86_64.egg/pynbody/snapshot/__init__.py", line 1942, in _derive_array
    self.base._derive_array(array_name, self._unifamily)
  File "/Users/mrey/anaconda3/lib/python3.6/site-packages/pynbody-0.43-py3.6-macosx-10.7-x86_64.egg/pynbody/snapshot/__init__.py", line 1774, in _derive_array
    self.base._derive_array(array_name, fam)
  File "/Users/mrey/anaconda3/lib/python3.6/site-packages/pynbody-0.43-py3.6-macosx-10.7-x86_64.egg/pynbody/snapshot/__init__.py", line 1471, in _derive_array
    result = fn(self[fam])
  File "/Users/mrey/anaconda3/lib/python3.6/site-packages/pynbody-0.43-py3.6-macosx-10.7-x86_64.egg/pynbody/sph/__init__.py", line 137, in smooth
    self.kdtree.populate('hsm', config['sph']['smooth-particles'])
  File "/Users/mrey/anaconda3/lib/python3.6/site-packages/pynbody-0.43-py3.6-macosx-10.7-x86_64.egg/pynbody/sph/kdtree.py", line 114, in populate
    smx = kdmain.nn_start(self.kdtree, int(nn), self.boxsize)
ValueError: The particles span a region larger than the specified boxsize

I tracked it down to a problem when calculating the smooth pynbody array rather than the actual plotting itself. I think pynbody or tangos tries to calculate the smooth array on the ancestor snapshot rather than the particle data. Any insight would be much appreciated.

[apontzen]

• As you say, plot_extent() is a class method so at present there is no way for it to vary according to other properties. Removing this limitation in a backward-compatible way would be worth considering. In the meantime, worth remembering that plot_extent literally just affects plotting in the browser interface so if you change the axis labels to be "x/virial_radius" etc, the situation can still be made clear.
• Calculating the smooth array should be fine so this seems like a bug in the pynbody server implementation. However I have not been able to reproduce this as yet.

[Martin-Rey]

After a bit of investigation:

1. The problem is generated by this check in the C routines of pynbody: https://github.com/pynbody/pynbody/blob/cc70f0dd5cd7649ee977254c2e223f8271a671c5/pynbody/sph/smooth.cpp#L11-L27

2. For some reason, when given full access to the ancestor snapshot, the smooth array is always calculated across the entire snapshot, not just on a sub snap. I verified this in an independent ipython notebook:

   
   In [1]: s51 = pynbody.load("/Users/mrey/Documents/simulations/ramses/DM_only/ramses_dmo_L50_N256_ZOOM4_Halo740/output_00051/")
   ...: s51.physical_units()
   ...: 
   /Users/mrey/anaconda3/lib/python3.6/site-packages/pynbody-0.44-py3.6-macosx-10.7-x86_64.egg/pynbody/snapshot/ramses.py:722: UserWarning: Using field at offset 9 to distinguish stars. If this is wrong, try editing your config.ini, section [ramses], entry particle-distinguisher.
   warnings.warn("Using field %s to distinguish stars. If this is wrong, try editing your config.ini, section [ramses], entry particle-distinguisher."%pb_name)

In [2]: pynbody.analysis.halo.center(s51.halos()[10], vel=False) Out[2]: <pynbody.transformation.GenericTranslation at 0x18173b2f98>

In [3]: region = s51[pynbody.filt.Sphere(350)]

In [4]: s51.keys() Out[4]: ['hop_grp', 'pos', 'x', 'y', 'z', 'vel', 'vx', 'vy', 'vz', 'mass']

In [5]: region.keys() Out[5]: ['hop_grp', 'pos', 'x', 'y', 'z', 'vel', 'vx', 'vy', 'vz', 'mass']

In [6]: region.d['smooth'] Out[6]: SimArray([ 8.26333618, 17.87956619, 20.67194176, ..., 22.76054001, 17.42362976, 24.66303253], 'kpc')

In [7]: region.keys() Out[7]: ['hop_grp', 'pos', 'x', 'y', 'z', 'vel', 'vx', 'vy', 'vz', 'mass', 'smooth']

In [8]: s51.keys() Out[8]: ['hop_grp', 'pos', 'x', 'y', 'z', 'vel', 'vx', 'vy', 'vz', 'mass', 'smooth']

The smooth array has been calculated across the full s51 snapshot even though I only requested it on the spherical region. I don't know if this intentional behaviour or not.

3. The check fails not because of a boxsize problem but because of the state of the position arrays:
- The @centred calculation centres and wraps only the selected particle of the region specification meaning that particle_data['pos'] spawns [-sthing, +sthing]
- Pynbody then tries to calculate the smooth array on the full snapshot. The previously wrapped particles spawn [-sthing, +sthing], while the remaining spawn [0, boxsize]. Hence the full snapshot positions spawn [-sthing, boxsize] failing the test to fit in a boxsize.

4. One funny side effect is that removing access to the full ancestor by running calculations in server mode removes that behaviour and produces sensible images, which exactly what I want. 

[apontzen]

Thanks, this is a really helpful analysis. (2) is actually intentional behaviour; for example if you are creating an image of a subregion you need smooth the particles on the boundary to get the strictly correct internal smoothing properties. But in this case it would be better if the behaviour were different and perhaps there is a way tangos can actually prevent access to the ancestor.

[apontzen]

It turns out that fixing this at the level of changing the behaviour of .ancestor or .base in a sensible, consistent way is extremely hard and likely to have side-effects that are undesirable.

I think, certainly for now, we have to deal with the conceptual wrinkle. On a practical front, PR #51 should fix the images problem - please let me know.

[Martin-Rey]

I can confirm that this get ride of the bug in default load-mode and still works in server mode.

On a conceptual level, what does halo.ancestor represents in server/partial loading mode since the ancestor snapshot is not meant to be fully available ?

[Martin-Rey]

Closing this issue following PR #51