Astropy vs CosmoloPy

[mjlarson]

While checking out the #25, I decided to give Astropy a try for the cosmology backend for Firesong. After a little bit of playing with it, there's now a branch designed to allow us to directly compare back and forth between astropy and cosmolopy. It introduces a switch like this:

x = firesong_simulation(".", 
                        filename=None, 
                        bins=100,
                        use_astropy=True, 
                        seed=1)

There's some unexplained oddness with the two packages. Using the Astropy Planck15 cosmology parameters in cosmolopy leads to an error while writing a custom LambdaCDM cosmology for Astropy using the existing default Firesong values in Evolution.py leads the process to hang indefinitely. I'm not sure why either is failing. Maybe there are different assumptions in the meaning of the parameters?

Even using the slightly different cosmologies, we get extremely similar output:

Firesong with cosmolopy!
##### FIRESONG initializing  - Calculating Neutrino CDFs #####
Standard candle sources+
Source evolution assumed: HB2006SFR
Local density of neutrino sources: 1e-09/Mpc^3
Total number of neutrinos sources in the Universe: 18460
Redshift range: 0 - 10.0
Standard Candle Luminosity: 3.6471e+52 erg/yr
##### FIRESONG initialization done #####

Actual diffuse flux simulated :
E^2 dNdE = 1.2047888712065207e-08 (E/100 TeV)^(-0.1299999999999999) [GeV/cm^2.s.sr]

vs

Firesong with Astropy!
##### FIRESONG initializing  - Calculating Neutrino CDFs #####
Standard candle sources+
Source evolution assumed: HB2006SFR
Local density of neutrino sources: 1e-09/Mpc^3
Total number of neutrinos sources in the Universe: 18438
Redshift range: 0 - 10.0
Standard Candle Luminosity: 3.6478e+52 erg/yr
##### FIRESONG initialization done #####

Actual diffuse flux simulated :
E^2 dNdE = 1.2059045448963603e-08 (E/100 TeV)^(-0.1299999999999999) [GeV/cm^2.s.sr]

Astropy seems to be consistently slower than cosmolpy by anywhere from 10-30%. That may be less of an issue now that things are running quickly, though.

%timeit firesong_simulation(".", filename=None, bins=100, density=1e-7, use_astropy=False, seed=1, verbose=False)
1.57 s ± 141 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit firesong_simulation(".", filename=None, bins=100, density=1e-7, use_astropy=True, seed=1, verbose=False)
1.97 s ± 363 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

That said, cosmolopy's original author stopped updates in 2012 and there have only been rare community updates since then. I also couldn't get it to install without forking my own version with edits in python 3.8.

So do we want to consider switching to Astropy instead of cosmolopy? Or do we try to adopt cosmolopy somehow for the extra speed and to ensure it's maintained for the future?

[ChrisCFTung]

I wrote some functions to calculate the distance. The version is available in the branch no_cosmolopy The good thing is the execution time is reduced by about quite a bit. The bad thing is that it won't support non LambdaCDM universe.

In [2]: %timeit firesong_simulation(".", filename=None, bins=100, density=1e-7, use_cosmolopy=False, seed=1, verbose=False)
457 ms ± 3.15 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [3]: %timeit firesong_simulation(".", filename=None, bins=100, density=1e-7, use_cosmolopy=True, seed=1, verbose=False)
689 ms ± 19.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

The output are also the same.

[itaboada]

Neither! We are going with a custom calculation of the cosmological functions we need.