Merge pull request #45 from mirochaj/docs

documentation overhaul

README.md

@@ -10,9 +10,10 @@ A few papers on how it works:
 
 - 1-D radiative transfer: [Mirocha et al. (2012)](http://adsabs.harvard.edu/abs/2012ApJ...756...94M)
 - Uniform backgrounds \& global 21-cm signal: [Mirocha (2014)](http://adsabs.harvard.edu/abs/2014MNRAS.443.1211M)
-- Galaxy luminosity functions: [Mirocha, Furlanetto, & Sun (2017)](http://adsabs.harvard.edu/abs/2016arXiv160700386M)
+- Galaxy luminosity functions: [Mirocha, Furlanetto, & Sun (2017)](https://ui.adsabs.harvard.edu/abs/2017MNRAS.464.1365M/abstract)
 - Population III star formation: [Mirocha et al. (2018)](http://adsabs.harvard.edu/abs/2018MNRAS.478.5591M)
 - Rest-ultraviolet colours at high-z: [Mirocha, Mason, & Stark (2020)](https://ui.adsabs.harvard.edu/abs/2020arXiv200507208M/abstract)
+- Near-infrared background and nebular emission: [Sun et al. (2021)](https://ui.adsabs.harvard.edu/abs/2021MNRAS.508.1954S/abstract)
 
 Plus some more applications:
 
@@ -26,15 +27,22 @@ risk! Correctness of results is not guaranteed.
 
 ## Citation
 
-If you use ARES in paper please reference [Mirocha (2014)](http://adsabs.harvard.edu/abs/2014MNRAS.443.1211M) if it's an application of the global 21-cm modeling machinery and [Mirocha et al. (2012)](http://adsabs.harvard.edu/abs/2012ApJ...756...94M) if you use the 1-D radiative transfer. Note that for some applications, ARES relies heavily on lookup tables and publicly-available software packages that should be referenced as well. These include:
-
-- [CAMB](https://camb.readthedocs.io/en/latest/)
-- [hmf](https://hmf.readthedocs.io/en/latest/)
-- Collisional coupling coefficients from [Zygelman (2005)](https://ui.adsabs.harvard.edu/abs/2005ApJ...622.1356Z/abstract).
-- Lyman-alpha transition probabilities from [Pritchard \& Furlanetto](https://ui.adsabs.harvard.edu/abs/2006MNRAS.367.1057P/abstract).
+If you use ARES in paper please reference [Mirocha (2014)](http://adsabs.harvard.edu/abs/2014MNRAS.443.1211M) if it's an application of the global 21-cm modeling machinery and [Mirocha et al. (2012)](http://adsabs.harvard.edu/abs/2012ApJ...756...94M) if you use the 1-D radiative transfer and/or SED optimization. For galaxy semi-analytic modeling, please have a look at [Mirocha, Furlanetto, & Sun (2017)](http://adsabs.harvard.edu/abs/2016arXiv160700386M), [Mirocha, Mason, & Stark (2020)](https://ui.adsabs.harvard.edu/abs/2020arXiv200507208M/abstract), and [Mirocha (2020)](https://ui.adsabs.harvard.edu/abs/2020MNRAS.499.4534M/abstract), and for PopIII star modeling, see [Mirocha et al. (2018)](http://adsabs.harvard.edu/abs/2018MNRAS.478.5591M).
 
 Please also provide a link to [this page](https://github.com/mirochaj/ares) as a footnote.
 
+Note that for some applications, ARES relies heavily on lookup tables and publicly-available software packages that should be referenced as well. These include:
+
+- Code for Anisotropies in the Microwave Background ([CAMB](https://camb.readthedocs.io/en/latest/)).
+- The Halo Mass Function ([hmf](https://hmf.readthedocs.io/en/latest/)) package (see [Murray et al.(2013)](https://arxiv.org/abs/1306.6721)).
+- Lookup tables and fitting formulae for the fraction of photo-electron energy deposited in heat, ionization, excitation from [Shull \& van Steenberg (1985)](https://ui.adsabs.harvard.edu/abs/1985ApJ...298..268S/abstract), [Ricotti, Gnedin, \& Shull (2002)](https://ui.adsabs.harvard.edu/abs/2002ApJ...575...33R/abstract), and [Furlanetto \& Stoever (2010)](https://ui.adsabs.harvard.edu/abs/2010MNRAS.404.1869F/abstract) (see `secondary_ionization` parameter, values of 2, 3, and 4, respectively).
+- Collisional coupling coefficients for the 21-cm line from [Zygelman (2005)](https://ui.adsabs.harvard.edu/abs/2005ApJ...622.1356Z/abstract).
+- Wouthuysen-Field coupling coefficients for the 21-cm line from [Chuzhoy, Alvarez, & Shapiro (2006)](https://ui.adsabs.harvard.edu/abs/2006ApJ...651....1C/abstract), [Furlanetto \& Pritchard (2006)](https://ui.adsabs.harvard.edu/abs/2006MNRAS.372.1093F/abstract), and [Hirata (2006)](https://ui.adsabs.harvard.edu/abs/2006MNRAS.367..259H/abstract) (see `approx_Salpha` parameter, values of 2, 3, and 4 respectively).
+- Lyman-alpha transition probabilities from [Pritchard \& Furlanetto (2006)](https://ui.adsabs.harvard.edu/abs/2006MNRAS.367.1057P/abstract).
+- Stellar population synthesis model options include starburst99 ([Leitherer et al. (1999)](https://ui.adsabs.harvard.edu/abs/1999ApJS..123....3L/abstract)) and BPASS versions 1 ([Eldridge \& Stanway (2009)](https://ui.adsabs.harvard.edu/abs/2009MNRAS.400.1019E/abstract)) and 2 ([Eldridge et al. (2017)](https://ui.adsabs.harvard.edu/abs/2017PASA...34...58E/abstract),[Stanway \& Eldridge (2018)](https://ui.adsabs.harvard.edu/abs/2018MNRAS.479...75S/abstract)) (via `pop_sed` parameter, values `'starburst99'`, `'bpass_v1'`, and `'bpass_v2'`, respectively).
+
+Feel free to get in touch if you are unsure of whether any of these tools are being used under the hood for your application.
+
 ## Dependencies
 
 You will need:
@@ -46,6 +54,7 @@ You will need:
 
 and optionally,
 
+- [camb](https://camb.readthedocs.io/en/latest/)
 - [hmf](https://github.com/steven-murray/hmf)
 - [mpi4py](http://mpi4py.scipy.org)
 - [pymp](https://github.com/classner/pymp)
@@ -57,6 +66,10 @@ and optionally,
 - [shapely](https://pypi.python.org/pypi/Shapely)
 - [descartes](https://pypi.python.org/pypi/descartes)
 
+If you'd like to build the documentation locally, you'll need:
+
+- [numpydoc](https://numpydoc.readthedocs.io/en/latest/)
+- [nbsphinx](https://nbsphinx.readthedocs.io/en/0.8.8/)
 
 Note: **ares** has been tested only with Python 2.7.x and Python 3.7.x.
 
@@ -70,7 +83,7 @@ cd ares
 python setup.py install
 ```
 
-**ares** will look in ``$ARES/input`` for lookup tables of various kinds. To download said lookup tables, run:
+**ares** will look in ``ares/input`` for lookup tables of various kinds. To download said lookup tables, run:
 
 ```
 python remote.py
@@ -143,4 +156,6 @@ Additional contributions / corrections / suggestions from:
 - Matteo Leo
 - Emma Klemets
 - Felix Bilodeau-Chagnon
+- Venno Vipp
+- Oscar Hernandez
 - Joshua Hibbard

ares/analysis/GalaxyPopulation.py

@@ -562,8 +562,8 @@ def PlotColors(self, pop, axes=None, fig=1, z_uvlf=[4,6,8,10],
         ##
         # Plot models
         ##
-        Ms = np.arange(6, 13.+dMst, dMst)
-        mags = np.arange(-25, -12-dmag, dmag)
+        Ms = np.arange(2, 13.+dMst, dMst)
+        mags = np.arange(-30, 10-dmag, dmag)
         mags_cr = np.arange(-25, -10, dmag)
         hst_shallow = b14.filt_shallow
         hst_deep = b14.filt_deep
@@ -603,7 +603,7 @@ def PlotColors(self, pop, axes=None, fig=1, z_uvlf=[4,6,8,10],
                         _ax = ax_uvlf
 
                     if show_MUV:
-                        phi = pop.LuminosityFunction(z, mags)
+                        _mags, phi = pop.get_lf(z, mags)
                         uvlf_by_pop[h][z] = phi
 
                         if not fill:
@@ -617,7 +617,7 @@ def PlotColors(self, pop, axes=None, fig=1, z_uvlf=[4,6,8,10],
                         else:
                             _ax2 = ax_smf
 
-                        phi = pop.StellarMassFunction(z, bins=Ms)
+                        _bins, phi = pop.StellarMassFunction(z, bins=Ms)
                         smf_by_pop[h][z] = phi
 
                         if not fill:
@@ -1754,7 +1754,7 @@ def _MegaPlotPop(self, kw, pop, redshifts=None, include_colors=True,
         for j, z in enumerate(redshifts):
 
             # UVLF
-            phi = pop.LuminosityFunction(z, _mags)
+            _mags_, phi = pop.LuminosityFunction(z, _mags)
             ax_phi.semilogy(_mags, phi, color=colors[j], drawstyle='steps-mid')
 
             # Binned version

ares/physics/Cosmology.py

@@ -31,7 +31,7 @@
  'primordial_index': 'ns',
 }
 
-class Cosmology(InitialConditions):
+class Cosmology(object):
     def __init__(self, pf=None, **kwargs):
         if pf is not None:
             self.pf = pf
@@ -111,23 +111,24 @@ def __init__(self, pf=None, **kwargs):
         # 'sigma_8':self.sigma8,
         # 'n': self.primordial_index}
 
+    @property
+    def _ics(self):
+        if not hasattr(self, '_ics_'):
+            self._ics_ = InitialConditions(pf=self.pf)
+            self._ics_.prefix = self.get_prefix()
+        return self._ics_
+
     def nH(self, z):
         return self.nH0 * (1. + z)**3
 
     def nHe(self, z):
         return self.nHe0 * (1. + z)**3
 
-    @property
-    def path_ARES(self):
-        if not hasattr(self, '_path_ARES'):
-            self._path_ARES = ARES
-        return self._path_ARES
-
     @property
     def path_Planck(self):
         if not hasattr(self, '_path_Planck'):
             name = self.pf['cosmology_name'].replace('planck_', '')
-            self._path_Planck = self.path_ARES \
+            self._path_Planck = ARES \
                  + '/input/planck/base/plikHM_{}'.format(name)
         return self._path_Planck
 
@@ -250,7 +251,7 @@ def _load_planck(self):
                 s = "# Set cosmological parameters to values in {}th element of".format(num)
                 s += " concatenated array made from the following files:"
                 print(s)
-                path_str = path.replace(self.path_ARES, '$ARES')
+                path_str = path.replace(ARES, '$ARES')
                 print("# {}_{}_?.txt".format(path_str, prefix))
 
         return
@@ -276,7 +277,7 @@ def get_prefix(self):
     @property
     def inits(self):
         if not hasattr(self, '_inits'):
-            self._inits = self.get_inits_rec()
+            self._inits = self._ics.get_inits_rec()
         return self._inits
 
     def TimeToRedshiftConverter(self, t_i, t_f, z_i):

ares/physics/HaloMassFunction.py

@@ -538,7 +538,7 @@ def _load_hmf(self):
 
         if self.pf['verbose'] and rank == 0:
             name = self.tab_name
-            print("# Loaded {}.".format(name.replace(self.cosm.path_ARES, '$ARES')))
+            print("# Loaded {}.".format(name.replace(ARES, '$ARES')))
 
         if self.pf['hmf_func'] is not None:
             if self.pf['verbose']:

ares/physics/Hydrogen.py

@@ -585,7 +585,7 @@ def DifferentialBrightnessTemperature(self, z, xavg, Ts, Tr=0.0):
     @property
     def inits(self):
         if not hasattr(self, '_inits'):
-            self._inits = self.cosm.get_inits_rec()
+            self._inits = self.cosm._ics.get_inits_rec()
         return self._inits
 
     def saturated_limit(self, z):

ares/physics/InitialConditions.py

@@ -15,6 +15,7 @@
 import os
 import re
 import numpy as np
+from ..data import ARES
 
 try:
     import camb
@@ -34,19 +35,29 @@ class InitialConditions(object):
     """
     This should be inherited by Cosmology.
     """
+    def __init__(self, pf):
+        self.pf = pf
+
+    @property
+    def prefix(self):
+        return self._prefix
+
+    @prefix.setter
+    def prefix(self, value):
+        self._prefix = value
+
     def get_inits_rec(self):
         """
         Get recombination history from file or directly from CosmoRec.
         """
-        fn = '{}/input/inits/inits_{}.txt'.format(self.path_ARES,
-            self.get_prefix())
+        fn = '{}/input/inits/inits_{}.txt'.format(ARES, self.prefix)
 
         # Look for table first, then run if we don't find it.
         if os.path.exists(fn):
             z, xe, Tk = np.loadtxt(fn, unpack=True)
             if self.pf['verbose']:
                 name = fn
-                print("# Loaded {}.".format(fn.replace(self.path_ARES, '$ARES')))
+                print("# Loaded {}.".format(fn.replace(ARES, '$ARES')))
             return {'z': z, 'xe': xe, 'Tk': Tk}
         else:
             if self.pf['verbose']:
@@ -86,9 +97,9 @@ def _run_CosmoRec(self, save=True): # pragma: no cover
         CR_pars = [self.pf[par] for par in _pars_CosmoRec]
 
         # Correct output dir. Just add provided path on top of $ARES
-        CR_pars[-2] = '{}/{}/'.format(self.path_ARES, CR_pars[-2])
+        CR_pars[-2] = '{}/{}/'.format(ARES, CR_pars[-2])
 
-        fn_pars = 'cosmorec_{}.dat'.format(self.get_prefix())
+        fn_pars = 'cosmorec_{}.dat'.format(self.prefix)
 
         # Create parameter file for reference
         to_outputs = CR_pars[-2]
@@ -119,8 +130,8 @@ def _run_CosmoRec(self, save=True): # pragma: no cover
          'Tk': data[:,2][-1::-1],
         }
 
-        fn_out = '{}/input/inits/inits_{}.txt'.format(self.path_ARES,
-            self.get_prefix())
+        fn_out = '{}/input/inits/inits_{}.txt'.format(ARES,
+            self.prefix)
 
         np.savetxt(fn_out, data[-1::-1,0:3], header='z; xe; Te')
 

ares/populations/GalaxyEnsemble.py

@@ -2030,7 +2030,7 @@ def get_mags(self, z, MUV=None, wave=1600., cam=None, filters=None,
         # Either load previous result or compute from scratch
         fil = filters
         if cached_result is not None:
-            M, mags = cached_result
+            M, mags, xph = cached_result
         else:
             # Take monochromatic (or within some window) MUV
             L = self.get_lum(z, wave=wave, window=window, load=load)
@@ -2064,10 +2064,11 @@ def get_mags(self, z, MUV=None, wave=1600., cam=None, filters=None,
 
                 mags = np.array(mags)
             else:
+                xph = None
                 mags = M + magcorr
 
             if hasattr(self, '_cache_mags_'):
-                self._cache_mags_[kw_tup] = M, mags
+                self._cache_mags_[kw_tup] = M, mags, xph
 
         ##
         # Interpolate etc.
@@ -3601,7 +3602,7 @@ def load(self):
 
                 f.close()
                 if self.pf['verbose']:
-                    print("# Loaded {}.".format(fn_hist.replace(self.cosm.path_ARES, '$ARES')))
+                    print("# Loaded {}.".format(fn_hist.replace(ARES, '$ARES')))
 
             else:
                 # Assume pickle?
@@ -3611,7 +3612,7 @@ def load(self):
                 f.close()
                 if self.pf['verbose']:
                     name = fn_hist + '.pkl'
-                    print("# Loaded {}.".format(name.replace(self.cosm.path_ARES, '$ARES')))
+                    print("# Loaded {}.".format(name.replace(ARES, '$ARES')))
 
                 if self.pf['verbose']:
                     print("# Read `pop_histories` as dictionary")

ares/populations/Population.py

@@ -420,7 +420,8 @@ def is_emissivity_scalable(self):
 
             if (self.affects_cgm) and (not self.affects_igm):
                 if self.pf['pop_fesc'] != self.pf['pop_fesc_LW']:
-                    print("# WARNING: revisit scalability wrt fesc.")
+                    if self.pf['verbose']:
+                        print("# WARNING: revisit scalability wrt fesc.")
                     #print("Not scalable cuz fesc pop={}".format(self.id_num))
             #        self._is_emissivity_scalable = False
             #        return False

ares/sources/SynthesisModel.py

@@ -707,7 +707,7 @@ def data(self):
                         self._litinst._load(**self.pf)
 
                     if self.pf['verbose']:
-                        print("# Loaded {}".format(_fn.replace(self.cosm.path_ARES,
+                        print("# Loaded {}".format(_fn.replace(ARES,
                             '$ARES')))
             else:
                 if self.pf['source_sed_by_Z'] is not None:
@@ -720,7 +720,7 @@ def data(self):
 
                     if self.pf['verbose']:
                         for _fn_ in _fn:
-                            print("# Loaded {}".format(_fn_.replace(self.cosm.path_ARES, '$ARES')))
+                            print("# Loaded {}".format(_fn_.replace(ARES, '$ARES')))
 
                 # Shape is (Z, wavelength, time)?
                 to_interp = np.array(_tmp)

ares/util/ParameterBundles.py

@@ -583,18 +583,19 @@ def __add__(self, other):
 
                 if self.verbose:
                     if type(other[key]) in [int, float, np.float64]:
-                        msg1 = "UPDATE: Setting {0} -> {1:9.2e}".format(key.ljust(20),
+                        msg1 = "UPDATE: Setting {0} -> {1:9.2e}".format(key,
                             other[key])
                     else:
-                        msg1 = "UPDATE: Setting {0} -> {1}".format(key.ljust(20),
+                        msg1 = "UPDATE: Setting {0} -> {1}".format(key,
                             str(other[key]).ljust(12))
                     if type(tmp[key]) in [int, float, np.float64]:
-                        msg2 = "previously {0} = {1:9.2e}".format(str(key).ljust(20),
+                        msg2 = "       previously {0} = {1:9.2e}".format(str(key),
                             tmp[key])
                     else:
-                        msg2 = "previously {0} = {1}".format(str(key).ljust(20), tmp[key])
+                        msg2 = "       previously {0} = {1}".format(str(key), tmp[key])
 
-                    print(line('{0} [{1}]'.format(msg1, msg2)))
+                    print(line('{0}'.format(msg1)))
+                    print(line('[{0}]'.format(msg2)))
 
             tmp[key] = other[key]