tests for emissivity from GalaxyEnsemble

Author: mirochaj

.coveragerc

@@ -3,6 +3,8 @@
 branch = False
 
 [report]
+precision = 2
+
 # Regexes for lines to exclude from consideration
 exclude_lines =
     # Have to re-enable the standard pragma

ares/populations/GalaxyEnsemble.py

@@ -76,6 +76,7 @@ def __dict__(self, name):
     def tab_z(self):
         if not hasattr(self, '_tab_z'):
             h = self._gen_halo_histories()
+            self._tab_z = h['z']
         return self._tab_z
 
     @tab_z.setter
@@ -89,16 +90,6 @@ def tab_t(self):
             self._tab_t = self.cosm.t_of_z(self.tab_z) / s_per_yr
         return self._tab_t
 
-    @property
-    def tab_dz(self):
-        if not hasattr(self, '_tab_dz'):
-            dz = np.diff(self.tab_z)
-            if np.allclose(np.diff(dz), 0):
-                dz = dz[0]
-            self._tab_dz = dz
-
-        return self._tab_dz
-
     @property
     def _b14(self):
         if not hasattr(self, '_b14_'):
@@ -132,7 +123,7 @@ def _roman(self): # pragma: no cover
     def run(self):
         return
 
-    def cSFRD(self, z, Mh):
+    def get_sfrd_in_mass_range(self, z, Mlo, Mhi=None):
         """
         Compute cumulative SFRD as a function of lower-mass bound.
 
@@ -141,26 +132,23 @@ def cSFRD(self, z, Mh):
         Cumulative *FRACTION* of SFRD in halos above Mh.
         """
 
-        if type(Mh) not in [list, np.ndarray]:
-            Mh = np.array([Mh])
-
         iz = np.argmin(np.abs(z - self.histories['z']))
         _Mh = self.histories['Mh'][:,iz]
         _sfr = self.histories['SFR'][:,iz]
         _w = self.histories['nh'][:,iz]
 
-        # Really this is the number of galaxies that formed in a given
-        # differential redshift slice.
-        SFRD = np.zeros_like(Mh)
-        for i, _mass in enumerate(Mh):
-            ok = _Mh >= _mass
-            SFRD[i] = np.sum(_sfr[ok==1] * _w[ok==1]) / rhodot_cgs
+        if Mhi is None:
+            Mhi = _Mh.max()
 
-        SFRD_tot = self.SFRD(z)
+        ok = np.logical_and(_Mh >= Mlo, _Mh <= Mhi)
+        SFRD = np.sum(_sfr[ok==1] * _w[ok==1]) / rhodot_cgs
 
-        return SFRD / SFRD_tot
+        return SFRD
 
     def SFRD(self, z, Mmin=None):
+        return self.get_sfrd(z, Mmin=Mmin)
+
+    def get_sfrd(self, z, Mmin=None):
         """
         Will convert to internal cgs units.
         """
@@ -241,11 +229,11 @@ def tile(self, arr, thin, renorm=False):
         else:
             return new
 
-    def noise_normal(self, arr, sigma):
+    def get_noise_normal(self, arr, sigma):
         noise = np.random.normal(scale=sigma, size=arr.size)
         return np.reshape(noise, arr.shape)
 
-    def noise_lognormal(self, arr, sigma):
+    def get_noise_lognormal(self, arr, sigma):
         lognoise = np.random.normal(scale=sigma, size=arr.size)
         #noise = 10**(np.log10(arr) + np.reshape(lognoise, arr.shape)) - arr
         noise = np.power(10, np.log10(arr) + np.reshape(lognoise, arr.shape)) \
@@ -409,11 +397,11 @@ def _gen_halo_histories(self):
         # Two potential kinds of scatter in MAR
         mar = self.tile(mar_raw, thin)
         if sigma_env > 0:
-            mar *= (1. + self.noise_normal(mar, sigma_env))
+            mar *= (1. + self.get_noise_normal(mar, sigma_env))
 
         if sigma_mar > 0:
             np.random.seed(self.pf['pop_scatter_mar_seed'])
-            noise = self.noise_lognormal(mar, sigma_mar)
+            noise = self.get_noise_lognormal(mar, sigma_mar)
             mar += noise
             # Normalize by mean of log-normal to preserve mean MAR?
             mar /= np.exp(0.5 * sigma_mar**2)
@@ -668,18 +656,23 @@ def _TabulateEmissivity(self, E=None, Emin=None, Emax=None, wave=None):
 
         return zarr, tab / cm_per_mpc**3
 
-    def Emissivity(self, z, E=None, Emin=None, Emax=None):
+    def get_emissivity(self, z, E=None, Emin=None, Emax=None):
         """
         Compute the emissivity of this population as a function of redshift
-        and rest-frame photon energy [eV].
+        and (potentially) rest-frame photon energy [eV].
+
+        .. note :: If Emin and Emax are supplied, this is a luminosity density,
+            and will have units of erg/s/(co-moving cm)^3. If `E` is supplied,
+            will also carry units of eV^-1.
 
         Parameters
         ----------
         z : int, float
+            Redshift.
 
         Returns
         -------
-        Emissivity in units of erg / s / c-cm**3 [/ eV]
+        Emissivity in units of erg / s / c-cm**3 [/ eV].
 
         """
 
@@ -707,9 +700,9 @@ def Emissivity(self, z, E=None, Emin=None, Emax=None):
         return 10**func(z)
         #return self._cache_ehat_[(E, Emin, Emax)](z)
 
-    def PhotonLuminosityDensity(self, z, E=None, Emin=None, Emax=None):
+    def get_photon_density(self, z, E=None, Emin=None, Emax=None):
         # erg / s / cm**3
-        rhoL = self.Emissivity(z, E=E, Emin=Emin, Emax=Emax)
+        rhoL = self.get_emissivity(z, E=E, Emin=Emin, Emax=Emax)
         erg_per_phot = self._get_energy_per_photon(Emin, Emax) * erg_per_ev
 
         return rhoL / np.mean(erg_per_phot)
@@ -1243,7 +1236,7 @@ def _gen_prescribed_galaxy_histories(self, zstop=0):
                 noise = np.zeros_like(Sd)
                 np.random.seed(self.pf['pop_dust_scatter_seed'])
                 for _i, _z in enumerate(z):
-                    noise[:,_i] = self.noise_lognormal(Sd[:,_i], sigma[:,_i])
+                    noise[:,_i] = self.get_noise_lognormal(Sd[:,_i], sigma[:,_i])
 
                 Sd += noise
 
@@ -1578,32 +1571,6 @@ def _gen_active_galaxy_histories(self):
 
         return results
 
-
-    def Slice(self, z, slc):
-        """
-        slice format = {'field': (lo, hi)}
-        """
-
-        iz = np.argmin(np.abs(z - self.tab_z))
-        hist = self.histories
-
-        c = np.ones(hist['Mh'].shape[0], dtype=int)
-        for key in slc:
-            lo, hi = slc[key]
-
-            ok = np.logical_and(hist[key][:,iz] >= lo, hist[key][:,iz] <= hi)
-            c *= ok
-
-        # Build output
-        to_return = {}
-        for key in self.histories:
-            if self.histories[key].ndim == 1:
-                to_return[key] = self.histories[key][c==1]
-            else:
-                to_return[key] = self.histories[key][c==1,iz]
-
-        return to_return
-
     def get_field(self, z, field):
         iz = np.argmin(np.abs(z - self.histories['z']))
         return self.histories[field][:,iz]

tests/test_populations_ensemble.py

@@ -12,7 +12,7 @@
 
 import ares
 import numpy as np
-from ares.physics.Constants import rhodot_cgs
+from ares.physics.Constants import rhodot_cgs, E_LL, cm_per_mpc, ev_per_hz
 
 def test():
     pars = ares.util.ParameterBundle('mirocha2020:univ')
@@ -26,14 +26,17 @@ def test():
     # Test I/O. Should add more here eventually.
     pop.save('test_ensemble', clobber=True)
 
+    z = pop.tab_z
+    t = pop.tab_t
+
     # Test SFRD
-    sfrd = pop.SFRD(6.) * rhodot_cgs
+    sfrd = pop.get_sfrd(6.) * rhodot_cgs
 
     assert 1e-3 <= sfrd <= 1, "SFRD unreasonable"
 
-    csfrd = pop.cSFRD(6., Mh=1e10)
+    sfrd10 = pop.get_sfrd_in_mass_range(6., Mlo=1e10)
 
-    assert csfrd < 1., "cSFRD >= SFRD!?"
+    assert sfrd10 < sfrd, "sfrd(Mh>1e10) >= sfrd(all Mh)!?"
 
     logMh, logf_stell, std = pop.get_smhm(6.)
     assert np.all(logf_stell < 1)
@@ -144,5 +147,20 @@ def test():
     # Test routines to retrieve MUV-Beta, AUV, etc.
     AUV = pop.get_AUV(6.)
 
+    # Emissivity stuff: just OOM check at the moment.
+    zarr = np.arange(6, 30)
+    e_ion = np.array([pop.get_emissivity(z, Emin=E_LL, Emax=1e2) \
+        for z in zarr]) * cm_per_mpc**3
+    e_ion2 = np.array([pop.get_emissivity(z, Emin=E_LL, Emax=1e2) \
+        for z in zarr]) * cm_per_mpc**3     # check caching
+
+    assert 1e37 <= np.mean(e_ion) <= 1e41
+    assert np.allclose(e_ion, e_ion2)
+
+    n_ion = np.array([pop.get_photon_density(z, Emin=E_LL, Emax=1e2) \
+        for z in zarr]) * cm_per_mpc**3
+
+    assert 1e47 <= np.mean(n_ion) <= 1e51
+
 if __name__ == '__main__':
     test()