Speed up perform_pop_syn by ~50%

[MichaelMedford]

This branch implements two changes:

• How read_ind works in ebf.py
• How cKDTrees are constructed in perform_pop_syn for finding the extinction to compact objects

read_ind The current version of perform_pop_syn spends about half of its time reading in stars from the ebf file with the read_ind command, which asks for specific indices for a field in the ebf file. Our loop structure means that we are always reading in stars with identical popid and similar ages (and in the future similar metallicities). After some investigating it turns out that these elements are stored in memory in similar locations.

read_ind currently reads in each index one at a time, executing a copy command for every element in the index array. This branch groups the indices into blocks of contiguous elements before reading them in. Then each block is loaded in all at once, greatly reducing the number of copy commands. This simple change dropped the execution time of 'perform_pop_syn` by about ~30% and makes no changes at all to what data is read in. Identical behavior but with much better performance!

Compact Object Extinction KDTree Compact objects are assigned extinction by giving them the same extinction as the nearest star in 3D space. In order to find this star, a KDTree is constructed and the location of the compact object is queried against the tree. Currently a new KDTree is constructed for each batch of random stars in the last sub-bin loop. This branch creates a single KDTree before starting the loops over sub-bins. The stars in this KDTree are made up of the same number of elements as go into each sub-bin, but randomly drawn from across the population. Then each sub-bin of compact objects is queried against this tree. Speed up results vary, as the number of sub-bins required greatly depends on the stellar density and area of the PopSyCLE run. However I've been seeing speed ups of perform_pop_syn execution time of 10%-20% with compact objects having nearly identical extinctions as before.

For a field of (l,b,A) = (2, 1, 0.33):

N_comp_objs = 1584
median(exbv_master - exbv_pps) = 0.0
std(exbv_master - exbv_pps) = 0.165
std(sigmaclip(exbv_master - exbv_pps, low=5, high=5)[0]) = 0.0305

Plots showing the effects on the extinction to compact objects shown here.

As shown by both the plots and the sigma-clipped standard deviation, there are a few large outliers. However I have no way of knowing at this time whether that is error introduced by this method, or the noise that was present in the previous method. I am open to suggestions on how to tease apart that difference.

[MichaelMedford]

@jluastro @caseylam Could I please get some feedback on these improvements? I think that this branch as written is ready to be merged into master, and will help greatly improve the execution time of perform_pop_syn.

[caseylam]

Nice!

However I have no way of knowing at this time whether that is error introduced by this method, or the noise that was present in the previous method. I am open to suggestions on how to tease apart that difference. I'd suggest looking at the WD-star pairs, and then looking at how the extinction differs between separation in 3-D space (and also broken into perpendicular/parallel components to LOS), and galactic latitude/longitude. Then maybe some pattern will pop out about the really big outliers. My guess is that those nearby stars happened to have a very different extinction. The method of getting the extinction is in Section 2.6: https://iopscience.iop.org/article/10.1088/0004-637X/730/1/3/pdf

[MichaelMedford]

Doh, I meant the previous comment to be part of this review.

I also meant to ask, are the plots you showed in the confluence page sigma clipped or not?

The median difference is centered at zero which is good. If the differences are symmetric, then I think it is okay.

The plots are not sigma clipped. The differences are symmetric. I have found that they have a spatial orientation, but I have not yet determined why. Here is a plot of the log difference in extinctions, as well as the log of the difference.

[MichaelMedford]

@caseylam @jluastro Can we pull the trigger on this PR? The differences between compact object extinction are both small and symmetric, so I think we're not loosing any precision by adopting this faster method. I would like to have it in here before I start the PR for the new galaxia_params branch the supports the new galaxia installation.

[jluastro]

Doh, I meant the previous comment to be part of this review. I also meant to ask, are the plots you showed in the confluence page sigma clipped or not? The median difference is centered at zero which is good. If the differences are symmetric, then I think it is okay.

The plots are not sigma clipped. The differences are symmetric. I have found that they have a spatial orientation, but I have not yet determined why. Here is a plot of the log difference in extinctions, as well as the log of the difference.

Can you remind me what the top and bottom are for these figures? The bottom has some spurious swaths.

[MichaelMedford]

The top and bottom rows are two different ways to representing the difference between the old and new methods. In the top plot, the color corresponds to log(exbv_master / exbv_pps). Values close to 0 show that the two methods are converging to the same result. I think this is a better indicator of how the method is doing. The color in the bottom plot is log(exbv_master - exbv_pps), which is a bit more contrived. Either way, the histograms of the differences show that these are very much outliers and that the std(exbv_master - exbv_pps) = 0.165.