Split `positions_from_delta()` into individual functions

[ntessore]

Add your issue here

The current positions_from_delta() function carries out a large number of steps:

• Apply bias model to given delta map.
• Convert biased delta to mean number count in each map pixel.
• Apply a visibility map.
• Compute observed number counts from a distribution (e.g., Poisson).
• Accumulate the galaxies in the number count map into batches of a given size (e.g., 1_000_000)
• Sample random positions for each galaxy within its pixel.

This structure is very rigid, and makes it difficult to implement things such as #136.

We should split this up into individual functions, so that they can be mixed and matched with different implementations:

for i, delta in enumerate(matter):
    # mean number count in each pixel
    nbar = glass.nbar_from_density(arcmin2=ngal[i], nside=nside)
    # galaxy density from linear bias model, could be a complicated model
    gal = nbar * (1 + beff[i] * delta)
    # apply visibility map for this shell, could be a complicated model
    gal *= vmap[i]
    # Poisson-sample galaxy indices from map in batches
    for ipix in glass.batched_poisson(gal, rng=rng):
        # randomly sample galaxy positions within each pixel
        lon, lat = glass.random_healpix_positions(nside, ipix)

This needs a bit of experimentation for multi-dimensional arrays, but in principle seems feasible.

[paddyroddy]

We should split this up into individual functions, so that they can be mixed and matched with different implementations:

for i, delta in enumerate(matter): # mean number count in each pixel nbar = glass.nbar_from_density(arcmin2=ngal[i], nside=nside) # galaxy density from linear bias model, could be a complicated model gal = nbar * (1 + beff[i] * delta) # apply visibility map for this shell, could be a complicated model gal *= vmap[i] # Poisson-sample galaxy indices from map in batches for ipix in glass.batched_poisson(gal, rng=rng): # randomly sample galaxy positions within each pixel lon, lat = glass.random_healpix_positions(nside, ipix)

This needs a bit of experimentation for multi-dimensional arrays, but in principle seems feasible.

Just to clarify on this, are you saying this is how you envisage the function being split up? Been trying to map the current function to this in my head, and it's not entirely clear to me.

[ntessore]

Just to clarify on this, are you saying this is how you envisage the function being split up? Been trying to map the current function to this in my head, and it's not entirely clear to me.

Yes, something along those lines. The functional split in the current positions_from_delta() is roughly this:

https://github.com/glass-dev/glass/blob/fa012ac54454b19dc4f83bf1a8b758896c308555/glass/points.py#L229-L240

Figure out how many "populations" of objects we are dealing with, by broadcasting all inputs together. This is probably going to be the most tricky part to generalise neatly when not doing all subsequent steps in one function.

https://github.com/glass-dev/glass/blob/fa012ac54454b19dc4f83bf1a8b758896c308555/glass/points.py#L242-L243

This loops over the leading dimensions ("populations") and processes each "population" independently. This will need some thought on how to do it without having this big loop over everything. For the time being, imagine we are dealing with just one population of objects.

https://github.com/glass-dev/glass/blob/fa012ac54454b19dc4f83bf1a8b758896c308555/glass/points.py#L244-L249

Applies the bias model to delta.

https://github.com/glass-dev/glass/blob/fa012ac54454b19dc4f83bf1a8b758896c308555/glass/points.py#L251-L257

Computes the expected number of objects per pixel nbar = ARCMIN2_SPHERE / n.size * ngal[k]. Then turns the biased delta_b in to an expected number count n = nbar * (1 + delta_b).

https://github.com/glass-dev/glass/blob/fa012ac54454b19dc4f83bf1a8b758896c308555/glass/points.py#L259-L261

Apply a visibility by multiplying it with the expected number count n = vis * n.

https://github.com/glass-dev/glass/blob/fa012ac54454b19dc4f83bf1a8b758896c308555/glass/points.py#L263-L267

Sample the actual number of galaxies in each pixel from the Poisson distribution.

https://github.com/glass-dev/glass/blob/fa012ac54454b19dc4f83bf1a8b758896c308555/glass/points.py#L269-L315

Sample the individual galaxies in each pixel, randomly distributed over sub pixels, in batches.

[mwiet]

I have a similar local implementation of this which creates a new classes within glass.observations that act like a visibility mask for a given combination of tomographic bins and redshift shells: glass.observations.angular_variable_depth_mask and glass.observations.angular_los_variable_depth_mask. The class constructs a lookup table based on a model of variable depth akin to SALMO. It is designed to output the ratio between the galaxy counts when variable depth occurs and the galaxy counts, so it can be applied as a factor to the visibility mask when sampling galaxies from a matter field with glass.points.positions_from_delta.

The implementation adds two novel features:

• The interpolation between the variable depth tracer variable and the change in galaxy density can be interpolated with any function (previously, in SALMO only linear interpolation was possible).
• The effect of variable depth along the line-of-sight and in the angular direction can be modelled as fully independent (previously, in SALMO the selection of galaxy shapes and galaxy redshifts was always assumed to be the same).

I have documented and tested the code (incl. tests on the KiDS-1000). It seems stable and consistent with previous results, so I have created a pull request. Please provide any feedback on any additional functionalities and/or changes to make it fit within the rest of the GLASS code base.