Improve SOAP runtime

[robjmcgibbon]

SOAP is taking a long time to run on COLIBRE. This PR adds some recording of how long we are spending processing each subhalo, and which properties are dominating the runtime. I will try and add comments for any speedups I add.

For satellites (for which we skip SO properties) we were still checking if the load region around them dropped below the target density of the lowest density SO. This had a significant effect on runtime, especially for high resolution COLIBRE runs where we identify lots of substructure.

EDIT There's a bunch of other changes here too, but they're all described by the comments. I'd like this branch to be the one that we use when we publish the SOAP paper. TODO before then:

• Refactor so everything is packaged more cleanly
• Update main README & parameter file README
• Add mpi-based compression scripts
• Subsampled inertia tensors
• Clean up github issues/PRs
• Test fully against master

[robjmcgibbon]

We were calling astropy every time we wanted the age of a stellar particle. I've move to using a lookup table. I verified it should be accurate using the following script

import matplotlib.pyplot as plt
import numpy as np
import unyt

z_now = 0
t_now = unyt.unyt_quantity.from_astropy(
    cosmology.lookback_time(z_now)
).to("Myr")

# Create a lookup table for z=49 to z=0
a_lookup = np.linspace(1/50, 1, 1000)
z = (1.0 / a_lookup) - 1.0
# Remember that lookback time runs backwards
age_lookup = unyt.unyt_array.from_astropy(
    cosmology.lookback_time(z)
).to('Myr') - t_now

# Compare astropy vs interpolated ages from z=29 to z=0
a_check = np.linspace(1/30, 1, 3000)
z_check = (1.0 / a_check) - 1.0
age_check = unyt.unyt_array.from_astropy(
    cosmology.lookback_time(z_check)
).to('Myr') - t_now
age_interpolate = np.interp(a_check, a_lookup, age_lookup)

i_max = np.argmax(np.abs(age_check - age_interpolate))
max_delta = np.max(np.abs(age_check - age_interpolate))
print(f'Max delta: {max_delta:.3g} occurs at age {age_check[i_max]:.5g}')

# Plot results
# This has a strange beat frequency effect
fig, ax = plt.subplots(1)
ax.scatter(age_check, age_check-age_interpolate, s=0.1, lw=0)
ax.plot(age_lookup, np.zeros_like(age_lookup), '.', color='C1')
ax.set_xlabel('True age [Myr]')
ax.set_ylabel('True age - Interpolated age [Myr]')
plt.savefig('test_age_interpolation.png', dpi=200)
plt.close()

[robjmcgibbon]

I changed the documentation so that properties now link to any footnotes they have

[robjmcgibbon]

I removed a redundant sort from SO_properties which I had introduced for the concentration calculation.

[robjmcgibbon]

I've switched off iterative calculation of inertia tensors for the time being, although we may instead use a subsample of the particles to determine the inertia tensor.

[robjmcgibbon]

I have merged in #126, which allows for SOAP to restart.

I have updated the property table and parameter files so that filters (e.g. basic, general, stars) are now specified in the parameter file.

I'm removed some leftover VR FOFSubhalo code

[robjmcgibbon]

I'm now sorting halos on each chunk so that the largest objects are processed first

[robjmcgibbon]

For projected apertures and exclusive spheres we were often recomputing properties using the same particles each time. This occurred for small subhalos, where their EncloseRadius (the distance from the halo centre to the furthest particle) was smaller than multiple apertures (e.g. if EncloseRadius=30kpc, then the 50kpc & 100kpc ExlcusiveSpheres would be the same). In these cases we now just copy across the values from the previous aperture, as they would be the same.

There is one wrinkle, which is that the area of the circle we use for calculating the ProjectedInertiaTensor depends on the size of the apertures. This means that for the later iterations of the inertia tensor calculation there can be particles excluded from the smaller apertures but not the larger ones, even if the same particles are being considered in both cases. This should cause very minor differences to the output value. However, I have added a new parameter called strict_halo_copy (default False) which will prevent copying occurring for any properties like this.

COLIBRE has quite a large number of high z snapshots. This means we were spending a lot of time on InclusiveSpheres. I've added a parameter that allows halo to be skipped if the InclusiveSphere is larger than the EncloseRadius. This parameter defaults to False, but we will set it to True for COLIBRE.

[robjmcgibbon]

I've added some scripts that can be used to analyse the runtime of SOAP. There's lots of plots, such as

[robjmcgibbon]

I added a modulus to the COM calculation, since we were getting some values outside the box.

I merged Victor's script for creating symmetric matrices from the flattened ones we store in SOAP (#103). I'm about to reorganise the repo, and want to group all the auxiliary scripts together

[robjmcgibbon]

I've made the min_read_radius parameter optional. We read in a value for each halo finder from the input halo catalogues, so we don't really need this parameter anymore. It is still useful if we're calculating large SOs (e.g. 50_crit), so I have left the functionality in place.

I've added a FLAMINGO VR test

[robjmcgibbon]

I had to spread out the files more evenly when loading the HBT catalogues. Load imbalance was causing the group membership script to crash when running on some of the large COLIBRE runs.

[robjmcgibbon]

I've changed the compression script to use MPI instead of multiprocessing. This was required for the 2p8 flamingo catalogues. I've added more information to the README about how to run the compression.