Bug in position assignment of synthetic centrals

[aphearin]

Line 278-280 in extend_subhalo_mpeak_range.py assign random xyz spatial positions to synthetic galaxies:

gals_x = np.random.uniform(np.min(healpix_mock['target_halo_x']), np.max(healpix_mock['target_halo_x']), num_needed)

gals_y = np.random.uniform(np.min(healpix_mock['target_halo_y']), np.max(healpix_mock['target_halo_y']), num_needed)

gals_z = np.random.uniform(np.min(healpix_mock['target_halo_z']), np.max(healpix_mock['target_halo_z']), num_needed)

These lines of code should be assigning a random xyz-position to a galaxy within a rectangle that encloses the healpixel, so that the subsequent call to mask_galaxies_outside_healpix on line 281 throws out those proposed points that happen to lie outside the healpixel, and the enclosing while loop continues in this fashion until all points lie within the healpixel.

This algorithm has two deficiencies.

1. By using the xyz min/max of halos within the healpixel, the enclosing rectangle is actually not enclosing, and so the resulting points do not actually uniformly span the healpixel, as is intended.
2. More critically, for healpixels with a single halo there is a critical failure: all synthetics are assigned to the exact same xyz position.

This algorithm should be improved prior to running SkySim5000.

Proposed patch:

A crude and simple solution would be to use (x.min() - 2, x.max()+2) as the enclosing rectangle: this will solve the critical issue of position multiplicity immediately, with the drawback that the while loop will typically need to churn through many iterations before termination.

Proposed solution:

We should really be using a smarter algorithm to randomly draw positions that fall within the healpixel.

CC @evevkovacs @patricialarsen @bleeml

[evevkovacs]

@aphearin If the healpixel has lots of halos, I don't think that issue 1 is a big problem. The healpixel has a diamond shape, and the rectangular box constructed from the halo min/max positions will lie outside the healpix boundary for most of the area of the healpixel (but it is true that the very corners of the diamond may not be uniformly filled) But for the case where there are very few halos in the healpixel, then indeed, the rectangular box will be too small. However these are edge cases, where the octant overlaps a tiny portion of a healpixel. So if we randomize the synthetic galaxies throughout the healpixel in these cases, then they will actually be outside of the octant (otherwise there would have been more halos in the first place). We are assigning the number of synthetics to the healpixel based on the area of the healpixel (see lines 189-192), so that number is wrong in these edge cases too.

[evevkovacs]

I will submit a PR for the skysim_5000 branch which contains the new code to deal with this problem. I am computing the boundaries of the healpixel and using that, together with the redshifts of the lightcone shell to define a volume in which to generate trial synthetic galaxy positions. The code then checks to see if the trial position is within the healpixel. The only complication occurs for healpixels overlapping the octant. There is an additional mask to make sure that the coordinates of the galaxy also lie within the octant. The number of synthetics generated is reduced for edge cases.