Chengzhang's first HW

.gitignore

@@ -102,3 +102,7 @@ venv.bak/
 
 # mypy
 .mypy_cache/
+.DS_Store
+concentration.pdf
+halomergerrate.pdf
+volumemergerrate.pdf

HW2/Schwarzchild_radius.py

@@ -0,0 +1,12 @@
+import pint
+
+def calculate_schwarzchild_radius(mass):
+    ureg = pint.UnitRegistry()
+    ureg.define("Msun = 1.98855*10**30 * kg")
+
+    schwarzchild_radius = 2 * ureg.newtonian_constant_of_gravitation * ureg.Msun * mass / ureg.c**2
+    return schwarzchild_radius.to_base_units()
+
+if __name__ == "__main__":
+    mass = 1 # in solar mass unit
+    print(calculate_schwarzchild_radius(mass))

HW2/matrix_multiplication.py

@@ -0,0 +1,51 @@
+import time
+
+import numpy as np
+
+np.random.seed(2046)
+
+
+
+def loop_matrix_multiply(m1, m2):
+    result = np.zeros((m1.shape[0], m2.shape[1]))
+    for i in range(m1.shape[0]):
+        for j in range(m2.shape[1]):
+            for k in range(m1.shape[1]):
+                result[i, j] += m1[i, k] * m2[k, j]
+    return result
+
+def test_loop(m1, m2):
+    start = time.time()
+    assert np.allclose(loop_matrix_multiply(m1, m2), np.eye(m1.shape[0]))
+    end = time.time()
+
+    print(f"Loop Method: {end - start:.6f} seconds")
+
+def list_matrix_multiply(m1, m2):
+    return np.array([[sum(m1[i, :] * m2[:, j]) for j in range(m2.shape[1])] for i in range(m1.shape[0])])
+
+def test_list(m1, m2):
+    start = time.time()
+    assert np.allclose(list_matrix_multiply(m1, m2), np.eye(m1.shape[0]))
+    end = time.time()
+
+    print(f"List Method: {end - start:.6f} seconds")
+
+def numpy_matrix_multiply(m1, m2):
+    return np.dot(m1, m2)
+
+def test_numpy(m1, m2):
+    start = time.time()
+    assert np.allclose(numpy_matrix_multiply(m1, m2), np.eye(m1.shape[0]))
+    end = time.time()
+
+    print(f"Numpy Method: {end - start:.6f} seconds")
+
+if __name__ == "__main__":
+    m1 = np.random.rand(500, 500)
+    m2 = np.linalg.inv(m1)
+
+    test_loop(m1, m2)
+    test_list(m1, m2)
+    test_numpy(m1, m2)
+

pbhmergers.py

@@ -25,9 +25,11 @@ def __init__(self,*args,conc_model="ludlow", conc_value=1.,hubble=0.67, **kwargs
         #Hubble constant and related objects!
         #The try except is because the method names change between pint 0.6 and pint 0.8
         try:
-            self.ureg.define(pint.unit.UnitDefinition('hub', '', (),pint.unit.ScaleConverter(hubble)))
-        except AttributeError:
+        #     self.ureg.define(pint.unit.UnitDefinition('hub', '', (),pint.unit.ScaleConverter(hubble)))
+        # except AttributeError:
             self.ureg.define(pint.definitions.UnitDefinition('hub', '', (),pint.converters.ScaleConverter(hubble)))
+        except AttributeError:
+            self.ureg.define('hub = %g' % hubble)
         self.ureg.define("Msolarh = Msolar / hub")
         self.ureg.define("Mpch = Mpc / hub")
         #Factor of R_s at which the maximum circular velocity of the halo is reached.
@@ -63,7 +65,7 @@ def rhocrit(self):
         hubz2 = (self.overden.omega_matter0/aa**3 + self.overden.omega_lambda0) * hubble**2
         #Critical density at redshift in units of kg m^-3
         rhocrit = 3 * hubz2 / (8*math.pi* self.ureg.newtonian_constant_of_gravitation)
-        print "rhocrit = ", rhocrit
+        print("rhocrit = ", rhocrit)
         return rhocrit.to_base_units()
 
     def R200(self, mass):
@@ -116,7 +118,7 @@ def cross_section(self, mass):
         gammaint = sigma**(10/7)*scipy.special.gammainc(5/7,vbysig**2)* scipy.special.gamma(5/7)/2
         #We also need to normalise the probability function for v:
         #Integrate[4*Pi*v^2*P[v, sigma, Vvir], {v, 0, Vvir}]
-        probnorm = math.pi**(3/2)*sigma**3*scipy.special.erf(vbysig) - 2*math.pi/3*np.exp(-(vvir**2/sigma**2))*(3*sigma**2*vvir + 2*vvir**3) *0
+        probnorm = math.pi**(3/2)*sigma**3*scipy.special.erf(vbysig) - 2*math.pi/3*np.exp(-(vvir**2/sigma**2))*(3*sigma**2*vvir + 2*vvir**3)
         assert np.all(probnorm.magnitude > 0)
         cross_section = prefac*(gammaint + cutoff)/probnorm
         assert self.ureg.get_dimensionality('[length]**3 [time]**(-1) [mass]**(-2)') == self.ureg.get_dimensionality(cross_section)
@@ -224,7 +226,7 @@ def mergerhalflife(self,mass,threefac=True, bhmass=None):
             threefac = self.threebodyratio(mass)
             threefac = np.max([threefac, np.ones_like(threefac)],axis=0)
             rate *= threefac
-        return 0.5*(mass/bhmass)/rat
+        return 0.5*(mass/bhmass)/rate
 
     def bias(self,mass):
         """The formula for halo bias in EPS theory (Mo & White 1996), eq. 13"""

power_specs.py

@@ -638,7 +638,7 @@ def plot_z(self,Knot,redshift,title="",ylabel=""):
         """Get a power spectrum in the flat format we use"""
         """for inputting some cosmomc tables"""
         def GetFlat(self,dir):
-                Pk_sdss=np.empty([11, 12])
+                # Pk_sdss=np.empty([11, 12])
                 #For z=2.07 we need to average snap_011 and snap_010
                 z=2.07
                 (k,PF_a)=self.loadpk(dir+self.suf+"snapshot_011"+self.ext, self.box)
@@ -660,12 +660,12 @@ def GetFlat(self,dir):
         flux=flux_pow()
         matter=matter_pow()
         fpdf=flux_pdf()
-        A_knot=knot(("A0.54/","A0.74/","A0.84/","A1.04/", "A1.14/","A1.34/"), (0.54,0.74,0.84,1.04,1.14,1.34),0.94,"best-fit/", 60)
-        AA_knot=knot(("AA0.54/","AA0.74/","AA1.14/","AA1.34/"), (0.54,0.74,1.14,1.34),0.94,"boxcorr400/", 120)
-        B_knot=knot(("B0.33/","B0.53/","B0.73/", "B0.83/", "B1.03/","B1.13/", "B1.33/"), (0.33,0.53,0.73,0.83,1.03, 1.13,1.33),0.93,"best-fit/", 60)
-        C_knot=knot(("C0.11/", "C0.31/","C0.51/","C0.71/","C1.11/","C1.31/","C1.51/"),(0.11, 0.31,0.51,0.71, 1.11,1.31,1.51),0.91,"bf2/", 60)
-        D_knot=knot(("D0.50/","D0.70/","D1.10/","D1.20/", "D1.30/", "D1.50/", "D1.70/"),(0.50, 0.70,1.10,1.20,1.30, 1.50, 1.70),0.90,"bfD/", 48)
-        interp=flux_interp(flux, (AA_knot, B_knot, C_knot, D_knot))
+        # A_knot=knot(("A0.54/","A0.74/","A0.84/","A1.04/", "A1.14/","A1.34/"), (0.54,0.74,0.84,1.04,1.14,1.34),0.94,"best-fit/", 60)
+        # AA_knot=knot(("AA0.54/","AA0.74/","AA1.14/","AA1.34/"), (0.54,0.74,1.14,1.34),0.94,"boxcorr400/", 120)
+        # B_knot=knot(("B0.33/","B0.53/","B0.73/", "B0.83/", "B1.03/","B1.13/", "B1.33/"), (0.33,0.53,0.73,0.83,1.03, 1.13,1.33),0.93,"best-fit/", 60)
+        # C_knot=knot(("C0.11/", "C0.31/","C0.51/","C0.71/","C1.11/","C1.31/","C1.51/"),(0.11, 0.31,0.51,0.71, 1.11,1.31,1.51),0.91,"bf2/", 60)
+        # D_knot=knot(("D0.50/","D0.70/","D1.10/","D1.20/", "D1.30/", "D1.50/", "D1.70/"),(0.50, 0.70,1.10,1.20,1.30, 1.50, 1.70),0.90,"bfD/", 48)
+        # interp=flux_interp(flux, (AA_knot, B_knot, C_knot, D_knot))
 #Thermal parameters for models
 
 #Models where gamma is varied
@@ -692,6 +692,6 @@ def GetFlat(self,dir):
         bf2ag=np.array([1.0184905 ,1.0332849 ,1.0404092 ,1.0538235 ,1.0598905 ,1.0695964 ,1.0771414 ,1.0827464 ,1.0904784 ,1.0964669 ,1.1009428,1.1068357])
         bf2at=np.array([26914.856, 27497.555, 28407.030, 28357.785, 28919.376, 28610.038, 28478.749, 28484.856, 27841.160, 27335.328, 26933.453,26099.643])/1e3
 
-        G_knot=knot(("bf2z/","bf2b/","bf2c/","bf2d/", "bf2T15/","bf2T35/","bf2T45/","bf2a/"),(bf2zg,bf2bg,bf2cg,bf2dg,T15g,T35g,T45g,bf2ag),bf2g,"bf2/",60, (bf2zt,bf2bt,bf2ct,bf2dt,T15t,T35t,T45t,bf2at),bf2t)
-        g_int=flux_interp(flux,(G_knot))
+        # G_knot=knot(("bf2z/","bf2b/","bf2c/","bf2d/", "bf2T15/","bf2T35/","bf2T45/","bf2a/"),(bf2zg,bf2bg,bf2cg,bf2dg,T15g,T35g,T45g,bf2ag),bf2g,"bf2/",60, (bf2zt,bf2bt,bf2ct,bf2dt,T15t,T35t,T45t,bf2at),bf2t)
+        # g_int=flux_interp(flux,(G_knot))