diff --git a/README.rst b/README.rst index a7134f83..105fe4b1 100644 --- a/README.rst +++ b/README.rst @@ -10,7 +10,10 @@ PopSyCLE Dependencies ------------ -`Galaxia `_ +`galaxia `_ +PopSyCLE requires a custom version of galaxia in order to support +user selected galaxy models. Please follow the installation instructions +found at our galaxia GitHub repo: https://github.com/jluastro/galaxia_. `PyPopStar `_ diff --git a/popsycle/add_pbh_plots.py b/popsycle/add_pbh_plots.py new file mode 100644 index 00000000..67b68f30 --- /dev/null +++ b/popsycle/add_pbh_plots.py @@ -0,0 +1,189 @@ +import matplotlib +matplotlib.use('PDF') +import matplotlib.pyplot as plt +import numpy as np +from matplotlib.backends.backend_pdf import PdfPages + +def print_plots(output_root, galactic_lin_distance, galactic_lin_b, galactic_lin_l, galactocen_lin_spherical_distance, + galactocen_lin_spherical_b, galactocen_lin_spherical_l, rho_lin, r_max, n_lin, cdf_los, + x_cyl, y_cyl, r_cyl, r_proj_los_cyl, n_pbh, d_galac, b_galac, l_galac, area_proj_los_cyl, + mask_obs_cone, field_of_view_radius, l_radian, b_radian, f_cdf_d, pbh_mass): + + ################################################################################## + #Line of Sight - Galactic Coordinates + ################################################################################## + fig_arr = [] + fig1, axs = plt.subplots(3, sharex=True, gridspec_kw={'hspace': 0}, figsize=(8,8)) + fig_arr.append(fig1) + axs[0].plot(galactic_lin_distance, '.', label='Distance', c='C0') + axs[1].plot(galactic_lin_b, '.', label='Latitude',c='C1') + axs[2].plot(galactic_lin_l, '.', label='Longitude', c='C2') + + + plt.legend(loc='best') + # Hide x labels and tick labels for all but bottom plot. + for ax in axs: + ax.label_outer() + ax.legend(loc='upper left') + + axs[0].set_title('Line of Sight Linear Space Galactic Coordinates') + axs[0].set_ylabel('[kpc]') + + axs[1].set_ylabel('[deg]') + + axs[2].set_ylabel('[deg]') + axs[2].set_xlabel('Line-of-Sight Linear Space Coordinate Index') + + axs[0].set_ylim(0, axs[0].set_ylim()[1]) + + ################################################################################## + #Line of Sight - Galactocentric Coordinates + ################################################################################## + fig2, axs = plt.subplots(3, sharex=True, gridspec_kw={'hspace': 0}, figsize=(8,8)) + fig_arr.append(fig2) + axs[0].plot(galactic_lin_distance, + galactocen_lin_spherical_distance, '.', label='Distance', c='C0') + axs[1].plot(galactic_lin_distance, + galactocen_lin_spherical_b, '.', label='Latitude', c='C1') + axs[2].plot(galactic_lin_distance, + galactocen_lin_spherical_l, '.', label='Longitude',c='C2') + + plt.legend(loc='best') + # Hide x labels and tick labels for all but bottom plot. + for ax in axs: + ax.label_outer() + ax.legend(loc='upper left') + + axs[0].set_title('Line of Sight Linear Space\nGalactocentric Coordinates') + axs[0].set_ylabel('[kpc]') + + axs[1].set_ylabel('[deg]') + + axs[2].set_ylabel('[deg]') + axs[2].set_xlabel('Galactic Coordinate Distance [kpc]') + + axs[0].set_ylim(0, axs[0].set_ylim()[1]) + + ################################################################################## + #Denisty along Line of Sight + ################################################################################## + fig3, (ax1, ax2) = plt.subplots(1, 2, figsize=(9,4), + sharey=True, gridspec_kw={'wspace': 0}) + fig_arr.append(fig3) + # + ax1.semilogy(galactic_lin_distance, + rho_lin, '.',c='k') + ax1.set_ylabel('DENSITY Msun*pc^-3') + ax1.set_xlabel('Galactic Coordinate Distance [kpc]') + ax1.set_xlim(0, ax1.set_xlim()[1]) + ax2.semilogy(galactocen_lin_spherical_distance, + rho_lin, '.',c='k') + ax2.set_xlabel('Galactocentric Distance [kpc]') + ax2.set_xlim(0, ax2.set_xlim()[1]) + + plt.suptitle('Milky Way Density along Line-of-Sight') + + ################################################################################## + #Cumulative Projected Desnity + ################################################################################## + fig4 = plt.figure(figsize=(10,7)) + fig_arr.append(fig4) + + plt.plot(galactic_lin_distance, + np.cumsum(rho_lin) * r_max / n_lin, + '.', label='Density',c='k') + plt.legend(loc='best') + plt.ylabel('Cumulative Projected Density} [M_sun*pc^{-2}]$') + plt.xlabel('Galactic Distance [kpc]') + plt.xlim(0, plt.xlim()[1]) + + ################################################################################## + #Discrete CDF + ################################################################################## + fig5 = plt.figure(figsize=(10,7)) + fig_arr.append(fig5) + plt.plot(cdf_los, lw=3) + plt.ylabel('CDF') + plt.xlabel('Line-of-Sight Linear Space Coordinate Index') + + ################################################################################## + #Interpolated CDF + ################################################################################## + fig6 = plt.figure(figsize=(10,7)) + fig_arr.append(fig6) + x = np.linspace(0,1,100) + y = f_cdf_d(x) + plt.plot(y, x, lw=3) + plt.ylabel('CDF') + plt.xlabel('Galactic Coordinate Distance along LOS [kpc]') + + ################################################################################## + #Cylindrical distribution of PBHs + #Infered Radii Distribution + #Galactic Distance Distribution + ################################################################################## + fig7, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(14,4)) + fig_arr.append(fig7) + # Create a plot to look at the projected cylindrical distribution of PBH + ax1.plot(x_cyl * 1000, y_cyl * 1000, ',') + ax1.set_aspect('equal', 'datalim') + ax1.set_ylabel('Cylindrical LOS Samples l-axis [pc]') + ax1.set_xlabel('Cylindrical LOS Samples b-axis [pc]') + # Create a histogram of the infered radii to make sure that it shows the expected linear trend + bins = 20 + ax2.hist(r_cyl * 1000, bins=bins, label='Sampled Number') + ax2.plot((0, r_proj_los_cyl * 1000), (0, 2 * n_pbh / bins), c='C2', label='Expected Number') + ax2.legend(loc='upper left') + ax2.set_ylabel('Number of PBH per LOS Cylindrical Annuli') + ax2.set_xlabel('Cylindrical Radius [pc]') + + # Create a histogram of the galactic distance distribution + bins=100 + n_per_bin, _, _ = ax3.hist(d_galac, bins=bins, label='Sampled Number') + ax3.plot(galactic_lin_distance, + rho_lin * area_proj_los_cyl / pbh_mass * r_max / bins * 1000**3,c='C2', label='Approximate Expectation') + ax3.legend(loc='best') + ax3.set_ylabel('Number of PBH') + ax3.set_xlabel('Galactic Distance [kpc]') + + fig7.tight_layout(pad=2.0) + + plt.suptitle('Distribution of PBHs in Cylindrical Line-of-Sight Tube (i.e. not light cone)') + + ################################################################################## + #Checking Light Cone Boundaries + ################################################################################## + fig8, (ax1, ax2) = plt.subplots(2, 1, figsize=(5,8), + sharex=True, gridspec_kw={'hspace': 0}) + fig_arr.append(fig8) + ax1.plot(d_galac[~mask_obs_cone], r_cyl[~mask_obs_cone] * 1000, + '.', alpha=0.1, label='Outside Light Cone') + ax1.plot(d_galac[mask_obs_cone], r_cyl[mask_obs_cone] * 1000, + '.', alpha=0.1, label='Inside Light Cone') + ax1.plot((0,r_max), + (0, field_of_view_radius * np.pi / 180 * r_max *1000), + label='Light Cone Boundry') + ax1.legend(loc='upper left') + ax1.set_ylabel('Radius on LOS Cylinder [pc]') + + ax2.plot(d_galac[~mask_obs_cone], np.sqrt((l_galac[~mask_obs_cone] - l_radian)**2 + + (b_galac[~mask_obs_cone] - b_radian)**2) * 180/np.pi, + '.', alpha=0.1, label='Outside Light Cone' ) + ax2.plot(d_galac[mask_obs_cone], np.sqrt((l_galac[mask_obs_cone] - l_radian)**2 + + (b_galac[mask_obs_cone] - b_radian)**2) * 180/np.pi, + '.', alpha=0.1, label='Inside Light Cone') + ax2.plot((0,r_max), + (field_of_view_radius, field_of_view_radius), label='Light Cone Boundry') + ax2.set_ylim(0,5*field_of_view_radius) + ax2.legend(loc='best') + ax2.set_ylabel('Field of View Radius [deg]') + ax2.set_xlabel('Galactic Coordinate Distance [kpc]') + + ################################################################################## + #Saving all plots to PDF + ################################################################################## + pp = PdfPages(output_root+'_pbh_diagnostic_plots.pdf') + for fig in fig_arr: + pp.savefig(fig) + plt.close(fig) + pp.close() diff --git a/popsycle/converter.py b/popsycle/converter.py index 4d63fcdf..9684846b 100644 --- a/popsycle/converter.py +++ b/popsycle/converter.py @@ -47,6 +47,7 @@ def convert_h5_array_dtype_to_compound_dtype(hdf5_file): hdf5_file_new = hdf5_file.replace('.h5', '.compound_dtype.h5') f_in = h5py.File(hdf5_file, 'r') f_out = h5py.File(hdf5_file_new, 'w') + f_out['add_pbh'] = False # Looping over all of the datasets for key in f_in: diff --git a/popsycle/data/pbh_config.yaml b/popsycle/data/pbh_config.yaml new file mode 100644 index 00000000..cb551090 --- /dev/null +++ b/popsycle/data/pbh_config.yaml @@ -0,0 +1,9 @@ +# Example configuration parameters for PBHs + +fdm: 1 +pbh_mass: 40 +r_max: 16.6 +r_s: 18.6 +gamma: 1 +v_esc: 550 +rho_0 : 0.0093 diff --git a/popsycle/data/popsycle_config.yaml b/popsycle/data/popsycle_config.yaml index 5106aba6..92ed5bd0 100644 --- a/popsycle/data/popsycle_config.yaml +++ b/popsycle/data/popsycle_config.yaml @@ -6,6 +6,7 @@ n_obs: 101 theta_frac: 2 blend_rad: 0.75 isochrones_dir: /Users/myself/popsycle_isochrones +galaxia_galaxy_model_filename: /Users/myself/galaxia_galaxy_model_filename bin_edges_number: None BH_kick_speed_mean: 50 NS_kick_speed_mean: 400 diff --git a/popsycle/data/radial_velocity_profile_middle.csv b/popsycle/data/radial_velocity_profile_middle.csv new file mode 100644 index 00000000..2319354c --- /dev/null +++ b/popsycle/data/radial_velocity_profile_middle.csv @@ -0,0 +1,62 @@ +r,v +0.010792642776929309, 111.71180517031452 +0.012590321940393834, 118.88165941281024 +0.014687431970009386, 126.54067684534101 +0.01713384764066691, 134.0005466632754 +0.01993184060733583, 142.11168012051704 +0.023448011680987736, 150.78286434254127 +0.027455988964315874, 159.3571491494456 +0.03213404568533922, 168.4661816124074 +0.03764430047370784, 177.73468233692108 +0.043859894720648336, 186.54417763189218 +0.05108586704256386, 196.01349797292718 +0.05934511962521002, 204.67635482054385 +0.07083192093472244, 215.1824749096273 +0.08253367166054373, 223.90379647052103 +0.09684442649610563, 232.571184061413 +0.11350401149180162, 242.54850927105662 +0.13427738741960332, 251.89354878059038 +0.1586056689369005, 261.0497880274381 +0.186179188876288, 269.5935014641335 +0.21617856025114549, 276.95879153868964 +0.2526914812832889, 283.38932791498735 +0.2962593087128887, 290.36094394136614 +0.35792137580603484, 297.77197612125576 +0.4309060654759373, 304.89151975943804 +0.5059594624841087, 309.244695083941 +0.5778483816986392, 312.37804867591933 +0.6816933312568133, 314.7048135621079 +0.927699047785235, 320.24095532420716 +1.0671697449159139, 319.8132280123578 +1.4727675574625494, 314.96530138192946 +1.7718422284532012, 311.9846498190986 +2.1197411537756987, 307.83085241102907 +2.4598627413134806, 301.78571980162417 +2.8847016095688685, 295.5628794127563 +3.369909155131431, 289.9952727173605 +4.037236963747592, 283.1343596020251 +4.820932000343341, 279.00220082094125 +5.650248971118326, 275.5083826876096 +6.703101332022676, 271.60158332331133 +7.770476112914143, 265.4265828672343 +9.096562802640236, 259.22858778705637 +10.611733993556404, 252.9216664462353 +12.3792800416127, 246.90957929498413 +14.44123782613885, 240.62655413150787 +16.846646109472406, 235.17599223163336 +19.652711807301024, 229.3228663966612 +22.766181747138315, 223.97434631744974 +26.744873675700347, 218.17494681231219 +31.068822348129252, 212.83026716581725 +36.396415667350134, 206.61996451113927 +43.05764653615527, 199.86831942509124 +50.22955388134484, 193.3372953786207 +58.59605171407531, 186.69858418697999 +68.54786206569099, 179.3094384837924 +79.24703527090092, 171.75301799473309 +93.02409457669958, 163.93495203048826 +108.51867546639197, 155.9656915244779 +126.5941149824405, 147.0301382654459 +147.68029446830633, 138.0518741736095 +172.27869855775518, 128.43932707311797 +192.7048730125345, 121.55334325539347 diff --git a/popsycle/data/radial_velocity_profile_shallow.csv b/popsycle/data/radial_velocity_profile_shallow.csv new file mode 100644 index 00000000..8e412d1d --- /dev/null +++ b/popsycle/data/radial_velocity_profile_shallow.csv @@ -0,0 +1,66 @@ +r,v +0.010792642776929309, 204.26643023097807 +0.012590321940393834, 210.92487127516114 +0.014687431970009386, 218.17494681231219 +0.01713384764066691, 225.4158429807519 +0.019987751131241467, 232.89705353032554 +0.023317015749352955, 240.21341604790845 +0.02720082013657757, 247.6177427638289 +0.031731531344141896, 255.39654951292243 +0.03701690156357718, 262.5159561489671 +0.04318263075634799, 269.98842786464206 +0.05037535612850712, 277.3556810971669 +0.05876613954792163, 284.43477324360026 +0.0685545358439959, 291.3605743829068 +0.07997333874472501, 298.11329956878745 +0.0932941173219056, 304.6732953099007 +0.10883367461568658, 311.0211336366553 +0.1269615820414861, 316.5932048851325 +0.14810896876722457, 322.0805612681069 +0.17277877509531006, 326.913144325661 +0.20155771370168304, 331.43832299459757 +0.235130223201544, 335.44920768526305 +0.2742947458941446, 338.73163728292803 +0.3199827168140894, 342.046186032515 +0.3732807156985778, 343.8140575474295 +0.4354563087024429, 346.3837933152392 +0.5079881944447303, 346.3837933152392 +0.5926013713388398, 346.3837933152392 +0.6913081625775506, 346.3837933152392 +0.8064560744546262, 343.61717701711524 +0.9407836262193813, 341.654561829232 +1.0974854792445174, 338.15005999198166 +1.2802884144496132, 334.10688186636213 +1.4935399649225864, 329.73408677356304 +1.7423118116552716, 324.4878564913184 +2.03252040141478, 318.41185653421763 +2.3710676553598873, 312.4496290180272 +2.7660051148222395, 306.07263557365945 +3.2267254280694635, 299.1396191737312 +3.764185732107431, 292.02890586026 +4.391168242126657, 285.9048791942064 +5.122584246093038, 278.7891986399261 +5.975828734271417, 271.383867634092 +6.971194097701775, 265.692780180841 +8.132352734462563, 258.3391991378944 +9.48692004136333, 253.79240406567155 +11.067111180485346, 247.90157653062977 +12.970922127990187, 241.5318118152786 +15.060949011985775, 236.1208479753841 +17.618863818498795, 229.54628085055066 +20.496060936083154, 223.2314835842758 +23.909992300202614, 217.3019029730845 +27.892565970532186, 211.4086969647845 +32.53849799081271, 206.147098218244 +37.958280805597894, 199.9828369269898 +44.28081106028432, 193.7807781680185 +51.65645510129202, 186.37803645873726 +60.26062508202531, 179.36080012164197 +70.29795072379233, 172.80561961785477 +82.0071459470937, 165.2548659942107 +95.6666861145325, 158.0340430777453 +111.60143090506897, 149.06535054405288 +130.19035032893146, 139.90937104204863 +151.8755376280755, 130.51589873409813 +177.1727234125981, 121.49765534647908 +195.4228471568372, 116.70016866901982 diff --git a/popsycle/data/radial_velocity_profile_steep.csv b/popsycle/data/radial_velocity_profile_steep.csv new file mode 100644 index 00000000..dd5109fd --- /dev/null +++ b/popsycle/data/radial_velocity_profile_steep.csv @@ -0,0 +1,64 @@ +r,v +0.010792642776929309, 51.189721234746486 +0.012590321940393834, 56.2637730438137 +0.014687431970009386, 61.73460135767326 +0.01713384764066691, 67.58236521833585 +0.019987751131241467, 73.81473571195039 +0.023317015749352955, 80.43733893511352 +0.02720082013657757, 87.40343466110791 +0.031731531344141896, 94.80975218695623 +0.03701690156357718, 102.66708492081379 +0.04318263075634799, 110.8576402505859 +0.05037535612850712, 119.0179284411037 +0.05876613954792163, 127.63259980070472 +0.0685545358439959, 136.4012214489796 +0.07997333874472501, 145.5219851562832 +0.0932941173219056, 155.0748727428563 +0.10883367461568658, 164.21690036223077 +0.1269615820414861, 173.59930211335876 +0.14810896876722457, 183.3076420995835 +0.17277877509531006, 193.00534947186358 +0.20155771370168304, 202.28701613851908 +0.235130223201544, 211.89363246804209 +0.2742947458941446, 221.19495574296548 +0.3199827168140894, 230.37612812090563 +0.3732807156985778, 239.11517900904093 +0.4354563087024429, 247.7596190021874 +0.5079881944447303, 255.5428827333543 +0.5926013713388398, 263.2688804002836 +0.6913081625775506, 269.98842786464206 +0.8064560744546262, 275.92953461668526 +0.9407836262193813, 282.1629525216277 +1.0974854792445174, 284.43477324360026 +1.2802884144496132, 287.7119979375535 +1.4935399649225864, 289.6964297280955 +1.7423118116552716, 290.02849651765854 +2.03252040141478, 292.53116033258596 +2.3710676553598873, 288.7025088087155 +2.7660051148222395, 288.5371871710757 +3.2267254280694635, 285.9048791942064 +3.764185732107431, 282.1629525216277 +4.340214672149436, 278.4168360406884 +6.424289627691652, 275.5083826876096 +7.582335742589992, 269.21628633919846 +8.942171302597005, 263.2521253122571 +10.318607100674031, 258.0434148632248 +12.037328396676813, 252.19833425421467 +14.042328922474109, 246.06245599323356 +16.381292847455335, 240.62655413150787 +19.109846866257, 236.1208479753841 +22.29288314740797, 229.84889395002654 +26.00610263923727, 224.89996081862702 +30.337815436904652, 218.4250315448684 +35.39104102800242, 213.35512334606335 +41.285958366076336, 207.68788330277627 +48.16276404123605, 200.44156345555936 +56.185006522648195, 193.66981219227398 +65.54347576993862, 187.55607596607996 +76.46074071866175, 180.18459004092205 +89.19644255313948, 172.31141178905165 +104.05346965457406, 163.6534877832804 +121.38516108089176, 156.32345022691644 +141.603709895957, 146.35792791879018 +165.18996620135306, 135.46686778654006 +187.38180510452284, 127.19470448022776 diff --git a/popsycle/data/slurm_config.yaml b/popsycle/data/slurm_config.yaml index 83fbb641..ae66bada 100644 --- a/popsycle/data/slurm_config.yaml +++ b/popsycle/data/slurm_config.yaml @@ -9,6 +9,11 @@ account: ulens queue: regular # Name of the resource that will be used for the run resource: haswell +# Set to true if slurm scheduler requires a CONSTRAINT specification +# Set to false if slurm scheduler does not +include_constraint: true +# Options for slurm's srun command +srun_options: # Additional lines to be run before executing run.py additional_lines: - module load cray-hdf5/1.10.5.2 diff --git a/popsycle/ebf.py b/popsycle/ebf.py old mode 100644 new mode 100755 index e16c2922..25874bfa --- a/popsycle/ebf.py +++ b/popsycle/ebf.py @@ -259,6 +259,9 @@ import sys import time import os +from itertools import groupby +from operator import itemgetter + __version__ = "0.0.14" @@ -2749,13 +2752,39 @@ def read(self,i,nsize=1): else: return None - def read_ind(self,ind): - if numpy.max(ind)__ + pbh_config_filename : str + Name of pbh_config.yaml file containing the PBH parameters + that will be passed along to the run_on_slurm.py command in the + slurm script. If set to None, `add_pbh` is skipped over. + Default None. + seed : int If non-None, all random sampling will be seeded with the specified seed, forcing identical output for PyPopStar and PopSyCLE. @@ -524,23 +615,48 @@ def generate_slurm_script(slurm_config_filename, popsycle_config_filename, """ # Check for files + # Enforce slurm_config_filename is an absolute path + slurm_config_filename = os.path.abspath(slurm_config_filename) + # Check for slurm config file. Exit if not present. if not os.path.exists(slurm_config_filename): raise Exception('Slurm configuration file {0} does not exist. ' 'Write out file using ' 'run.generate_slurm_config_file ' 'before proceeding.'.format(slurm_config_filename)) + + # Enforce popsycle_config_filename is an absolute path + popsycle_config_filename = os.path.abspath(popsycle_config_filename) + # Check for popsycle config file. Exit if not present. if not os.path.exists(popsycle_config_filename): raise Exception('PopSyCLE configuration file {0} does not exist. ' 'Write out file using ' 'run.generate_popsycle_config_file ' 'before proceeding.'.format(popsycle_config_filename)) - # Enforce popsycle_config_filename is an absolute path - popsycle_config_filename = os.path.abspath(popsycle_config_filename) - popsycle_config = load_config_file(popsycle_config_filename) - - # Load the slurm configuration file + if pbh_config_filename is not None: + # Enforce pbh_config_filename is an absolute path + pbh_config_filename = os.path.abspath(pbh_config_filename) + # Check for pbh config file, if provided. Exit if not present. + if not os.path.exists(pbh_config_filename): + raise Exception('PBH configuration file {0} does not exist. ' + 'Write out file using ' + 'run.generate_pbh_config_file ' + 'before proceeding.'.format(pbh_config_filename)) + + # Load the configuration files + # Load slurm config slurm_config = load_config_file(slurm_config_filename) + slurm_config['include_constraint'] = bool(slurm_config['include_constraint']) # Enforce boolean for 'include_constraint' + if slurm_config['srun_options'] is None: + slurm_config['srun_options'] = '' + # Load popsycle config + popsycle_config = load_config_file(popsycle_config_filename) + if popsycle_config['bin_edges_number'] == 'None': # Enforce None for 'bin_edges_number' + popsycle_config['bin_edges_number'] = None + # Load pbh config, if provided. + if pbh_config_filename is not None: + pbh_config = load_config_file(pbh_config_filename) + pbh_config['diagnostic_plots'] = bool(pbh_config['diagnostic_plots']) # Check pipeline stages for valid inputs _check_slurm_config(slurm_config, walltime) @@ -549,10 +665,9 @@ def generate_slurm_script(slurm_config_filename, popsycle_config_filename, longitude=longitude, latitude=latitude, area=area, + galaxia_galaxy_model_filename=popsycle_config['galaxia_galaxy_model_filename'], seed=seed) if not skip_perform_pop_syn: - if popsycle_config['bin_edges_number'] == 'None': - popsycle_config['bin_edges_number'] = None _check_perform_pop_syn(ebf_file='test.ebf', output_root=output_root, iso_dir=popsycle_config['isochrones_dir'], @@ -562,6 +677,19 @@ def generate_slurm_script(slurm_config_filename, popsycle_config_filename, additional_photometric_systems=[popsycle_config['photometric_system']], overwrite=overwrite, seed=seed) + if pbh_config_filename is not None: + _check_add_pbh(hdf5_file='test.h5', + ebf_file='test.ebf', + fdm=pbh_config['fdm'], + pbh_mass=pbh_config['pbh_mass'], + r_max=pbh_config['r_max'], + r_s=pbh_config['r_s'], + gamma=pbh_config['gamma'], + v_esc=pbh_config['v_esc'], + rho_0=pbh_config['rho_0'], + diagnostic_plots=pbh_config['diagnostic_plots'], + new_output_root=None, + seed=seed) if not skip_calc_events: _check_calc_events(hdf5_file='test.h5', output_root2=output_root, @@ -580,7 +708,6 @@ def generate_slurm_script(slurm_config_filename, popsycle_config_filename, overwrite=overwrite, output_file='default') - # Make a run directory for the PopSyCLE output path_run = os.path.abspath(path_run) if not os.path.exists(path_run): @@ -607,6 +734,8 @@ def generate_slurm_script(slurm_config_filename, popsycle_config_filename, account = slurm_config['account'] # Queue queue = slurm_config['queue'] + # Options for slurm's srun command + srun_options = slurm_config['srun_options'] # Name of the resource that will be ussed for the run resource = slurm_config['resource'] # Maximum number of cores per node @@ -640,8 +769,10 @@ def generate_slurm_script(slurm_config_filename, popsycle_config_filename, # Job name #SBATCH --account={account} #SBATCH --qos={queue} -#SBATCH --constraint={resource} -#SBATCH --nodes=1 +""" + if slurm_config['include_constraint']: + slurm_template += '#SBATCH --constraint={resource}\n' + slurm_template += """#SBATCH --nodes=1 #SBATCH --time={walltime} #SBATCH --job-name={jobname} #SBATCH --output={jobname}-%j.out @@ -660,11 +791,11 @@ def generate_slurm_script(slurm_config_filename, popsycle_config_filename, cd {path_run} """ - for line in slurm_config['additional_lines']: - slurm_template += '%s\n' % line + if slurm_config['additional_lines'] is not None: + for line in slurm_config['additional_lines']: + slurm_template += '%s\n' % line slurm_template += """ -srun -N 1 -n 1 {path_python} {run_filepath}/run.py --output-root={output_root} --field-config-filename={field_config_filename} --popsycle-config-filename={popsycle_config_filename} --n-cores-calc-events={n_cores_calc_events} {optional_cmds} - +srun -N 1 -n 1 {srun_options} {path_python} {run_filepath}/run.py --output-root={output_root} --field-config-filename={field_config_filename} --popsycle-config-filename={popsycle_config_filename} --n-cores-calc-events={n_cores_calc_events} {optional_cmds} date echo "All done!" """ @@ -690,6 +821,9 @@ def generate_slurm_script(slurm_config_filename, popsycle_config_filename, if skip_refine_events: optional_cmds += '--skip-refine-events ' + if pbh_config_filename: + optional_cmds += '--pbh-config-filename={0}'.format(pbh_config_filename) + # Populate the mpi_template specified inputs job_script = slurm_template.format(**locals()) @@ -790,6 +924,11 @@ def run(): optional.add_argument('--skip-refine-events', help="Skip running refine_events.", action='store_true') + optional.add_argument('--pbh-config-filename', type=str, + help='Name of configuration file containing ' + 'pbh inputs. Default None.', + default=None) + args = parser.parse_args() t0 = time.time() @@ -812,6 +951,16 @@ def run(): Exiting...""".format(args.popsycle_config_filename)) sys.exit(1) + # Check for pbh config file, if provided. Exit if not present. + if args.pbh_config_filename is not None: + if not os.path.exists(args.pbh_config_filename): + print("""Error: PBH configuration file {0} missing, + cannot continue. In order to execute run.py with 'add_pbh', + generate a PBH configuration file using + popsycle.run.generate_pbh_config_file. + Exiting...""".format(args.pbh_config_filename)) + sys.exit(1) + # Load the config files for field parameters field_config = load_config_file(args.field_config_filename) @@ -840,6 +989,7 @@ def run(): longitude=field_config['longitude'], latitude=field_config['latitude'], area=field_config['area'], + galaxia_galaxy_model_filename=popsycle_config['galaxia_galaxy_model_filename'], seed=args.seed) if not args.skip_perform_pop_syn: _check_perform_pop_syn(ebf_file=filename_dict['ebf_filename'], @@ -851,6 +1001,21 @@ def run(): additional_photometric_systems=additional_photometric_systems, overwrite=args.overwrite, seed=args.seed) + if args.pbh_config_filename is not None: + pbh_config = load_config_file(args.pbh_config_filename) + pbh_config['diagnostic_plots'] = bool(pbh_config['diagnostic_plots']) + _check_add_pbh(hdf5_file=filename_dict['hdf5_filename'], + ebf_file=filename_dict['ebf_filename'], + fdm=pbh_config['fdm'], + pbh_mass=pbh_config['pbh_mass'], + r_max=pbh_config['r_max'], + r_s=pbh_config['r_s'], + gamma=pbh_config['gamma'], + v_esc=pbh_config['v_esc'], + rho_0=pbh_config['rho_0'], + diagnostic_plots=pbh_config['diagnostic_plots'], + new_output_root=None, + seed=args.seed) if not args.skip_calc_events: _check_calc_events(hdf5_file=filename_dict['hdf5_filename'], output_root2=args.output_root, @@ -881,6 +1046,7 @@ def run(): longitude=field_config['longitude'], latitude=field_config['latitude'], area=field_config['area'], + galaxia_galaxy_model_filename=popsycle_config['galaxia_galaxy_model_filename'], seed=args.seed) if not args.skip_perform_pop_syn: @@ -902,6 +1068,33 @@ def run(): overwrite=args.overwrite, seed=args.seed) + # If optional pbh_config_filename is provided: + if args.pbh_config_filename is not None: + # Check if .h5 file exists from perform popsyn, + # use as input for following function + if not os.path.exists(filename_dict['hdf5_filename']): + print("""Error: hdf5 file was not created properly by + synthetic.perform_pop_syn. + Exiting....""") + sys.exit(1) + + print("** COMMENT ON WARNING **") + print(" run.py executes 'add_pbh' without ") + print(" using the 'new_output_root' argument") + print(" and instead replaces %s." % filename_dict['hdf5_filename']) + print(" Therefore ignore the following warning:") + synthetic.add_pbh(hdf5_file=filename_dict['hdf5_filename'], + ebf_file=filename_dict['ebf_filename'], + fdm=pbh_config['fdm'], + pbh_mass=pbh_config['pbh_mass'], + r_max=pbh_config['r_max'], + r_s=pbh_config['r_s'], + gamma=pbh_config['gamma'], + v_esc=pbh_config['v_esc'], + rho_0=pbh_config['rho_0'], + diagnostic_plots=pbh_config['diagnostic_plots'], + seed=args.seed) + if not args.skip_calc_events: # Remove calc_events output if already exists and overwrite=True if _check_for_output(filename_dict['events_filename'], diff --git a/popsycle/synthetic.py b/popsycle/synthetic.py index 5bcb6f0a..cdef14e1 100755 --- a/popsycle/synthetic.py +++ b/popsycle/synthetic.py @@ -7,6 +7,7 @@ - perform_pop_syn - calc_events - refine_events +- add_pbh """ import numpy as np import h5py @@ -23,6 +24,7 @@ from popstar import synthetic, evolution, reddening, ifmr from scipy.interpolate import interp1d from scipy.spatial import cKDTree +from scipy import special, integrate, interpolate import time import datetime import gc @@ -34,8 +36,10 @@ import inspect import numpy.lib.recfunctions as rfn import copy +from distutils import spawn from popsycle import ebf from popsycle.filters import transform_ubv_to_ztf +import shutil from popsycle import utils @@ -172,7 +176,7 @@ def write_galaxia_params(output_root, def _check_run_galaxia(output_root, longitude, latitude, area, - seed): + galaxia_galaxy_model_filename, seed): """ Check that the inputs to run_galaxia are valid @@ -194,6 +198,9 @@ def _check_run_galaxia(output_root, longitude, latitude, area, area : float Area of the sky that will be generated, in square degrees + galaxia_galaxy_model_filename : str + Name of the galaxia galaxy model, as outlined at https://github.com/jluastro/galaxia + seed : int Seed Galaxia will use to generate objects. If not set, script will generate a seed from the current time. Setting this seed guarantees @@ -215,12 +222,45 @@ def _check_run_galaxia(output_root, longitude, latitude, area, if type(area) != float: raise Exception('area (%s) must be an integer or a float.' % str(area)) + # Check that galaxia is in the executable PATH + if spawn.find_executable('galaxia') is None: + raise Exception('galaxia is not an executable currently in $PATH') + + # Confrim that galaxia is the PopSyCLE compliant version + stdout, _ = utils.execute('galaxia --version') + galaxia_version = stdout.replace('\n', '').split()[1] + if galaxia_version != '0.7.2.1': + raise Exception('galaxia must be version 0.7.2.1 installed from https://github.com/jluastro/galaxia') + + # Check the galaxia_galaxy_model_filename for correct type and existence + if type(galaxia_galaxy_model_filename) != str: + raise Exception('galaxia_galaxy_model_filename (%s) must be a string.' % str(galaxia_galaxy_model_filename)) + + if not os.path.exists(galaxia_galaxy_model_filename): + raise Exception('galaxia_galaxy_model_filename (%s) does not exist' % galaxia_galaxy_model_filename) + + # Check that GalaxiaData is a valid galaxia folder + GalaxiaData = None + for line in open(galaxia_galaxy_model_filename, 'r'): + if 'GalaxiaData' in line: + GalaxiaData = line.replace('\n','').split()[1] + break + + if GalaxiaData is None: + raise Exception('GalaxiaData missing from ' + 'galaxia_galaxy_model_filename (%s)' % galaxia_galaxy_model_filename) + + if not os.path.exists(GalaxiaData): + raise Exception('GalaxiaData (%s) in galaxia_galaxy_model_filename ' + '(%s) does not exist' % (GalaxiaData, galaxia_galaxy_model_filename)) + if seed is not None: if type(seed) != int: raise Exception('seed (%s) must be None or an integer.' % str(seed)) def run_galaxia(output_root, longitude, latitude, area, + galaxia_galaxy_model_filename, seed=None): """ Given an object root, sky location and area, creates the parameter @@ -245,6 +285,9 @@ def run_galaxia(output_root, longitude, latitude, area, area : float Area of the sky that will be generated, in square degrees + galaxia_galaxy_model_filename : str + Name of the galaxia galaxy model, as outlined at https://github.com/jluastro/galaxia + Optional Parameters ------------------- seed : int @@ -255,7 +298,7 @@ def run_galaxia(output_root, longitude, latitude, area, # Error handling/complaining if input types are not right. _check_run_galaxia(output_root, longitude, latitude, area, - seed) + galaxia_galaxy_model_filename, seed) # Writes out galaxia params to disk write_galaxia_params(output_root=output_root, @@ -265,8 +308,9 @@ def run_galaxia(output_root, longitude, latitude, area, seed=seed) # Execute Galaxia - cmd = 'galaxia -r galaxia_params.%s.txt' % output_root - print('Executing with galaxia_params.%s.txt' % output_root) + cmd = 'galaxia -r galaxia_params.%s.txt %s' % (output_root, galaxia_galaxy_model_filename) + print('** Executing Galaxia with galaxia_params.%s.txt ' + 'and %s **' % (output_root, galaxia_galaxy_model_filename)) t0 = time.time() _ = utils.execute(cmd) t1 = time.time() @@ -291,6 +335,7 @@ def run_galaxia(output_root, longitude, latitude, area, line1 = 'longitude , ' + str(longitude) + '\n' line2 = 'latitude , ' + str(latitude) + '\n' line3 = 'area , ' + str(area) + '\n' + line3b = 'galaxia_galaxy_model_filename , ' + galaxia_galaxy_model_filename + '\n' line4 = 'seed , ' + str(seed) + '\n' line8 = 'VERSION INFORMATION' + '\n' @@ -305,7 +350,7 @@ def run_galaxia(output_root, longitude, latitude, area, line18 = output_root + '.ebf : ebf file' + '\n' with open(output_root + '_galaxia.log', 'w') as out: - out.writelines([line0, dash_line, line1, line2, line3, line4, + out.writelines([line0, dash_line, line1, line2, line3, line3b, line4, empty_line, line8, dash_line, line9, line10, line11, empty_line, line12, dash_line, line13, empty_line, line17, dash_line, line18]) @@ -570,11 +615,13 @@ def perform_pop_syn(ebf_file, output_root, iso_dir, long_bin_edges[wrap_id] -= 360 ########## - # Create h5py file to store lat/long binned output + # Create h5py file to store lat/long binned output and + # that this hdf5 file was created without / before add_pbh ########## h5file = h5py.File(output_root + '.h5', 'w') h5file['lat_bin_edges'] = lat_bin_edges h5file['long_bin_edges'] = long_bin_edges + h5file['add_pbh'] = False h5file.close() ########## @@ -644,6 +691,29 @@ def perform_pop_syn(ebf_file, output_root, iso_dir, num_stars_in_bin = 2e6 num_bins = int(math.ceil(len_adx / num_stars_in_bin)) + # Create a KDTree from randomly selected stars in the + # pop_id / age_bin used for calculating extinction to luminous + # white dwarfs. Because the same KDTree is used for each sub-bin, + # two compact objects randomly selected to have nearly identical + # positions would have identical extinctions. This low + # probability event is a reasonable trade-off for the reduced + # compute time gained by only constructing the KDTree once. + kdt_star_p = None + exbv_arr4kdt = None + if len_adx > 0: + num_kdtree_samples = int(min(len_adx, 2e6)) + kdt_idx = np.random.choice(np.arange(len_adx), + size=num_kdtree_samples, + replace=False) + bin_idx = popid_idx[age_idx[kdt_idx]] + star_px = ebf.read_ind(ebf_file, '/px', bin_idx) + star_py = ebf.read_ind(ebf_file, '/py', bin_idx) + star_pz = ebf.read_ind(ebf_file, '/pz', bin_idx) + star_xyz = np.array([star_px, star_py, star_pz]).T + kdt_star_p = cKDTree(star_xyz) + exbv_arr4kdt = ebf.read_ind(ebf_file, '/exbv_schlegel', bin_idx) + del bin_idx, star_px, star_py, star_pz + ########## # Loop through bins of 2 million stars at a time. ########## @@ -746,6 +816,7 @@ def perform_pop_syn(ebf_file, output_root, iso_dir, comp_dict, next_id = _make_comp_dict(iso_dir, age_of_bin, mass_in_bin, stars_in_bin, next_id, + kdt_star_p, exbv_arr4kdt, BH_kick_speed_mean=BH_kick_speed_mean, NS_kick_speed_mean=NS_kick_speed_mean, additional_photometric_systems=additional_photometric_systems, @@ -767,6 +838,7 @@ def perform_pop_syn(ebf_file, output_root, iso_dir, ########## del star_dict gc.collect() + del kdt_star_p, exbv_arr4kdt t1 = time.time() print('perform_pop_syn runtime : {0:f} s'.format(t1 - t0)) @@ -777,7 +849,7 @@ def perform_pop_syn(ebf_file, output_root, iso_dir, binned_counter = 0 hf = h5py.File(output_root + '.h5', 'r') for key in hf: - if 'bin_edges' not in key: + if 'bin_edges' not in key and 'add_pbh' not in key: binned_counter += len(hf[key]) ########## @@ -978,6 +1050,7 @@ def current_initial_ratio(logage, ratio_file, iso_dir, seed=None): def _make_comp_dict(iso_dir, log_age, currentClusterMass, star_dict, next_id, + kdt_star_p, exbv_arr4kdt, BH_kick_speed_mean=50, NS_kick_speed_mean=400, additional_photometric_systems=None, seed=None): @@ -998,8 +1071,16 @@ def _make_comp_dict(iso_dir, log_age, currentClusterMass, star_dict : dictionary The number of entries for each key is the number of stars. - next_id : The next unique ID number (int) that will be assigned to - the new compact objects created. + next_id : int + The next unique ID number (int) that will be assigned to + the new compact objects created. + + kdt_star_p : scipy cKDTree + KDTree constructed from the positions of randomly selected stars + that all share the same popid and similar log_age. + + exbv_arr4kdt : numpy + Array of galactic extinctions for the stars in kdt_star_p Optional Parameters ------------------- @@ -1233,16 +1314,14 @@ def _make_comp_dict(iso_dir, log_age, currentClusterMass, lum_WD_idx = np.argwhere(~np.isnan(comp_table['m_ubv_I'])) if len(lum_WD_idx) > 0: - star_xyz = np.array([star_dict['px'], - star_dict['py'], - star_dict['pz']]).T + # The extinction to the luminous white dwarfs is calculated + # by finding the nearest star in the pop_id / age_bin KDTree + # to the compact object and copying that star's extinction. comp_xyz = np.array([comp_dict['px'][lum_WD_idx], comp_dict['py'][lum_WD_idx], comp_dict['pz'][lum_WD_idx]]).T - kdt = cKDTree(star_xyz) - dist, indices = kdt.query(comp_xyz) - - comp_dict['exbv'][lum_WD_idx] = star_dict['exbv'][indices.T] + dist, indices = kdt_star_p.query(comp_xyz) + comp_dict['exbv'][lum_WD_idx] = exbv_arr4kdt[indices.T] comp_dict['ubv_I'][lum_WD_idx] = comp_table['m_ubv_I'][lum_WD_idx].data comp_dict['ubv_K'][lum_WD_idx] = comp_table['m_ukirt_K'][lum_WD_idx].data @@ -1394,6 +1473,610 @@ def _bin_lb_hdf5(lat_bin_edges, long_bin_edges, obj_arr, output_root): return +############################################################################ +########### Primordial black hole injection and associated functions ####### +############################################################################ + + +def rho_dmhalo(r, rho_0=.0093, r_s=18.6, gamma=1): + """ + Density profile of the dark matter halo. + We are using the parametrization from McMillan (2017) Equation 5, + with defaults based on the mean values in Table 2. + r: galactocentric radius [units: kpc] + rho_0: characteristic density in [units: m_sun / pc**3] + r_s: scale radius in [units: kpc]. + gamma: gamma=1 for NFW, gamma > 1 cuspy, gamma < 1 cored + + returns: density at r [units: m_sun / pc**3] + """ + x = r / r_s + rho = rho_0 / (x ** gamma * (1 + x) ** (3 - gamma)) + return rho + + +def _check_add_pbh(hdf5_file, ebf_file, + fdm, pbh_mass, + r_max, r_s, gamma, v_esc, rho_0, + diagnostic_plots, new_output_root, seed): + """ + Checks that the inputs of add_pbj are valid + + Parameters + ---------- + hdf5_file : str or hdf5 file + str : name of the hdf5 file from the output of perform_pop_syn + + ebf_file : str or ebf file + str : name of the ebf file from Galaxia + ebf file : actually the ebf file from Galaxia + + fdm : float + Fraction of dark matter. + The fraction of dark matter that you want to consist of PBHs. + Defaults to 1. + + pbh_mass : int + The single mass that all PBHs will have (in units of Msun). + Defaults to 40 Msun (from LIGO detections thought to be primordial) + + r_max : float + The maximum radius from the Earth that you want to find PBHs. + Defaults to 16.6 kpc. (2 * distance to galactic center) + + r_s: float + The scale radius of the Milky Way (in kpc). r_s = r_vir / c (virial radius / concentration index) + Defaults to 18.6 kpc. The median value given in McMillan 2017. + + rho_0: float + The initial density that will be used in the NFW profile equations (in units of Msun/pc^3). + Defaults to .0093 [Msun / pc^3]. The median value given in McMillan 2017. + + gamma: float + The inner slope of the MW dark matter halo as described in LaCroix 2018. + Gamma goes into the determination of the velocities and each value returns a slightly different distribution. + The default value is 1, corresponding to an NFW profile. + + v_esc: int + The escape velocity of the Milky Way (in km/s). + v_esc is used in calculating the velocities. + Default is 550 km/s. Most papers cite values of 515-575. + + diagnostic_plots: bool + If set to True, pbh_diagnostic_plots.py is run, and diagnostic plots are saved into a pdf file. + Default False. + + new_output_root : str + If set to None, 'add_pbh' overwrites the original hdf5 file with a + new hdf5 file of the same name. If set to a string, this string is the + prefix of the new hdf5 file. + Default None. + + seed : int + If set to non-None, all random sampling will be seeded with the + specified seed, forcing identical output for PyPopStar and PopSyCLE. + Default None. + """ + if hdf5_file[-3:] != '.h5': + raise Exception('hdf5_file (%s) must be an ebf file.' % str(hdf5_file)) + + if ebf_file[-4:] != '.ebf': + raise Exception('ebf_file (%s) must be an ebf file.' % str(ebf_file)) + + if type(fdm) != int: + if type(fdm) != float: + raise Exception('fdm (%s) must be an integer or a float.' % str(fdm)) + + if type(pbh_mass) != int: + if type(pbh_mass) != float: + raise Exception('pbh_mass (%s) must be an integer or a float.' % str(pbh_mass)) + + if type(r_max) != int: + if type(r_max) != float: + raise Exception('r_max (%s) must be an integer or a float.' % str(r_max)) + + if type(r_s) != int: + if type(r_s) != float: + raise Exception('r_s (%s) must be an integer or a float.' % str(r_s)) + + if gamma not in [.25, .5, 1]: + raise Exception('gamma (%s) must be either .25, .5, 1' % str(gamma)) + + if type(v_esc) != int: + if type(v_esc) != float: + raise Exception('v_esc (%s) must be an integer or a float.' % str(v_esc)) + + if type(rho_0) != int: + if type(rho_0) != float: + raise Exception('rho_0 (%s) must be an integer or a float.' % str(rho_0)) + + if type(diagnostic_plots) != bool: + raise Exception('diagnostic_plots (%s) must be a boolean.' % str(diagnostic_plots)) + + if new_output_root is not None: + if type(new_output_root) != str: + raise Exception('new_output_root (%s) must be None or a string.' % str(new_output_root)) + + if seed is not None: + if type(seed) != int: + raise Exception('seed (%s) must be None or an integer.' % str(seed)) + + +def add_pbh(hdf5_file, ebf_file, fdm=1, pbh_mass=40, + r_max=16.6, r_s=18.6, gamma=1, v_esc=550, + rho_0=0.0093, diagnostic_plots=False, + new_output_root=None, seed=None): + """ + Given some hdf5 file from perform_pop_syn output, creates PBH positions, velocities, etc, + and saves them in a new HDF5 file with the PBHs added. + + Parameters + ---------- + hdf5_file : str or hdf5 file + str : name of the hdf5 file from the output of perform_pop_syn + + ebf_file : str or ebf file + str : name of the ebf file from Galaxia + ebf file : actually the ebf file from Galaxia + + fdm : float + Fraction of dark matter. + The fraction of dark matter that you want to consist of PBHs. + Defaults to 1. + + pbh_mass : int + The single mass that all PBHs will have (in units of Msun). + Defaults to 40 Msun (from LIGO detections thought to be primordial) + + r_max : float + The maximum radius from the Earth that you want to find PBHs. + Defaults to 16.6 kpc. (2 * distance to galactic center) + + r_s: float + The scale radius of the Milky Way (in kpc). r_s = r_vir / c (virial radius / concentration index) + Defaults to 18.6 kpc. The median value given in McMillan 2017. + + gamma: float + The inner slope of the MW dark matter halo as described in LaCroix 2018. + Gamma goes into the determination of the velocities and each value returns a slightly different distribution. + The default value is 1, corresponding to an NFW profile. + + v_esc: int + The escape velocity of the Milky Way (in km/s). + v_esc is used in calculating the velocities. + Default is 550 km/s. Most papers cite values of 515-575. + + rho_0: float + The initial density that will be used in the NFW profile equations (in units of Msun/pc^3). + Defaults to .0093 [Msun / pc^3]. The median value given in McMillan 2017. + + Optional Parameters + ------------------- + diagnostic_plots: bool + If set to True, pbh_diagnostic_plots.py is run, and diagnostic plots are saved into a pdf file. + Default False. + + new_output_root : str + If set to None, 'add_pbh' overwrites the original hdf5 file with a + new hdf5 file of the same name. If set to a string, this string is the + prefix of the new hdf5 file. + Default None. + + seed : int + If set to non-None, all random sampling will be seeded with the + specified seed, forcing identical output for PyPopStar and PopSyCLE. + Default None. + + Outputs + ------- + .h5 : hdf5 file + The new .h5 file with PBHs injected in. + """ + ########## + # Error handling: check whether files exist and + # whether input types are correct. + ########## + + _check_add_pbh(hdf5_file=hdf5_file, ebf_file=ebf_file, + fdm=fdm, pbh_mass=pbh_mass, + r_max=r_max, r_s=r_s, gamma=gamma, v_esc=v_esc, + rho_0=rho_0, diagnostic_plots=diagnostic_plots, + new_output_root=new_output_root, seed=seed) + + if new_output_root is None: + output_root = hdf5_file.replace('.h5', '') + output_hdf5_file = hdf5_file.replace('.h5', '_pbh_tmp.h5') + print('** WARNING **') + print(" 'add_pbh' will overwrite %s, with PBHs appended to each key." % hdf5_file) + print(" To generate a new hdf5 file, rerun 'add_pbh' with the 'new_output_root' argument.") + else: + output_root = new_output_root + output_hdf5_file = '%s.h5' % new_output_root + + ########## + # Start of code + ######### + + # Set random seed + np.random.seed(seed) + + t0 = time.time() + + # Read in the hdf5 file that doesn't have PBHs (product of perform_pop_syn) + no_pbh_hdf5_file = h5py.File(hdf5_file, 'r') + key_list = list(no_pbh_hdf5_file) + # Delete lat_bin_edges, long_bin_edges, and add_pbh from key_list + key_list = [key for key in key_list if 'bin_edges' not in key] + key_list = [key for key in key_list if 'add_pbh' not in key] + + # Get data from lat_bin_edges and long_bin_edges + lat_bin = no_pbh_hdf5_file['lat_bin_edges'][:] + long_bin = no_pbh_hdf5_file['long_bin_edges'][:] + bin_edges_number = len(long_bin) + + # Get the maximum ID from all of the stars and compact objects + # (Later used to set the IDs of the PBHs) + max_id_no_pbh = [] + for key in key_list: + if len(no_pbh_hdf5_file[key]) == 0: + continue + max_id_no_pbh.append(np.max(no_pbh_hdf5_file[key]['obj_id'])) + max_id = np.amax(max_id_no_pbh) + + hdf5_dset_names = no_pbh_hdf5_file[key_list[0]][:].dtype.names + + no_pbh_hdf5_file.close() + + # Read in ebf file + t = ebf.read_ind(ebf_file, '/log', 0) + # Convert log to useful dictionary + ebf_log = make_ebf_log(t) + + # Obtain survey area, center latitude, and center longitude + b = float(ebf_log['latitude']) # deg + b_radian = b * np.pi / 180 # rad + l = float(ebf_log['longitude']) # deg + l_radian = l * np.pi / 180 # rad + surveyArea = float(ebf_log['surveyArea']) # deg^2 + + # Calculate the field of view for the current field + field_of_view_radius = (surveyArea / np.pi) ** (1 / 2) + + # Generate an array of heliocentric radii + # (Just used to numerically integrate the line-of-sight density) + n_lin = 100000 + r_h_linspace = np.linspace(0, r_max, num=n_lin) + + # Represent the line-of-sight line as galactic coordinates + galactic_lin = coord.Galactic(l=l_radian * units.rad, + b=b_radian * units.rad, + distance=r_h_linspace * units.kpc) + + # Transform the line-of-sight into to galactocentric coordinates + # (Outputs l, b, and distance [units: deg, deg, kpc]) + galacto_lin = galactic_lin.transform_to( + coord.Galactocentric(representation_type='spherical')) + + # Determine dark matter density at all galactocentric radii along the line-of-sight + rho_lin = rho_dmhalo(galacto_lin.spherical.distance.value, + rho_0=rho_0, r_s=r_s, gamma=gamma) + + # Estimate the total mass within the line-of-sight cylinder [units: M_sun kpc**-2] + # Total mass = projected line-of-sight density x projected line-of-sight area + rho_marg_r = np.trapz(rho_lin, dx=(r_max) / n_lin) * 1000 ** 3 + print("Projected density along line-of-sight = {0:0.2e} [M_sun kpc**-2]".format( + rho_marg_r)) + # Determine line-of-sight cylinder radius, assuming small angle approximation [units: kpc] + r_proj_los_cyl = field_of_view_radius * np.pi / 180 * (r_max) + # Get projected area of the LOS cylinder [units: kpc**2] + area_proj_los_cyl = np.pi * r_proj_los_cyl ** 2 + # Get the total mass within the line-of-sight cylinder + mass_los_cyl = rho_marg_r * area_proj_los_cyl + print("Mass within line-of-sight cylinder = {0:0.2e} [M_sun]".format( + mass_los_cyl)) + + # Get the total number of black holes to randomly draw + n_pbh = int(np.round(fdm * mass_los_cyl / pbh_mass)) + + # Estimate the discrete CDF based on the discrete PDF + # (Switch into galactic coordinates) + rho_marg_r_cum = integrate.cumtrapz(y=rho_lin, + x=galactic_lin.distance.kpc, + dx=(r_max) / n_lin) + cdf_los = rho_marg_r_cum / rho_marg_r_cum[-1] + # Since cumtrapz does not include zero for the first element, insert it + cdf_los = np.insert(cdf_los, 0, 0) + + # Create a function to interpolate the CDF so that we can randomly sample from it + f_cdf_d = interpolate.interp1d(cdf_los, galactic_lin.distance.kpc) + + # Randomly sample galactic coordinates for the PBHs based on the CDF + d_galactic = f_cdf_d(np.random.uniform(size=n_pbh)) + + # Randomly assign an l & b galactic coordinate to each PBH within the + # line-of-sight cone (sample the angle from 0 to 2pi [units: radians]) + theta = np.random.uniform(size=n_pbh) * 2 * np.pi + # Make uniform by sampling a radius and correcting for annular area + r_cyl = r_proj_los_cyl * np.sqrt(np.random.uniform(size=n_pbh)) # kpc + y_cyl = r_cyl * np.sin(theta) # kpc + x_cyl = r_cyl * np.cos(theta) # kpc + + # Create mask for masking out PBHs falling outside of the observation cone + mask_obs_cone = r_cyl <= r_proj_los_cyl * d_galactic / (r_max) + print( + 'Number of PBH before and after light cone masking: {0} and {1}, respectively'.format( + n_pbh, np.sum(mask_obs_cone))) + + # Get galactic l and b, assuming the small angle approximation + b_galactic = r_cyl * np.sin(theta) / d_galactic + b_radian # rad + l_galactic = r_cyl * np.cos(theta) / np.cos(b_radian) / d_galactic + l_radian # rad + + if diagnostic_plots: + print('Saving diagnostic plots') + from popsycle.add_pbh_plots import print_plots + print_plots(output_root=output_root, galactic_lin_distance=galactic_lin.distance.kpc, + galactic_lin_b=galactic_lin.b.deg, galactic_lin_l=galactic_lin.l.deg, + galactocen_lin_spherical_distance=galacto_lin.spherical.distance.kpc, + galactocen_lin_spherical_b=galacto_lin.spherical.lat.deg, + galactocen_lin_spherical_l=galacto_lin.spherical.lon.deg, rho_lin=rho_lin, + r_max=r_max, n_lin=n_lin, cdf_los=cdf_los, + x_cyl=x_cyl, y_cyl=y_cyl, r_cyl=r_cyl, r_proj_los_cyl=r_proj_los_cyl, n_pbh=n_pbh, + d_galac=d_galactic, b_galac=b_galactic, l_galac=l_galactic, area_proj_los_cyl=area_proj_los_cyl, + mask_obs_cone=mask_obs_cone, field_of_view_radius=field_of_view_radius, l_radian=l_radian, + b_radian=b_radian, f_cdf_d=f_cdf_d, pbh_mass=pbh_mass) + + # Mask out PBHs outside of the observation cone + d_galactic = d_galactic[mask_obs_cone] + b_galactic = b_galactic[mask_obs_cone] + l_galactic = l_galactic[mask_obs_cone] + + latitude = b_galactic * (180 / np.pi) # degrees + longitude = l_galactic * (180 / np.pi) # degrees + + # Determine how many PBHs are in the current field + N_PBHs_in_field = len(d_galactic) + print('%i PBHs in the field' % N_PBHs_in_field) + + if N_PBHs_in_field == 0: + print('-- No PBHs in the field') + print('-- Copying %s to %s' % (hdf5_file, output_hdf5_file)) + shutil.copy(hdf5_file, output_hdf5_file) + return + + # Convert PBH positions back to galactocentric for determining velocities + galactic_pbh = coord.Galactic(l=longitude * units.deg, + b=latitude * units.deg, + distance=d_galactic * units.kpc) + galacto_pbh = galactic_pbh.transform_to( + coord.Galactocentric(representation_type='cartesian')) + # Get the radial galactocentric PBH distances + pbh_r_galacto = galacto_pbh.spherical.distance.value + + # Inner slope of the MW halo: + # From Lacroix et al 2018, Figure 11 (top left panel) + data_dir = '%s/data' % os.path.dirname(inspect.getfile(add_pbh)) + if gamma == 1: + vel_data = np.genfromtxt('%s/radial_velocity_profile_steep.csv' % data_dir, + names=True, delimiter=',') + elif gamma == .25: + vel_data = np.genfromtxt('%s/radial_velocity_profile_shallow.csv' % data_dir, + names=True, delimiter=',') + elif gamma == .5: + vel_data = np.genfromtxt('%s/radial_velocity_profile_middle.csv' % data_dir, + names=True, delimiter=',') + else: + raise Exception('gamma (%s) must be either .25, .5, 1' % str(gamma)) + + # Interpolate v values from the above data, given PBH radii + pbh_vrms = np.interp(pbh_r_galacto, vel_data['r'], vel_data['v']) + v_vals = np.arange(0, v_esc) + a = (1 / 2) * pbh_vrms * ((np.pi / 2) ** (1 / 2)) + + # Calculate RMS velocities for each PBH by randomly sampling the Maxwellian CDF + rms_velocities = [] + + for a_val in a: + cdf = special.erf(v_vals / (a_val * 2 ** (1 / 2))) - \ + (((2 / np.pi) ** (1 / 2)) * ((v_vals * np.exp(-v_vals ** 2 / (2 * a_val ** 2)) / a_val))) + rand = np.random.uniform(0,1) + rms_velocities.append(np.interp(rand, cdf, v_vals)) + + # Complete PBH spherical velocities by sampling a latitude and longitude + # (arcsin is used so that PBHs are placed in more realistic locations) + sin_lat_vel = np.random.uniform(-1, 1, len(d_galactic)) + lat_vel = np.arcsin(sin_lat_vel) + long_vel = np.random.uniform(0, 2 * np.pi, len(d_galactic)) + + # Transform velocities to cartesian to get vx, vy, and vz + cart_vel = coord.spherical_to_cartesian(rms_velocities, lat_vel, long_vel) + + # Load up a numpy array + comp_dtype = _generate_comp_dtype(hdf5_dset_names) + pbh_data = np.empty(len(d_galactic), dtype=comp_dtype) + + # Transform latitude and longitude into the right format for PopSyCLE + longitude = np.where(longitude > 180, longitude - 360, longitude) + + # Set the galactic coordinates + pbh_data['rad'] = d_galactic + pbh_data['glon'] = longitude + pbh_data['glat'] = latitude + + # Set heliocentric velocities equal to galactocentric velocities + pbh_data['vx'] = cart_vel[0] + pbh_data['vy'] = cart_vel[1] + pbh_data['vz'] = cart_vel[2] + + # Get the rest of the PBH data for the combined .h5 file + pbh_data['mass'] = np.full(len(d_galactic), pbh_mass) + pbh_data['zams_mass'] = np.full(len(d_galactic), pbh_mass) + pbh_data['age'] = np.full(len(d_galactic), np.nan) + pbh_data['popid'] = np.full(len(d_galactic), 10) + pbh_data['rem_id'] = np.full(len(d_galactic), 104) + + # Set the heliocentric PBH positions + cart_helio = coord.spherical_to_cartesian(d_galactic, b_galactic, l_galactic) + pbh_data['px'] = cart_helio[0] + pbh_data['py'] = cart_helio[1] + pbh_data['pz'] = cart_helio[2] + + # Set the galactic velocity and proper motions + vr, mu_b, mu_lcosb = calc_sph_motion(cart_vel[0], + cart_vel[1], + cart_vel[2], + d_galactic, b_galactic, l_galactic) + pbh_data['vr'] = vr + pbh_data['mu_b'] = mu_b + pbh_data['mu_lcosb'] = mu_lcosb + + pbh_data['obj_id'] = np.arange((max_id + 1), (max_id + len(d_galactic) + 1)) + + pbh_data['exbv'] = np.full(len(d_galactic), np.nan) + pbh_data['ubv_K'] = np.full(len(d_galactic), np.nan) + pbh_data['ubv_J'] = np.full(len(d_galactic), np.nan) + pbh_data['ubv_I'] = np.full(len(d_galactic), np.nan) + pbh_data['ubv_U'] = np.full(len(d_galactic), np.nan) + pbh_data['ubv_R'] = np.full(len(d_galactic), np.nan) + pbh_data['ubv_B'] = np.full(len(d_galactic), np.nan) + pbh_data['ubv_H'] = np.full(len(d_galactic), np.nan) + pbh_data['ubv_V'] = np.full(len(d_galactic), np.nan) + + # Note: The first four elements may not be in the compound datatype + # due to legacy files. The ztf filters will only be present if + # the hdf5 file was created with additional_photometric_systems = ['ztf'] + for field in ['teff', 'grav', 'mbol', 'feh', + 'ztf_g', 'ztf_r', 'ztf_i']: + if field in hdf5_dset_names: + pbh_data[field] = np.full(len(d_galactic), np.nan) + + # Calculate the maximum and minimum l and b values for each dataset in the + # file with no PBHs, to determine the datasets to correctly append the PBHs to + lat_long_list = [] + for key in key_list: + idx_l = int(key.split('b')[0].replace('l', '')) + max_l = long_bin[idx_l + 1] + min_l = long_bin[idx_l] + + idx_b = int(key.split('b')[1]) + max_b = lat_bin[idx_b + 1] + min_b = lat_bin[idx_b] + + lat_long_list.append((min_l, max_l, min_b, max_b)) + + # Open the file with no PBHs and create a new file for data with PBHs added + no_pbh_hdf5_file = h5py.File(hdf5_file, 'r') + pbh_hdf5_file = h5py.File(output_hdf5_file, 'w') + pbh_hdf5_file['add_pbh'] = True + + # Append the PBH data to the data with no PBHs, and write to an .h5 file + N_objs_no_pbh = 0 + N_objs_pbh = 0 + N_pbhs_removed = 0 + for idx, key in enumerate(key_list): + key_data = no_pbh_hdf5_file[key][:] + + # Remove any PBHs that are already in the h5 file + cond = key_data['rem_id'] != 104 + key_data = key_data[cond] + N_pbhs_removed += np.sum(~cond) + + # Count the number of objects in the key + N_objs_no_pbh += key_data.shape[0] + + # Build a mask on the PBHs that fits the bounds of the key + min_l, max_l, min_b, max_b = lat_long_list[idx] + mask = (pbh_data['glon'] >= min_l) & \ + (pbh_data['glon'] <= max_l) & \ + (pbh_data['glat'] >= min_b) & \ + (pbh_data['glat'] <= max_b) + + # If there are no PBHs in the key, copy over the original key + if np.sum(mask) == 0: + combined_data = key_data + # If there are PBHs in the key, append them to the original key + else: + pbh_data_in_key = pbh_data[mask] + combined_data = np.hstack((key_data, pbh_data_in_key)) + N_objs_pbh += combined_data.shape[0] + _ = pbh_hdf5_file.create_dataset(key, + shape=(combined_data.shape[0],), + dtype=comp_dtype, + data=combined_data) + _ = pbh_hdf5_file.create_dataset('lat_bin_edges', (len(lat_bin), 1), + data=lat_bin) + _ = pbh_hdf5_file.create_dataset('long_bin_edges', (len(lat_bin), 1), + data=long_bin) + no_pbh_hdf5_file.close() + pbh_hdf5_file.close() + + # If 'new_output_root' is None, replace temporary file with original + if new_output_root is None: + os.remove(hdf5_file) + os.rename(output_hdf5_file, hdf5_file) + + t1 = time.time() + + ########## + # Make log file + ########## + now = datetime.datetime.now() + popsycle_path = os.path.dirname(inspect.getfile(perform_pop_syn)) + popsycle_hash = subprocess.check_output(['git', 'rev-parse', 'HEAD'], + cwd=popsycle_path).decode('ascii').strip() + dash_line = '-----------------------------' + '\n' + empty_line = '\n' + + line0 = 'FUNCTION INPUT PARAMETERS' + '\n' + line1 = 'hdf5_file , ' + hdf5_file + '\n' + line2 = 'ebf_file , ' + ebf_file + '\n' + line3 = 'fdm , ' + str(fdm) + '\n' + line4 = 'pbh_mass , ' + str(pbh_mass) + ' , (Msun)' + '\n' + line5 = 'r_max , ' + str(r_max) + ' , (kpc)' + '\n' + line6 = 'r_s , ' + str(r_s) + ' , (kpc)' + '\n' + line7 = 'gamma , ' + str(gamma) + '\n' + line8 = 'v_esc , ' + str(v_esc) + ' , (km/s)' + '\n' + line9 = 'rho_0 , ' + str(rho_0) + ' , (Msun / pc^3)' + '\n' + line10 = 'diagnostic_plots , ' + str(diagnostic_plots) + '\n' + line11 = 'new_output_root , ' + str(new_output_root) + '\n' + line12 = 'seed , ' + str(seed) + '\n' + + line13 = 'VERSION INFORMATION' + '\n' + line14 = str(now) + ' : creation date' + '\n' + line15 = popsycle_hash + ' : PopSyCLE commit' + '\n' + + line16 = 'OTHER INFORMATION' + '\n' + line17 = str(t1 - t0) + ' : total runtime (s)' + '\n' + line18 = str(N_objs_no_pbh) + ' : original objects' + '\n' + line19 = str(N_PBHs_in_field) + ' : PBHs in the field' + '\n' + line20 = str(N_objs_pbh) + ' : new total objects' + '\n' + + line21 = 'FILES CREATED' + '\n' + line22 = output_hdf5_file + ' : HDF5 file' + '\n' + + with open(output_root + '_add_pbh.log', 'w') as out: + out.writelines([line0, dash_line, line1, line2, line3, line4, line5, + line6, line7, line8, line9, line10, line11, line12, + empty_line, line13, dash_line, line14, line15, + empty_line, line16, dash_line, line17, line18, line19, + line20, empty_line, line21, dash_line, line22]) + + ########## + # Informative print statements. + ########## + print('- %i PBHS removed from %s' % (N_pbhs_removed, hdf5_file)) + print('--- %i original objects' % N_objs_no_pbh) + print('--- %i PBHs in the field' % N_PBHs_in_field) + print('--- %i new total objects' % N_objs_pbh) + + if N_objs_pbh == N_objs_no_pbh + N_PBHs_in_field: + print('Binned PBHs equals total PBHs') + else: + print('** MISSING PBHs!! **') + + print('add_pbh runtime: {0:f} s'.format(t1 - t0)) + return + + ############################################################################ ########### Candidate event calculation and associated functions ########### ############################################################################ @@ -1608,31 +2291,30 @@ def calc_events(hdf5_file, output_root2, if results[ii][1] is not None: results_bl.append(results[ii][1]) - if len(results_ev) == 0: - print('No events!') - return - else: + if len(results_ev) > 0: events_tmp = np.concatenate(results_ev, axis=0) if len(results_bl) == 0: blends_tmp = np.array([]) else: blends_tmp = np.concatenate(results_bl, axis=0) - # Convert the events numpy recarray into an - # Astropy Table for easier consumption. - events_tmp = unique_events(events_tmp) - events_final = Table(events_tmp) - N_events = len(events_final) - print('Candidate events detected: ', N_events) - - if len(results_bl) != 0: - blends_tmp = unique_blends(blends_tmp) - blends_final = Table(blends_tmp) + # Convert the events numpy recarray into an + # Astropy Table for easier consumption. + events_tmp = unique_events(events_tmp) + events_final = Table(events_tmp) + N_events = len(events_final) + print('Candidate events detected: ', N_events) - # Save out file - events_final.write(output_root2 + '_events.fits', overwrite=overwrite) - blends_final.write(output_root2 + '_blends.fits', overwrite=overwrite) + if len(results_bl) != 0: + blends_tmp = unique_blends(blends_tmp) + blends_final = Table(blends_tmp) + # Save out file + events_final.write(output_root2 + '_events.fits', overwrite=overwrite) + blends_final.write(output_root2 + '_blends.fits', overwrite=overwrite) + else: + N_events = 0 + print('No events!') t1 = time.time() ########## @@ -1642,7 +2324,7 @@ def calc_events(hdf5_file, output_root2, radius_cut = radius_cut / 1000.0 # back to arcsec popsycle_path = os.path.dirname(inspect.getfile(perform_pop_syn)) popsycle_hash = subprocess.check_output(['git', 'rev-parse', 'HEAD'], - cwd=popsycle_path).decode('ascii').strip() + cwd=popsycle_path).decode('ascii').strip() dash_line = '-----------------------------' + '\n' empty_line = '\n' line0 = 'FUNCTION INPUT PARAMETERS' + '\n' @@ -1654,6 +2336,7 @@ def calc_events(hdf5_file, output_root2, line6 = 'theta_frac , ' + str(theta_frac) + ' , (thetaE)' + '\n' line7 = 'blend_rad , ' + str(blend_rad) + ' , (arcsec)' + '\n' line8 = 'n_proc , ' + str(n_proc) + '\n' + line9 = 'VERSION INFORMATION' + '\n' line10 = str(now) + ' : creation date' + '\n' line11 = popsycle_hash + ' : PopSyCLE commit' + '\n' @@ -1662,9 +2345,14 @@ def calc_events(hdf5_file, output_root2, line13 = str(t1 - t0) + ' : total runtime (s)' + '\n' line14 = str(N_events) + ' : total number of events' + '\n' - line15 = 'FILES CREATED' + '\n' - line16 = output_root2 + '_events.fits : events file' + '\n' - line17 = output_root2 + '_blends.fits : blends file' + '\n' + if N_events > 0: + line15 = 'FILES CREATED' + '\n' + line16 = output_root2 + '_events.fits : events file' + '\n' + line17 = output_root2 + '_blends.fits : blends file' + '\n' + else: + line15 = 'NO FILES CREATED' + '\n' + line16 = '\n' + line17 = '\n' with open(output_root2 + '_calc_events.log', 'w') as out: out.writelines([line0, dash_line, line1, line2, line3, @@ -2492,30 +3180,30 @@ def refine_events(input_root, filter_name, photometric_system, red_law, event_tab['u0'] = u0 if len(event_tab) == 0: print('No events!') - return - - ########## - # Calculate apparent magnitudes - ########## + output_file = 'NO FILE CREATED' + else: + ########## + # Calculate apparent magnitudes + ########## - # Einstein crossing time - # THIS HAS TO GO BEFORE _CALC_OBSERVABLES - t_E = event_tab['theta_E'] / event_tab['mu_rel'] # yr - t_E *= 365.25 # days - event_tab['t_E'] = t_E # days + # Einstein crossing time + # THIS HAS TO GO BEFORE _CALC_OBSERVABLES + t_E = event_tab['theta_E'] / event_tab['mu_rel'] # yr + t_E *= 365.25 # days + event_tab['t_E'] = t_E # days - # Add stuff to event_tab... shouldn't have any direct outputs - _calc_observables(filter_name, red_law, event_tab, blend_tab, photometric_system) + # Add stuff to event_tab... shouldn't have any direct outputs + _calc_observables(filter_name, red_law, event_tab, blend_tab, photometric_system) - # Relative parallax - pi_rel = event_tab['rad_L'] ** -1 - event_tab['rad_S'] ** -1 - event_tab['pi_rel'] = pi_rel # mas + # Relative parallax + pi_rel = event_tab['rad_L'] ** -1 - event_tab['rad_S'] ** -1 + event_tab['pi_rel'] = pi_rel # mas - # Microlensing parallax - pi_E = pi_rel / event_tab['theta_E'] - event_tab['pi_E'] = pi_E # dim'less + # Microlensing parallax + pi_E = pi_rel / event_tab['theta_E'] + event_tab['pi_E'] = pi_E # dim'less - event_tab.write(output_file, overwrite=overwrite) + event_tab.write(output_file, overwrite=overwrite) t_1 = time.time() @@ -2911,6 +3599,7 @@ def make_label_file(h5file_name, overwrite=False): return + ########################################################################### ############ General formulas, conversions, and calculations ############## ########################################################################### @@ -3423,3 +4112,5 @@ def calc_f(lambda_eff): f = L * (B - V) ** -1 return f + +