Commit ce0cb480 by Uddeepta Deka

updated waveform generation script

parent 6401add5
This source diff could not be displayed because it is too large. You can view the blob instead.
#!/home1/uddeepta.deka/softwares/miniconda3/envs/igwn-py39/bin/python
'''
Use this script to generate parameters to run
lensed_waveform_exec.py on a cluster for different parameter values.
'''
import numpy as np
# Define your range for y_value and Ml_value
y_values = [0.1, 0.5, 1.0]
q_values = [1, 1e-3, 1e-6, 1e-9, 1e-11, 0]
# Write to a parameter file
with open('params.txt', 'w') as f:
for y in y_values:
for qeff in q_values:
f.write(f"{y} {qeff}\n")
# We use this script to generate a .sub file for ht_condor
import os
import numpy as np
## The following config. is used to generate F(t), F(w), h_U, h_L sample waveforms
# y_arr = [0.5, 1.0, 2.0]
# q_arr = [-0.5, 0.1, 0.5]
# t_max = 200
## The following config. is used to generate data for Q=0 mismatch checks
y_arr = np.geomspace(0.1, 3, 50)
Ml_arr = np.geomspace(1, 1e5, 50)
t_max = 200
# y_arr = np.logspace(-1, np.log10(2.0), 3000)
dir_path = "/home1/uddeepta.deka/lensing/charged_lens_updated/data/condor_output/"
if not os.path.exists(dir_path+"outfiles/"):
os.makedirs(dir_path+"outfiles/")
if not os.path.exists(dir_path+"logs/"):
os.makedirs(dir_path+"logs/")
if not os.path.exists(dir_path+"errors/"):
os.makedirs(dir_path+"errors/")
if not os.path.exists(dir_path+"results/"):
os.makedirs(dir_path+"results/")
condor_filename = "condor_sub.sub"
vars = {
"Universe": "vanilla",
"Executable": "/home1/uddeepta.deka/lensing/charged_lens_updated/scripts/lensed_waveform_exec.py",
"Getenv": True,
"Output": dir_path+"outfiles/output.out.$(Cluster).$(Process).txt",
"Error": dir_path+"errors/error.out.$(Cluster).$(Process).txt",
"Log": dir_path+"logs/log.log.$(Cluster).$(Process).txt",
"notify_user": "uddeepta.deka@icts.res.in",
"Notification": "Never",
"requirements" : '(Machine != "node72.alice.icts.res.in")',
"request_cpus": 4,
"request_memory": 1024,
}
with open(dir_path + condor_filename, 'w') as fp:
for key, val in vars.items():
fp.write(f"{key}\t = {val}\n")
for y in y_arr:
# for q in q_arr:
for Ml in Ml_arr:
fp.write(f"\narguments = y {y} Ml {Ml} t_max {t_max} folder_path {dir_path}results/\nQueue\n")
\ No newline at end of file
......@@ -20,7 +20,12 @@ import lal
import pycbc.psd as psd
from pycbc import detector
from pycbc.waveform import get_fd_waveform
import matplotlib.pyplot as plt
from plotsettings import *
from cycler import cycler
style = 'Notebook'
plotsettings(style)
plt.rcParams['axes.prop_cycle'] = cycler(color=dark2)
#===========================================================================
#============================ lensed waveform ==============================
......@@ -210,6 +215,84 @@ class MicrolensedWaveform:
hh_l_recons = self.f_plus * lensed_hp_recons + self.f_cross * lensed_hc_recons
return hh_l, hh_l_recons
#===========================================================================
#================================ Plotting =================================
#===========================================================================
def make_plots(wf_obj, outdir=None):
max_td = np.max(wf_obj.lens.img_td)
x_max = 1.2 * np.sqrt(2 * max_td)
ngrid = 1000
x1_ = np.linspace(-x_max + wf_obj.lens.y1, x_max + wf_obj.lens.y1, ngrid, endpoint=True)
x2_ = np.linspace(-x_max + wf_obj.lens.y2, x_max + wf_obj.lens.y2, ngrid, endpoint=True)
X1, X2 = np.meshgrid(x1_, x2_)
timedelay_map = wf_obj.lens.timeDelay(X1, X2)
gradtd_map_X1, gradtd_map_X2 = wf_obj.lens.gradTimeDelay(X1, X2)
gradtd_map_norm = np.sqrt(gradtd_map_X1 ** 2 + gradtd_map_X2 ** 2)
fig, ax = plt.subplots(figsize=(figWidthsOneColDict[style],figHeightsDict[style]))
cs1 = ax.contourf(X1, X2, timedelay_map, 20, cmap='RdPu')
cs2 = ax.contour(X1, X2, gradtd_map_norm, levels=[0.0, 0.01, 0.1], colors=['r', 'r', 'r'])
ax.scatter(wf_obj.lens.y1, wf_obj.lens.y2, marker='^', label='source')
cbar = fig.colorbar(cs1, ax=ax)
ax.set_xlabel(r'$x_1$')
ax.set_ylabel(r'$x_2$')
for img in wf_obj.lens.img:
ax.scatter(img, 0, c='k', marker='+', s=30)
ax.set_xlim(np.min(wf_obj.lens.img)-0.5, np.max(wf_obj.lens.img)+0.5)
ax.set_ylim(-0.5, +0.5)
ax.plot([], [],ls='', marker='+', c='k', label='image')
ax.legend(loc='best', frameon=True)
plt.savefig(outdir+"verif_plot.pdf", format='pdf', bbox_inches='tight')
plt.close()
plt.figure(figsize=(figWidthsOneColDict[style], figHeightsDict[style]))
plt.plot(wf_obj.t, wf_obj.Ft, label='vanilla')
plt.plot(wf_obj.t, wf_obj.Ft_recons, label='reconstructed')
plt.axhline(1, ls='--', c='magenta')
image_time_delays = wf_obj.lens.img_td - np.min(wf_obj.lens.img_td)
for td in image_time_delays:
plt.axvline(td, ls='-', c='grey')
plt.plot([], [], ls='-', c='grey', label='images')
plt.xlabel(r'$t$')
plt.ylabel(r'$\widetilde{F}(t)$')
plt.xscale('log')
plt.legend()
plt.xlim(left=1e-3, right=wf_obj.t_max)
plt.savefig(outdir+"Ft_plot.pdf", format='pdf', bbox_inches='tight')
plt.close()
plt.figure(figsize=(figWidthsTwoColDict[style]*1.5, figHeightsDict[style]))
plt.subplot(131)
plt.plot(wf_obj.f_UL, np.abs(wf_obj.Ff), label='vanilla')
plt.plot(wf_obj.f_UL, np.abs(wf_obj.Ff_recons), label='reconstructed')
plt.xlabel(r'$f$[Hz]')
plt.ylabel(r'$|F(f)|$')
plt.xscale('log')
plt.legend()
plt.subplot(132)
plt.plot(wf_obj.f_UL, np.angle(wf_obj.Ff), label='vanilla')
plt.plot(wf_obj.f_UL, np.angle(wf_obj.Ff_recons), label='reconstructed')
plt.xlabel(r'$f$[Hz]')
plt.ylabel(r'arg($F(f)$)')
plt.xscale('log')
plt.legend()
plt.subplot(133)
plt.plot(wf_obj.f_UL, np.abs(wf_obj.hh_UL), lw=2, alpha=0.6, c='grey', label='unlensed')
plt.plot(wf_obj.f_UL, np.abs(wf_obj.hh_L), label='lensed')
plt.plot(wf_obj.f_UL, np.abs(wf_obj.hh_L_recons), label='lensed reconstructed')
plt.xlabel(r'$f$[Hz]')
plt.ylabel(r'$|h(f)|$')
plt.xlim(left=wf_obj.f_lo, right=min(max(wf_obj.f_UL), max(wf_obj.f_Sh)))
plt.xscale('log')
plt.yscale('log')
plt.legend()
plt.tight_layout()
plt.savefig(outdir+"strain_plot.pdf", format='pdf', bbox_inches='tight')
plt.close()
return
#===========================================================================
#=============================== Driver code ===============================
#===========================================================================
......@@ -234,7 +317,7 @@ def main():
"f_hi":1024.0, "zs":2.0, "m1":20.0, "m2":20.0, "chi1":0.0, "chi2":0.0,
"iota":np.pi/8., "phi0":np.pi/4., "ra":1.7, "dec":0.6, "pol":4.8,
"approx":"IMRPhenomXP", "delta_f":1./32, "t_max":30., "dt_fac":7.5,}
folder_path = os.getcwd() + "/lensed_waveform_data/"
folder_path = "/home1/uddeepta.deka/lensing/charged_lens_updated/data/lensed_waveform_data/"
for i in range(len(sys.argv)):
if sys.argv[i] == "-h":
......@@ -266,6 +349,7 @@ def main():
np.savetxt(f, np.array([mw.t, mw.Ft, mw.Ft_recons]).T, header=header)
np.savez(waveform_datafile, w=mw.w, f=mw.f_UL, Ff=mw.Ff, Ff_recons=mw.Ff_recons,
h_ul=mw.hh_UL, h_l=mw.hh_L, h_l_recons=mw.hh_L_recons)
make_plots(mw, outdir=outdir)
tEnd = time.perf_counter()
print("Data stored in : ", outdir)
print("Time taken to generate data : %.2f s."%(tEnd-tStart))
......
Universe = vanilla
Executable = /home1/uddeepta.deka/lensing/charged_lens_updated/scripts/lensed_waveform_exec.py
Arguments = y $(y_value) q $(q_value) t_max 200
Getenv = True
Output = /home1/uddeepta.deka/lensing/charged_lens_updated/data/lensed_waveform_data/outfiles/output_$(y_value)_$(q_value).out
Error = /home1/uddeepta.deka/lensing/charged_lens_updated/data/lensed_waveform_data/errors/error_$(y_value)_$(q_value).err
Log = /home1/uddeepta.deka/lensing/charged_lens_updated/data/lensed_waveform_data/logs/log_$(y_value)_$(q_value).log
notify_user = uddeepta.deka@icts.res.in
Notification = Never
requirements = (Machine != "node72.alice.icts.res.in")
request_cpus = 4
request_memory = 1024
queue y_value q_value from params.txt
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or sign in to comment