icecube/firesong

FIRst Extragalactic Simulation Of Neutrino and Gamma-ray

Firesong

FIRst Extragalactic Simulation Of Neutrinos and Gamma-rays.

See the Docs

See the FIRESONG paper

Documentation for the astrophysics and cosmology of neutrino sources can be found on:

• René Rieman's PhD thesis, section 10.
• Theo Glauch's Master thesis, chapter 6.
• Chun Fai Tung's PhD thesis, chapter 4.

Set up

This package is developed for Python 3, and can be installed via pip:

pip install firesong

Or by downloading the repository and running:

python setup.py install

If you want to execute scripts of FIRESONG directly via shell command, you can specify where you would like the output of simulations to go, by default. In bash: export FIRESONG=/location/of/FIRESONG/.

Basic usage

Several scripts are provided:

• firesong/Firesong.py - Generates an instance of all neutrino sources in the Universe according to the parameters provided (e.g. local neutrino source density). Both the flux and the redshift of neutrino sources are calculated.
• firesong/FluxPDF.py - Generates the flux probability density distribution of a source. It complements Firesong.py because is it much faster for large source densities. Only the flux of neutrino sources is calculated.
• firesong/Legend.py - Generates an instance of gamma-ray sources in the universe according to a luminosity dependent density evolution (LDDE). Both the redshift and the gamma-ray flux (without attenuation) are calculated.

Examples:

• A muon neutrino diffuse flux saturation example:

If installed via pip, in in the python console,

from firesong.Firesong import firesong_simulation
firesong_simulation('./', density=1e-6, Evolution='CC2015SNR', zmax=4.0,
                    fluxnorm=1.44e-8, index=2.28, LF='SC')

or with the repository downloaded

python Firesong.py -d 1e-6 --evolution CC2015SNR --zmax 4.0
--fluxnorm 1.44e-8 --index 2.28 --LF SC

wlll simulate neutrino sources with a local density of 10^-6 Mpc^-3 with a source density evolution that follows the Clash and Candels 2015 Supernova Rate (CC2015SNR). The simulation will be done up to a redshift of 4.0. The neutrino luminosity, because it is not specified as an option, will be calculated internally to saturate a muon neutrino diffuse flux with a normalization, at 100 TeV, of E^2d\phi/dE = 1.44 x 10^-8 GeV.cm^-2.s^-1.sr^-1 and with a spectral index of -2.28. Neutrino luminosity is distributed as a delta function, i.e., standard candle (SC). The result will be output as a text file firesong.out in the current directoy, or in the directory of environment variable FIRESONG if it is set.

• An exploration of the luminosity vs. local density plane (aka Kowalski plot) example: in the console

firesong_simulation('./', density=1e-6, Evolution='MD2016SFR', zmax=8.0,
                    index=2.28, LF='SC', luminosity=1e51)

or

python Firesong.py -d 1e-6 --evolution MD2016SFR --zmax 8.0
--index 2.28 --LF SC -L 1e51

will simulate neutrino sources with a local density of 10^-6 Mpc^-3 with a source density evolution that follows the Madau and Dickinson 2016 Star Formation Rate History (MD2014SFR). The simulation will be done up to a redshift of 8.0. The neutrino power law spectral index is -2.28. Neutrino luminosity is distributed as a delta function, i.e., standard candle (SC) and is set to be 10^51 erg/year. Note that muon neutrino diffuse flux normalization is ignored when a luminosity is specified, but the spectral index should still be provided. This mode of operation allows the exploration of the luminosity vs. local density plane (aka Kowalski plot).

More examples are included in the notebooks directory. Jupyter notebook and matplotlib are required to run the examples.

Tests

All unittest could be run by

python -m unittest discover tests/

If you would like to suppress the printed output, you may add a -b flag to this command. If you want to run a test for a certain file separately use either

python -m unittest tests/test_<...>

or

python tests/test_<...>.py

How to request support

Questions about support for FIRESONG can be sent to one of the development leads for FIRESONG. See AUTHORS.md

How to contribute

Community contributions to Firesong are accepted and welcome. Issues can be reported by any user, even if not a member of IceCube. Pull requests can be requested by any user, even if not a member of IceCube.

firesong's Issues

Provide option to calculate Luminosity as done currently or as Luminosity per log neutrino energy

"Waxman and Murase (https://arxiv.org/abs/1607.01601) use effective luminosity. It’d be nice to have the option to use one or the other (N.B. no change needed for density vs. effective density)"

Brainstorm on how to improve speed for high density/density rate

I believe there are a few ideas out there right now, including:

• Making an analysis dependent cut on redshift (Ignacio)
• For large densities, return a histogram instead of a list. This will result in a small loss of per-source precision (Alex)

Remove non-FIRESONG code

Remove the CTA directory and contents. Remove CTAsensitivity_LiMa.py. Remove (empty) directories Results/HAWC and Results/CTA. This is because these items are not really the objective of Firesong, but a byproduct. But also, the performance files for CTA are outdated.

Add Ephemerids

A helper tool that randomizes the time, and calculates distance to Sun and Moon would be useful for non-IceCube users. This is only interesting once the response of IceCube has been included.

Installation instructions

[As part of the JOSS review.]

I ran into several issues while trying to install FIRESONG.

• The installation instructions in the README file are currently incomplete: If you recommend installing FIRESONG as an editable install from a local folder (pip install --editable FIRESONG), you need to state explicitly that the user should download the repo first.
  However, note that edible installs are mainly for developers of a package. For users, publishing FIRESONG to PyPI and installing via pip install FIRESONG (i.e. without downloading the repo) may be preferable; both for simplicity and to ensure they install a tagged release (as opposed to the most recent git commit on the base branch, which is more prone to cause confusion).
• Right now, installing FIRESONG does not install any dependencies; so after installing it into an empty conda environment, I get an error message ("ModuleNotFoundError: No module named 'numpy'") when trying to run e.g. from firesong.Firesong import firesong_simulation. You can use setup.py’s install_requires keyword to let pip install these dependencies automatically.
• Since FIRESONG documentation is mainly in the form of Jupyter notebooks, it would be useful to add jupyter and matplotlib as dependencies.

Add utility to relate g-ray / nu fluxes via p-p and p-gamma processes

The utility can be useful to, e.g. produce the high-energy flux of gamma rays associated with neutrinos simulated with Firesong.py (but note that cascading/EBL should not be part of FIRESONG) or to simulate neutrinos from simulated gamma ray sources, such as is done in Legend.py.

This is probably a very simple script

Fix test or revert to prior defaults

Do we fix tests to match the results of the current defaults (SFR evolution MD2014, ICRC 2019 numu diffuse best fit, Planck 2018 cosmological parameters) or revert to prior defaults?

Interaction between fluxnorm and zmin

[As part of the JOSS review.]

Here’s a usage example:

from firesong.Firesong import firesong_simulation

sim = firesong_simulation(None, filename=None,
                                density=1e-10,
                                Evolution='MD2014SFR',
                                zmin=0.0005,  # <-- Default value
                                verbose=False)
print("=== zmin = 0.0005 ===")
print(f"Expected flux: {sim['header']['fluxnorm']}")
print(f"Actual flux of {len(sim['sources']['flux'])} sources: {sim['total_flux']}")

sim = firesong_simulation(None, filename=None,
                                density=1e-10,
                                Evolution='MD2014SFR',
                                zmin=4,  # <-- Note this non-default value!
                                verbose=False)
print("=== zmin = 4 ===")
print(f"Expected flux: {sim['header']['fluxnorm']}")
print(f"Actual flux of {len(sim['sources']['flux'])} sources: {sim['total_flux']}")

And this is the output I got:

=== zmin = 0.0005 ===
Expected flux: 1.44e-08
Actual flux of 1329 sources: 1.528777753275981e-08
=== zmin = 4 ===
Expected flux: 1.44e-08
Actual flux of 1329 sources: 5.704090191660857e-10

The actual flux will vary randomly but is usually between 1e-8 and 2e-8 for zmin=-0.0005 (comparable to the expected flux normalization) while it is between 5.5e-10 and 6e-10 for zmin=4 (i.e. more than 20 times smaller than the expected flux normalization). The number of sources is constant across runs.

I would have expected one of two things:

1. The flux normalization describes the flux from sources in the selected redshift range, between zmin and zmax. (This would imply that the actual flux should always be close to the flux normalization.)
2. The flux normalization describes the flux integrated over all redshifts (starting at 0.0), but the code outputs only the fraction of sources that are between zmin and zmax. (This would imply that selecting a smaller z range should result in fewer sources being generated.)

Neither option appears to be the case. Am I misunderstanding something or is this a bug?

FluxPDF.py apparently not working

But Ignacio may be clueless

command line "-h" does not provide a correct description of luminosity options

This needs to be made consistent with the rewriting of the Luminosity class.

Clean up the units

"erg/year? Yuck! this includes providing more clear units for the case of transient simulation."

Get rid of FIRESONG environmental variable

It's not really needed. Continue having the directory Results as the default output dir.

InverseCDF Evaluates at redshift bin centers

The current InverseCDF function looks like this:

        cdf = np.cumsum(pdf)/np.sum(pdf)
        mask = np.diff(cdf) > 0
        self.invCDF = UnivariateSpline(cdf[1:][mask], x[1:][mask] + (x[1]-x[0])/2., k=1, s=0)

Here, cdf is the value of the redshift evolution PDF evaluated at a redshift of x. The code, however, interprets cdf as the evolution PDF evaluated at the midpoint between given redshifts. For the default binning of 10000, this is likely to be a very small difference, but it can be a noticeable source of error if the binning is reduced. Since FIRESONG isn't doing anything binned at the moment, we should fix these lines.

Add luminosity function for a parametrization of the observed blazar population.

This is a paper for blazars:
https://iopscience.iop.org/article/10.1088/2041-8205/800/2/L27

Warning: The Ajello 2012 and 2014 on this topic had a wrong sign in the function.
It's described in Formulas 1-10

Implement tests

There are various test scripts, but most of them don't do anything. It would be good to actually add tests for test_Firesong at the very least, potentially checking things in the header as well as average fluxes / redshifts.
Similar tests should be implemented for the other modules

support for broken powerlaw luminosity functions

the local luminosity function of x-ray AGN is well approximated by a broken powerlaw
(e.g. eq. 30 of https://arxiv.org/pdf/1503.01120.pdf) with different slopes for the faint and bright ends of the population.

I'd volunteer to contribute an implementation after #70 is resolved.

Rebuilding documentation

The published documentation is not aligned with the docstrings (notably UNITS placeholders are still there).

• documentation in docs needs to be rebuilt and pushed to the repository;
• It could be nice to add instructions to the README on how to build documentation. Installing and running pdoc seems straightforward but LaTeX does not render (what dependencies should be satisfied?);
• can this be automatised with a GitHub Action?

Put the G back in FIRESONG: Add utility/module to calculate g-ray fluxes given nu fluxes

This should be done for pp and p-gamma interactions.
It could be a module that is called as an option itself in FIRESONG, which would output gamma-ray flux or as a standalone program that converts numu fluxes to gamma ray fluxes. Should be easy to code.

documentation and argument names for flux_pdf

[As part of the JOSS review.]

I’m overall impressed by the good documentation across the whole code base; even for internal functions that users of FIRESONG will never interact with directly. Well done!

I just noticed two very minor issues:

• In the docstring for flux_pdf, the description of LumMin/LumMax and logFMin/logFMax says things like "Max luminosity to consider, in UNITS". Could you replace that with the correct units?
• In the docstrings for flux_pdf and firesong_simulation, the argument timescale is described as "timescale (bool, optional, default=1000): Timescale in seconds" — I suspect that should be an int, not a bool?

---

While I’m talking about minor issues regarding the arguments to flux_pdf, I may as well expand on this point that @rafaelab made:

Maybe this is me being picky, but I found it annoying that the keyword "Evolution" in firesong_simulation() is capitalised when all the other keywords are not. I got this wrong twice when trying to run a simple test without copy-pasting. The same goes for "Transient".

There are a few additional inconsistent argument names in flux_pdf:

• zmin/zmax and emin/emax, but LumMin/LumMax
• nFluxBins vs. nLbins ("bins" being uppercase in the first, but lowercase in the second one)
• "flux" is written out in nFluxBins and fluxnorm but abbreviated in logFMin/logFMax

This obviously shouldn’t block acceptance of the JOSS paper, but it would remove a tiny pain point when using the code. (Like Rafael, I too have been tripped up by this once or twice during my testing.)
Of course, this is effectively a backwards-incompatible change to the API; so if you do decide do rename these arguments, I completely understand if you do it slowly and carefully so you don’t break all analysis scripts that use FIRESONG across multiple experiments. 😉

Astropy vs CosmoloPy

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?

problems with $FIRESONG environment variable

[As part of the JOSS review.]

The two code samples in the README, which assume that the user has downloaded the repository, result in a FileNotFoundError if the $FIRESONG environment variable is not set to the main FIRESONG directory:

> pwd
/Users/jost/Documents/FIRESONG/firesong/

> python Firesong.py -d 1e-6 --evolution CC2015SNR --zmax 4.0 --fluxnorm 1.44e-8 --index 2.28 --LF SC
Enviromental variable FIRESONG not set
set to ./
##############################################################################
##### FIRESONG initializing  - Calculating Neutrino CDFs #####
Standard candle sources+
Source evolution assumed: CC2015SNR
Local density of neutrino sources: 1e-06/Mpc^3
Total number of neutrinos sources in the Universe: 5566172
Redshift range: 0 - 4.0
Standard Candle Luminosity: 8.2520e+49 erg/yr
##### FIRESONG initialization done #####

Traceback (most recent call last):
  File "Firesong.py", line 214, in <module>
    firesong_simulation(outputdir,
  File "Firesong.py", line 126, in firesong_simulation
    out = output_writer(outputdir, filename)
  File "/Users/jost/opt/anaconda3/envs/firesong/lib/python3.8/site-packages/firesong/input_output.py", line 30, in __init__
    self.output = self.open_file(outputdir, filename)
  File "/Users/jost/opt/anaconda3/envs/firesong/lib/python3.8/site-packages/firesong/input_output.py", line 49, in open_file
    output = open(outputdir+str(filename), "w")
FileNotFoundError: [Errno 2] No such file or directory: './Results/Firesong.out'

Since the README file implies that setting FIRESONG is optional, this example should also work when it is not set.

On a related note, I found it surprising that the output files are written to ${FIRESONG}/Results/ and not to ${FIRESONG}/. While it’s good practice to have a separate folder for output files, I think it would be more intuitive to point this environment variable directly to that output folder. (Or if there is a good reason to point it to the main FIRESONG folder, at least document clearly in the README that this subfolder name is hardcoded.)

Negative distances for small redshifts

Using UnivariateSpline for the calculation of the luminosity in def LuminosityDistance(self, z): is pretty dangerous. The spline doesn't necessary behave very well for small redshifts and can even produce negative values.

For example

from firesong.Evolution import get_evolution, SourcePopulation from firesong.Evolution import cosmology population = SourcePopulation(cosmology, get_evolution('NoEvolution')) population.LuminosityDistance(np.full(1001, 0.0017507872720130618))

produces array([-1.04664586, -1.04664586, -1.04664586, ..., -1.04664586, -1.04664586, -1.04664586])

a safer alternative would it be to simply use interp1d or another more stable interpolator

Remove NeutrinoAlert.py

It doesn't work in realistic scenarios.
The code is correct when the neutrino sources are 'weak', i.e. a signal expectation of <<1 event in IceCube. But in practice, in interesting simulations, there's typically at least a source (or many!) in which the signal expection is not <<1.
This script is probably confusing. It's better for users to run Firesong.py instead.

Look into options for continuous integration

Would be nice to have automated checks of the tests with every new commit

Remove CTA/HAWC/etc code related and references

Alternatively, we can branch this out.

FluxPDF for non-factorable luminosity functions

The way that FIRESONG is currently implemented this is not super important, but in general, it would be nice if FluxPDF could also deal with non-factorable luminosity functions, i.e., those that do not fulfill p(z, L) = p(z) * p(L). Technically, this simply requires us to sum over all combinations of luminosity and redshift. We did have a working implementation for it at some point, but I fear we have lost it during migration so FIRESONG. Hence, a re-implementation would be necessary.

Make cosmological parameters configurable

Currently Omega_M, Omega_Lambda and h_0 are hardcoded on Evolution.py

This is not high priority for neutrino astrophysics, to be honest. Many theorists simply do 0.3, 0.7 and 0.7

Simulate transient sources, NOT TESTED YET!

In https://github.com/icecube/FIRESONG/blob/base/Firesong.py#L147 the helper string says "NOT TESTED YET!".
We should make sure that the functionality is tested and remove the warning.

Same in https://github.com/icecube/FIRESONG/blob/base/FluxPDF.py#L138

allow to parameterize luminosity functions in "natural parameters" instead of mean_luminosity

The mean_luminosity is not a good property to characterize powerlaws or broken powerlaws. A more natural choice would be Lmin, Lmax, gamma (or Lmin, Lmax, Lbreak, gamma1, gamma2).

The current design puts a lot of emphasis on the mean_luminosity

FIRESONG/firesong/Luminosity.py

Line 232 in 9955d4c

def get_LuminosityFunction(mean_luminosity, LF, **kwargs):

and performs the transformation to the natural parameters implicitly e.g. here (powerlaw)

FIRESONG/firesong/Luminosity.py

Line 182 in 9955d4c

self.Fmin = (self.index+2.)/(self.index+1.)*self.mean

Suggest to rework the interface to the underlying PDFs such that the user can specify the natural parameters directly.

sample code in README file creates "hidden" output files

[As part of the JOSS review.]

When using this sample code from the README file

from firesong.Firesong import firesong_simulation
firesong_simulation('.', density=1e-6, Evolution='CC2015SNR', zmax=4.0,
                    fluxnorm=1.44e-8, index=2.28, LF='SC')

Firesong will generate a file called .Firesong.out (note the period at the start of the file name!), which is hidden by default in most GUI-based file managers. This is potentially very confusing to users.

Instead of simply changing the initial argument in this example (as well as the second example in the README), I would recommend improving the FIRESONG code itself to handle these file paths better (perhaps by using os.path.join() instead of string concatenation). That’s because '.' is such a common shorthand that people would likely use it even if it wasn’t given in an example.

Add option to skip file output

file-I/O can be one of the slowest parts of FIRESONG. Would be nice to have an option to allow the user to choose to just return a dictionary with all of the relevant simulation results instead of writing to a file

Review documentation

"I [Ignacio] will take care of this over the next few months actually, and we can steal material from Chris’ PhD thesis"