Example: Fitting your own waveform¶
In [1]:
Copied!
from jaxqualin.waveforms import delayed_QNM, waveform
from jaxqualin.qnmode import mode, mode_list
from jaxqualin.fit import QNMFitVaryingStartingTime
from jaxqualin.selection import ModeSearchAllFreeVaryingN
from jaxqualin.plot import (plot_amplitudes, plot_phases,
plot_omega_free, plot_predicted_qnms,
plot_mode_searcher_results)
import numpy as np
import matplotlib.pyplot as plt
from jaxqualin.waveforms import delayed_QNM, waveform
from jaxqualin.qnmode import mode, mode_list
from jaxqualin.fit import QNMFitVaryingStartingTime
from jaxqualin.selection import ModeSearchAllFreeVaryingN
from jaxqualin.plot import (plot_amplitudes, plot_phases,
plot_omega_free, plot_predicted_qnms,
plot_mode_searcher_results)
import numpy as np
import matplotlib.pyplot as plt
PyCBC.libutils: pkg-config call failed, setting NO_PKGCONFIG=1
Make waveform¶
As example, we will make a waveform that contains three QNMs that are distorted close to the merger peak. You can replace this with whatever time domain waveform you would like to fit.
In [2]:
Copied!
Mf = 1
af = 0.7
modes = mode_list(['2.2.0', '2.2.1', '3.2.0'], Mf, af)
A_phi_dict = {'2.2.0': dict(A = 1., phi = 0.),
'2.2.1': dict(A = 3., phi = np.pi/2),
'3.2.0': dict(A = 1e-2, phi = np.pi)}
t_arr = np.linspace(0, 120, 1000)
h_arr = np.zeros(t_arr.shape, dtype = np.complex128)
for i, mode in enumerate(modes):
if i == 0:
A_delay = 0
A_sig = 10
phi_sig = 5
else:
A_delay = 5
A_sig = 2
phi_sig = 2
h_arr += delayed_QNM(mode, t_arr,
A_phi_dict[mode.string()]['A'],
A_phi_dict[mode.string()]['phi'],
A_delay = A_delay, A_sig = A_sig, phi_sig = phi_sig)
Mf = 1
af = 0.7
modes = mode_list(['2.2.0', '2.2.1', '3.2.0'], Mf, af)
A_phi_dict = {'2.2.0': dict(A = 1., phi = 0.),
'2.2.1': dict(A = 3., phi = np.pi/2),
'3.2.0': dict(A = 1e-2, phi = np.pi)}
t_arr = np.linspace(0, 120, 1000)
h_arr = np.zeros(t_arr.shape, dtype = np.complex128)
for i, mode in enumerate(modes):
if i == 0:
A_delay = 0
A_sig = 10
phi_sig = 5
else:
A_delay = 5
A_sig = 2
phi_sig = 2
h_arr += delayed_QNM(mode, t_arr,
A_phi_dict[mode.string()]['A'],
A_phi_dict[mode.string()]['phi'],
A_delay = A_delay, A_sig = A_sig, phi_sig = phi_sig)
Make waveform object¶
You can pass in whatever time domain waveform.
If you know that the peak strain is located at t = 0, then pass in t_peak = 0
, or else the waveform object detects the peak on its own.
If t_peak
is not passed, the first remove_num = 500
data points are removed by default because BBH merger simulations often contain junk radiation in the beginning.
You can set remove_num = 0
if you do not want to remove any data points.
In [3]:
Copied!
h = waveform(t_arr, h_arr, t_peak = 0)
h = waveform(t_arr, h_arr, t_peak = 0)
In [4]:
Copied!
fig, ax = plt.subplots()
ax.semilogy(h.time, np.abs(h.hr))
ax.set_xlabel(r'$t$')
ax.set_ylabel(r'$|h_r|$')
fig, ax = plt.subplots()
ax.semilogy(h.time, np.abs(h.hr))
ax.set_xlabel(r'$t$')
ax.set_ylabel(r'$|h_r|$')
Out[4]:
Text(0, 0.5, '$|h_r|$')
Free QNMs (unfixed frequencies)¶
In [5]:
Copied!
t0_arr = np.linspace(0, 50, num = 101) # array of starting times to fit for
# t0 = 0 is the peak of the strain
qnm_fixed_list = [] # list of QNMs with fixed frequencies in the fit model
run_string_prefix = f"custom_example_lm_2.2" # prefix of pickle file for saving the results
N_free = 3 # number of free modes to use
# fitter object
fitter = QNMFitVaryingStartingTime(
h, t0_arr, N_free = N_free,
qnm_fixed_list = qnm_fixed_list, load_pickle = True,
run_string_prefix = run_string_prefix)
t0_arr = np.linspace(0, 50, num = 101) # array of starting times to fit for
# t0 = 0 is the peak of the strain
qnm_fixed_list = [] # list of QNMs with fixed frequencies in the fit model
run_string_prefix = f"custom_example_lm_2.2" # prefix of pickle file for saving the results
N_free = 3 # number of free modes to use
# fitter object
fitter = QNMFitVaryingStartingTime(
h, t0_arr, N_free = N_free,
qnm_fixed_list = qnm_fixed_list, load_pickle = True,
run_string_prefix = run_string_prefix)
In [6]:
Copied!
# do fits, in series from lowest to highest starting time
# This takes a little longer to run the first time
fitter.do_fits()
# do fits, in series from lowest to highest starting time
# This takes a little longer to run the first time
fitter.do_fits()
Runname: custom_example_lm_2.2, fitting for N_free = 3. Status: 0%| | 0/101 [00:00<?, ?it/s]
In [7]:
Copied!
# fitter results object
result = fitter.result_full
# fitter results object
result = fitter.result_full
In [8]:
Copied!
fig, ax = plt.subplots()
# mode locations to visualize on the plot
predicted_qnms = mode_list(['2.2.0', '2.2.1', '3.2.0'], Mf, af)
plot_omega_free(result, ax)
plot_predicted_qnms(ax, predicted_qnms)
fig, ax = plt.subplots()
# mode locations to visualize on the plot
predicted_qnms = mode_list(['2.2.0', '2.2.1', '3.2.0'], Mf, af)
plot_omega_free(result, ax)
plot_predicted_qnms(ax, predicted_qnms)
Fixed QNMs (fixed frequencies)¶
In [9]:
Copied!
t0_arr = np.linspace(0, 50, num = 101) # array of starting times to fit for
# t0 = 0 is the peak of the straisn
qnm_fixed_list = mode_list(['2.2.0', '2.2.1', '3.2.0'],
Mf, af) # list of QNMs with fixed frequencies in the fit model
run_string_prefix = f"custom_example_lm_2.2" # prefix of pickle file for saving the results
N_free = 0 # number of free modes to use
# fitter object
fitter = QNMFitVaryingStartingTime(
h, t0_arr, N_free = N_free,
qnm_fixed_list = qnm_fixed_list, load_pickle = True,
run_string_prefix = run_string_prefix)
t0_arr = np.linspace(0, 50, num = 101) # array of starting times to fit for
# t0 = 0 is the peak of the straisn
qnm_fixed_list = mode_list(['2.2.0', '2.2.1', '3.2.0'],
Mf, af) # list of QNMs with fixed frequencies in the fit model
run_string_prefix = f"custom_example_lm_2.2" # prefix of pickle file for saving the results
N_free = 0 # number of free modes to use
# fitter object
fitter = QNMFitVaryingStartingTime(
h, t0_arr, N_free = N_free,
qnm_fixed_list = qnm_fixed_list, load_pickle = True,
run_string_prefix = run_string_prefix)
In [10]:
Copied!
fitter.do_fits()
result = fitter.result_full
fitter.do_fits()
result = fitter.result_full
Runname: custom_example_lm_2.2, fitting with the following modes: 2.2.0, 2.2.1, 3.2.0. Status: 0%| …
In [11]:
Copied!
fig, axs = plt.subplots(1, 2, figsize = (12, 5))
plot_amplitudes(result, fixed_modes = qnm_fixed_list, ax = axs[0])
plot_phases(result, fixed_modes = qnm_fixed_list, ax = axs[1], legend = False)
fig, axs = plt.subplots(1, 2, figsize = (12, 5))
plot_amplitudes(result, fixed_modes = qnm_fixed_list, ax = axs[0])
plot_phases(result, fixed_modes = qnm_fixed_list, ax = axs[1], legend = False)
Mode search¶
In [12]:
Copied!
h.set_lm(2, 2)
relevant_lm_list = [(2,2)]
run_string_prefix = f"custom_example_lm_2.2" # prefix of pickle file for saving the results
mode_searcher = ModeSearchAllFreeVaryingN(h, Mf, af,
relevant_lm_list = relevant_lm_list,
N_list = [3],
t0_arr = t0_arr,
run_string_prefix = run_string_prefix
)
h.set_lm(2, 2)
relevant_lm_list = [(2,2)]
run_string_prefix = f"custom_example_lm_2.2" # prefix of pickle file for saving the results
mode_searcher = ModeSearchAllFreeVaryingN(h, Mf, af,
relevant_lm_list = relevant_lm_list,
N_list = [3],
t0_arr = t0_arr,
run_string_prefix = run_string_prefix
)
In [13]:
Copied!
mode_searcher.do_mode_searches()
mode_searcher.do_mode_searches()
Loaded fit custom_example_lm_2.2_N_3_t0_0.0000_50.0000_101 from an old run. Runname: custom_example_lm_2.2, N_free = 3, potential modes: 2.2.1, 3.2.0, 2.2.0 Performing amplitude and phase flatness check for N_free = 3 Loaded fit custom_example_lm_2.2_N_0_fix_2.2.0_2.2.1_3.2.0_t0_0.0000_50.0000_101 from an old run. Runname: custom_example_lm_2.2, N_free = 3, found the following 3 modes: 2.2.1, 3.2.0, 2.2.0
In [14]:
Copied!
plot_mode_searcher_results(mode_searcher)
plot_mode_searcher_results(mode_searcher)