Currently, the project plan involves:

1. assign the particles to a regular gridded cube of size N
2. smooth the field
3. apply the NEXUS method to assign each grid cell to one of the following structure classes: 0-voids, 1-walls, 2-filaments, 3-halos/nodes
4. compute the power spectrum (and PDF) of the DM-only grid, but with some cells swapped out by the corresponding hydro cells.

During today's telecon, I talked about an alternative to step-3; instead of applying NEXUS, identify the regions that are in the 4 different quartiles of the density fields (top 25% densest, 25-50% quantile density,.. and 25% emptiest regions); which can be done by rank order all the pixels, and take the top 1/4, next 1/4 values, etc.

This issue will attempt to address this alternative method (for step 3)

pseudo code:

rho = rho_DM ### where rho_DM is a grid of shape N^3

structure_class = NEXUS(rho_DM)

## structure_class is also a grid of shape N^3, where each cell has value=0,1,2,or 3, for either (1) the 4 different structures (void, walls, filaments, halos) or (2) different quantile values, (0-top 25% densest regions, 1-[25%-50%] quantile... etc)

### if we want to swap out filaments cells, i.e. only look at baryonic effects on filaments, we do
rho[np.where (structure_class==2)] = rho_hydro[np.where (structure_class==2)]

## finally, compare the resulting power spectrum to 
ratio = power_spectrum(rho)/power_spectrum(rho_DM)