diverse

Library/plots.py

@@ -98,7 +98,7 @@ def plot_regression_on_ax(data_true, data_predicted, ax=None, label='', data_lab
     ax.plot([x_min, x_max], [x_min, x_max], 'r--')
     R2 = metrics.r2_score(data_true, data_predicted)
     # mse = metrics.mean_squared_error(data_true, data_predicted)
-    ax.text(0.05, 0.9, f'{label}$R^2 = {R2:.2f}$', transform=ax.transAxes, fontsize=28, **kwargs)
+    ax.text(0.05, 0.9, f'{label}$R^2 = {R2:.2f}$', transform=ax.transAxes, fontsize=40, color='white', bbox={'facecolor':'black', 'alpha':0.5, 'pad':5})
     # ax.text(0.05, 0.8, f'{label}$MSE = {mse:.2f}$', transform=ax.transAxes, fontsize=28, **kwargs)
     return ax
 

Technical_notes/timing_dlogL/notes.md

@@ -230,7 +230,7 @@ timing = n_swg*swg + n_exp2pifdt*exp2pifdt + n_stimesexp*stimesexp + n_nwip*nwip
 
 
 
-## Analysis for IMRPhenomD_NRTidal with 12 parameters
+## Analysis for IMRPhenomD_NRTidal with 12 parameters - 59sec
 
 ### Caracteristics
 - IMRPhenomD_NRTidal
@@ -365,6 +365,114 @@ phase1_timing = traj_timing * 1500 / 3600
 
 
 
+
+
+
+## Analysis for IMRPhenomD_NRTidal with 12 parameters - 64sec
+
+### Caracteristics
+- IMRPhenomD_NRTidal
+- 12 parameters: aligned spins, tides but phase marginalization
+- flow = 30 Hz => duration = 64 sec
+```py
+len(frequency_domain_strain) = 120833
+ifo.frequency_mask.sum() = 119063
+```
+
+### Timing at the end of phase1
+
+### Timing of functions directly
+```py
+%timeit likelihood_gradient.calculate_dlogL_s_minus_h_dh():
+  900 ms ± 14.5 ms ??
+%timeit waveform_generator.frequency_domain_strain(self.likelihood.parameters):
+  31.3 ms ± 283 µs
+%timeit ifo.get_fd_time_translation(self.likelihood.parameters):
+  4.05 ms ± 29.1 µs
+%timeit ifo.get_detector_response(waveform_polarizations, self.likelihood.parameters, exp_two_pi_f_dt=ifo.exp_two_pi_f_dt_at_point):
+  2.19 ms ± 34.1 µs
+%timeit ifo.inner_product(template_plus):
+  2.4 ms ± 34.3 µs
+%timeit nwip(ifo.fd_strain, template_plus, ifo.psd_array, ifo.strain_data.duration):
+  1.15 ms ± 93.9 µs
+```
+
+### Why `ifo.inner_product()` takes longer than `nwip()`
+
+```py
+def inner_product(self, signal):
+return gwutils.noise_weighted_inner_product(
+    aa=signal[self.strain_data.frequency_mask],
+    bb=self.strain_data.frequency_domain_strain[self.strain_data.frequency_mask],
+    power_spectral_density=self.power_spectral_density_array[self.strain_data.frequency_mask],
+    duration=self.strain_data.duration)
+
+%timeit signal[self.strain_data.frequency_mask]:
+  209 µs ± 27.5 µs
+%timeit self.strain_data.frequency_domain_strain[self.strain_data.frequency_mask]:
+  504 µs ± 1.62 µs
+%timeit self.power_spectral_density_array[self.strain_data.frequency_mask]:
+  305 µs ± 3.52 µs
+```
+
+- So "masking" an array from 131072 elements to 129153 takes 209 µs.
+- Calling the `ifo.strain_data.frequency_domain_strain` is longer because the following is done internally each time:
+```py
+ifo.strain_data._frequency_domain_strain * ifo.strain_data.frequency_mask
+```
+- Calling the psd is longer because when `ifo.power_spectral_density_array` is called, the psd array is multiplied each time by the `ifo.strain_data.window_factor`.
+
+I could create a new nwip function for my Interferometer Class which would be like:
+```py
+def inner_product_mask_fixed(self, signal_masked):
+  return gwutils.noise_weighted_inner_product(
+      aa=signal_masked,
+      bb=ifo.strain_data.frequency_domain_strain_masked,
+      power_spectral_density=ifo.power_spectral_density_array_windowed_masked,
+      duration=ifo.strain_data.duration)
+```
+
+### Timing analysis
+- Parameters contributions to `dlogL`
+
+
+|  | Diff method | h_source | exp(i * 2pi * f * dt) | s = s * exp(i * 2pi * f * dt) | < h1 \| h2 > |
+|:----|:----|:----:|:----:|:----:|:----:|
+| At point | - | 1 | 3 | 3 | 6 |
+| cos(iota) | Forward | 1 | 0 | 3 | 6 |
+| psi | Forward | 0 | 0 | 3 | 6 |
+| ln(D) | Forward | 1 | 0 | 3 | 6 |
+| ln(Mc) | Central | 2 | 0 | 6 | 9 |
+| ln(mu) | Central | 2 | 0 | 6 | 9 |
+| sin(dec) | Forward | 0 | 3 | 3 | 6 |
+| ra | Forward | 0 | 3 | 3 | 6 |
+| ln(tc) | Central | 0 | 6 | 6 | 9 |
+| chi_1 | Central | 2 | 0 | 6 | 9 |
+| chi_2 | Central | 2 | 0 | 6 | 9 |
+| lambda_1 | Forward | 1 | 0 | 3 | 6 |
+| lambda_2 | Forward | 1 | 0 | 3 | 6 |
+| **TOTAL** | - | 15 | 15 | 54 | 93 |
+| **IMRPhenomD_NRTidal, 30Hz** | **909ms** | 15x31.3ms | 15x4.05ms | 54x2.2ms | 93x2.8ms |
+
+```py
+swg = 31.3
+exp2pifdt = 4.05
+stimesexp = 2.2
+nwip = 2.4
+
+n_swg = 15
+n_exp2pifdt = 15
+n_stimesexp = 54
+n_nwip = 93
+
+timing = n_swg*swg + n_exp2pifdt*exp2pifdt + n_stimesexp*stimesexp + n_nwip*nwip
+```
+
+
+
+
+
+
 ## Analysis forecasts IMRPhenomPv2_NRTidal with 16 parameters
 
 |  | Diff method | h_source | exp(i * 2pi * f * dt) | s = s * exp(i * 2pi * f * dt) | < h1 \| h2 > |

Tests/evaluate_dnn_perf.py

+import sys
+sys.path.append('../')
+import numpy as np
+import Library.Fit_NR as fnr
+import Library.plots as plots
+import Library.CONST as CONST
+
+
+if __name__ == '__main__':
+
+    global_fit_dir = '../__output_data/GW190412real/0_1_2_3_4_5_6_7_8_15_16/PSD_2/501500_1500_2000_200_200_0.0123/IMRPhenomHM/sampler_state/global_fit'
+
+    train_size = 286200
+    train_size = 504297
+    train_size_max = 543185
+    input_dir = global_fit_dir + f'/train_size_{train_size}'
+    # input_dir = global_fit_dir + '/train_size_286200'
+    # input_dir = global_fit_dir + '/train_size_424623'
+    # input_dir = global_fit_dir + '/train_size_504297'
+    # input_dir = global_fit_dir + '/train_size_543185'
+    outdir_plots = input_dir + '/plots'
+
+
+    x_train = np.loadtxt(input_dir + '/qpos_trained.dat')
+    y_train = np.loadtxt(input_dir + '/dlogL_trained.dat')
+
+    if train_size == train_size_max:
+    # if True:
+        x_val = None
+        y_val = None
+        validation_data = None
+    else:
+        x_val = np.loadtxt(global_fit_dir + f'/train_size_{train_size_max}/qpos_trained.dat')
+        y_val = np.loadtxt(global_fit_dir + f'/train_size_{train_size_max}/dlogL_trained.dat')
+        val_size = train_size_max - 504297
+        x_val = x_val[-val_size:]
+        y_val = y_val[-val_size:]
+        validation_data = (x_val, y_val)
+
+
+    epochs = 60
+    batch_size = 128
+    fit_method = fnr.DNNRegressorMultipleGradients(
+        n_dim=x_train.shape[1],
+        batch_size=batch_size,
+        epochs=epochs
+    )
+
+    print(f'Loading the {fit_method.nick_name} fit from saved directory {input_dir}...')
+    fit_method.load(input_dir)
+
+    search_trajectory_parameters_keys = ['cos_theta_jn', 'phase', 'psi', 'log_luminosity_distance', 'log_chirp_mass', 'log_reduced_mass', 'sin_dec', 'ra', 'log_geocent_duration', 'chi_1', 'chi_2']
+    dlogL_latex_labels = CONST.get_dlogL_latex_labels(search_trajectory_parameters_keys)
+    plots.make_regression_plots(fit_method, x_train, y_train, outdir_plots, validation_data=validation_data, y_labels=dlogL_latex_labels)

Tests/time_waveform_gen.py

@@ -109,12 +109,12 @@ if __name__ == '__main__':
     #   print(lalsim.GetStringFromApproximant(i))
     approximant = 'TaylorF2'
     approximant = 'IMRPhenomD'
-    # approximant = 'IMRPhenomD_NRTidal'
+    approximant = 'IMRPhenomD_NRTidal'
     # approximant = 'IMRPhenomPv2'
-    approximant = 'IMRPhenomPv2_NRTidal'
-    approximant = 'SEOBNRv4'
+    # approximant = 'IMRPhenomPv2_NRTidal'
+    # approximant = 'SEOBNRv4'
     # ROM file can be found here:https://git.ligo.org/lscsoft/lalsuite-extra/-/tree/master/data/lalsimulation
-    approximant = 'SEOBNRv4_ROM'
+    # approximant = 'SEOBNRv4_ROM'
     ifos = ['H1', 'L1', 'V1']
     sampling_frequency = 4096
     maximum_frequency = sampling_frequency/2

examples/config_files/config_GW170817_IMRPhenomD_NRTidal.ini

@@ -19,7 +19,7 @@ prior_file=../gw/prior_files/GW170817_LALInf_IMRPhenomD_NRTidal.prior
 
 [hmc]
 # Total number of trajectories of HMC
-n_traj_hmc_tot = 1500
+n_traj_hmc_tot = 10
 # Number of trajectories using  numericalgradients in phase 1.
 n_traj_fit = 1500
 # Number of iterations if phase1 loop to go over for this run. This permits to split phase1 into several chunks. If set to 0, then the value will be set to that of n_traj_fit
@@ -31,10 +31,10 @@ n_fit_2 = 200
 # Length of numerical trajectories during phase1
 length_num_traj = 200
 # Stepsize between steps of trajectories used for phase1.
-epsilon0 = 1.12e-2
+epsilon0 = 5e-3
 # Python file to configure your search parameter space
 ; config_hmc_parameters_file = ../examples/config_params/config_hmc_parameters.json
 config_hmc_parameters_file = ../examples/config_params/GW170817_IMRPhenomD_NRTidal.json
 ; config_hmc_parameters_file = ../examples/config_params/config_hmc_parameters_spins_P.json
 # Number of statistically independant samples desired.
-n_sis = 5000
+n_sis = 6622

examples/config_files/config_GW190412_IMRPhenomHM.ini

@@ -18,7 +18,7 @@ prior_file=../gw/prior_files/GW190412_alignedspins_withphic.prior
 
 [hmc]
 # Total number of trajectories of HMC
-n_traj_hmc_tot = 1500
+n_traj_hmc_tot = 501500
 # Number of trajectories using  numericalgradients in phase 1.
 n_traj_fit = 1500
 # Number of iterations if phase1 loop to go over for this run. This permits to split phase1 into several chunks. If set to 0, then the value will be set to that of n_traj_fit
@@ -30,7 +30,7 @@ n_fit_2 = 200
 # Length of numerical trajectories during phase1
 length_num_traj = 200
 # Stepsize between steps of trajectories used for phase1.
-epsilon0 = 5e-3
+epsilon0 = 1.23e-2
 # Python file to configure your search parameter space
 ; config_hmc_parameters_file = ../examples/config_params/config_hmc_parameters.json
 config_hmc_parameters_file = ../examples/config_params/GW190412_IMRPhenomHM.json

examples/config_files/config_GW190814_IMRPhenomHM.ini

@@ -18,9 +18,9 @@ prior_file=../gw/prior_files/GW190814_alignedspins_withphic.prior
 
 [hmc]
 # Total number of trajectories of HMC
-n_traj_hmc_tot = 1500
+n_traj_hmc_tot = 15
 # Number of trajectories using  numericalgradients in phase 1.
-n_traj_fit = 1500
+n_traj_fit = 15
 # Number of iterations if phase1 loop to go over for this run. This permits to split phase1 into several chunks. If set to 0, then the value will be set to that of n_traj_fit
 n_traj_for_this_run = 1000000
 # Number of points used to sub-select points in Ordered Look-Up Tables for local fit bimodal parameters like (cos(inc), psi, logD)