Merge pull request #94 from nanograv/dev-qlf

Evolving MMBulge Amplitude

holodeck/population.py

@@ -634,7 +634,7 @@ def modify(self, pop):
         # Store old version
         pop._mass = pop.mass
         # if `scatter` is `True`, then it is set to the value in `mhost.SCATTER_DEX`
-        pop.mass = self.mhost.mbh_from_host(pop, scatter)
+        pop.mass = self.mhost.mbh_from_host(pop, scatter=scatter)
         return
 
 

holodeck/relations.py

@@ -266,12 +266,12 @@ def __init__(self, mamp=None, mamp_log10=None, mplaw=None, mref=None, bulge_mfra
     def bulge_mass_frac(self, mstar):
         return self._bulge_mfrac
 
-    def mbh_from_host(self, pop, scatter) -> ArrayLike:
+    def mbh_from_host(self, pop, scatter=None) -> ArrayLike:
         host = self.get_host_properties(pop)
         mbulge = host['mbulge']
-        return self.mbh_from_mbulge(mbulge, scatter=scatter)
+        return self.mbh_from_mbulge(mbulge, redz=None, scatter=scatter)
 
-    def mbh_from_mbulge(self, mbulge, scatter):
+    def mbh_from_mbulge(self, mbulge, redz=None, scatter=None):
         """Convert from stellar-bulge mass to black-hole mass.
 
         Parameters
@@ -292,7 +292,7 @@ def mbh_from_mbulge(self, mbulge, scatter):
         mbh = _log10_relation(mbulge, self._mamp, self._mplaw, scatter_dex, x0=self._mref)
         return mbh
 
-    def mbulge_from_mbh(self, mbh, scatter):
+    def mbulge_from_mbh(self, mbh, redz=None, scatter=None):
         """Convert from black-hole mass to stellar-bulge mass.
 
         Parameters
@@ -312,7 +312,7 @@ def mbulge_from_mbh(self, mbh, scatter):
         mbulge = _log10_relation_reverse(mbh, self._mamp, self._mplaw, scatter_dex, x0=self._mref)
         return mbulge
 
-    def mstar_from_mbulge(self, mbulge):
+    def mstar_from_mbulge(self, mbulge, redz=None):
         """Convert from stellar bulge-mass to black-hole mass.
 
         Parameters
@@ -331,7 +331,7 @@ def mstar_from_mbulge(self, mbulge):
         """
         return mbulge / self._bulge_mfrac
 
-    def mbh_from_mstar(self, mstar, scatter):
+    def mbh_from_mstar(self, mstar, redz=None, scatter=None):
         """Convert from total stellar mass to black-hole mass.
 
         Parameters
@@ -351,7 +351,7 @@ def mbh_from_mstar(self, mstar, scatter):
         mbulge = self.mbulge_from_mstar(mstar)
         return self.mbh_from_mbulge(mbulge, scatter)
 
-    def mstar_from_mbh(self, mbh, scatter):
+    def mstar_from_mbh(self, mbh, redz=None, scatter=None):
         """Convert from black-hole mass to total stellar mass.
 
         Parameters
@@ -368,10 +368,10 @@ def mstar_from_mbh(self, mbh, scatter):
             Total stellar mass of host galaxy.  [grams]
 
         """
-        mbulge = self.mbulge_from_mbh(mbh, scatter)
+        mbulge = self.mbulge_from_mbh(mbh, redz=redz, scatter=scatter)
         return self.mstar_from_mbulge(mbulge)
 
-    def dmstar_dmbh(self, mstar):
+    def dmstar_dmbh(self, mstar, redz=None):
         """Calculate the partial derivative of stellar mass versus BH mass :math:`d M_star / d M_bh`.
 
         .. math::
@@ -392,7 +392,7 @@ def dmstar_dmbh(self, mstar):
         plaw = self._mplaw
         fbulge = self._bulge_mfrac
         mbulge = mstar * fbulge
-        mbh = self.mbh_from_mbulge(mbulge, scatter=False)
+        mbh = self.mbh_from_mbulge(mbulge, scatter=False, redz=redz)
         deriv = mstar / (plaw * mbh)
         return deriv
 
@@ -435,6 +435,7 @@ class MMBulge_Redshift(MMBulge_Standard):
 
     TODO: make sure all of the inherited methods from `MMBulge_Standard` are appropriate for
           redshift dependencies!!  In particular, check `dmstar_dmbh`
+          check which redshifts need to be passed into this function. does not pass all cases as is
 
     """
 
@@ -460,13 +461,16 @@ def mbh_from_host(self, pop, scatter):
         host = self.get_host_properties(pop, copy=False)
         mbulge = host['mbulge']    # shape (N, 2)
         redz = host['redz']        # shape (N,)
-        return self.mbh_from_mbulge(mbulge, redz, scatter=scatter)
+        return self.mbh_from_mbulge(mbulge, redz=redz, scatter=scatter)
 
     def mbh_from_mbulge(self, mbulge, redz, scatter):
         scatter_dex = self._scatter_dex if scatter else None
         # Broadcast `redz` to match shape of `mbulge`, if needed
         # NOTE: this will work for (N,) ==> (N,)    or   (N,) ==> (N,X)
-        redz = np.broadcast_to(redz, mbulge.T.shape).T
+        try:
+            redz = np.broadcast_to(redz, mbulge.T.shape).T
+        except:
+            redz = redz
         zmamp = self._mamp * (1.0 + redz)**self._zplaw
         mbh = _log10_relation(mbulge, zmamp, self._mplaw, scatter_dex, x0=self._mref)
         return mbh

holodeck/sams/sam.py

@@ -200,11 +200,26 @@ def mass_stellar(self):
             Galaxy total stellar masses for all MBH. [0, :] is primary, [1, :] is secondary [grams].
 
         """
+        redz = self.redz[np.newaxis, np.newaxis, :]
         # total-mass, mass-ratio ==> (M1, M2)
         masses = utils.m1m2_from_mtmr(self.mtot[:, np.newaxis], self.mrat[np.newaxis, :])
         # BH-masses to stellar-masses
-        masses = self._mmbulge.mstar_from_mbh(masses, scatter=False)
-        return masses
+        mbh_pri = masses[0]
+        mbh_sec = masses[1]
+        args = [mbh_pri[..., np.newaxis], mbh_sec[..., np.newaxis], redz]
+        # Convert to shape (M, Q, Z)
+        mbh_pri, mbh_sec, redz = np.broadcast_arrays(*args)
+        mstar_pri = self._mmbulge.mstar_from_mbh(mbh_pri, redz=redz, scatter=False)
+        mstar_sec = self._mmbulge.mstar_from_mbh(mbh_sec, redz=redz, scatter=False)
+
+        # q = m2 / m1
+        mstar_rat = mstar_sec / mstar_pri
+        # M = m1 + m2
+        mstar_tot = mstar_pri + mstar_sec
+        # args = [mstar_rat[..., np.newaxis]]
+        # # Convert to shape (M, Q, Z)
+        # mstar_rat = np.broadcast_arrays(*args)
+        return mstar_pri, mstar_rat, mstar_tot, redz
 
     @property
     def static_binary_density(self):
@@ -227,17 +242,9 @@ def static_binary_density(self):
             log = self._log
 
             # ---- convert from MBH ===> mstar
-
-            # `mstar_tot` starts as the secondary mass, sorry
-            mstar_pri, mstar_tot = self.mass_stellar()
-            # q = m2 / m1
-            mstar_rat = mstar_tot / mstar_pri
-            # M = m1 + m2
-            mstar_tot = mstar_pri + mstar_tot
+
             redz = self.redz[np.newaxis, np.newaxis, :]
-            args = [mstar_pri[..., np.newaxis], mstar_rat[..., np.newaxis], mstar_tot[..., np.newaxis], redz]
-            # Convert to shape (M, Q, Z)
-            mstar_pri, mstar_rat, mstar_tot, redz = np.broadcast_arrays(*args)
+            mstar_pri, mstar_rat, mstar_tot, redz = self.mass_stellar()
 
             # choose whether the primary mass, or total mass, is used in different calculations
             mass_gsmf = mstar_tot if GSMF_USES_MTOT else mstar_pri
@@ -650,11 +657,8 @@ def gwb_ideal(self, fobs_gw, sum=True, redz_prime=None):
         * There is no coalescence of binaries cutting them off at high-frequencies.
 
         """
-        mstar_pri, mstar_tot = self.mass_stellar()
-        # q = m2 / m1
-        mstar_rat = mstar_tot / mstar_pri
-        # M = m1 + m2
-        mstar_tot = mstar_pri + mstar_tot
+        redz = self.redz[np.newaxis, np.newaxis, :]
+        mstar_pri, mstar_rat, mstar_tot, redz = self.mass_stellar()
 
         # default to using `redz_prime` values if a GMT instance is stored
         if redz_prime is None:
@@ -665,12 +669,7 @@ def gwb_ideal(self, fobs_gw, sum=True, redz_prime=None):
             raise AttributeError(err)
 
         rz = self.redz
-        rz = rz[np.newaxis, np.newaxis, :]
         if redz_prime:
-            args = [mstar_pri[..., np.newaxis], mstar_rat[..., np.newaxis], mstar_tot[..., np.newaxis], rz]
-            # Convert to shape (M, Q, Z)
-            mstar_pri, mstar_rat, mstar_tot, rz = np.broadcast_arrays(*args)
-
             gmt_mass = mstar_tot if GMT_USES_MTOT else mstar_pri
             rz, _ = self._gmt.zprime(gmt_mass, mstar_rat, rz)
             print(f"{self} :: {utils.stats(rz)=}")