From 65fb8aae0e457af952e9e444a22e778028fabf02 Mon Sep 17 00:00:00 2001 From: sugar-joh Date: Mon, 30 Sep 2024 10:33:53 -0700 Subject: [PATCH 1/6] fix grammar --- pbhmergers.py | 4 ++-- power_specs.py | 18 +++++++++--------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/pbhmergers.py b/pbhmergers.py index 95b3d81..0599a0f 100644 --- a/pbhmergers.py +++ b/pbhmergers.py @@ -63,7 +63,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): @@ -224,7 +224,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""" diff --git a/power_specs.py b/power_specs.py index 1e00346..571aabd 100644 --- a/power_specs.py +++ b/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)) From f3e079c150b2fb6c4bdaae9942899f38c9064b46 Mon Sep 17 00:00:00 2001 From: sugar-joh Date: Mon, 30 Sep 2024 10:52:19 -0700 Subject: [PATCH 2/6] HW1 Problem 3 --- HW2/Schwarzchild_radius.py | 12 +++++++++ HW2/matrix_multiplication.py | 51 ++++++++++++++++++++++++++++++++++++ 2 files changed, 63 insertions(+) create mode 100644 HW2/Schwarzchild_radius.py create mode 100644 HW2/matrix_multiplication.py diff --git a/HW2/Schwarzchild_radius.py b/HW2/Schwarzchild_radius.py new file mode 100644 index 0000000..881b0c8 --- /dev/null +++ b/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)) \ No newline at end of file diff --git a/HW2/matrix_multiplication.py b/HW2/matrix_multiplication.py new file mode 100644 index 0000000..53b21b9 --- /dev/null +++ b/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) + From 1a401f96b8cf1a9f61d16edce43cc1344c961884 Mon Sep 17 00:00:00 2001 From: sugar-joh Date: Mon, 30 Sep 2024 11:15:18 -0700 Subject: [PATCH 3/6] fix version problem --- pbhmergers.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/pbhmergers.py b/pbhmergers.py index 0599a0f..b749a90 100644 --- a/pbhmergers.py +++ b/pbhmergers.py @@ -24,10 +24,11 @@ def __init__(self,*args,conc_model="ludlow", conc_value=1.,hubble=0.67, **kwargs #Mpc newton's constant and light speed are already defined. #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.definitions.UnitDefinition('hub', '', (),pint.converters.ScaleConverter(hubble))) + # try: + # self.ureg.define(pint.unit.UnitDefinition('hub', '', (),pint.unit.ScaleConverter(hubble))) + # except AttributeError: + # self.ureg.define(pint.definitions.UnitDefinition('hub', '', (),pint.converters.ScaleConverter(hubble))) + self.ureg.define(f'hub={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. From 97664aa43225100ece1fe8581e13482497da067e Mon Sep 17 00:00:00 2001 From: sugar-joh Date: Mon, 30 Sep 2024 11:24:21 -0700 Subject: [PATCH 4/6] gitignore --- .gitignore | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.gitignore b/.gitignore index 894a44c..a9c339a 100644 --- a/.gitignore +++ b/.gitignore @@ -102,3 +102,7 @@ venv.bak/ # mypy .mypy_cache/ +.DS_Store +concentration.pdf +halomergerrate.pdf +volumemergerrate.pdf From 1865478a75230132dd3b4e7ba8edfc06d2d6beca Mon Sep 17 00:00:00 2001 From: sugar-joh Date: Tue, 1 Oct 2024 19:09:25 -0700 Subject: [PATCH 5/6] remove multipy by 0 --- pbhmergers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pbhmergers.py b/pbhmergers.py index b749a90..d4bb5da 100644 --- a/pbhmergers.py +++ b/pbhmergers.py @@ -117,7 +117,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) From df9a2eecf43c27e072cac414ee430435e5538a75 Mon Sep 17 00:00:00 2001 From: sugar-joh Date: Tue, 1 Oct 2024 19:11:56 -0700 Subject: [PATCH 6/6] Fix AttributeError with modern pint version --- pbhmergers.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pbhmergers.py b/pbhmergers.py index d4bb5da..c72dbb6 100644 --- a/pbhmergers.py +++ b/pbhmergers.py @@ -24,11 +24,12 @@ def __init__(self,*args,conc_model="ludlow", conc_value=1.,hubble=0.67, **kwargs #Mpc newton's constant and light speed are already defined. #Hubble constant and related objects! #The try except is because the method names change between pint 0.6 and pint 0.8 - # try: + try: # self.ureg.define(pint.unit.UnitDefinition('hub', '', (),pint.unit.ScaleConverter(hubble))) # except AttributeError: - # self.ureg.define(pint.definitions.UnitDefinition('hub', '', (),pint.converters.ScaleConverter(hubble))) - self.ureg.define(f'hub={hubble}') + 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.