Merge pull request #40 from dangilman/galacticus

Galacticus

pyHalo/Halos/HaloModels/NFW.py

@@ -46,14 +46,11 @@ def lenstronomy_params(self):
 
         (concentration) = self.profile_args
         Rs_angle, theta_Rs = self._lens_cosmo.nfw_physical2angle(self.mass, concentration, self.z)
-
         x, y = np.round(self.x, 4), np.round(self.y, 4)
         Rs_angle = np.round(Rs_angle, 10)
         theta_Rs = np.round(theta_Rs, 10)
-
         kwargs = [{'alpha_Rs': self._rescale_norm * theta_Rs, 'Rs': Rs_angle,
                   'center_x': x, 'center_y': y}]
-
         return kwargs, None
 
     @property

pyHalo/Halos/HaloModels/TNFW.py

@@ -1,6 +1,7 @@
 from pyHalo.Halos.halo_base import Halo
 from lenstronomy.LensModel.Profiles.tnfw import TNFW as TNFWLenstronomy
 import numpy as np
+from pyHalo.Halos.util import tnfw_mass_fraction
 
 class TNFWFieldHalo(Halo):
 
@@ -94,23 +95,6 @@ def profile_args(self):
 
         return self._profile_args
 
-    @property
-    def bound_mass(self):
-        """
-        Computes the mass inside the virial radius (with truncation effects included)
-        :return: the mass inside r = c * r_s
-        """
-        prof = TNFWLenstronomy()
-        kwargs_profile = self.lenstronomy_params[0][0]
-        alpha_rs = kwargs_profile['alpha_Rs']
-        rs = kwargs_profile['Rs']
-        r_trunc = kwargs_profile['r_trunc']
-        r = self.c * rs
-        rho0 = prof.alpha2rho0(alpha_rs, rs)
-        mass_3d = prof.mass_3d(r, rs, rho0, r_trunc)
-        mass_3d_infall = prof.mass_3d(r, rs, rho0, 1000*rs)
-        return (mass_3d / mass_3d_infall) * self.mass
-
 
 class TNFWSubhalo(TNFWFieldHalo):
     """
@@ -142,3 +126,17 @@ def params_physical(self):
             self._params_physical = {'rhos': rhos, 'rs': rs, 'r200': r200, 'r_trunc_kpc': rt}
 
         return self._params_physical
+
+    @property
+    def bound_mass(self):
+        """
+        Computes the mass inside the virial radius (with truncation effects included)
+        :return: the mass inside r = c * r_s
+        """
+        if hasattr(self, '_kwargs_lenstronomy'):
+            tau = self._kwargs_lenstronomy[0]['r_trunc'] / self._kwargs_lenstronomy[0]['Rs']
+        else:
+            params_physical = self.params_physical
+            tau = params_physical['r_trunc_kpc'] / params_physical['rs']
+        f = tnfw_mass_fraction(tau, self.c)
+        return f * self.mass

pyHalo/Halos/galacticus_truncation/interp_mass_loss.py

@@ -0,0 +1,82 @@
+import numpy as np
+from scipy.interpolate import RegularGridInterpolator
+from pyHalo.Halos.galacticus_truncation.tabulated_mass_loss import _log10_mbound_over_minfall, \
+    _log10_mbound_over_minfall_scatter_dex
+
+class InterpGalacticus(object):
+    """
+    This class interpolates output from the semi-analytic model galacticus to predict the bound mass of a subhalo
+    as a function of its infall mass, concentration, host halo concentration, and the time since infall.
+    """
+    def __init__(self):
+        _t_ellapsed_min, _t_ellapsed_max = 0.0, 7.658449372663094
+        _chostmin, _chostmax = 2.0, 8.0
+        _log10cmin, _log10cmax = 0.3010299956639812, 2.3010299956639813
+        _n = 20
+        _t_coords = np.linspace(_t_ellapsed_min, _t_ellapsed_max, _n)
+        _chost_coords = np.linspace(_chostmin, _chostmax, _n)
+        _log10c_coords = np.linspace(_log10cmin, _log10cmax, _n)
+        _values = np.array(_log10_mbound_over_minfall).reshape(_n, _n, _n)
+        _values_scatter = np.array(_log10_mbound_over_minfall_scatter_dex).reshape(_n, _n, _n)
+        _points = (_t_coords, _chost_coords, _log10c_coords)
+        self._mfrac_interp = RegularGridInterpolator(_points, _values, bounds_error=False, fill_value=None)
+        self._mfrac_scatter_dex_interp = RegularGridInterpolator(_points, _values_scatter, bounds_error=False,
+                                                                 fill_value=None)
+
+    def evaluate_mean_mass_loss(self, log10_concentration_infall, time_since_infall, host_concentration):
+        """
+
+        :param log10_concentration_infall: log10(c) where c is the halo concentration at infall
+        :param time_since_infall: the time ellapsed since infall and the deflector redshift
+        :param host_concentration: the concentration of the host halo
+        :return: bound mass divided by the infall mass
+        """
+        point = (time_since_infall, host_concentration, log10_concentration_infall)
+        y = self._mfrac_interp(point)
+        if isinstance(log10_concentration_infall, float) and \
+            isinstance(time_since_infall, float) and \
+            isinstance(host_concentration, float):
+            return float(y)
+        else:
+            return np.squeeze(y)
+
+    def evaluate_scatter_dex(self, log10_concentration_infall, time_since_infall, host_concentration):
+        """
+
+        :param log10_concentration_infall: log10(c) where c is the halo concentration at infall
+        :param time_since_infall: the time ellapsed since infall and the deflector redshift
+        :param host_concentration: the concentration of the host halo
+        :return: the scatter in dex of the bound mass divided by infall mass
+        """
+        point = (time_since_infall, host_concentration, log10_concentration_infall)
+        y = self._mfrac_scatter_dex_interp(point)
+        if isinstance(log10_concentration_infall, float) and \
+            isinstance(time_since_infall, float) and \
+            isinstance(host_concentration, float):
+            return float(y)
+        else:
+            return np.squeeze(y)
+
+    def __call__(self, log10_concentration_infall, time_since_infall, host_concentration):
+        """
+        Evaluates the prediction from galacticus for subhalo bound mass
+        :param log10_concentration_infall: log10(c) where c is the halo concentration at infall
+        :param time_since_infall: the time ellapsed since infall and the deflector redshift
+        :param host_concentration: the concentration of the host halo
+        :return: the log10(bound mass divided by the infall mass), plus scatter
+        """
+        mean = self.evaluate_mean_mass_loss(log10_concentration_infall, time_since_infall, host_concentration)
+        scatter_dex = self.evaluate_scatter_dex(log10_concentration_infall, time_since_infall, host_concentration)
+        output = np.random.normal(mean, scatter_dex)
+
+        if isinstance(log10_concentration_infall, float) and \
+            isinstance(time_since_infall, float) and \
+            isinstance(host_concentration, float):
+            output = min(0.0, max(-4.0, output))
+        else:
+            inds_high = np.where(output > 0.0)
+            output[inds_high] = 0.0
+            inds_low = np.where(output < -4.0)
+            output[inds_low] = -3.0
+        return output
+

pyHalo/Halos/lens_cosmo.py

@@ -7,10 +7,18 @@
 from colossus.lss.bias import twoHaloTerm
 from scipy.integrate import quad
 
+
 class LensCosmo(object):
 
     def __init__(self, z_lens=None, z_source=None, cosmology=None):
+        """
 
+        This class performs calcuations relevant for certain halo mass profiles and lensing-related quantities for a
+        given lens/source redshift and cosmology
+        :param z_lens: deflector redshift
+        :param z_source: source redshift
+        :param cosmology: and instance of the Cosmology class (see pyhalo.Cosmology.cosmology.py)
+        """
         if cosmology is None:
             from pyHalo.Cosmology.cosmology import Cosmology
             cosmology = Cosmology()
@@ -74,13 +82,14 @@ def twohaloterm(self, r, M, z, mdef='200c'):
         rho_2h = twoHaloTerm(r_h, M_h, z, mdef=mdef) / self.cosmo._colossus_cosmo.rho_m(z)
         return rho_2h
 
-    def nfw_physical2angle(self, m, c, z):
+    def nfw_physical2angle(self, m, c, z, no_interp=False):
         """
         converts the physical mass and concentration parameter of an NFW profile into the lensing quantities
         updates lenstronomy implementation with arbitrary redshift
 
         :param m: mass enclosed 200 rho_crit in units of M_sun (physical units, meaning no little h)
         :param c: NFW concentration parameter (r200/r_s)
+        :param no_interp: bool; compute NFW params with interpolation
         :return: Rs_angle (angle at scale radius) (in units of arcsec), alpha_Rs (observed bending angle at the scale radius
         """
         dd = self.cosmo.D_A_z(z)

pyHalo/Halos/tidal_truncation.py

@@ -1,11 +1,9 @@
 import numpy as np
-import pickle
 from pyHalo.Halos.concentration import ConcentrationDiemerJoyce
-from scipy.interpolate import RegularGridInterpolator
 from pyHalo.Halos.util import tau_mf_interpolation
 from colossus.lss import peaks
 from colossus.halo import splashback
-
+from pyHalo.Halos.galacticus_truncation.interp_mass_loss import InterpGalacticus
 
 class TruncationSplashBack(object):
 
@@ -122,7 +120,7 @@ class AdiabaticTidesTruncation(object):
     def __init__(self, lens_cosmo, log_m_host, z_host, mass_loss_interp):
         """
 
-        :param lens_cosmo: an instacee of the LensCosmo class
+        :param lens_cosmo: an instance of the LensCosmo class
         :param log_m_host: the host halo mass
         :param z_host: the redshift of the host halo
         :param mass_loss_interp: an instance of RegularGridInterpolator
@@ -194,8 +192,8 @@ def _make_params_in_bounds_tau_evaluate(self, point):
         This routine makes sure the arguments for the initerpolation are inside the domain of the function.
         """
         (log10c, log10mass_loss_fraction) = point
-        log10c = max(self._min_c, log10c)
-        log10c = min(self._max_c, log10c)
+        log10c = max(np.log10(self._min_c), log10c)
+        log10c = min(np.log10(self._max_c), log10c)
         log10mass_loss_fraction = max(-1.5, log10mass_loss_fraction)
         log10mass_loss_fraction = min(-0.01, log10mass_loss_fraction)
         return (log10c, log10mass_loss_fraction)
@@ -241,3 +239,62 @@ def truncation_radius(self, halo_mass, halo_redshift, c_median, c_actual, r_peri
         rt_over_rs = self._norm * (c_actual / c_median) ** self._cpower * (r_peri / 0.5)
         _, rs, _ = self.lens_cosmo.NFW_params_physical(halo_mass, c_actual, halo_redshift)
         return rs * rt_over_rs
+
+class TruncationGalacticus(object):
+
+    def __init__(self, lens_cosmo, c_host):
+        """
+
+        :param lens_cosmo:
+        """
+
+        self._chost = c_host
+        self._lens_cosmo = lens_cosmo
+        self._mass_loss_interp = InterpGalacticus()
+        self._tau_mf_interpolation = tau_mf_interpolation()
+
+    @staticmethod
+    def _make_params_in_bounds_tau_evaluate(point):
+        """
+        This routine makes sure the arguments for the interpolation are inside the domain of the function.
+        """
+        min_c, max_c = np.log10(1.0), np.log10(10 ** 2.7)
+        (log10c, log10mass_loss_fraction) = point
+        log10c = max(min_c, log10c)
+        log10c = min(max_c, log10c)
+        log10mass_loss_fraction = max(-3.0, log10mass_loss_fraction)
+        log10mass_loss_fraction = min(-0.001, log10mass_loss_fraction)
+        return (log10c, log10mass_loss_fraction)
+
+    def truncation_radius_halo(self, halo):
+
+        """
+        Thiis method computess the truncation radius using the class attributes of an instance of Halo
+        :param halo: an instance of halo
+        :return: the truncation radius in physical kpc
+        """
+
+        return self.truncation_radius(halo.mass, halo.c,
+                                      halo.time_since_infall, self._chost, halo.z)
+
+    def truncation_radius(self, halo_mass, infall_concentration,
+                          time_since_infall, chost, z_eval_angles):
+        """
+
+        :param halo_mass:
+        :param infall_concentration:
+        :param time_since_infall:
+        :param chost:
+        :param z_eval_angles:
+        :return:
+        """
+        log10c = np.log10(infall_concentration)
+        log10mbound_over_minfall = self._mass_loss_interp(log10c, time_since_infall, chost)
+        point = self._make_params_in_bounds_tau_evaluate((log10c, log10mbound_over_minfall))
+        log10tau = float(self._tau_mf_interpolation(point))
+        tau = 10 ** log10tau
+        _, rs, _ = self._lens_cosmo.NFW_params_physical(halo_mass,
+                                                        infall_concentration,
+                                                        z_eval_angles)
+        return tau * rs
+

pyHalo/Halos/util.py

@@ -40,11 +40,11 @@ def tau_mf_interpolation():
     and final bound mass
     :return: an instance of RegularGridInterpolator that returns (r_t / r_s) given a final mass and concentration
     """
-    N = 50
-    tau = np.logspace(-1., 2., N)
+    N = 80
+    tau = np.logspace(-2., 2.3, N)
     log10_c = np.linspace(0, 2.7, N)
     # mass_fraction_1d = np.logspace(-1.45, -0.02, N)
-    mass_fraction_1d = np.logspace(-1.5, np.log10(0.999), N)
+    mass_fraction_1d = np.logspace(-3.0, np.log10(0.9999), N)
     log10tau_2d = np.zeros((N, N))
 
     # This computes the value of tau that correponds to each pair of (concentration, mass_loss)

pyHalo/Rendering/SpatialDistributions/nfw.py

@@ -35,19 +35,22 @@ def __init__(self, rmax2d_arcsec, rs_arcsec, r_core_arcsec, r_200_arcsec, arcsec
         super(ProjectedNFW, self).__init__()
 
     @classmethod
-    def from_Mhost(cls, m_host, zlens, rmax2d_arcsec, r_core_units_rs, lens_cosmo):
+    def from_Mhost(cls, m_host, zlens, rmax2d_arcsec, r_core_units_rs, lens_cosmo, c_host=None):
         """
 
-        :param m_host:
-        :param zlens:
-        :param rmax2d_arcsec:
-        :param r_core_units_rs:
-        :param lens_cosmo:
+        :param m_host: host halo mass
+        :param zlens: the deflector redshift
+        :param rmax2d_arcsec: the maximum radius in projection to rendering subhalos
+        :param r_core_units_rs: the core size in units rs of a cored nfw profile (subhalos are distributed
+        following a cored nfw profile)
+        :param lens_cosmo: an instance of LensCosmo
+        :param c_host: host halo concentration
         :return:
         """
-        c_model = ConcentrationDiemerJoyce(lens_cosmo)
-        c = c_model.nfw_concentration(m_host, zlens)
-        _, rs_mpc, r_200_mpc = lens_cosmo.nfwParam_physical(m_host, c, zlens)
+        if c_host is None:
+            c_model = ConcentrationDiemerJoyce(lens_cosmo)
+            c_host = c_model.nfw_concentration(m_host, zlens)
+        _, rs_mpc, r_200_mpc = lens_cosmo.nfwParam_physical(m_host, c_host, zlens)
         arcsec_to_kpc = lens_cosmo.cosmo.kpc_proper_per_asec(zlens)
         rs_arcsec = rs_mpc * 1000 / arcsec_to_kpc
         r_200_arcsec = r_200_mpc * 1000 / arcsec_to_kpc