diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index dc7639680..b724642bd 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -134,6 +134,27 @@ jobs: - name: Run simulation_manager smoke test run: bash .travis/test-simulation-manager.sh + lisa-check: + needs: install + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: '3.10' + cache: 'pip' + cache-dependency-path: requirements.txt + - name: Enable symlink + run: sudo ln -sf $(which python3) /usr/bin/python + - name: Install dependencies + run: | + python -m pip install --upgrade pip --break-system-packages + python -m pip install -r requirements.txt --break-system-packages + python -m pip install coverage pytest --break-system-packages + python -m pip install --editable . --break-system-packages + - name: Run LISA smoke and contract tests + run: bash .travis/test-lisa.sh + integration-check: needs: install runs-on: ubuntu-latest diff --git a/.travis/test-lisa.sh b/.travis/test-lisa.sh new file mode 100644 index 000000000..e821b5739 --- /dev/null +++ b/.travis/test-lisa.sh @@ -0,0 +1,18 @@ +#!/usr/bin/env bash +set -euo pipefail + +PYTHON_BIN="${RIFT_LISA_PYTHON:-${PYTHON:-python}}" +if ! command -v "${PYTHON_BIN}" >/dev/null 2>&1; then + PYTHON_BIN="$(command -v python3)" +fi + +"${PYTHON_BIN}" -m pytest -q \ + MonteCarloMarginalizeCode/Code/test/test_lisa_auxiliary_imports.py \ + MonteCarloMarginalizeCode/Code/test/test_lisa_response_import.py \ + MonteCarloMarginalizeCode/Code/test/test_lisa_lalsimutils_compat.py \ + MonteCarloMarginalizeCode/Code/test/test_lisa_run_checks.py \ + MonteCarloMarginalizeCode/Code/test/test_lisa_demo_contract.py \ + MonteCarloMarginalizeCode/Code/test/test_lisa_helper_contract.py \ + MonteCarloMarginalizeCode/Code/test/test_lisa_pseudo_pipe_contract.py \ + MonteCarloMarginalizeCode/Code/test/test_lisa_pp_surface.py \ + MonteCarloMarginalizeCode/Code/test/test_lisa_synthetic_demo.py diff --git a/MonteCarloMarginalizeCode/Code/RIFT/LISA/initial_grid/__init__.py b/MonteCarloMarginalizeCode/Code/RIFT/LISA/initial_grid/__init__.py new file mode 100644 index 000000000..b54333b16 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/LISA/initial_grid/__init__.py @@ -0,0 +1 @@ +"""Initial-grid helpers for LISA analyses.""" diff --git a/MonteCarloMarginalizeCode/Code/RIFT/LISA/initial_grid/fisher_errors.py b/MonteCarloMarginalizeCode/Code/RIFT/LISA/initial_grid/fisher_errors.py new file mode 100644 index 000000000..463fe8fca --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/LISA/initial_grid/fisher_errors.py @@ -0,0 +1,264 @@ +#!/usr/bin/env python + +"""The purpose of this code is to give rough estimates of width of posteriors for Mc, eta, spin and skylocation for initial grid generation.""" +import numpy as np +import RIFT.lalsimutils as lsu +from RIFT.LISA.response.LISA_response import * +from argparse import ArgumentParser +from scipy.interpolate import interp1d +import os + +__author__ = "A. Jan" +########################################################################################### +# Functions to generate 2-PN waveforms as per http://arxiv.org/abs/gr-qc/9502040. +########################################################################################### +def amplitude_2PN(fvals, Mc): + """Returns amplitude in frequency domains for a 2-PN waveform""" +# print(f"amplitude_2PN: {locals()}") + return 1 * Mc**(5/6) * fvals**(-7/6) + +def phase_2PN(fvals, Mc, eta, sigma, beta, coa_phase=0, coa_time=0): + """Returns phase in frequency domains for a 2-PN waveform""" +# print(f"phase_2PN: {locals()}") + M = Mc/eta**(3/5) + fac = np.pi*M*fvals + newtonian = 3/128 * (np.pi*Mc*fvals)**(-5/3) + one_PN = 20/9 * (743/336 + 11/4* eta)*(fac)**(2/3) + one_five_PN = 4 * (4*np.pi - beta)*fac + two_PN = 10*(3058673/1016064 + 5429/1008*eta + 617/144*eta**2 - sigma)*(fac)**(4/3) + phase_vals = 2*np.pi*fvals*coa_time - coa_phase - np.pi/4 + newtonian * (1 + one_PN - one_five_PN + two_PN) + return phase_vals + +def time_2PN(fvals, Mc, eta, sigma, beta, coa_time=0): + """Returns time corresponding to input frequency for a 2-PN waveform""" + M = Mc/eta**(3/5) + fac = np.pi*M*fvals + time_vals = coa_time - 5/256* Mc*(np.pi*Mc*fvals)**(-8/3) * (1 + 4/3 * (743/336 + 11/4 * eta)*(fac)**(2/3) - 8/5 *(4*np.pi - beta)*fac + 2*(3058673/1016064 + 5429/1008 * eta + 617/144*eta**2 - sigma)*(fac)**(4/3)) + return time_vals + +def get_sigma_beta(Mc, eta, a1z, a2z): + """Returns sigma and beta parameters that are used in construction of the waveform. These parameters contain information about spin.""" + alpha = Mc / eta**(3/5) + beta = Mc**2 / eta**(1/5) + m1 = 0.5 * (alpha + np.sqrt(alpha**2 - 4*beta)) + m2 = 0.5 * (alpha - np.sqrt(alpha**2 - 4*beta)) + M = m1 + m2 + sigma_val = eta/48 * (-247 * a1z*a1z + 721 * a1z*a2z) + beta_val = 1/12 * ((113 * (m1/M)**2 + 75 * eta)*a1z + (113 * (m2/M)**2 + 75 * eta)*a2z ) + return sigma_val, beta_val + +def get_derivative(fvals, Mc, eta, sigma, beta, wf): + """Derivates of the waveform with respect to coalescence time, coalescence phase, Mc, eta, sigma, beta""" +# print(f"get_derivative: {locals()}") + M = Mc/eta**(3/5) + v = (np.pi*M*fvals)**(1/3) + + A4 = 4/3 * (743/336 + 11/4 * eta) + B4 = 8/5 * (4*np.pi - beta) + C4 = 2 * (3058673/1016064 + 5429/1008 * eta + 617/144 * eta**2 - sigma) + + A5 = 743/168 - 33/4 * eta + B5 = 27/5 * (4*np.pi - beta) + C5 = 18 * (3058673/1016064 - 5429/4032 * eta - 617/96*eta**2 - sigma) + + d_tc = 2*np.pi*1j*(fvals) * wf + d_phi = -1j * wf + + d_log_mc = -1j*(5/128 * (np.pi*Mc*fvals)**(-5/3) * (1 + A4*v**2 - B4*v**3 + C4*v**4)) * wf + d_log_eta = -1j*(1/96 * (np.pi*Mc*fvals)**(-5/3) * ( A5*v**2 - B5*v**3 + C5*v**4)) * wf + + d_beta = 1j * 3/32 * eta**(-3/5) * (np.pi*Mc*fvals)**(-2/3) * wf + d_sigma = -1j * 15/64 * eta**(-4/5) * (np.pi*Mc*fvals)**(-1/3) * wf + return np.array([d_tc, d_phi, d_log_mc, d_log_eta, d_beta, d_sigma]) + +def get_wf(fvals, Mc, eta, sigma, beta, psd_vals, coa_phase=0, coa_time=0, snr=None, LISA_response=False, skylocation = None): + """Generate a 2-PN waveform""" + phase = phase_2PN(fvals, Mc, eta, sigma, beta, coa_phase, coa_time) + amp = amplitude_2PN(fvals, Mc) + wf = amp * np.exp(1j*phase) + if LISA_response: + H = transformed_Hplus_Hcross(skylocation[0], skylocation[1], 0.0, 0.0, 0.0, 2, 2) + time = time_2PN(fvals, Mc, eta, sigma, beta, coa_time) + A, E, T = Evaluate_Gslr(time + coa_time, fvals, H, skylocation[0], skylocation[1]) + wf = wf * A + if snr: + # bring the source closer or further, depending on SNR + deltaF = np.diff(fvals)[0] + snr_fiducial = np.sqrt(get_inner_product(wf, wf, psd_vals, deltaF)) + correction = snr / snr_fiducial + wf = correction * wf + return wf + +########################################################################################### +# Utilities +########################################################################################### +def load_psd(psdf, fvals): + """Loads in PSD""" + # load in psd + psd_dict = {} + inst = "A" + print( "Reading PSD for instrument %s from %s" % (inst, psdf)) + psd_dict[inst] = lsu.get_psd_series_from_xmldoc(psdf, inst) + psd_fvals = psd_dict[inst].f0 + psd_dict[inst].deltaF*np.arange(psd_dict[inst].data.length) + interp_func = interp1d(psd_fvals, psd_dict[inst].data.data) + return interp_func(fvals) + +def get_mass_from_mc_eta(mc, eta): + """Returns m1, m2 from mc and eta.""" + alpha = mc / eta**(3/5) + beta = mc**2 / eta**(1/5) + m1 = 0.5 * (alpha + np.sqrt(alpha**2 - 4*beta)) + m2 = 0.5 * (alpha - np.sqrt(alpha**2 - 4*beta)) + return m1, m2 + +def get_mc_eta_from_mass(m1, m2): + """Returns m1, m2 from mc and eta.""" + mc = (m1*m2)**(3/5) / (m1+m2)**(1/5) + eta = (m1*m2) / (m1+m2)**(2) + if eta==0.25: + eta=0.24999 + return mc, eta + +def get_inner_product(wf1, wf2, psd_vals, deltaF): + """Calculate inner product""" + assert len(wf1) == len(wf2) == len(psd_vals) + weight = 1/psd_vals + intgd = np.sum(np.conj(wf1) * wf2 * weight) * deltaF + return 4 * np.real(intgd) + +def get_massratio_error(eta, q, eta_error): + return eta_error * (1 + q)**3 / (1 - q) + +def get_spin_error(eta, q, a1z, a2z, eta_error, beta_error, sigma_error): + q_error = get_massratio_error(eta, q, eta_error) + if round(a1z,4)==round(a2z,4): + a1z = a1z + 0.001 * a1z + a2z = a2z - 0.001 * a2z + + c1 = sigma_error - eta_error/48 * (474*a1z*a2z) + a1 = eta/48 * 474 * a2z + b1 = eta/48 * 474 * a1z + + c2 = beta_error - eta_error/12 * ( (113/q + 75)*a1z + (113*q + 75)*a2z) - eta/12 * ( (-113/q**2)*a1z*q_error + 113*a2z*q_error) + a2 = eta/12 * (113/q + 75) + b2 = eta/12 * (113*q + 75) + + coefficients = np.array([ [a1, a2], [b1, b2] ]) + dependents = np.array( [c1, c2] ) + answers = np.abs(np.linalg.solve(coefficients, dependents)) + return [np.min(answers), np.max(answers)] + +########################################################################################### +# Fisher matrix +########################################################################################### +def get_fisher_matrix(Mc, eta, sigma, beta, fvals, psd_vals, deltaF, wf): + """Get fisher information matrix""" + derivatives = get_derivative(fvals, Mc, eta, sigma, beta, wf) + N = 6 + tau_ij = np.zeros((N,N)) + for i in np.arange(0, N): + for j in np.arange(0, N): + tau_ij[i,j] = get_inner_product(derivatives[i], derivatives[j], psd_vals, deltaF) + inv_tau_ij = (np.linalg.inv(tau_ij)) + return tau_ij, inv_tau_ij + +def get_error_bounds(P_inj, snr, psd_path, snr_fmin=0.0001): + """Get error bounds on parameters using fisher information matrix.""" + response=True # use LISA response + deltaF = 0.00001 # hardcoded deltaF + mc, eta = get_mc_eta_from_mass(P_inj.m1/lsu.lsu_MSUN, P_inj.m2/lsu.lsu_MSUN) + q = P_inj.m2/P_inj.m1 + + if q == 1: + q = 0.9 + + # convert chirp mass to seconds + Mc = mc * 5 * 10**(-6) + M = Mc/eta**(3/5) + sigma, beta = get_sigma_beta(Mc, eta, P_inj.s1z, P_inj.s2z) + fmax = 6**(-3/2) / np.pi/ M + print(f"Fmax is = {fmax} Hz") + fvals = np.arange(float(snr_fmin), fmax, deltaF) + + # Load psd + psd_vals = load_psd(psd_path, fvals) + + # generate waveform + wf = get_wf(fvals, Mc, eta, sigma, beta, psd_vals, 0, float(P_inj.tref), snr=float(snr), LISA_response=response, skylocation=[P_inj.theta, P_inj.phi]) + print(f"SNR of generated waveform is = {np.sqrt(get_inner_product(wf,wf,psd_vals, deltaF))}") + + # Calculate fisher matrix + tau_ij, inv_tau_ij = get_fisher_matrix(Mc, eta, sigma, beta, fvals, psd_vals, deltaF, wf) + if eta < 0.24: + factor_eta = 12.5 + if eta >=0.24: # the errors estimates seem to be large for q~1 case. + factor_eta = 1.0 + factor_mc = 50 + factor_spin1 = 60 + factor_spin2 = 60 + spin_bounds = get_spin_error(eta, q, P_inj.s1z, P_inj.s2z, (np.sqrt(1/tau_ij[3,3]))*eta, np.sqrt(1/tau_ij[4,4]), np.sqrt(1/tau_ij[5,5])) + + print(f"Mc span = {2*factor_mc*np.sqrt(1/tau_ij[2,2])*mc}, eta span = {2*np.sqrt(1/tau_ij[3,3])*eta*factor_eta}, s1z span = {2*factor_spin1*spin_bounds[0]}, s2z span = {2*factor_spin2*spin_bounds[1]}, beta span = {0.036*(210/snr)**2}, lambda span = {0.044*(210/snr)**2}") + + + return np.array([ mc - factor_mc*np.sqrt(1/tau_ij[2,2])*mc, mc + factor_mc*np.sqrt(1/tau_ij[2,2])*mc, eta-(np.sqrt(1/tau_ij[3,3]))*eta*factor_eta, eta+(np.sqrt(1/tau_ij[3,3]))*eta*factor_eta, P_inj.s1z - factor_spin1*spin_bounds[0], P_inj.s1z + factor_spin1*spin_bounds[0], P_inj.s2z - factor_spin2*spin_bounds[1], P_inj.s2z + factor_spin2*spin_bounds[1], P_inj.theta - 0.018*(210/snr), P_inj.theta + 0.018*(210/snr), P_inj.phi - 0.022*(210/snr), P_inj.phi + 0.022*(210/snr)]) + + +########################################################################################### + +if __name__ =='__main__': + ########################################################################################### + parser=ArgumentParser() + parser.add_argument("--inj", help="Full path to mdc.xml.gz") + parser.add_argument("--psd-path", help="Full path to A-psd.xm.gz") + parser.add_argument("--snr", help="SNR of the signal") + parser.add_argument("--snr-fmin", help="fmin used in snr calculations", default=0.0001) + parser.add_argument("--generate-grid", help="Use the fisherbounds to generate grid", default=True) + parser.add_argument("--points", help="number of points in the grid", default=25000) + opts = parser.parse_args() + print(f"Loading file:\n {opts.inj}") + P_inj_list = lsu.xml_to_ChooseWaveformParams_array(opts.inj) + P_inj = P_inj_list[0] + print("######") + print(f"m1 = {P_inj.m1/lsu.lsu_MSUN}, m2 = {P_inj.m2/lsu.lsu_MSUN}, s1z = {P_inj.s1z}, s2z = {P_inj.s2z}, beta = {P_inj.theta}, lambda = {P_inj.phi}, tref = {P_inj.tref} s") + print("######") + error_bounds = get_error_bounds(P_inj, float(opts.snr), opts.psd_path, snr_fmin=opts.snr_fmin) + print(f"Mc bounds = [{error_bounds[0]:0.2f}, {error_bounds[1]:0.2f}]") + print(f"eta bounds = [{error_bounds[2]:0.8f}, {error_bounds[3]:0.8f}]") + print(f"s1z bounds = [{error_bounds[4]:0.6f}, {error_bounds[5]:0.6f}]") + print(f"s2z bounds = [{error_bounds[6]:0.6f}, {error_bounds[7]:0.6f}]") + print(f"beta bounds = [{error_bounds[8]:0.6f}, {error_bounds[9]:0.6f}]") + print(f"lambda bounds = [{error_bounds[10]:0.6f}, {error_bounds[11]:0.6f}]") + mc_span, eta_span = (error_bounds[1] - error_bounds[0]), (error_bounds[3] - error_bounds[2]) + s1z_span, s2z_span = (error_bounds[5] - error_bounds[4]), (error_bounds[7] - error_bounds[6]) + beta_span, lambda_span = (error_bounds[9] - error_bounds[8]), (error_bounds[11] - error_bounds[10]) + if opts.generate_grid: + import os + from RIFT.LISA.utils.utils import * + cmd = f"util_ManualOverlapGrid.py --inj {opts.inj} " + cmd += f"--parameter mc --parameter-range '[{error_bounds[0]:0.2f}, {error_bounds[1]:0.2f}]' " + cmd += f"--parameter eta --parameter-range '[{error_bounds[2]:0.5f}, {error_bounds[3]:0.5f}]' " + cmd += f"--random-parameter s1z --random-parameter-range '[{error_bounds[4]:0.6f}, {error_bounds[5]:0.6f}]' " + cmd += f"--random-parameter s2z --random-parameter-range '[{error_bounds[6]:0.6f}, {error_bounds[7]:0.6f}]' " + cmd += f"--random-parameter theta --random-parameter-range '[{error_bounds[8]:0.6f}, {error_bounds[9]:0.6f}]' " + cmd += f"--random-parameter phi --random-parameter-range '[{error_bounds[10]:0.6f}, {error_bounds[11]:0.6f}]' " + cmd += f"--grid-cartesian-npts {int(opts.points)} --skip-overlap" + print(f"\t Generating grid\n{cmd}") + os.system(cmd) + os.system('mv overlap-grid.xml.gz overlap-grid-primary.xml.gz') + + secondary_peak = get_secondary_mode_for_skylocation(float(P_inj.tref), P_inj.phi, P_inj.theta) + beta_sec, lambda_sec = secondary_peak[0,2], secondary_peak[0,1] + print(f"Secondary peak: lambda {lambda_sec}, beta {beta_sec}") + cmd = f"util_ManualOverlapGrid.py --inj {opts.inj} " + cmd += f"--parameter mc --parameter-range '[{error_bounds[0]:0.2f}, {error_bounds[1]:0.2f}]' " + cmd += f"--parameter eta --parameter-range '[{error_bounds[2]:0.5f}, {error_bounds[3]:0.5f}]' " + cmd += f"--random-parameter s1z --random-parameter-range '[{error_bounds[4]:0.6f}, {error_bounds[5]:0.6f}]' " + cmd += f"--random-parameter s2z --random-parameter-range '[{error_bounds[6]:0.6f}, {error_bounds[7]:0.6f}]' " + cmd += f"--random-parameter theta --random-parameter-range '[{beta_sec-0.5*beta_span}, {beta_sec+0.5*beta_span}]' " + cmd += f"--random-parameter phi --random-parameter-range '[{lambda_sec-0.5*lambda_span}, {lambda_sec+0.5*lambda_span}]' " + cmd += f"--grid-cartesian-npts {int(opts.points)} --skip-overlap" + print(f"\t Generating grid\n{cmd}") + os.system(cmd) + os.system('mv overlap-grid.xml.gz overlap-grid-secondary.xml.gz') + os.system('ligolw_add overlap-grid-primary.xml.gz overlap-grid-secondary.xml.gz -o overlap-grid.xml.gz') diff --git a/MonteCarloMarginalizeCode/Code/RIFT/LISA/injections/LISA_injections.py b/MonteCarloMarginalizeCode/Code/RIFT/LISA/injections/LISA_injections.py new file mode 100644 index 000000000..d55f2512e --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/LISA/injections/LISA_injections.py @@ -0,0 +1,183 @@ +import numpy as np +import RIFT.lalsimutils as lsu +from RIFT.LISA.response.LISA_response import * +import lal +import lalsimulation +import matplotlib.pyplot as plt +import os + +__author__ = "A. Jan" + + +def load_psd(param_dict): + """ + Load Power Spectral Density (PSD) data for the LISA instrument. + + Parameters: + param_dict (dict): A dictionary containing parameters for loading the PSD, + should contain 'psd_path', 'deltaF', and 'snr_fmin'. + + Returns: + dict: A dictionary containing PSD for A, E, and T channels. + """ + print(f"Reading PSD to calculate SNR for LISA instrument from {param_dict['psd_path']}.") + psd = {} + psd["A"] = lsu.get_psd_series_from_xmldoc(param_dict["psd_path"] + "/A-psd.xml.gz", "A") + psd["A"] = lsu.resample_psd_series(psd["A"], param_dict['deltaF']) + psd_fvals = psd["A"].f0 + param_dict['deltaF']*np.arange(psd["A"].data.length) + psd["A"].data.data[ psd_fvals < param_dict['snr_fmin']] = 0 + + psd["E"] = lsu.get_psd_series_from_xmldoc(param_dict["psd_path"]+ "/E-psd.xml.gz", "E") + psd["E"] = lsu.resample_psd_series(psd["E"], param_dict['deltaF']) + psd_fvals = psd["E"].f0 + param_dict['deltaF']*np.arange(psd["E"].data.length) + psd["E"].data.data[ psd_fvals < param_dict['snr_fmin']] = 0 + + psd["T"] = lsu.get_psd_series_from_xmldoc(param_dict["psd_path"]+ "/T-psd.xml.gz", "T") + psd["T"] = lsu.resample_psd_series(psd["T"], param_dict['deltaF']) + psd_fvals = psd["T"].f0 + param_dict['deltaF']*np.arange(psd["T"].data.length) + psd["T"].data.data[ psd_fvals < param_dict['snr_fmin']] = 0 + return psd + +def calculate_snr(data_dict, fmin, fmax, fNyq, psd, only_positive_modes=True): + """ + Calculate the zero-noise Signal-to-Noise Ratio (SNR) for LISA signals. + + Parameters: + data_dict (dict): A dictionary containing the A, E, and T signal data, + fmin (float): The minimum frequency for integration in Hz, + fmax (float): The maximum frequency for integration in Hz, + fNyq (float): The Nyquist frequency in Hz, + psd (dict): A dictionary containing the PSD for A, E, and T channels. + + Returns: + float: The total zero-noise SNR calculated across all channels. + """ + + assert data_dict["A"].deltaF == data_dict["E"].deltaF == data_dict["T"].deltaF + print(f"Integrating from {fmin} to {fmax} Hz.") + + # create instance of inner product + IP_A = lsu.ComplexIP(fmin, fmax, fNyq, psd["A"].deltaF, psd["A"], False, False, 0.0,) + IP_E = lsu.ComplexIP(fmin, fmax, fNyq, psd["A"].deltaF, psd["E"], False, False, 0.0,) + IP_T = lsu.ComplexIP(fmin, fmax, fNyq, psd["A"].deltaF, psd["T"], False, False, 0.0,) + + IP_factor = 1 + if only_positive_modes: + IP_factor = 2 + + # calculate SNR of each channel + A_snr, E_snr, T_snr = np.sqrt(IP_factor*IP_A.ip(data_dict["A"], data_dict["A"])), np.sqrt(IP_factor*IP_E.ip(data_dict["E"], data_dict["E"])), np.sqrt(IP_factor*IP_T.ip(data_dict["T"], data_dict["T"])) + + # combine SNR + snr = np.real(np.sqrt(A_snr**2 + E_snr**2 + T_snr**2)) # SNR (zero noise) = sqrt() + + print(f"A-channel snr = {A_snr.real:0.3f}, E-channel snr = {E_snr.real:0.3f}, T-channel snr = {T_snr.real:0.3f},\n\tTotal SNR = {snr:0.3f}.") + return snr + +def create_PSD_injection_figure(data_dict, psd, injection_save_path, snr): + """ + Create a frequency-domain injection figure with PSD plotted against A, E and T data.. + + Parameters: + data_dict (dict): A dictionary containing signal data for each channel, + psd (dict): A dictionary containing the PSD for A, E, and T channels, + injection_save_path (str): The file path where the generated figure will be saved. + snr (float): The SNR to display in the figure title. + + Returns: + None: This function saves the figure to the specified path. + """ + channels = list(data_dict.keys()) + fvals = get_fvals(data_dict[channels[0]]) + + # plot data + plt.title(f"Injection vs PSD (SNR = {snr:0.2f})") + plt.xlabel("Frequency [Hz]") + plt.ylabel("Characterstic strain") + psd_fvals = psd[channels[0]].f0 + data_dict[channels[0]].deltaF*np.arange(psd[channels[0]].data.length) + + for channel in channels: + # For m > 0, hlm is define for f < 0 in lalsimulation. That's why abs is over fvals too. + data = np.abs(2*fvals*data_dict[channel].data.data) # we need get both -m and m modes, right now this only has positive modes present. + plt.loglog(-fvals, data, label = channel, linewidth = 1.2) + plt.loglog(psd_fvals, np.sqrt(psd_fvals * psd[channel].data.data), label = channel + "-psd", linewidth = 0.8) + + plt.legend(loc="upper right") + + # place x-y limits + plt.gca().set_ylim([10**(-24), 10**(-17)]) + plt.gca().set_xlim([10**(-4), 1]) + plt.grid(alpha = 0.5) + + # save + plt.savefig(injection_save_path + "/injection-psd.png", bbox_inches = "tight") + + +def generate_lisa_TDI_dict(param_dict): + print(param_dict) + P = lsu.ChooseWaveformParams() + P.m1 = param_dict["m1"] * lal.MSUN_SI + P.m2 = param_dict["m2"] * lal.MSUN_SI + P.s1x, P.s1y, P.s1z = 0.0, 0.0, param_dict["s1z"] + P.s2x, P.s2y, P.s2z = 0.0, 0.0, param_dict["s2z"] + P.dist = param_dict["dist"] * 1e6 * lal.PC_SI + + P.deltaT, P.deltaF = param_dict["deltaT"], param_dict["deltaF"] + P.fref = param_dict["wf-fref"] + P.approx = lalsimulation.GetApproximantFromString(param_dict["approx"]) + P.fmin, P.fmax = param_dict["fmin"], 0.5/P.deltaT + P.psi, P.phiref, P.inclination, P.tref, P.theta, P.phi = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + P.eccentricity, P.meanPerAno = 0.0, 0.0 + + modes = np.array(param_dict["modes"]) + lmax = np.max(modes[:,0]) + + path_to_NR_hdf5 = param_dict["path_to_NR_hdf5"] if 'path_to_NR_hdf5' in param_dict else None + + number_of_bins = 1/(P.deltaF*P.deltaT) + power_of_number_of_bins = np.log2(number_of_bins) + assert power_of_number_of_bins == np.ceil(power_of_number_of_bins), f'Number of bins needs to be a power of 2, increase 1/deltaF from {1/P.deltaF} to {1/ (2**np.ceil(power_of_number_of_bins)*P.deltaT)}.' + + print("###############") + if 1/P.deltaF/60/60/24 >0.5: + print(f"Data length = {1/P.deltaF/60/60/24:2f} days.") + else: + print(f"Data length = {1/P.deltaF/60/60:2f} hrs.") + + + print(f"\nWaveform is being generated with m1 = {P.m1/lsu.lsu_MSUN}, m2 = {P.m2/lsu.lsu_MSUN}, s1z = {P.s1z}, s2z = {P.s2z}, distance = {P.dist/1e6/lal.PC_SI}") + print(f"deltaF = {P.deltaF}, fmin = {P.fmin}, fmax = {P.fmax}, deltaT = {P.deltaT}, modes = {list(modes)}, lmax = {lmax}, tref = {param_dict['tref']}") + print(f"phiref = {param_dict['phi_ref']}, psi = {param_dict['psi']}, inclination = {param_dict['inclination']}, beta = {param_dict['beta']}, lambda = {param_dict['lambda']}") + print(f"path_to_NR_hdf5 = {path_to_NR_hdf5}, approx = {lalsimulation.GetStringFromApproximant(P.approx)}\n") + print("###############") + + hlmf = lsu.hlmoff_for_LISA(P, Lmax=lmax, modes=modes, path_to_NR_hdf5=path_to_NR_hdf5) + modes = list(hlmf.keys()) + + # create injections + data_dict = create_lisa_injections(hlmf, P.fmax, param_dict["fref"], param_dict["beta"], param_dict["lambda"], param_dict["psi"], param_dict["inclination"], param_dict["phi_ref"], param_dict["tref"]) + return data_dict + +def generate_lisa_injections(data_dict, param_dict, get_snr = True): + if not(os.path.exists(param_dict['save_path'])): + print(f"Provided path doesn't exist {param_dict['save_path']}, creating it.") + os.mkdir(param_dict["save_path"]) + create_h5_files_from_data_dict(data_dict, param_dict["save_path"]) + cmd = f"util_WriteInjectionFile.py --parameter m1 --parameter-value {param_dict['m1']} \ + --parameter m2 --parameter-value {param_dict['m2']} \ + --parameter s1x --parameter-value 0.0 --parameter s1y --parameter-value 0.0 --parameter s1z --parameter-value {param_dict['s1z']} \ + --parameter s2x --parameter-value 0.0 --parameter s2y --parameter-value 0.0 --parameter s2z --parameter-value {param_dict['s2z']} \ + --parameter eccentricity --parameter-value 0 --approx {param_dict['approx']} --parameter dist --parameter-value {param_dict['dist']} \ + --parameter fmin --parameter-value {param_dict['fmin']} --parameter incl --parameter-value {param_dict['inclination']} \ + --parameter tref --parameter-value {param_dict['tref']} --parameter phiref --parameter-value {param_dict['phi_ref']} \ + --parameter theta --parameter-value {param_dict['beta']} --parameter phi --parameter-value {param_dict['lambda']} \ + --parameter psi --parameter-value {param_dict['psi']} " + print(f"Executing command to create mdc.xml.gz\n{cmd}") + os.system(cmd) + os.system(f"mv mdc.xml.gz {param_dict['save_path']}/mdc.xml.gz") + os.system(f"ls {param_dict['save_path']}/*h5 | lal_path2cache > {param_dict['save_path']}/local.cache") + os.system(f" util_SimInspiralToCoinc.py --sim-xml {param_dict['save_path']}/mdc.xml.gz --event 0 --ifo A --ifo E --ifo T ; mv coinc.xml {param_dict['save_path']}/coinc.xml") + if get_snr and 'psd_path' in param_dict: + psd = load_psd(param_dict) + snr = calculate_snr(data_dict, param_dict['snr_fmin'], param_dict['snr_fmax'], 0.5/param_dict['deltaT'], psd) + create_PSD_injection_figure(data_dict, psd, param_dict["save_path"], snr) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/LISA/injections/__init__.py b/MonteCarloMarginalizeCode/Code/RIFT/LISA/injections/__init__.py new file mode 100644 index 000000000..0c84a5988 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/LISA/injections/__init__.py @@ -0,0 +1 @@ +"""Helpers for generating LISA injection data products.""" diff --git a/MonteCarloMarginalizeCode/Code/RIFT/LISA/injections/create_injections.py b/MonteCarloMarginalizeCode/Code/RIFT/LISA/injections/create_injections.py new file mode 100644 index 000000000..d6a80f87c --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/LISA/injections/create_injections.py @@ -0,0 +1,104 @@ +#!/usr/bin/env python + +"""Create LISA A/E/T injection frames from a RIFT injection XML file.""" + +from argparse import ArgumentParser +import ast +import os + +import lal +import lalsimulation +import numpy as np + +import RIFT.lalsimutils as lalsimutils +from RIFT.LISA.injections.LISA_injections import ( + generate_lisa_TDI_dict, + generate_lisa_injections, +) + +__author__ = "A. Jan" + + +def parse_args(argv=None): + parser = ArgumentParser() + parser.add_argument("--save-path", default=os.getcwd(), help="Path where h5 files should be written.") + parser.add_argument("--psd-path", default=None, help="Directory containing A/E/T PSD XML files for SNR.") + parser.add_argument("--inj", required=True, help="Inspiral XML file containing injection information.") + parser.add_argument("--fNyq", default=0.125, type=float, help="Nyquist frequency for generated waveforms.") + parser.add_argument("--deltaF", default=1 / (64 * 32768), type=float, help="Injection deltaF.") + parser.add_argument( + "--modes", + default="[(2,2),(2,1),(3,3),(3,2),(3,1),(4,4),(4,3),(4,2),(5,5)]", + help="List of modes to use in injection.", + ) + parser.add_argument("--path-to-NR-hdf5", default=None, help="NRHDF5 path when using NR injection data.") + parser.add_argument("--snr-fmin", default=0.0001, type=float, help="fmin while calculating SNR.") + parser.add_argument("--skip-snr", action="store_true", help="Write frames without PSD/SNR products.") + return parser.parse_args(argv) + + +def parameter_dict_from_xml(opts): + P_inj = lalsimutils.xml_to_ChooseWaveformParams_array(str(opts.inj))[0] + modes = np.array(ast.literal_eval(opts.modes)) + param_dict = { + "m1": P_inj.m1 / lal.MSUN_SI, + "m2": P_inj.m2 / lal.MSUN_SI, + "s1z": P_inj.s1z, + "s2z": P_inj.s2z, + "dist": P_inj.dist / (1e6 * lal.PC_SI), + "fmin": P_inj.fmin, + "fmax": opts.fNyq, + "deltaF": opts.deltaF, + "deltaT": 0.5 / opts.fNyq, + "fref": None, + "wf-fref": P_inj.fref, + "tref": float(P_inj.tref), + "beta": P_inj.theta, + "lambda": P_inj.phi, + "psi": P_inj.psi, + "phi_ref": P_inj.phiref, + "inclination": P_inj.incl, + "approx": lalsimulation.GetStringFromApproximant(P_inj.approx), + "modes": modes, + "save_path": opts.save_path or os.getcwd(), + "path_to_NR_hdf5": opts.path_to_NR_hdf5, + "snr_fmin": opts.snr_fmin, + "snr_fmax": opts.fNyq, + } + if opts.psd_path: + param_dict["psd_path"] = opts.psd_path + return param_dict + + +def main(argv=None): + opts = parse_args(argv) + param_dict = parameter_dict_from_xml(opts) + print(f"Saving frames in {param_dict['save_path']}") + print("###############") + if 1 / param_dict["deltaF"] / 60 / 60 / 24 > 0.5: + print(f"Data length = {1 / param_dict['deltaF'] / 60 / 60 / 24} days.") + else: + print(f"Data length = {1 / param_dict['deltaF'] / 60 / 60} hrs.") + print( + f"\nWaveform is being generated with m1 = {param_dict['m1']}, " + f"m2 = {param_dict['m2']}, s1z = {param_dict['s1z']}, s2z = {param_dict['s2z']}" + ) + print( + f"deltaF = {param_dict['deltaF']}, fmin = {param_dict['fmin']}, " + f"fmax = {param_dict['fmax']}, deltaT = {param_dict['deltaT']}, " + f"modes = {list(param_dict['modes'])}, tref = {param_dict['tref']}" + ) + print( + f"phiref = {param_dict['phi_ref']}, psi = {param_dict['psi']}, " + f"inclination = {param_dict['inclination']}, beta = {param_dict['beta']}, " + f"lambda = {param_dict['lambda']}" + ) + print(f"path_to_NR_hdf5 = {param_dict['path_to_NR_hdf5']}, approx = {param_dict['approx']}\n") + print("###############") + + data_dict = generate_lisa_TDI_dict(param_dict) + generate_lisa_injections(data_dict, param_dict, get_snr=not opts.skip_snr and opts.psd_path is not None) + + +if __name__ == "__main__": + main() diff --git a/MonteCarloMarginalizeCode/Code/RIFT/LISA/psd_generation/R.txt b/MonteCarloMarginalizeCode/Code/RIFT/LISA/psd_generation/R.txt new file mode 100644 index 000000000..dfbbf60b4 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/LISA/psd_generation/R.txt @@ -0,0 +1,899 @@ +1.023300e-05 1.500000e-01 +1.047143e-05 1.500000e-01 +1.071541e-05 1.500000e-01 +1.096508e-05 1.500000e-01 +1.122057e-05 1.500000e-01 +1.148201e-05 1.500000e-01 +1.174954e-05 1.500000e-01 +1.202330e-05 1.500000e-01 +1.230345e-05 1.500000e-01 +1.259012e-05 1.500000e-01 +1.288347e-05 1.500000e-01 +1.318365e-05 1.500000e-01 +1.349083e-05 1.500000e-01 +1.380517e-05 1.500000e-01 +1.412683e-05 1.500000e-01 +1.445598e-05 1.500000e-01 +1.479281e-05 1.500000e-01 +1.513748e-05 1.500000e-01 +1.549018e-05 1.500000e-01 +1.585110e-05 1.500000e-01 +1.622043e-05 1.500000e-01 +1.659837e-05 1.500000e-01 +1.698511e-05 1.500000e-01 +1.738086e-05 1.500000e-01 +1.778584e-05 1.500000e-01 +1.820025e-05 1.500000e-01 +1.862431e-05 1.500000e-01 +1.905826e-05 1.500000e-01 +1.950232e-05 1.500000e-01 +1.995672e-05 1.500000e-01 +2.042171e-05 1.500000e-01 +2.089754e-05 1.500000e-01 +2.138445e-05 1.500000e-01 +2.188271e-05 1.500000e-01 +2.239258e-05 1.500000e-01 +2.291433e-05 1.500000e-01 +2.344823e-05 1.500000e-01 +2.399457e-05 1.500000e-01 +2.455365e-05 1.500000e-01 +2.512575e-05 1.500000e-01 +2.571118e-05 1.500000e-01 +2.631025e-05 1.500000e-01 +2.692328e-05 1.500000e-01 +2.755059e-05 1.500000e-01 +2.819252e-05 1.500000e-01 +2.884940e-05 1.500000e-01 +2.952159e-05 1.500000e-01 +3.020945e-05 1.500000e-01 +3.091333e-05 1.500000e-01 +3.163361e-05 1.500000e-01 +3.237067e-05 1.500000e-01 +3.312491e-05 1.500000e-01 +3.389672e-05 1.500000e-01 +3.468651e-05 1.500000e-01 +3.549471e-05 1.500000e-01 +3.632173e-05 1.500000e-01 +3.716803e-05 1.500000e-01 +3.803404e-05 1.500000e-01 +3.892024e-05 1.500000e-01 +3.982708e-05 1.500000e-01 +4.075505e-05 1.500000e-01 +4.170464e-05 1.500000e-01 +4.267636e-05 1.500000e-01 +4.367072e-05 1.500000e-01 +4.468825e-05 1.500000e-01 +4.572948e-05 1.500000e-01 +4.679498e-05 1.500000e-01 +4.788530e-05 1.500000e-01 +4.900103e-05 1.500000e-01 +5.014275e-05 1.500000e-01 +5.131108e-05 1.500000e-01 +5.250663e-05 1.500000e-01 +5.373003e-05 1.500000e-01 +5.498193e-05 1.500000e-01 +5.626301e-05 1.500000e-01 +5.757396e-05 1.500000e-01 +5.891541e-05 1.500000e-01 +6.028813e-05 1.500000e-01 +6.169287e-05 1.500000e-01 +6.313031e-05 1.500000e-01 +6.460123e-05 1.500000e-01 +6.610646e-05 1.500000e-01 +6.764675e-05 1.500000e-01 +6.922293e-05 1.500000e-01 +7.083579e-05 1.500000e-01 +7.248627e-05 1.500000e-01 +7.417520e-05 1.500000e-01 +7.590349e-05 1.500000e-01 +7.767202e-05 1.500000e-01 +7.948178e-05 1.500000e-01 +8.133372e-05 1.500000e-01 +8.322878e-05 1.500000e-01 +8.516802e-05 1.500000e-01 +8.715241e-05 1.500000e-01 +8.918308e-05 1.500000e-01 +9.126106e-05 1.500000e-01 +9.338745e-05 1.500000e-01 +9.556335e-05 1.500000e-01 +9.778998e-05 1.500000e-01 +1.000685e-04 1.500000e-01 +1.024001e-04 1.500000e-01 +1.047860e-04 1.500000e-01 +1.072275e-04 1.500000e-01 +1.097259e-04 1.500000e-01 +1.122826e-04 1.500000e-01 +1.148987e-04 1.500000e-01 +1.175759e-04 1.500000e-01 +1.203154e-04 1.500000e-01 +1.231187e-04 1.500000e-01 +1.259874e-04 1.500000e-01 +1.289229e-04 1.500000e-01 +1.319268e-04 1.500000e-01 +1.350007e-04 1.500000e-01 +1.381462e-04 1.500000e-01 +1.413650e-04 1.500000e-01 +1.446588e-04 1.500000e-01 +1.480294e-04 1.500000e-01 +1.514785e-04 1.500000e-01 +1.550079e-04 1.500000e-01 +1.586196e-04 1.500000e-01 +1.623155e-04 1.500000e-01 +1.660974e-04 1.500000e-01 +1.699675e-04 1.500000e-01 +1.739277e-04 1.500000e-01 +1.779802e-04 1.500000e-01 +1.821272e-04 1.500000e-01 +1.863707e-04 1.500000e-01 +1.907132e-04 1.500000e-01 +1.951568e-04 1.500000e-01 +1.997039e-04 1.500000e-01 +2.043570e-04 1.500000e-01 +2.091186e-04 1.500000e-01 +2.139910e-04 1.500000e-01 +2.189770e-04 1.500000e-01 +2.240792e-04 1.500000e-01 +2.293002e-04 1.500000e-01 +2.346429e-04 1.500000e-01 +2.401101e-04 1.500000e-01 +2.457047e-04 1.500000e-01 +2.514296e-04 1.500000e-01 +2.572879e-04 1.500000e-01 +2.632827e-04 1.500000e-01 +2.694172e-04 1.500000e-01 +2.756946e-04 1.500000e-01 +2.821183e-04 1.500000e-01 +2.886916e-04 1.500000e-01 +2.954182e-04 1.500000e-01 +3.023014e-04 1.500000e-01 +3.093450e-04 1.500000e-01 +3.165528e-04 1.500000e-01 +3.239285e-04 1.500000e-01 +3.314760e-04 1.500000e-01 +3.391994e-04 1.500000e-01 +3.471027e-04 1.500000e-01 +3.551902e-04 1.500000e-01 +3.634661e-04 1.500000e-01 +3.719349e-04 1.500000e-01 +3.806010e-04 1.500000e-01 +3.894690e-04 1.500000e-01 +3.985436e-04 1.500000e-01 +4.078297e-04 1.500000e-01 +4.173321e-04 1.500000e-01 +4.270560e-04 1.500000e-01 +4.370064e-04 1.500000e-01 +4.471886e-04 1.500000e-01 +4.576081e-04 1.500000e-01 +4.682704e-04 1.500000e-01 +4.791811e-04 1.500000e-01 +4.903460e-04 1.500000e-01 +5.017711e-04 1.500000e-01 +5.134623e-04 1.500000e-01 +5.254263e-04 1.500000e-01 +5.376686e-04 1.500000e-01 +5.501960e-04 1.500000e-01 +5.630158e-04 1.500000e-01 +5.761336e-04 1.500000e-01 +5.895581e-04 1.500000e-01 +6.032948e-04 1.500000e-01 +6.173511e-04 1.500000e-01 +6.317354e-04 1.500000e-01 +6.464550e-04 1.500000e-01 +6.615173e-04 1.500000e-01 +6.769307e-04 1.500000e-01 +6.927035e-04 1.500000e-01 +7.088431e-04 1.500000e-01 +7.253594e-04 1.500000e-01 +7.422603e-04 1.500000e-01 +7.595547e-04 1.500000e-01 +7.772525e-04 1.500000e-01 +7.953627e-04 1.500000e-01 +8.138942e-04 1.500000e-01 +8.328579e-04 1.500000e-01 +8.522639e-04 1.500000e-01 +8.721215e-04 1.500000e-01 +8.924417e-04 1.500000e-01 +9.132357e-04 1.500000e-01 +9.345142e-04 1.500000e-01 +9.562885e-04 1.500000e-01 +9.785700e-04 1.500000e-01 +1.001371e-03 1.499459e-01 +1.024703e-03 1.501543e-01 +1.048578e-03 1.498621e-01 +1.073010e-03 1.499323e-01 +1.098011e-03 1.499193e-01 +1.123595e-03 1.500063e-01 +1.149774e-03 1.499914e-01 +1.176564e-03 1.498757e-01 +1.203978e-03 1.500625e-01 +1.232031e-03 1.500393e-01 +1.260737e-03 1.498941e-01 +1.290112e-03 1.499813e-01 +1.320172e-03 1.499910e-01 +1.350932e-03 1.499849e-01 +1.382409e-03 1.500440e-01 +1.414619e-03 1.499732e-01 +1.447579e-03 1.500271e-01 +1.481308e-03 1.500182e-01 +1.515823e-03 1.500439e-01 +1.551141e-03 1.499906e-01 +1.587283e-03 1.500412e-01 +1.624266e-03 1.499937e-01 +1.662112e-03 1.499962e-01 +1.700839e-03 1.500501e-01 +1.740469e-03 1.500018e-01 +1.781022e-03 1.500179e-01 +1.822519e-03 1.500075e-01 +1.864984e-03 1.500124e-01 +1.908438e-03 1.500182e-01 +1.952905e-03 1.499975e-01 +1.998407e-03 1.499926e-01 +2.044970e-03 1.499962e-01 +2.092618e-03 1.500017e-01 +2.141376e-03 1.500031e-01 +2.191270e-03 1.500008e-01 +2.242327e-03 1.500006e-01 +2.294573e-03 1.499953e-01 +2.348037e-03 1.499949e-01 +2.402746e-03 1.499914e-01 +2.458730e-03 1.499924e-01 +2.516018e-03 1.499944e-01 +2.574641e-03 1.500025e-01 +2.634631e-03 1.499996e-01 +2.696018e-03 1.499951e-01 +2.758835e-03 1.499999e-01 +2.823116e-03 1.500005e-01 +2.888894e-03 1.499993e-01 +2.956205e-03 1.499994e-01 +3.025085e-03 1.500028e-01 +3.095569e-03 1.499984e-01 +3.167696e-03 1.500010e-01 +3.241504e-03 1.499992e-01 +3.317031e-03 1.500009e-01 +3.394317e-03 1.500014e-01 +3.473405e-03 1.499993e-01 +3.554335e-03 1.500013e-01 +3.637151e-03 1.499991e-01 +3.721897e-03 1.499988e-01 +3.808617e-03 1.500003e-01 +3.897358e-03 1.499996e-01 +3.988166e-03 1.499991e-01 +4.081091e-03 1.499988e-01 +4.176180e-03 1.499991e-01 +4.273485e-03 1.499991e-01 +4.373057e-03 1.499994e-01 +4.474950e-03 1.499987e-01 +4.579216e-03 1.499992e-01 +4.685912e-03 1.499992e-01 +4.795093e-03 1.499987e-01 +4.906819e-03 1.499988e-01 +5.021148e-03 1.499986e-01 +5.138141e-03 1.499987e-01 +5.257857e-03 1.499981e-01 +5.380370e-03 1.499988e-01 +5.505732e-03 1.499982e-01 +5.634014e-03 1.499985e-01 +5.765287e-03 1.499985e-01 +5.899615e-03 1.499981e-01 +6.037076e-03 1.499981e-01 +6.177744e-03 1.499981e-01 +6.321682e-03 1.499980e-01 +6.468978e-03 1.499979e-01 +6.619706e-03 1.499979e-01 +6.773944e-03 1.499977e-01 +6.931777e-03 1.499977e-01 +7.093288e-03 1.499975e-01 +7.258561e-03 1.499974e-01 +7.427685e-03 1.499972e-01 +7.600750e-03 1.499972e-01 +7.777849e-03 1.499970e-01 +7.959071e-03 1.499969e-01 +8.144522e-03 1.499967e-01 +8.334285e-03 1.499965e-01 +8.528475e-03 1.499964e-01 +8.727188e-03 1.499962e-01 +8.930532e-03 1.499959e-01 +9.138613e-03 1.499958e-01 +9.351545e-03 1.499956e-01 +9.569435e-03 1.499953e-01 +9.792401e-03 1.499952e-01 +1.002057e-02 1.499949e-01 +1.025404e-02 1.499946e-01 +1.049296e-02 1.499944e-01 +1.073745e-02 1.499942e-01 +1.098763e-02 1.499939e-01 +1.124364e-02 1.499936e-01 +1.150562e-02 1.499934e-01 +1.177370e-02 1.499930e-01 +1.204803e-02 1.499927e-01 +1.232875e-02 1.499923e-01 +1.261601e-02 1.499920e-01 +1.290996e-02 1.499916e-01 +1.321076e-02 1.499913e-01 +1.351857e-02 1.499908e-01 +1.383356e-02 1.499903e-01 +1.415588e-02 1.499899e-01 +1.448571e-02 1.499894e-01 +1.482323e-02 1.499890e-01 +1.516861e-02 1.499884e-01 +1.552204e-02 1.499879e-01 +1.588370e-02 1.499873e-01 +1.625379e-02 1.499867e-01 +1.663251e-02 1.499860e-01 +1.702004e-02 1.499855e-01 +1.741661e-02 1.499848e-01 +1.782242e-02 1.499841e-01 +1.823768e-02 1.499833e-01 +1.866262e-02 1.499824e-01 +1.909746e-02 1.499816e-01 +1.954243e-02 1.499808e-01 +1.999777e-02 1.499799e-01 +2.046371e-02 1.499790e-01 +2.094052e-02 1.499779e-01 +2.142843e-02 1.499769e-01 +2.192771e-02 1.499758e-01 +2.243863e-02 1.499747e-01 +2.296145e-02 1.499735e-01 +2.349645e-02 1.499722e-01 +2.404392e-02 1.499709e-01 +2.460414e-02 1.499696e-01 +2.517742e-02 1.499682e-01 +2.576405e-02 1.499666e-01 +2.636436e-02 1.499650e-01 +2.697864e-02 1.499634e-01 +2.760725e-02 1.499617e-01 +2.825050e-02 1.499598e-01 +2.890873e-02 1.499579e-01 +2.958231e-02 1.499560e-01 +3.027157e-02 1.499539e-01 +3.097690e-02 1.499518e-01 +3.169866e-02 1.499495e-01 +3.243724e-02 1.499470e-01 +3.319303e-02 1.499446e-01 +3.396643e-02 1.499420e-01 +3.475784e-02 1.499392e-01 +3.556770e-02 1.499363e-01 +3.639643e-02 1.499333e-01 +3.724447e-02 1.499303e-01 +3.811226e-02 1.499269e-01 +3.900028e-02 1.499234e-01 +3.990899e-02 1.499198e-01 +4.083886e-02 1.499161e-01 +4.179041e-02 1.499122e-01 +4.276413e-02 1.499080e-01 +4.376053e-02 1.499037e-01 +4.478015e-02 1.498992e-01 +4.582353e-02 1.498944e-01 +4.689122e-02 1.498894e-01 +4.798378e-02 1.498842e-01 +4.910180e-02 1.498788e-01 +5.024588e-02 1.498731e-01 +5.141661e-02 1.498670e-01 +5.261462e-02 1.498608e-01 +5.384053e-02 1.498543e-01 +5.509500e-02 1.498474e-01 +5.637876e-02 1.498402e-01 +5.769238e-02 1.498327e-01 +5.903660e-02 1.498248e-01 +6.041216e-02 1.498166e-01 +6.181973e-02 1.498079e-01 +6.326015e-02 1.497988e-01 +6.473410e-02 1.497893e-01 +6.624243e-02 1.497795e-01 +6.778586e-02 1.497690e-01 +6.936529e-02 1.497581e-01 +7.098150e-02 1.497468e-01 +7.263533e-02 1.497348e-01 +7.432773e-02 1.497223e-01 +7.605958e-02 1.497092e-01 +7.783177e-02 1.496956e-01 +7.964525e-02 1.496812e-01 +8.150097e-02 1.496662e-01 +8.339996e-02 1.496505e-01 +8.534318e-02 1.496341e-01 +8.733166e-02 1.496168e-01 +8.936652e-02 1.495988e-01 +9.144874e-02 1.495798e-01 +9.357948e-02 1.495600e-01 +9.575989e-02 1.495394e-01 +9.799113e-02 1.495177e-01 +1.002743e-01 1.494950e-01 +1.026107e-01 1.494712e-01 +1.050015e-01 1.494463e-01 +1.074481e-01 1.494203e-01 +1.099516e-01 1.493930e-01 +1.125135e-01 1.493644e-01 +1.151350e-01 1.493345e-01 +1.178177e-01 1.493032e-01 +1.205628e-01 1.492703e-01 +1.233719e-01 1.492361e-01 +1.262465e-01 1.492001e-01 +1.291881e-01 1.491625e-01 +1.321981e-01 1.491231e-01 +1.352784e-01 1.490819e-01 +1.384303e-01 1.490387e-01 +1.416558e-01 1.489935e-01 +1.449563e-01 1.489463e-01 +1.483338e-01 1.488967e-01 +1.517900e-01 1.488448e-01 +1.553267e-01 1.487906e-01 +1.589458e-01 1.487338e-01 +1.626493e-01 1.486743e-01 +1.664390e-01 1.486121e-01 +1.703170e-01 1.485469e-01 +1.742854e-01 1.484787e-01 +1.783463e-01 1.484073e-01 +1.825017e-01 1.483325e-01 +1.867540e-01 1.482543e-01 +1.911054e-01 1.481725e-01 +1.955581e-01 1.480867e-01 +2.001146e-01 1.479971e-01 +2.047773e-01 1.479032e-01 +2.095486e-01 1.478049e-01 +2.144311e-01 1.477021e-01 +2.194274e-01 1.475944e-01 +2.245400e-01 1.474819e-01 +2.297718e-01 1.473640e-01 +2.351255e-01 1.472407e-01 +2.406039e-01 1.471117e-01 +2.462100e-01 1.469767e-01 +2.519467e-01 1.468354e-01 +2.578170e-01 1.466877e-01 +2.638242e-01 1.465329e-01 +2.699713e-01 1.463712e-01 +2.762616e-01 1.462020e-01 +2.826985e-01 1.460250e-01 +2.892854e-01 1.458396e-01 +2.960257e-01 1.456460e-01 +3.029231e-01 1.454434e-01 +3.099812e-01 1.452314e-01 +3.172038e-01 1.450098e-01 +3.245946e-01 1.447779e-01 +3.321577e-01 1.445356e-01 +3.398970e-01 1.442821e-01 +3.478166e-01 1.440171e-01 +3.559207e-01 1.437401e-01 +3.642136e-01 1.434504e-01 +3.726998e-01 1.431476e-01 +3.813837e-01 1.428312e-01 +3.902700e-01 1.425004e-01 +3.993633e-01 1.421546e-01 +4.086684e-01 1.417934e-01 +4.181904e-01 1.414159e-01 +4.279342e-01 1.410216e-01 +4.379051e-01 1.406097e-01 +4.481083e-01 1.401794e-01 +4.585492e-01 1.397299e-01 +4.692334e-01 1.392607e-01 +4.801665e-01 1.387707e-01 +4.913544e-01 1.382591e-01 +5.028030e-01 1.377251e-01 +5.145183e-01 1.371677e-01 +5.265067e-01 1.365862e-01 +5.387742e-01 1.359793e-01 +5.513278e-01 1.353464e-01 +5.641737e-01 1.346862e-01 +5.773188e-01 1.339977e-01 +5.907705e-01 1.332800e-01 +6.045355e-01 1.325319e-01 +6.186211e-01 1.317522e-01 +6.330348e-01 1.309399e-01 +6.477843e-01 1.300938e-01 +6.628781e-01 1.292127e-01 +6.783229e-01 1.282953e-01 +6.941282e-01 1.273406e-01 +7.103007e-01 1.263473e-01 +7.268511e-01 1.253141e-01 +7.437866e-01 1.242397e-01 +7.611171e-01 1.231229e-01 +7.788511e-01 1.219626e-01 +7.969980e-01 1.207572e-01 +8.155682e-01 1.195058e-01 +8.345713e-01 1.182070e-01 +8.540165e-01 1.168597e-01 +8.739150e-01 1.154626e-01 +8.942772e-01 1.140147e-01 +9.151141e-01 1.125148e-01 +9.364361e-01 1.109619e-01 +9.582549e-01 1.093552e-01 +9.805825e-01 1.076937e-01 +1.003430e+00 1.059764e-01 +1.026810e+00 1.042030e-01 +1.050735e+00 1.023727e-01 +1.075217e+00 1.004851e-01 +1.100269e+00 9.853989e-02 +1.125905e+00 9.653700e-02 +1.152139e+00 9.447655e-02 +1.178984e+00 9.235875e-02 +1.206454e+00 9.018412e-02 +1.234565e+00 8.795341e-02 +1.263330e+00 8.566764e-02 +1.292766e+00 8.332812e-02 +1.322887e+00 8.093650e-02 +1.353710e+00 7.849472e-02 +1.385252e+00 7.600516e-02 +1.417528e+00 7.347058e-02 +1.450557e+00 7.089407e-02 +1.484354e+00 6.827922e-02 +1.518940e+00 6.563005e-02 +1.554331e+00 6.295115e-02 +1.590547e+00 6.024752e-02 +1.627607e+00 5.752470e-02 +1.665530e+00 5.478880e-02 +1.704337e+00 5.204640e-02 +1.744048e+00 4.930470e-02 +1.784684e+00 4.657143e-02 +1.826268e+00 4.385477e-02 +1.868820e+00 4.116355e-02 +1.912363e+00 3.850706e-02 +1.956921e+00 3.589496e-02 +2.002517e+00 3.333748e-02 +2.049176e+00 3.084513e-02 +2.096922e+00 2.842870e-02 +2.145780e+00 2.609922e-02 +2.195777e+00 2.386788e-02 +2.246938e+00 2.174577e-02 +2.299292e+00 1.974391e-02 +2.352866e+00 1.787301e-02 +2.407687e+00 1.614333e-02 +2.463786e+00 1.456451e-02 +2.521193e+00 1.314535e-02 +2.579936e+00 1.189361e-02 +2.640049e+00 1.081577e-02 +2.701562e+00 9.916862e-03 +2.764508e+00 9.200107e-03 +2.828922e+00 8.666817e-03 +2.894835e+00 8.316109e-03 +2.962285e+00 8.144684e-03 +3.031306e+00 8.146610e-03 +3.101936e+00 8.313258e-03 +3.174211e+00 8.633109e-03 +3.248170e+00 9.091687e-03 +3.323852e+00 9.671654e-03 +3.401298e+00 1.035279e-02 +3.480548e+00 1.111224e-02 +3.561645e+00 1.192476e-02 +3.644631e+00 1.276309e-02 +3.729551e+00 1.359853e-02 +3.816450e+00 1.440146e-02 +3.905373e+00 1.514214e-02 +3.996368e+00 1.579155e-02 +4.089484e+00 1.632234e-02 +4.184769e+00 1.670979e-02 +4.282274e+00 1.693294e-02 +4.382051e+00 1.697556e-02 +4.484153e+00 1.682720e-02 +4.588633e+00 1.648407e-02 +4.695548e+00 1.594971e-02 +4.804955e+00 1.523549e-02 +4.916910e+00 1.436075e-02 +5.031474e+00 1.335253e-02 +5.148708e+00 1.224491e-02 +5.268671e+00 1.107786e-02 +5.391436e+00 9.895565e-03 +5.517055e+00 8.744306e-03 +5.645599e+00 7.669990e-03 +5.777144e+00 6.715233e-03 +5.911750e+00 5.916486e-03 +6.049494e+00 5.301062e-03 +6.190445e+00 4.884625e-03 +6.334686e+00 4.669149e-03 +6.482281e+00 4.641825e-03 +6.633318e+00 4.775061e-03 +6.787876e+00 5.027824e-03 +6.946034e+00 5.348403e-03 +7.107875e+00 5.678547e-03 +7.273489e+00 5.958835e-03 +7.442964e+00 6.134753e-03 +7.616385e+00 6.162998e-03 +7.793845e+00 6.017228e-03 +7.975440e+00 5.692419e-03 +8.161268e+00 5.206925e-03 +8.351429e+00 4.601707e-03 +8.546012e+00 3.936231e-03 +8.745139e+00 3.281126e-03 +8.948897e+00 2.708240e-03 +9.157407e+00 2.279307e-03 +9.370775e+00 2.034981e-03 +9.589114e+00 1.986314e-03 +9.812542e+00 2.110701e-03 +1.004117e+01 2.353891e-03 +1.027317e+01 2.636361e-03 +1.051053e+01 2.876134e-03 +1.075338e+01 2.998850e-03 +1.100184e+01 2.956616e-03 +1.125603e+01 2.740923e-03 +1.151610e+01 2.387145e-03 +1.178218e+01 1.968075e-03 +1.205441e+01 1.576475e-03 +1.233292e+01 1.299934e-03 +1.261787e+01 1.194435e-03 +1.290941e+01 1.264840e-03 +1.320768e+01 1.459919e-03 +1.351284e+01 1.686029e-03 +1.382505e+01 1.837518e-03 +1.414448e+01 1.835049e-03 +1.447129e+01 1.657986e-03 +1.480564e+01 1.356479e-03 +1.514773e+01 1.034789e-03 +1.549771e+01 8.087459e-04 +1.585579e+01 7.531846e-04 +1.622213e+01 8.635563e-04 +1.659695e+01 1.053584e-03 +1.698042e+01 1.195855e-03 +1.737275e+01 1.189383e-03 +1.777414e+01 1.018974e-03 +1.818481e+01 7.690436e-04 +1.860497e+01 5.761498e-04 +1.903484e+01 5.423352e-04 +1.947463e+01 6.632118e-04 +1.992459e+01 8.248612e-04 +2.019389e+01 8.783485e-04 +2.046682e+01 8.745669e-04 +2.074344e+01 8.101724e-04 +2.102379e+01 6.997052e-04 +2.130794e+01 5.718728e-04 +2.159593e+01 4.611706e-04 +2.188781e+01 3.971086e-04 +2.218363e+01 3.943827e-04 +2.248346e+01 4.474117e-04 +2.278733e+01 5.315438e-04 +2.309531e+01 6.111163e-04 +2.340746e+01 6.520404e-04 +2.372382e+01 6.345889e-04 +2.404446e+01 5.614684e-04 +2.436943e+01 4.575385e-04 +2.469880e+01 3.605398e-04 +2.503262e+01 3.058727e-04 +2.537094e+01 3.114110e-04 +2.571385e+01 3.690265e-04 +2.606138e+01 4.472590e-04 +2.641361e+01 5.049100e-04 +2.677061e+01 5.101017e-04 +2.713242e+01 4.560236e-04 +2.749913e+01 3.652040e-04 +2.787080e+01 2.790962e-04 +2.824749e+01 2.371246e-04 +2.862926e+01 2.555412e-04 +2.901620e+01 3.177505e-04 +2.940837e+01 3.825940e-04 +2.980584e+01 4.073586e-04 +3.009166e+01 3.887271e-04 +3.038022e+01 3.438814e-04 +3.067155e+01 2.857917e-04 +3.096567e+01 2.319242e-04 +3.126261e+01 1.985148e-04 +3.156240e+01 1.949993e-04 +3.186507e+01 2.205519e-04 +3.217064e+01 2.640962e-04 +3.247913e+01 3.080108e-04 +3.279059e+01 3.344060e-04 +3.310503e+01 3.317895e-04 +3.342249e+01 2.995980e-04 +3.374299e+01 2.487022e-04 +3.406657e+01 1.974218e-04 +3.439324e+01 1.643695e-04 +3.472305e+01 1.608569e-04 +3.505603e+01 1.860106e-04 +3.539219e+01 2.268478e-04 +3.573158e+01 2.635958e-04 +3.607423e+01 2.782362e-04 +3.642016e+01 2.626335e-04 +3.676941e+01 2.224969e-04 +3.712200e+01 1.750612e-04 +3.747798e+01 1.411931e-04 +3.783737e+01 1.353931e-04 +3.820021e+01 1.584909e-04 +3.856653e+01 1.968728e-04 +3.893636e+01 2.290488e-04 +3.930973e+01 2.366041e-04 +3.968669e+01 2.140042e-04 +3.998188e+01 1.819763e-04 +4.027927e+01 1.479227e-04 +4.057888e+01 1.225620e-04 +4.088071e+01 1.136675e-04 +4.118478e+01 1.232799e-04 +4.149112e+01 1.468700e-04 +4.179974e+01 1.748518e-04 +4.211065e+01 1.960054e-04 +4.242387e+01 2.016272e-04 +4.273942e+01 1.888691e-04 +4.305733e+01 1.619464e-04 +4.337759e+01 1.306355e-04 +4.370024e+01 1.065223e-04 +4.402529e+01 9.839557e-05 +4.435275e+01 1.086158e-04 +4.468265e+01 1.320092e-04 +4.501501e+01 1.578725e-04 +4.534983e+01 1.743940e-04 +4.568715e+01 1.737045e-04 +4.602698e+01 1.553844e-04 +4.636933e+01 1.267961e-04 +4.671423e+01 9.993883e-05 +4.706170e+01 8.608552e-05 +4.741175e+01 9.057472e-05 +4.776440e+01 1.102055e-04 +4.811968e+01 1.345998e-04 +4.847760e+01 1.510683e-04 +4.883818e+01 1.507680e-04 +4.920145e+01 1.331575e-04 +4.956741e+01 1.064256e-04 +4.986865e+01 8.680727e-05 +5.017172e+01 7.617916e-05 +5.047664e+01 7.781201e-05 +5.078340e+01 9.059368e-05 +5.109203e+01 1.093575e-04 +5.140254e+01 1.267127e-04 +5.171493e+01 1.357541e-04 +5.202922e+01 1.326673e-04 +5.234543e+01 1.182317e-04 +5.266356e+01 9.759545e-05 +5.298359e+01 7.834212e-05 +5.330562e+01 6.754330e-05 +5.362958e+01 6.893133e-05 +5.395549e+01 8.134359e-05 +5.428338e+01 9.912163e-05 +5.461332e+01 1.143806e-04 +5.494520e+01 1.202724e-04 +5.527912e+01 1.139003e-04 +5.561508e+01 9.762130e-05 +5.595309e+01 7.810841e-05 +5.629314e+01 6.349379e-05 +5.663523e+01 5.978312e-05 +5.697942e+01 6.812534e-05 +5.732571e+01 8.420925e-05 +5.767409e+01 1.001511e-04 +5.802462e+01 1.080957e-04 +5.837725e+01 1.039042e-04 +5.873202e+01 8.917943e-05 +5.908900e+01 7.058682e-05 +5.944807e+01 5.669546e-05 +5.975354e+01 5.341806e-05 +6.006058e+01 5.892350e-05 +6.036919e+01 7.084252e-05 +6.067938e+01 8.444512e-05 +6.099119e+01 9.436675e-05 +6.130457e+01 9.661351e-05 +6.161957e+01 9.009443e-05 +6.193620e+01 7.707895e-05 +6.225446e+01 6.236464e-05 +6.257433e+01 5.143819e-05 +6.289589e+01 4.831915e-05 +6.321907e+01 5.393741e-05 +6.354393e+01 6.571449e-05 +6.387041e+01 7.856045e-05 +6.419861e+01 8.692548e-05 +6.452845e+01 8.709056e-05 +6.486006e+01 7.874823e-05 +6.519330e+01 6.517900e-05 +6.552832e+01 5.190111e-05 +6.586502e+01 4.432575e-05 +6.620345e+01 4.540303e-05 +6.654366e+01 5.429498e-05 +6.688554e+01 6.669922e-05 +6.722921e+01 7.673283e-05 +6.757465e+01 7.957552e-05 +6.792194e+01 7.368772e-05 +6.827089e+01 6.156549e-05 +6.862169e+01 4.862937e-05 +6.897431e+01 4.070960e-05 +6.932872e+01 4.127708e-05 +6.963733e+01 4.828538e-05 +6.994726e+01 5.850321e-05 +7.025859e+01 6.794498e-05 +7.057135e+01 7.289987e-05 +7.088546e+01 7.133557e-05 +7.120099e+01 6.369228e-05 +7.151788e+01 5.274286e-05 +7.183624e+01 4.254775e-05 +7.215596e+01 3.688476e-05 +7.247715e+01 3.775839e-05 +7.279975e+01 4.458373e-05 +7.312377e+01 5.439192e-05 +7.344926e+01 6.300346e-05 +7.377615e+01 6.671992e-05 +7.410457e+01 6.385066e-05 +7.443441e+01 5.543878e-05 +7.476571e+01 4.485563e-05 +7.509848e+01 3.640478e-05 +7.543276e+01 3.349968e-05 +7.576852e+01 3.718614e-05 +7.610574e+01 4.566012e-05 +7.644453e+01 5.502062e-05 +7.678479e+01 6.096179e-05 +7.712652e+01 6.068032e-05 +7.746982e+01 5.414837e-05 +7.781464e+01 4.415872e-05 +7.816103e+01 3.508549e-05 +7.850889e+01 3.089364e-05 +7.885837e+01 3.329783e-05 +7.920937e+01 4.093647e-05 +7.952034e+01 4.899495e-05 +7.983257e+01 5.500650e-05 +8.014600e+01 5.658556e-05 +8.046064e+01 5.302514e-05 +8.077654e+01 4.558480e-05 +8.109369e+01 3.702562e-05 +8.141205e+01 3.055404e-05 +8.173167e+01 2.857191e-05 +8.205259e+01 3.172097e-05 +8.237473e+01 3.860839e-05 +8.269811e+01 4.634181e-05 +8.302281e+01 5.167823e-05 +8.334877e+01 5.233652e-05 +8.367604e+01 4.794071e-05 +8.400451e+01 4.019433e-05 +8.433434e+01 3.219552e-05 +8.466543e+01 2.715530e-05 +8.499783e+01 2.704590e-05 +8.533154e+01 3.175122e-05 +8.566656e+01 3.909607e-05 +8.600290e+01 4.576396e-05 +8.634059e+01 4.872780e-05 +8.667954e+01 4.657830e-05 +8.701985e+01 4.016285e-05 +8.736153e+01 3.223154e-05 +8.770451e+01 2.622425e-05 +8.804886e+01 2.472457e-05 +8.839452e+01 2.827511e-05 +8.874159e+01 3.510678e-05 +8.908997e+01 4.193060e-05 +8.940288e+01 4.530638e-05 +8.971684e+01 4.465962e-05 +9.003195e+01 4.016015e-05 +9.034811e+01 3.345320e-05 +9.066542e+01 2.703868e-05 +9.098388e+01 2.330853e-05 +9.130339e+01 2.361072e-05 +9.162406e+01 2.770932e-05 +9.194582e+01 3.386416e-05 +9.226874e+01 3.951195e-05 +9.259281e+01 4.228944e-05 +9.291798e+01 4.099184e-05 +9.324430e+01 3.607440e-05 +9.357178e+01 2.948114e-05 +9.390041e+01 2.385534e-05 +9.423019e+01 2.144663e-05 +9.456112e+01 2.316504e-05 +9.489326e+01 2.818425e-05 +9.522650e+01 3.427528e-05 +9.556094e+01 3.874477e-05 +9.589654e+01 3.958841e-05 +9.623334e+01 3.636462e-05 +9.657130e+01 3.040146e-05 +9.691046e+01 2.423060e-05 +9.725082e+01 2.048053e-05 +9.759239e+01 2.070927e-05 +9.793512e+01 2.469782e-05 +9.827905e+01 3.053341e-05 +9.862423e+01 3.546296e-05 +9.897062e+01 3.714346e-05 +9.928505e+01 3.511593e-05 +9.960048e+01 3.042770e-05 +9.991690e+01 2.482488e-05 +1.002343e+02 2.041280e-05 +1.005528e+02 1.883992e-05 +1.008722e+02 2.064828e-05 +1.011927e+02 2.505010e-05 +1.015142e+02 3.022811e-05 +1.018367e+02 3.405068e-05 +1.021602e+02 3.492370e-05 +1.024848e+02 3.243802e-05 +1.028104e+02 2.754566e-05 +1.031370e+02 2.218252e-05 +1.034647e+02 1.848157e-05 +1.037934e+02 1.789295e-05 +1.041232e+02 2.057566e-05 +1.044540e+02 2.532073e-05 +1.047858e+02 3.004945e-05 +1.051187e+02 3.268731e-05 +1.054527e+02 3.204871e-05 +1.057877e+02 2.835301e-05 +1.061238e+02 2.314142e-05 +1.064609e+02 1.861870e-05 +1.067992e+02 1.669494e-05 +1.071385e+02 1.814313e-05 +1.074788e+02 2.224716e-05 +1.078203e+02 2.710748e-05 +1.081629e+02 3.048094e-05 +1.085065e+02 3.078786e-05 diff --git a/MonteCarloMarginalizeCode/Code/RIFT/LISA/psd_generation/__init__.py b/MonteCarloMarginalizeCode/Code/RIFT/LISA/psd_generation/__init__.py new file mode 100644 index 000000000..6126f61ad --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/LISA/psd_generation/__init__.py @@ -0,0 +1 @@ +"""Analytic PSD generation helpers for LISA.""" diff --git a/MonteCarloMarginalizeCode/Code/RIFT/LISA/psd_generation/generate_LISA_psd.py b/MonteCarloMarginalizeCode/Code/RIFT/LISA/psd_generation/generate_LISA_psd.py new file mode 100644 index 000000000..f573392ae --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/LISA/psd_generation/generate_LISA_psd.py @@ -0,0 +1,174 @@ +#!/usr/bin/env python + +"""Generate an analytic LISA PSD text file and optional diagnostic products.""" + +from argparse import ArgumentParser +import os +import subprocess + +import matplotlib.pyplot as plt +import numpy as np +from scipy import interpolate + +########################################################################################### +# CONSTANTS +########################################################################################### +fm = 3.168753575e-8 +YRSID_SI = 31558149.763545603 +C_SI = 299792458.0 +e = 0.004824185218078991 +a = 149597870700.0 +Larm = 2 * np.sqrt(3) * a * e +fstar = C_SI / (2 * np.pi * Larm) + +path_to_file = os.path.dirname(__file__) + + +########################################################################################### +# FUNCTIONS +# These functions were taken from LISA.py of LISA sensitivity +# (https://github.com/eXtremeGravityInstitute/LISA_Sensitivity) +########################################################################################### +def Pn(f): + """Calculate the strain power spectral density.""" + P_oms = (1.5e-11) ** 2 * (1.0 + (2.0e-3 / f) ** 4) + P_acc = (3.0e-15) ** 2 * (1.0 + (0.4e-3 / f) ** 2) * (1.0 + (f / (8.0e-3)) ** 4) + return (P_oms + 2.0 * (1.0 + np.cos(f / fstar) ** 2) * P_acc / (2.0 * np.pi * f) ** 4) / Larm**2 + + +def SnC(f, Tobs=0.5, NC=3): + """ + Estimate galactic binary confusion noise. + + Tobs is provided in seconds. Supported fitted regimes correspond to roughly + 0.5 yr, 1 yr, 2 yr, and 4 yr observations. + """ + if Tobs < 0.75 * YRSID_SI: + est = 1 + elif 0.75 * YRSID_SI < Tobs and Tobs < 1.5 * YRSID_SI: + est = 2 + elif 1.5 * YRSID_SI < Tobs and Tobs < 3.0 * YRSID_SI: + est = 3 + else: + est = 4 + + if est == 1: + alpha = 0.133 + beta = 243.0 + kappa = 482.0 + gamma = 917.0 + f_knee = 2.58e-3 + elif est == 2: + alpha = 0.171 + beta = 292.0 + kappa = 1020.0 + gamma = 1680.0 + f_knee = 2.15e-3 + elif est == 3: + alpha = 0.165 + beta = 299.0 + kappa = 611.0 + gamma = 1340.0 + f_knee = 1.73e-3 + else: + alpha = 0.138 + beta = -221.0 + kappa = 521.0 + gamma = 1680.0 + f_knee = 1.13e-3 + + A = 1.8e-44 / NC + Sc = 1.0 + np.tanh(gamma * (f_knee - f)) + Sc *= np.exp(-(f**alpha) + beta * f * np.sin(kappa * f)) + Sc *= A * f ** (-7.0 / 3.0) + return Sc + + +def Sn(f, Tobs=0.5, NC=3, R_exists=False, interp_func=None): + """Calculate the sensitivity curve.""" + if R_exists: + R = interpolate.splev(f, interp_func, der=0) + else: + R = 3.0 / 20.0 / (1.0 + 6.0 / 10.0 * (f / fstar) ** 2) * NC + + return Pn(f) / R + SnC(f, Tobs, NC) + + +def response_interpolant(NC): + if os.path.exists(f"{path_to_file}/R.txt"): + data = np.loadtxt(f"{path_to_file}/R.txt") + R = data[:, 1] * NC + f = data[:, 0] * fstar + return True, interpolate.splrep(f, R, s=0) + print("R.txt doesn't exist.") + return False, None + + +def generate_psd(fmin=5.0e-5, fmax=1.0, Tobs_years=0.5, NC=3, npts=500001): + """Return frequency and PSD arrays for the analytic LISA sensitivity curve.""" + R_exists, interp_func = response_interpolant(NC) + f = np.linspace(fmin, fmax, npts) + sens = Sn(f, Tobs_years * YRSID_SI, NC, R_exists, interp_func) + return f, sens + + +def write_lisa_psd(output_dir, fmin=5.0e-5, fmax=1.0, Tobs_years=0.5, NC=3, npts=500001, write_xml=True): + """Write LISA_psd.txt, LISA_psd_plot.png, and optionally A-psd.xml.gz.""" + f, sens = generate_psd(fmin=fmin, fmax=fmax, Tobs_years=Tobs_years, NC=NC, npts=npts) + os.makedirs(output_dir, exist_ok=True) + txt_path = os.path.join(output_dir, "LISA_psd.txt") + png_path = os.path.join(output_dir, "LISA_psd_plot.png") + + np.savetxt(txt_path, np.vstack([f, sens]).T) + + plt.figure() + plt.xlabel("Frequency [Hz]") + plt.ylabel("Characteristic strain") + plt.loglog(f, np.sqrt(f * sens)) + plt.savefig(png_path, bbox_inches="tight") + plt.close() + + if write_xml: + subprocess.run( + [ + "convert_psd_ascii2xml", + "--fname-psd-ascii", + txt_path, + "--conventional-postfix", + "--ifo", + "A", + ], + check=True, + cwd=output_dir, + ) + return txt_path, png_path + + +def parse_args(argv=None): + parser = ArgumentParser() + parser.add_argument("--NC", default=3, type=int, help="Number of channels.") + parser.add_argument("--Tobs", default=0.5, type=float, help="Observation time in years.") + parser.add_argument("--fmin", default=5.0e-5, type=float, help="Lowest PSD frequency.") + parser.add_argument("--fmax", default=1, type=float, help="Highest PSD frequency.") + parser.add_argument("--npts", default=500001, type=int, help="Number of PSD samples.") + parser.add_argument("--output-dir", default=os.getcwd(), help="Directory for generated PSD products.") + parser.add_argument("--skip-xml", action="store_true", help="Do not run convert_psd_ascii2xml.") + return parser.parse_args(argv) + + +def main(argv=None): + opts = parse_args(argv) + print(f"Argument parser has the following arguments:\n{vars(opts)}") + write_lisa_psd( + opts.output_dir, + fmin=opts.fmin, + fmax=opts.fmax, + Tobs_years=opts.Tobs, + NC=opts.NC, + npts=opts.npts, + write_xml=not opts.skip_xml, + ) + + +if __name__ == "__main__": + main() diff --git a/MonteCarloMarginalizeCode/Code/RIFT/LISA/run_checks/__init__.py b/MonteCarloMarginalizeCode/Code/RIFT/LISA/run_checks/__init__.py new file mode 100644 index 000000000..ff0ac4c01 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/LISA/run_checks/__init__.py @@ -0,0 +1 @@ +"""LISA run diagnostics.""" diff --git a/MonteCarloMarginalizeCode/Code/RIFT/LISA/run_checks/plot_RIFT.py b/MonteCarloMarginalizeCode/Code/RIFT/LISA/run_checks/plot_RIFT.py new file mode 100755 index 000000000..e6805c156 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/LISA/run_checks/plot_RIFT.py @@ -0,0 +1,127 @@ +#!/usr/bin/env python + +"""Lightweight diagnostics for LISA-RIFT run products.""" + +from __future__ import print_function + +from argparse import ArgumentParser +import json +import os + +import numpy as np + + +LISA_ILE_COLUMNS = ( + "event", + "m1", + "m2", + "s1x", + "s1y", + "s1z", + "s2x", + "s2y", + "s2z", + "ecliptic_longitude", + "ecliptic_latitude", + "lnL", + "sigma_lnL", + "n_total", + "n_eff", +) + + +def load_lisa_ile_table(path): + data = np.loadtxt(path, ndmin=2) + if data.shape[1] != len(LISA_ILE_COLUMNS): + raise ValueError( + "Expected {} LISA ILE columns in {}, found {}".format( + len(LISA_ILE_COLUMNS), path, data.shape[1] + ) + ) + return data + + +def component_masses_to_mc_eta(m1, m2): + total = m1 + m2 + eta = (m1 * m2) / total**2 + mc = (m1 * m2) ** (3.0 / 5.0) / total ** (1.0 / 5.0) + return mc, eta + + +def summarize_lisa_ile(path, lnL_window=15.0, error_threshold=0.4): + data = load_lisa_ile_table(path) + lnL = data[:, LISA_ILE_COLUMNS.index("lnL")] + sigma = data[:, LISA_ILE_COLUMNS.index("sigma_lnL")] + finite = np.isfinite(lnL) + if not np.any(finite): + raise ValueError("No finite lnL values in {}".format(path)) + + max_index = int(np.nanargmax(lnL)) + max_lnL = float(lnL[max_index]) + high = finite & (lnL >= max_lnL - lnL_window) + high_low_error = high & (sigma <= error_threshold) + m1 = data[:, LISA_ILE_COLUMNS.index("m1")] + m2 = data[:, LISA_ILE_COLUMNS.index("m2")] + mc, eta = component_masses_to_mc_eta(m1, m2) + + return { + "path": os.path.abspath(path), + "n_rows": int(data.shape[0]), + "max_lnL": max_lnL, + "max_index": max_index, + "high_lnL_points": int(np.count_nonzero(high)), + "high_lnL_low_error_points": int(np.count_nonzero(high_low_error)), + "best": { + "m1": float(m1[max_index]), + "m2": float(m2[max_index]), + "mc": float(mc[max_index]), + "eta": float(eta[max_index]), + "ecliptic_longitude": float(data[max_index, LISA_ILE_COLUMNS.index("ecliptic_longitude")]), + "ecliptic_latitude": float(data[max_index, LISA_ILE_COLUMNS.index("ecliptic_latitude")]), + "sigma_lnL": float(sigma[max_index]), + "n_eff": float(data[max_index, LISA_ILE_COLUMNS.index("n_eff")]), + }, + } + + +def build_parser(): + parser = ArgumentParser(description="Summarize LISA-RIFT ILE/all.net-style output.") + parser.add_argument("path", help="LISA ILE output, for example lisa_ile_0_.dat or all.net.") + parser.add_argument("--lnL-window", type=float, default=15.0) + parser.add_argument("--error-threshold", type=float, default=0.4) + parser.add_argument("--json", action="store_true", help="Emit machine-readable JSON.") + return parser + + +def main(argv=None): + opts = build_parser().parse_args(argv) + summary = summarize_lisa_ile( + opts.path, + lnL_window=opts.lnL_window, + error_threshold=opts.error_threshold, + ) + if opts.json: + print(json.dumps(summary, indent=2, sort_keys=True)) + return + + print("LISA RIFT diagnostic summary") + print(" file: {}".format(summary["path"])) + print(" rows: {}".format(summary["n_rows"])) + print(" max lnL: {:.6g}".format(summary["max_lnL"])) + print( + " high-lnL points within window: {} ({} with sigma <= {})".format( + summary["high_lnL_points"], + summary["high_lnL_low_error_points"], + opts.error_threshold, + ) + ) + best = summary["best"] + print( + " best: m1={m1:.6g} m2={m2:.6g} mc={mc:.6g} eta={eta:.6g} " + "lambda={ecliptic_longitude:.6g} beta={ecliptic_latitude:.6g} " + "sigma={sigma_lnL:.6g} neff={n_eff:.6g}".format(**best) + ) + + +if __name__ == "__main__": + main() diff --git a/MonteCarloMarginalizeCode/Code/RIFT/LISA/sangria_test/__init__.py b/MonteCarloMarginalizeCode/Code/RIFT/LISA/sangria_test/__init__.py new file mode 100644 index 000000000..fa004aa7c --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/LISA/sangria_test/__init__.py @@ -0,0 +1 @@ +"""Sangria conversion helpers for LISA/RIFT tests and data products.""" diff --git a/MonteCarloMarginalizeCode/Code/RIFT/LISA/sangria_test/pycbc_to_rift.py b/MonteCarloMarginalizeCode/Code/RIFT/LISA/sangria_test/pycbc_to_rift.py new file mode 100644 index 000000000..874f61513 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/LISA/sangria_test/pycbc_to_rift.py @@ -0,0 +1,161 @@ +"""Convert Sangria A/E/T time series into RIFT LISA HDF5 frames.""" + +from argparse import ArgumentParser +import os + +import h5py +import lal +import matplotlib.pyplot as plt +import numpy as np +from scipy.interpolate import interp1d + +import RIFT.lalsimutils as lsu + +__author__ = "A. Jan" + +DEFAULT_DELTA_T = 8 +DEFAULT_DURATION = 31536000 +DEFAULT_LENGTH = 4194304 + + +def _import_gwpy_timeseries(): + try: + from gwpy.timeseries import TimeSeries + except ImportError as exc: + raise ImportError("gwpy is required to read Sangria frame files") from exc + return TimeSeries + + +def _read_gwpy_frame(frame_path, channel): + TimeSeries = _import_gwpy_timeseries() + return TimeSeries.read(frame_path, channel=channel) + + +def _as_lal_time_series_input(time_series): + """Return data, delta_t, and epoch from either gwpy or PyCBC time series.""" + if hasattr(time_series, "value"): + data = np.asarray(time_series.value) + delta_t = float(time_series.dt.value) + epoch = float(time_series.t0.value) + return data, delta_t, epoch + data = np.asarray(time_series.data) + return data, float(time_series.delta_t), float(time_series._epoch) + + +def create_lal_COMPLEX16TimeSeries( + time_series, + delta_t=DEFAULT_DELTA_T, + duration=DEFAULT_DURATION, + target_length=DEFAULT_LENGTH, +): + """Resample a gwpy/PyCBC time series onto the LISA/RIFT cadence.""" + data, input_delta_t, epoch = _as_lal_time_series_input(time_series) + old_tvals = np.arange(0, input_delta_t * len(data), input_delta_t) + new_tvals = np.arange(0, duration, delta_t) + func = interp1d(old_tvals, data, fill_value=tuple([0, 0]), bounds_error=False) + new_data = func(new_tvals) + + ht_lal = lal.CreateCOMPLEX16TimeSeries( + "ht_lal", + epoch, + 0, + delta_t, + lal.DimensionlessUnit, + len(new_data), + ) + ht_lal.data.data = new_data + 0j + ht_lal = lal.ResizeCOMPLEX16TimeSeries(ht_lal, 0, target_length) + print( + f" Delta T = {ht_lal.deltaT} s, size = {ht_lal.data.length}, " + f"time = {ht_lal.data.length * ht_lal.deltaT / 3600 / 24:2f} days" + ) + return ht_lal + + +def _plot_time_series(tvals, series, channel, save_path): + plt.plot(tvals, series.data.data) + plt.xlabel("Time [s]") + plt.savefig(f"{save_path}/{channel}_time.png") + plt.cla() + + +def _plot_frequency_series(fvals, series, channel, save_path): + plt.loglog(fvals, 2 * fvals * np.abs(series.data.data)) + plt.xlabel("Frequency [Hz]") + plt.ylabel("Characteristic Strain") + plt.savefig(f"{save_path}/{channel}_frequency.png") + plt.cla() + + +def _write_h5_frame(channel, series, save_path): + frame_path = f"{save_path}/{channel}-fake_strain-1000000-10000.h5" + with h5py.File(frame_path, "w") as h5_file: + h5_file.create_dataset("data", data=series.data.data) + h5_file.attrs["deltaF"] = series.deltaF + h5_file.attrs["epoch"] = float(series.epoch) + h5_file.attrs["length"] = series.data.length + h5_file.attrs["f0"] = series.f0 + return frame_path + + +def create_injection_from_time_series(time_series, save_path): + """Write A/E/T diagnostic plots and RIFT HDF5 frames from time series objects.""" + os.makedirs(save_path, exist_ok=True) + time_domain = { + channel: create_lal_COMPLEX16TimeSeries(time_series[channel]) + for channel in ["A", "E", "T"] + } + + tvals = np.arange(0, time_domain["A"].data.length * time_domain["A"].deltaT, time_domain["A"].deltaT) + for channel, series in time_domain.items(): + _plot_time_series(tvals, series, channel, save_path) + + data_dict = { + channel: lsu.DataFourier(series) + for channel, series in time_domain.items() + } + fvals = -data_dict["A"].deltaF * np.arange( + data_dict["A"].data.length // 2, + -data_dict["A"].data.length // 2, + -1, + ) + + frame_paths = {} + for channel, series in data_dict.items(): + _plot_frequency_series(fvals, series, channel, save_path) + frame_paths[channel] = _write_h5_frame(channel, series, save_path) + return frame_paths + + +def create_injection_from_pycbc(pycbc_tseries, save_path): + """Backward-compatible alias for older scripts using PyCBC time series.""" + return create_injection_from_time_series(pycbc_tseries, save_path) + + +def read_gwpy_channels(frame_path, channels=("A", "E", "T")): + """Read named channels from a frame file using gwpy.""" + return {channel: _read_gwpy_frame(frame_path, channel) for channel in channels} + + +def read_pycbc_channels(frame_path, channels=("A", "E", "T")): + """Read named channels with gwpy; kept as a compatibility alias.""" + return read_gwpy_channels(frame_path, channels=channels) + + +def parse_args(argv=None): + parser = ArgumentParser() + parser.add_argument("frame_path", help="Input Sangria frame file.") + parser.add_argument("--save-path", default=os.getcwd(), help="Directory for generated RIFT products.") + parser.add_argument("--channels", default="A,E,T", help="Comma-separated channels to read from the frame.") + return parser.parse_args(argv) + + +def main(argv=None): + opts = parse_args(argv) + channels = tuple(channel.strip() for channel in opts.channels.split(",") if channel.strip()) + time_series = read_gwpy_channels(opts.frame_path, channels=channels) + create_injection_from_time_series(time_series, opts.save_path) + + +if __name__ == "__main__": + main() diff --git a/MonteCarloMarginalizeCode/Code/RIFT/LISA/template_ini/BBH_lisa.ini b/MonteCarloMarginalizeCode/Code/RIFT/LISA/template_ini/BBH_lisa.ini new file mode 100644 index 000000000..3c32ffe87 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/LISA/template_ini/BBH_lisa.ini @@ -0,0 +1,103 @@ +########################################################################################### +[data] +channels = {'A': 'A:FAKE-STRAIN', 'E': 'E:FAKE-STRAIN', 'T': 'T:FAKE-STRAIN'} + +[condor] +accounting_group=ligo.sim.o4.cbc.pe.rift +accounting_group_user=aasim.jan + +[analysis] +ifos=['A','E','T'] +singularity=False +osg=False + +[lalinference] +flow = {'A': 0.0001,'E': 0.0001,'T': 0.0001} +fhigh = {'A': 0.1, 'E': 0.1, 'T': 0.1} + +########################################################################################### +# MAIN ARGUMENTS +########################################################################################### +[rift-pseudo-pipe] +### LISA arguments ### +LISA=True +h5-frame-FD=True +lisa-reference-time=20900000.0 +lisa-reference-frequency=0.004078197479248047 +force-cip=True +data-integration-window-half=300 + +### Approximant arguments ### +approx="NRHybSur3dq8" +l-max=5 +modes="[(2,2),(2,1),(3,3),(3,2),(3,1),(4,4),(4,3),(4,2),(5,5)]" +fmin-template=0.00008 + +### ILE arguments ### +ile-sampler-method="AV" +ile-n-eff=40 +ile-copies=1 +ile-jobs-per-worker=18 +ile-runtime-max-minutes=700 +ile-memory=12144 +internal-ile-use-lnL=True +ile-no-gpu=True +ile-retires=3 +lisa-fixed-sky=False +ecliptic-latitude=0.31601 +ecliptic-longitude=3.48296 +use-gwsurrogate=False +internal-loud-signal-mitigation-suite=True +ile-distance-prior="uniform" + +### test, plotting, extrinsic arguements### +add-extrinsic=True +add-extrinsic-time-resampling=True +batch-extrinsic=True +archive-pesummary-label="run-summary" +internal-test-convergence-threshold=0.001 + +### CIP arguments ### +force-mc-range = "[1635000, 1660000]" +force-eta-range = "[0.18696172, 0.18803828]" +force-s1z-range = "[0.55, 0.65]" +force-s2z-range = "[0.202312, 0.297688]" +force-beta-range = "[0.4, 1.2]" +force-lambda-range = "[0.0, 1.0]" +internal-cip-use-lnL=True +n-output-samples=20000 +cip-sigma-cut=0.4 +fit-save-gp=True +cip-fit-method="rf" +cip-sampler-method="AV" +cip-explode-jobs=200 +cip-explode-jobs-last=300 +internal-use-aligned-phase-coordinates=True +internal-correlate-default=True +assume-nonprecessing=True +spin-magnitude-prior="uniform_aligned" +cip-request-disk="15M" +#internal-use-rescaled-transverse-spin-coordinates=True +#(for precession) + +### Algorithm arguments ### +#(forces subdag) +#internal-propose-converge-last-stage=True +#internal-n-iterations-subdag-max=20 +internal_n_evaluations_per_iteration=4000 +internal-force-iterations=5 +puff-iterations=5 +#fake-data-cache="/home/aasim.jan/LISA_PE/test-pipeline/local.cache" +#manual-initial-grid="/home/aasim.jan/LISA_PE/test-pipeline/overlap-grid.xml.gz" + +### misc arguments ### +use_osg=False +use_osg_cip=False +use_osg_file_transfer=False + +########################################################################################### +[engine] +fref=0.00008 +srate = 8 +distance-max = 200000 +distance-min = 1000 diff --git a/MonteCarloMarginalizeCode/Code/RIFT/LISA/utils/__init__.py b/MonteCarloMarginalizeCode/Code/RIFT/LISA/utils/__init__.py new file mode 100644 index 000000000..60611fbca --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/LISA/utils/__init__.py @@ -0,0 +1 @@ +"""Utilities for LISA-specific RIFT workflows.""" diff --git a/MonteCarloMarginalizeCode/Code/RIFT/LISA/utils/utils.py b/MonteCarloMarginalizeCode/Code/RIFT/LISA/utils/utils.py new file mode 100644 index 000000000..375ac40d8 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/LISA/utils/utils.py @@ -0,0 +1,571 @@ +import matplotlib.pyplot as plt +import lal +import numpy as np +import h5py +from scipy.interpolate import interp1d + +import RIFT.lalsimutils as lsu +from RIFT.LISA.response.LISA_response import * + +YRSID_SI = 31558149.763545603 +C_SI = lal.C_SI +ConstOmega = 1.99098659277e-7 +OrbitR =149597870700.0 + +__author__ = "A. Jan" +########################################################################################### +# Functions +########################################################################################### +def create_resampled_lal_COMPLEX16TimeSeries(tvals, data_dict, new_tvals=None): + """A helper function to create lal COMPLEX16TimeSeries. + Args: + tvals (numpy.array) : time values over which the data is defined, + data_dict (dictionary) : dictionary containing data stored in numpy array, + new_tvals (numpy.array): resampled time values. (set to None if resampling is not needed). + Returns: + data_dict : dictionary containing resampled data stored as lal.COMPLEX16TimeSeries objects.""" + data_dict_new = {} + for channel in data_dict.keys(): + print(f"Reading channel {channel}") + if not(new_tvals is None): + # new_tvals passed as arguments, check if interpolation is needed. + equal_length = np.equal(len(tvals), len(new_tvals)) + old_deltaT, new_deltaT = np.diff(tvals)[0], np.diff(new_tvals)[0] + equal_deltaT = np.equal(old_deltaT, new_deltaT) + else: + # new_tvals not passed as argument, set these to True to bypass interpolation. + equal_length=True + equal_deltaT=True + if equal_length and equal_deltaT: + print("Resampling not requested.") + new_deltaT = np.diff(tvals)[0] + new_data = data_dict[channel] + elif not(equal_length) and equal_deltaT: + print("Resampling not requested, but resizing requested.") + tmp = np.zeros(len(new_tvals)) + tmp[:len(tvals)] = data_dict[channel] + new_data = tmp + else: + new_deltaT = np.diff(new_tvals)[0] + print(f"Resampling from {old_deltaT} s to {new_deltaT} s.") + func = interp1d(tvals, data_dict[channel], fill_value=tuple([0,0]), bounds_error=False) + new_data = func(new_tvals) + ht_lal = lal.CreateCOMPLEX16TimeSeries("ht_lal", 0.0, 0, new_deltaT, lal.DimensionlessUnit, len(new_data)) + ht_lal.data.data = new_data + 0j + print(f" Delta T = {ht_lal.deltaT} s, size = {ht_lal.data.length}, time = {ht_lal.data.length*ht_lal.deltaT/3600/24:2f} days") + data_dict_new[channel] = ht_lal + + return data_dict_new + + +def get_ldc_psds(save_path=None, fvals=None, channels = ["A", "E", "T"], model = "SciRDv1"): + """This function generates LISA psds using the lisa data challenge package. + Args: + save_path (string): path where to save the psds as txt files, + fvals (boolean) : frequency values on which you want to evaluate the PSD, if None then it will generate frequency valuee, + channel (list) : list of channels ["A", "E", "T", "X", "XY"], + model (string) : "Proposal", "SciRDv1", "SciRDdeg1", "MRDv1","MRD_MFR","mldc", "newdrs", "LCESAcall", "redbook". + Returns: + psd_dictionary""" + try: + import ldc.lisa.noise as noise + except ImportError as exc: + raise ImportError("ldc is required to generate LDC LISA PSDs") from exc + if fvals is None: + fmin = 0.0 + deltaT = 8 + fNyq = 0.5/deltaT + deltaF = 1/(4194304*deltaT) + fvals = np.arange(fmin, fNyq, deltaF) + print(f"Generating psds using model = {model}, fmin = {fmin}Hz, fmax = {fNyq}Hz, 1/deltaF = {4194304*deltaT}s, deltaT = {deltaT}s.") + noise_model = noise.get_noise_model(model, fvals) + Sn = {} + Sn["fvals"] = fvals + for channel in channels: + print(f"Channel = {channel}") + Sn[channel] = noise_model.psd(fvals, channel) + if save_path: + np.savetxt(f"{save_path}/{channel}_psd.txt", np.vstack([Sn["fvals"], Sn[channel]]).T) + + return Sn + +def generate_data_from_radler(h5_path, output_as_AET = False, new_tvals = None, output_as_FD = False, condition=True, taper_percent=0.0001): + """This function takes in a radler h5 file and outputs a data dictionary. + Args: + h5_path (string) : path to radler h5 file, + output_as_AET (boolean): set to True if you want the data as A, E, T, + new_tvals (numpy array): pass the new tvals to truncate and/or resample the data, + output_as_FD (boolean) : set to True if you want the data in frequency domain, + Returns: + data_dictionary""" + # Load data + data = h5py.File(h5_path) + # Extract X, Y, Z data + XYZ_data = np.array(data["H5LISA"]['PreProcess']['TDIdata']) + # DeltaT + cadence = np.array(data["H5LISA"]["GWSources"]['MBHB-0']['Cadence']) + # Save as dictionary + data_dict = {} + data_dict["X"], data_dict["Y"], data_dict["Z"] = XYZ_data[:,1], XYZ_data[:,2], XYZ_data[:,3] + # tvals for this data + old_tvals = XYZ_data[:,0] + + # Convert into AET if requested + if output_as_AET: + tmp_dict = data_dict + data_dict = {} + data_dict["A"] = 1/np.sqrt(2) * (tmp_dict["Z"] - tmp_dict["X"]) + data_dict["E"] = 1/np.sqrt(6) * (tmp_dict["X"] - 2*tmp_dict["Y"] + tmp_dict["Z"]) + data_dict["T"] = 1/np.sqrt(3) * (tmp_dict["X"] + tmp_dict["Y"] + tmp_dict["Z"]) + + # new_tvals are none by default, so if they are not provided no interpolatin will occur and this function will just create lal tseries. + data_dict = create_resampled_lal_COMPLEX16TimeSeries(old_tvals, data_dict, new_tvals) + + # Condition if requested + if condition: + print("\tTapering requested") + for channel in data_dict: + TDlen = (data_dict[channel].data.length) + ntaper = int(taper_percent*TDlen) + # taper start of the time series + vectaper= 0.5 - 0.5*np.cos(np.pi*np.arange(ntaper)/(1.*ntaper)) + print(f"\t\t Tapering from index 0 ({vectaper[0]}) to {ntaper}.") + data_dict[channel].data.data[:ntaper] *= vectaper # 0 at 0 and slowly peak + # taper end of the time series + index_front = TDlen-ntaper + print(f"\t\t Tapering from index {index_front} ({vectaper[::-1][0]}) to {TDlen-1}.") + data_dict[channel].data.data[index_front:] *= vectaper[::-1] # slowly drop and then 0 at -1 + + # Convert into FD if requested + if output_as_FD: + if new_tvals is None: + power = np.log2(len(old_tvals)) + else: + power = np.log2(len(new_tvals)) + assert power == np.ceil(power), "The data bins need to be power of 2 for lal FFT routines, make sure len(new_tvals) is a power of 2." + tmp_dict = data_dict + data_dict = {} + for channel in tmp_dict: + data_dict[channel] = lsu.DataFourier(tmp_dict[channel]) + + return data_dict + + +def generate_data_from_sangria(h5_path, output_as_AET = False, new_tvals = None, output_as_FD = False, condition=True, resize = False, add_noise = False, taper_percent=0.0001): + """This function takes in a sangria h5 file and outputs a data dictionary. + Args: + h5_path (string) : path to radler h5 file, + output_as_AET (boolean): set to True if you want the data as A, E, T, + new_tvals (numpy array): pass the new tvals to truncate and/or resample the data, + output_as_FD (boolean) : set to True if you want the data in frequency domain, + Returns: + data_dictionary""" + # Load data + data = h5py.File(h5_path) + # Extract X, Y, Z data for mbhb + XYZ_data_mbhb = data['sky']['mbhb']['tdi'] + # Save as dictionary + data_dict_mbhb = {} + data_dict_mbhb.update({'X':np.array(XYZ_data_mbhb['X']).squeeze(1), + 'Y':np.array(XYZ_data_mbhb['Y']).squeeze(1), + 'Z':np.array(XYZ_data_mbhb['Z']).squeeze(1)}) + old_tvals = np.array(XYZ_data_mbhb['t']).squeeze(1) + + # data_dict is our return dictionary, save as no-noise right now, add noise based on user input. + data_dict = data_dict_mbhb + + # noise goes here + # full data with noise and gbs + if not(add_noise is False): + print(f"Adding noise {add_noise}") + full_data_dict = {} + XYZ_data_full = data['obs']['tdi'] + full_data_dict.update({"X":np.array(XYZ_data_full['X']).squeeze(1), + "Y":np.array(XYZ_data_full['Y']).squeeze(1), + "Z":np.array(XYZ_data_full['Z']).squeeze(1)}) + if add_noise=='with_gbs': + data_dict = full_data_dict + elif add_noise=='without_gbs': + # subtract mbhb signals + noise_with_gb = {} + for channel in ["X", "Y", "Z"]: + noise_with_gb[channel] = full_data_dict[channel] - data_dict_mbhb[channel] + # collect gbs + noise_just_gb = {} + noise_just_gb.update({"X":0, "Y":0, "Z":0}) + for i in ['v','d','i']: + print(i) + XYZ_data_gb = data['sky'][f'{i}gb']["tdi"] + for channel in ["X", "Y", "Z"]: + noise_just_gb[channel] += np.array(XYZ_data_gb[channel]).squeeze(1) + # just get noise + noise_without_gb = {} + for channel in ["X", "Y", "Z"]: + noise_without_gb[channel] = noise_with_gb[channel] - noise_just_gb[channel] + data_dict[channel] = data_dict_mbhb[channel] + noise_without_gb[channel] + + # Convert into AET if requested + if output_as_AET: + tmp_dict = data_dict + data_dict = {} + data_dict["A"] = 1/np.sqrt(2) * (tmp_dict["Z"] - tmp_dict["X"]) + data_dict["E"] = 1/np.sqrt(6) * (tmp_dict["X"] - 2*tmp_dict["Y"] + tmp_dict["Z"]) + data_dict["T"] = 1/np.sqrt(3) * (tmp_dict["X"] + tmp_dict["Y"] + tmp_dict["Z"]) + + # Condition if requested + if condition: + print("\tTapering requested") + for channel in data_dict: + TDlen = len(data_dict[channel]) + ntaper = int(taper_percent*TDlen) + # taper start of the time series + vectaper= 0.5 - 0.5*np.cos(np.pi*np.arange(ntaper)/(1.*ntaper)) + print(f"\t\t Tapering from index 0 ({vectaper[0]}) to {ntaper}.") + data_dict[channel][:ntaper] *= vectaper # 0 at 0 and slowly peak + # taper end of the time series + index_front = TDlen-ntaper + print(f"\t\t Tapering from index {index_front} ({vectaper[::-1][0]}) to {TDlen-1}.") + data_dict[channel][index_front:] *= vectaper[::-1] # slowly drop and then 0 at -1 + + # Resizing + if resize: + old_power = np.log2(len(old_tvals)) + new_power = np.ceil(old_power) + print(f"Resizing from {2**old_power*5/3600/24} to {2**new_power*5/3600/24}.") + for channel in data_dict: + tmp = np.zeros(int(2**new_power)) + tmp[:len(old_tvals)] = data_dict[channel] + data_dict[channel] = tmp + old_tvals = np.arange(0, len(data_dict[channel]), 1) * 5 + # new_tvals are none by default, so if they are not provided no interpolatin will occur and this function will just create lal tseries. + data_dict = create_resampled_lal_COMPLEX16TimeSeries(old_tvals, data_dict, new_tvals) + + # Convert into FD if requested + if output_as_FD: + if new_tvals is None: + power = np.log2(len(old_tvals)) + else: + power = np.log2(len(new_tvals)) + assert power == np.ceil(power), "The data bins need to be power of 2 for lal FFT routines, make sure len(new_tvals) is a power of 2." + tmp_dict = data_dict + data_dict = {} + for channel in tmp_dict: + data_dict[channel] = lsu.DataFourier(tmp_dict[channel]) + + return data_dict + +def get_ldc_mbhb_params(h5_path, dataset = "radler", sangria_signal=0): + """This function takes in a radler h5 file and outputs the parameter of the MBHB injection as a dictionary. + Args: + h5_path (string): path to radler h5 file, + Returns: + parameter dictionary""" + # Load data + data = h5py.File(h5_path) + if dataset == 'radler': + pGW = data["H5LISA"]["GWSources"]['MBHB-0'] + # Extract params + params = {} + params["m1"] = np.array(pGW.get('Mass1')) + params["m2"] = np.array(pGW.get('Mass2')) + params["chi1"] = np.array(pGW.get('Spin1')*np.cos(pGW.get('PolarAngleOfSpin1'))) + params["chi2"] = np.array(pGW.get('Spin2')*np.cos(pGW.get('PolarAngleOfSpin2'))) + + theL = np.array(pGW.get('InitialPolarAngleL')) + phiL = np.array(pGW.get('InitialAzimuthalAngleL')) + longt = np.array(pGW.get('EclipticLongitude')) + lat = np.array(pGW.get('EclipticLatitude')) + + params["tc"] = np.array(pGW.get('CoalescenceTime')) + params["phi0"] = np.array(pGW.get('PhaseAtCoalescence')) + params["DL"] = np.array(pGW.get('Distance')) + + dist = params["DL"] * 1.e6 * lal.PC_SI + # print ("DL = ", DL*1.e-3, "Gpc") + params["beta"] = np.array(lat) + params["lambda"] = np.array(longt) + params["incl"] = np.array(np.arccos( np.cos(theL)*np.sin(lat) + np.cos(lat)*np.sin(theL)*np.cos(longt-phiL))) + + up_psi = np.array(np.sin(params["beta"])*np.sin(theL)*np.cos(params["lambda"] - phiL) - np.cos(theL)*np.cos(params["beta"])) + down_psi = np.array(np.sin(theL)*np.sin(params["lambda"] - phiL)) + params["psi"] = np.array(np.arctan2(up_psi, down_psi)) + params["z"] = np.array(pGW.get("Redshift")) + if dataset == 'sangria': + pGW = data['sky']['mbhb']['cat'][sangria_signal] + # Extract params + params = {} + params["m1"] = np.array(pGW['Mass1']) + params["m2"] = np.array(pGW['Mass2']) + params["chi1"] = np.array(pGW['Spin1']*np.cos(pGW['PolarAngleOfSpin1'])) + params["chi2"] = np.array(pGW['Spin2']*np.cos(pGW['PolarAngleOfSpin2'])) + + theL = np.array(pGW['InitialPolarAngleL']) + phiL = np.array(pGW['InitialAzimuthalAngleL']) + longt = np.array(pGW['EclipticLongitude']) + lat = np.array(pGW['EclipticLatitude']) + + params["tc"] = np.array(pGW['CoalescenceTime']) + params["phi0"] = np.array(pGW['PhaseAtCoalescence']) + params["DL"] = np.array(pGW['Distance']) + + dist = params["DL"] * 1.e6 * lal.PC_SI + # print ("DL = ", DL*1.e-3, "Gpc") + params["beta"] = np.array(lat) + params["lambda"] = np.array(longt) + params["incl"] = np.array(np.arccos( np.cos(theL)*np.sin(lat) + np.cos(lat)*np.sin(theL)*np.cos(longt-phiL))) + + up_psi = np.array(np.sin(params["beta"])*np.sin(theL)*np.cos(params["lambda"] - phiL) - np.cos(theL)*np.cos(params["beta"])) + down_psi = np.array(np.sin(theL)*np.sin(params["lambda"] - phiL)) + params["psi"] = np.array(np.arctan2(up_psi, down_psi)) + params["z"] = np.array(pGW["Redshift"]) + + + + return params + +def modpi(phase): + """Modulus with pi as the period + + This function was originally in BBHx. + + Args: + phase (scalar or np.ndarray): Phase angle. + + Returns: + scalar or np.ndarray: Phase angle modulus by pi. + + """ + # from sylvain + return phase - np.floor(phase / np.pi) * np.pi + +def tSSBfromLframe(tL, lambdaSSB, betaSSB, t0=0.0): + """Get time in SSB frame from time in LISA-frame. + + Compute Solar System Barycenter time ``tSSB`` from retarded time at the center + of the LISA constellation ``tL``. **NOTE**: depends on the sky position + given in solar system barycenter (SSB) frame. + + This function was originally in BBHx. + + Args: + tL (scalar or np.ndarray): Time in LISA constellation reference frame. + lambdaSSB (scalar or np.ndarray): Ecliptic longitude in + SSB reference frame. + betaSSB (scalar or np.ndarray): Ecliptic latitude in SSB reference frame. + t0 (double, optional): Initial start time point away from zero. + (Default: ``0.0``) + + Returns: + scalar or np.ndarray: Time in the SSB frame. + + """ + ConstPhi0 = ConstOmega * t0 + phase = ConstOmega * tL + ConstPhi0 - lambdaSSB + RoC = OrbitR / C_SI + return ( + tL + + RoC * np.cos(betaSSB) * np.cos(phase) + - 1.0 / 2 * ConstOmega * pow(RoC * np.cos(betaSSB), 2) * np.sin(2.0 * phase) + ) + + +# Compute retarded time at the center of the LISA constellation tL from Solar System Barycenter time tSSB */ +def tLfromSSBframe(tSSB, lambdaSSB, betaSSB, t0=0.0): + """Get time in LISA frame from time in SSB-frame. + + Compute retarded time at the center of the LISA constellation frame ``tL`` from + the time in the SSB frame ``tSSB``. **NOTE**: depends on the sky position + given in solar system barycenter (SSB) frame. + + This function was originally in BBHx. + + Args: + tSSB (scalar or np.ndarray): Time in LISA constellation reference frame. + lambdaSSB (scalar or np.ndarray): Ecliptic longitude in + SSB reference frame. + betaSSB (scalar or np.ndarray): Time in LISA constellation reference frame. + t0 (double, optional): Initial start time point away from zero. + (Default: ``0.0``) + + Returns: + scalar or np.ndarray: Time in the LISA frame. + + """ + ConstPhi0 = ConstOmega * t0 + phase = ConstOmega * tSSB + ConstPhi0 - lambdaSSB + RoC = OrbitR / C_SI + return tSSB - RoC * np.cos(betaSSB) * np.cos(phase) + +def LISA_to_SSB(tL, lambdaL, betaL, psiL, t0=0.0): + """Convert sky/orientation from LISA frame to SSB frame. + + Convert the sky and orientation parameters from the center of the LISA + constellation reference to the SSB reference frame. + + The parameters that are converted are the reference time, ecliptic latitude, + ecliptic longitude, and polarization angle. + + This function was originally in BBHx. + + Args: + tL (scalar or np.ndarray): Time in LISA constellation reference frame. + lambdaL (scalar or np.ndarray): Ecliptic longitude in + LISA reference frame. + betaL (scalar or np.ndarray): Ecliptic latitude in LISA reference frame. + psiL (scalar or np.ndarray): Polarization angle in LISA reference frame. + t0 (double, optional): Initial start time point away from zero. + (Default: ``0.0``) + + Returns: + Tuple: (``tSSB``, ``lambdaSSB``, ``betaSSB``, ``psiSSB``) + + + """ + + t0 = t0 * YRSID_SI + + ConstPhi0 = ConstOmega * t0 + coszeta = np.cos(np.pi / 3.0) + sinzeta = np.sin(np.pi / 3.0) + coslambdaL = np.cos(lambdaL) + sinlambdaL = np.sin(lambdaL) + cosbetaL = np.cos(betaL) + sinbetaL = np.sin(betaL) + + lambdaSSB_approx = 0.0 + betaSSB_approx = 0.0 + # Initially, approximate alpha using tL instead of tSSB - then iterate */ + tSSB_approx = tL + for k in range(3): + alpha = ConstOmega * tSSB_approx + ConstPhi0 + cosalpha = np.cos(alpha) + sinalpha = np.sin(alpha) + lambdaSSB_approx = np.arctan2( + cosalpha * cosalpha * cosbetaL * sinlambdaL + - sinalpha * sinbetaL * sinzeta + + cosbetaL * coszeta * sinalpha * sinalpha * sinlambdaL + - cosalpha * cosbetaL * coslambdaL * sinalpha + + cosalpha * cosbetaL * coszeta * coslambdaL * sinalpha, + cosbetaL * coslambdaL * sinalpha * sinalpha + - cosalpha * sinbetaL * sinzeta + + cosalpha * cosalpha * cosbetaL * coszeta * coslambdaL + - cosalpha * cosbetaL * sinalpha * sinlambdaL + + cosalpha * cosbetaL * coszeta * sinalpha * sinlambdaL, + ) + betaSSB_approx = np.arcsin( + coszeta * sinbetaL + + cosalpha * cosbetaL * coslambdaL * sinzeta + + cosbetaL * sinalpha * sinzeta * sinlambdaL + ) + tSSB_approx = tSSBfromLframe(tL, lambdaSSB_approx, betaSSB_approx, t0) + + lambdaSSB_approx = lambdaSSB_approx % (2 * np.pi) + # /* Polarization */ + psiSSB = modpi( + psiL + + np.arctan2( + cosalpha * sinzeta * sinlambdaL - coslambdaL * sinalpha * sinzeta, + cosbetaL * coszeta + - cosalpha * coslambdaL * sinbetaL * sinzeta + - sinalpha * sinbetaL * sinzeta * sinlambdaL, + ) + ) + + return np.vstack([tSSB_approx, lambdaSSB_approx, betaSSB_approx, psiSSB]).T + +def SSB_to_LISA(tSSB, lambdaSSB, betaSSB, psiSSB, t0=0.0): + """Convert sky/orientation from SSB frame to LISA frame. + + Convert the sky and orientation parameters from the SSB reference frame to the center of the LISA + constellation reference frame. + + The parameters that are converted are the reference time, ecliptic latitude, + ecliptic longitude, and polarization angle. + + This function was originally in BBHx. + + Args: + tSSB (scalar or np.ndarray): Time in SSB reference frame. + lambdaSSB (scalar or np.ndarray): Ecliptic longitude in + SSB reference frame. + betaSSB (scalar or np.ndarray): Ecliptic latitude in SSB reference frame. + psiSSB (scalar or np.ndarray): Polarization angle in SSB reference frame. + t0 (double, optional): Initial start time point away from zero in years. + (Default: ``0.0``) + + Returns: + Tuple: (``tL``, ``lambdaL``, ``betaL``, ``psiL``) + + """ + t0 = t0 * YRSID_SI + + ConstPhi0 = ConstOmega * t0 + alpha = 0.0 + cosalpha = 0 + sinalpha = 0.0 + coslambda = 0 + sinlambda = 0.0 + cosbeta = 0.0 + sinbeta = 0.0 + + coszeta = np.cos(np.pi / 3.0) + sinzeta = np.sin(np.pi / 3.0) + coslambda = np.cos(lambdaSSB) + sinlambda = np.sin(lambdaSSB) + cosbeta = np.cos(betaSSB) + sinbeta = np.sin(betaSSB) + + alpha = ConstOmega * tSSB + ConstPhi0 + cosalpha = np.cos(alpha) + sinalpha = np.sin(alpha) + tL = tLfromSSBframe(tSSB, lambdaSSB, betaSSB, t0) + lambdaL = np.arctan2( + cosalpha * cosalpha * cosbeta * sinlambda + + sinalpha * sinbeta * sinzeta + + cosbeta * coszeta * sinalpha * sinalpha * sinlambda + - cosalpha * cosbeta * coslambda * sinalpha + + cosalpha * cosbeta * coszeta * coslambda * sinalpha, + cosalpha * sinbeta * sinzeta + + cosbeta * coslambda * sinalpha * sinalpha + + cosalpha * cosalpha * cosbeta * coszeta * coslambda + - cosalpha * cosbeta * sinalpha * sinlambda + + cosalpha * cosbeta * coszeta * sinalpha * sinlambda, + ) + betaL = np.arcsin( + coszeta * sinbeta + - cosalpha * cosbeta * coslambda * sinzeta + - cosbeta * sinalpha * sinzeta * sinlambda + ) + psiL = modpi( + psiSSB + + np.arctan2( + coslambda * sinalpha * sinzeta - cosalpha * sinzeta * sinlambda, + cosbeta * coszeta + + cosalpha * coslambda * sinbeta * sinzeta + + sinalpha * sinbeta * sinzeta * sinlambda, + ) + ) + + return np.vstack([tL, lambdaL, betaL, psiL]).T + + +def get_secondary_mode_for_skylocation(coalesence_time, lamda, beta, psi = 0.0, t0=0.0): + """ + Calculate the secondary mode location for a given sky location in the SSB frame. + + Parameters: + coalescence_time (float): The time of coalescence. + lambda_param (float): The lambda parameter for the SSB frame. + beta (float): The beta parameter, which is bimodal in the LISA frame. + psi (float, optional): The psi parameter (default is 0.0, not necessary). + t0 (float, optional): The initial time (default is 0.0). + + Returns: + numpy.ndarray: The secondary mode parameters in the SSB frame, needed for grid generation. + """ + lisa_params = SSB_to_LISA(coalesence_time, lamda, beta, psi, t0) + lisa_params_new = lisa_params + # in lisa frame only beta is bimodal (for full detector response) and its location is -beta_L + lisa_params_new[0, 2] = -1.0*lisa_params_new[0, 2] + SSB_params_new = LISA_to_SSB(lisa_params_new[0, 0], lisa_params_new[0, 1], lisa_params_new[0, 2], lisa_params_new[0, 3]) + return SSB_params_new diff --git a/MonteCarloMarginalizeCode/Code/RIFT/hyperpipe/drivers/base.py b/MonteCarloMarginalizeCode/Code/RIFT/hyperpipe/drivers/base.py index 6824a2c53..281d1e6e1 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/hyperpipe/drivers/base.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/hyperpipe/drivers/base.py @@ -238,11 +238,16 @@ def run(self, argv: Optional[Sequence[str]] = None) -> str: logger.info("Loaded grid %r with columns %r", opts.using_eos, column_names) start, stop = opts.eos_start_index, opts.eos_end_index - if start < 0 or stop > rows.shape[0]: - raise SystemExit( - f"marg driver: index range [{start},{stop}) exceeds grid " - f"size {rows.shape[0]}." - ) + if start < 0: + raise SystemExit(f"marg driver: negative start index {start}.") + # CLAMP the (fixed-chunk) range to the actual grid size. The DAG sizes + # the per-iteration MARG chunks from n-samples-per-job, but the grid an + # adaptive iteration actually places (e.g. a puffball) can be SMALLER, so + # the tail chunks run off the end. Process what exists; a fully + # out-of-range chunk just yields an empty (header-only) output that the + # consolidation step ignores. (Previously this raised SystemExit.) + n_rows = rows.shape[0] + start, stop = min(start, n_rows), min(stop, n_rows) for i in range(start, stop): row_values = rows[i, 2:] diff --git a/MonteCarloMarginalizeCode/Code/RIFT/simulation_manager/database.py b/MonteCarloMarginalizeCode/Code/RIFT/simulation_manager/database.py index 0e5266e9c..2fc7b79c2 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/simulation_manager/database.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/simulation_manager/database.py @@ -305,12 +305,31 @@ def needs_more_work(self) -> bool: return self.data.get("current_level", 0) < self.data.get("target_level", 0) def write(self, sim_dir: Union[str, Path]) -> None: - Path(sim_dir, self.FILENAME).write_text( - json.dumps(self.data, indent=2, sort_keys=True) + "\n") + # Atomic write: a plain write_text() truncates the file to 0 bytes + # before refilling it, so a concurrent reader (parallel marg jobs all + # share one archive) can observe an EMPTY status.json -> JSONDecodeError. + # Mirror IndexAppend._write_all: write a per-PID/thread temp then + # os.replace() it in (atomic on POSIX). Reader never sees a partial file. + path = Path(sim_dir, self.FILENAME) + tmp = path.with_name( + "{}.{}.{}.tmp".format(path.name, os.getpid(), threading.get_ident())) + tmp.write_text(json.dumps(self.data, indent=2, sort_keys=True) + "\n") + os.replace(tmp, path) # atomic on POSIX @classmethod def read(cls, sim_dir: Union[str, Path]) -> "StatusRecord": - return cls(json.loads(Path(sim_dir, cls.FILENAME).read_text())) + # Defensive: even with atomic writes, a networked fs can briefly expose + # an empty/partial read between create and rename. Retry a few times on + # a decode error before giving up. + path = Path(sim_dir, cls.FILENAME) + for attempt in range(5): + text = path.read_text() + try: + return cls(json.loads(text)) + except json.JSONDecodeError: + if attempt == 4: + raise + time.sleep(0.05) # --------------------------------------------------------------------------- diff --git a/MonteCarloMarginalizeCode/Code/bin/convergence_test_samples.py b/MonteCarloMarginalizeCode/Code/bin/convergence_test_samples.py index 90935d286..3a2c397ae 100755 --- a/MonteCarloMarginalizeCode/Code/bin/convergence_test_samples.py +++ b/MonteCarloMarginalizeCode/Code/bin/convergence_test_samples.py @@ -151,6 +151,23 @@ def read_samples(fname): if opts.verbose: print(" Samples 2 expanded fields ", samples2.dtype.names) +# Ensure mc/eta exist: hyperpipeline posteriors (RIFT_HYPERPIPELINE_FORMAT) carry +# only m1/m2, and standard_expand_samples does not always add the chirp-mass +# coordinates the convergence test fits -> "no field of name mc". Derive them. +def _ensure_mc_eta(samples): + if 'm1' not in samples.dtype.names or 'm2' not in samples.dtype.names: + return samples + m1 = samples['m1']; m2 = samples['m2'] + if 'mc' not in samples.dtype.names: + samples = add_field(samples, [('mc', float)]) + samples['mc'] = (m1 * m2) ** 0.6 / (m1 + m2) ** 0.2 + if 'eta' not in samples.dtype.names: + samples = add_field(samples, [('eta', float)]) + samples['eta'] = (m1 * m2) / (m1 + m2) ** 2 + return samples +samples1 = _ensure_mc_eta(samples1) +samples2 = _ensure_mc_eta(samples2) + # sanity check: is this a result for PE? if 'm1' in samples1.dtype.names: # Add missing fields needed for some tests diff --git a/MonteCarloMarginalizeCode/Code/bin/convert_primary_sky_mode_to_secondary b/MonteCarloMarginalizeCode/Code/bin/convert_primary_sky_mode_to_secondary new file mode 100755 index 000000000..9bd6fd7da --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/bin/convert_primary_sky_mode_to_secondary @@ -0,0 +1,54 @@ +#! /usr/bin/env python +"""Convert the PRIMARY LISA sky mode of a grid to the SECONDARY (reflected) mode. + +LISA's response makes the ecliptic latitude bimodal: a sky location and its +reflection fit the data nearly equally well. Once the primary mode is explored, +one iteration on the reflected mode suffices to map it. This tool rewrites each +grid point's ecliptic_longitude/latitude to the reflected location (the sky-mode +transform lives in RIFT.LISA.utils.utils.get_secondary_mode_for_skylocation). + +Hyperpipeline ASCII grids (named columns; ecliptic_longitude/ecliptic_latitude) +are handled by name -- no positional all.net hacking. Legacy XML grids are also +supported for back-compatibility. +""" +from argparse import ArgumentParser +import numpy as np + +import RIFT.lalsimutils as lsu +from RIFT.LISA.utils.utils import get_secondary_mode_for_skylocation +from RIFT.misc import hyperpipeline_io as hpio + +parser = ArgumentParser() +parser.add_argument("--lisa-reference-time", default=0.0, help="LISA coalescence time") +parser.add_argument("--fname", help="grid to reflect (hyperpipeline .dat or XML)") +parser.add_argument("--fname-out", default=None, help="output (default: overwrite --fname)") +opts = parser.parse_args() + +t_c = float(opts.lisa_reference_time) +out = opts.fname_out or opts.fname + + +def _reflect(lamda, beta): + sec = np.asarray(get_secondary_mode_for_skylocation(t_c, float(lamda), float(beta), 0.0)).reshape(-1) + # get_secondary_mode_for_skylocation returns [t, lambda, beta, psi] in the SSB frame + return sec[1], sec[2] + + +if hpio.sniff(opts.fname): + arr, cols = hpio.read_table(opts.fname) + arr = np.atleast_1d(arr) + if "ecliptic_longitude" not in cols: + raise SystemExit("convert_primary_sky_mode_to_secondary: grid has no sky columns") + mat = np.empty((len(arr), len(cols))) + for i in range(len(arr)): + lam_sec, bet_sec = _reflect(arr["ecliptic_longitude"][i], arr["ecliptic_latitude"][i]) + for j, c in enumerate(cols): + mat[i, j] = arr[c][i] + mat[i, cols.index("ecliptic_longitude")] = lam_sec + mat[i, cols.index("ecliptic_latitude")] = bet_sec + hpio.write_table(out, cols, mat) +else: + Plist = lsu.xml_to_ChooseWaveformParams_array(opts.fname) + for P in Plist: + P.phi, P.theta = _reflect(P.phi, P.theta) + lsu.ChooseWaveformParams_array_to_xml(Plist, out) diff --git a/MonteCarloMarginalizeCode/Code/bin/create_event_parameter_pipeline_BasicIteration b/MonteCarloMarginalizeCode/Code/bin/create_event_parameter_pipeline_BasicIteration index 0a6128457..2f6346bae 100755 --- a/MonteCarloMarginalizeCode/Code/bin/create_event_parameter_pipeline_BasicIteration +++ b/MonteCarloMarginalizeCode/Code/bin/create_event_parameter_pipeline_BasicIteration @@ -245,6 +245,10 @@ parser.add_argument("--puff-exe",default=None,help="util_ParameterPuffball.py") parser.add_argument("--puff-args",default=None,help="util_ParameterPuffball arguments. If not specified, puffball will not be performed ") parser.add_argument("--puff-cadence",default=None,type=int,help="Every n iterations (not including 0), the puffball code will be applied. Puffball points will be done *in addition* to the usual results from the DAG. (The puffball is based on perturbing points from that iteration, and this will roughly double that iteration in ILE job size). Proposed value 2 (i.e., puff overlap-grid-2, ...-4, ...-6. If not specified, puffball will not be performed ") parser.add_argument("--puff-max-it",default=-1,type=int,help="Maximum iteration number that puffball is applied. If negative, puffball is not applied ") +parser.add_argument("--search-reflected-sky-mode",action='store_true',help="LISA: at one iteration, reflect the CIP-proposed grid's sky to the secondary (reflected) sky mode (the LISA latitude is bimodal), so that iteration explores the reflected mode.") +parser.add_argument("--search-reflected-sky-mode-iteration",default=None,type=int,help="Iteration at which to reflect the sky grid. Default: n_iterations-2 (the primary mode is mapped by then).") +parser.add_argument("--reflected-sky-mode-exe",default=None,help="convert_primary_sky_mode_to_secondary executable.") +parser.add_argument("--lisa-reference-time",default=0.0,type=float,help="LISA coalescence time (for the reflected-sky transform).") parser.add_argument("--calmarg-pilot",action='store_true',help="Option C adaptive calibration: add per-iteration cal PILOT jobs (harvest top-lnL points from iteration N's composite, run ILE --calibration-dump-responsibilities, fit+consolidate a cal proposal that SEEDS wide_{N+1} via --calibration-proposal-breadcrumb). Requires the wide ILE args to already carry the calibration envelope + (per-iteration) proposal-breadcrumb path.") parser.add_argument("--calmarg-pilot-cadence",default=1,type=int,help="Run a cal pilot every n iterations. Default 1 (every iteration until the cap).") parser.add_argument("--calmarg-pilot-max-it",default=3,type=int,help="Stop launching cal pilots after this iteration (cal is boring -> freeze the proposal once learned). Default 3.") @@ -1245,6 +1249,38 @@ if puff_args and puff_cadence: _add_hpip_condor_env(puff_job) puff_job.write_sub_file() +# LISA reflected-sky-mode: a per-iteration node that reflects the CIP-proposed +# grid's sky to the secondary mode (convert_primary_sky_mode_to_secondary), +# overwriting overlap-grid-(it+1) in place so the NEXT iteration's ILE explores +# the reflected mode. arg_str ends with "--fname " so write_convert_sub's +# add_arg(file_input) becomes --fname . +reflected_mode_job = None +reflected_mode_iter = None +if opts.search_reflected_sky_mode and int(opts.n_iterations) > 1: + reflect_exe = opts.reflected_sky_mode_exe or dag_utils.which("convert_primary_sky_mode_to_secondary") + reflected_mode_iter = opts.search_reflected_sky_mode_iteration + if reflected_mode_iter is None: + reflected_mode_iter = int(opts.n_iterations) - 2 + if reflected_mode_iter <= 0: + reflected_mode_iter = int(opts.n_iterations) - 1 + reflected_mode_job, reflected_mode_job_name = dag_utils.write_convert_sub( + tag='reflected_sky_mode', log_dir=None, + arg_str="--lisa-reference-time {} --fname ".format(opts.lisa_reference_time), + file_input=opts.working_directory+'/overlap-grid-$(macroiterationnext).{}'.format(grid_suffix), + file_output=opts.working_directory+'/iteration_$(macroiteration)_cip/logs/reflected-sky-$(cluster)-$(process).out', + out_dir=opts.working_directory, exe=reflect_exe, + universe=local_worker_universe, no_grid=no_worker_grid) + reflected_mode_job.add_condor_cmd("initialdir", opts.working_directory+"/iteration_$(macroiteration)_cip") + reflected_mode_job.set_log_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/reflected-sky-$(cluster)-$(process).log") + reflected_mode_job.set_stderr_file(opts.working_directory+"/iteration_$(macroiteration)_cip/logs/reflected-sky-$(cluster)-$(process).err") + reflected_mode_job.add_condor_cmd('request_disk', opts.general_request_disk) + if opts.use_full_submit_paths: + reflected_mode_job.set_sub_file(opts.working_directory+"/"+reflected_mode_job.get_sub_file()) + if opts.condor_containerize_nonworker and singularity_image: + reflected_mode_job.add_condor_cmd("+SingularityImage", '"' + singularity_image + '"') + _add_hpip_condor_env(reflected_mode_job) + reflected_mode_job.write_sub_file() + ## calibration PILOT job (Option C adaptive cal; see RIFT/calmarg/DESIGN_adaptive_driver.md) calpilot_job = None @@ -2097,6 +2133,22 @@ for it in np.arange(it_start,opts.n_iterations): # separately in last_puff_node, applied to ilePuff jobs via extra_parent_nodes. last_puff_node = puff_node + # LISA reflected-sky-mode: at the chosen iteration, reflect the grid CIP just + # produced (overlap-grid-(it+1)) BEFORE the next iteration's ILE reads it. + # Unlike puffball this MUST gate the next ILE (it rewrites the grid in place), + # so fold it into parent_fit_node. + if (reflected_mode_job is not None) and (it == reflected_mode_iter): + print(" Reflected-sky-mode for iteration ", it) + reflect_node = pipeline.CondorDAGNode(reflected_mode_job) + reflect_node.add_macro("macroiteration", it) + reflect_node.add_macro("macroiterationnext", it+1) + reflect_node.set_category("CONVERT") + reflect_node.set_retry(opts.general_retries) + if not (parent_fit_node is None): + reflect_node.add_parent(parent_fit_node) + dag.add_node(reflect_node) + parent_fit_node = reflect_node + # Calibration PILOT for this iteration (Option C). Runs IN PARALLEL with CIP/puff -- # it does NOT gate them (parent_fit_node is left untouched) -- harvesting iteration # it's composite (via unify_node, which guarantees a non-empty composite) and emitting diff --git a/MonteCarloMarginalizeCode/Code/bin/helper_LISA_Events.py b/MonteCarloMarginalizeCode/Code/bin/helper_LISA_Events.py index 44afeda45..2802813fe 100755 --- a/MonteCarloMarginalizeCode/Code/bin/helper_LISA_Events.py +++ b/MonteCarloMarginalizeCode/Code/bin/helper_LISA_Events.py @@ -69,13 +69,28 @@ def _write_initial_grid(path, opts): ]) rows = [] - for idx in range(opts.grid_size): - row = base.copy() - if opts.grid_size > 1: - offset = (idx - (opts.grid_size - 1) / 2.0) * opts.grid_fractional_width - row[2] = opts.mass1 * (1.0 + offset) - row[3] = opts.mass2 * (1.0 - offset) - rows.append(row) + if opts.grid_size <= 1: + rows = [base.copy()] + else: + # Build a 2-D lattice that varies m1 and m2 INDEPENDENTLY, so the grid + # spans chirp mass AND symmetric mass ratio in two dimensions. The old + # 1-D line (m1*(1+off), m2*(1-off)) is collinear in (mc,eta) and is + # degenerate for CIP's 2-D (mc,eta) quadratic/rf fit (-> NaN/no output). + side = max(2, int(round(np.sqrt(opts.grid_size)))) + offs = np.linspace(-1.0, 1.0, side) * opts.grid_fractional_width + # vary-sky: jitter the ecliptic sky INDEPENDENTLY of the mass lattice, so + # the grid spans (mc, eta, phi, theta) in 4-D (a sky tied to the mass + # index would lie on a 2-D manifold -> degenerate for the 4-D CIP fit). + rng = np.random.RandomState(0) + for i, oi in enumerate(offs): + for j, oj in enumerate(offs): + row = base.copy() + row[2] = opts.mass1 * (1.0 + oi) + row[3] = opts.mass2 * (1.0 + oj) + if opts.vary_sky: + row[10] = opts.ecliptic_longitude + rng.uniform(-1.0, 1.0) * opts.sky_grid_width + row[11] = opts.ecliptic_latitude + rng.uniform(-1.0, 1.0) * opts.sky_grid_width + rows.append(row) hyperpipeline_io.write_table(path, columns, np.array(rows)) @@ -88,6 +103,9 @@ def build_parser(): parser.add_argument("--ile-args", default="args_ile.txt") parser.add_argument("--cip-args-list", default="args_cip_list.txt") parser.add_argument("--test-args", default="args_test.txt") + parser.add_argument("--puff-args", default="args_puff.txt") + parser.add_argument("--puff-factor", type=float, default=1.0, + help="util_ParameterPuffball --puff-factor: scale of the inter-iteration grid spread.") parser.add_argument("--transfer-file-list", default="helper_transfer_files.txt") parser.add_argument("--cepp-command-file", default="command-cepp-lisa.sh") @@ -95,13 +113,13 @@ def build_parser(): parser.add_argument( "--channel-name", action="append", - default=["A=fake_strain", "E=fake_strain", "T=fake_strain"], + default=None, help="LISA channel assignment, e.g. A=fake_strain.", ) parser.add_argument( "--psd-file", action="append", - default=["A=A_psd.xml.gz", "E=E_psd.xml.gz", "T=T_psd.xml.gz"], + default=None, help="PSD assignment, e.g. A=A_psd.xml.gz.", ) @@ -115,11 +133,13 @@ def build_parser(): parser.add_argument("--spin2z", type=float, default=0.0) parser.add_argument("--ecliptic-longitude", type=float, default=1.0) parser.add_argument("--ecliptic-latitude", type=float, default=0.3) + parser.add_argument("--vary-sky", action="store_true", help="Treat ecliptic sky location as an intrinsic grid parameter.") parser.add_argument("--grid-size", type=int, default=3) parser.add_argument("--grid-fractional-width", type=float, default=1.0e-3) + parser.add_argument("--sky-grid-width", type=float, default=0.02) # rad; vary-sky grid sky jitter parser.add_argument("--approximant", default="IMRPhenomD") - parser.add_argument("--fmin-template", type=float, default=1.0e-3) + parser.add_argument("--fmin-template", type=float, default=1.0e-4) # in-band for LISA (1e-3 starts near top of band) parser.add_argument("--fmax", type=float, default=0.125) parser.add_argument("--reference-freq", type=float, default=5.0e-3) parser.add_argument("--srate", type=float, default=0.25) @@ -127,20 +147,33 @@ def build_parser(): parser.add_argument("--modes", default="[(2,2)]") parser.add_argument("--lisa-reference-time", type=float, default=0.0) parser.add_argument("--lisa-reference-frequency", type=float, default=5.0e-3) - parser.add_argument("--d-max", type=float, default=5000.0) - parser.add_argument("--d-min", type=float, default=1.0) + parser.add_argument("--data-integration-window-half", type=float, default=300.0) # ~600s window: a 16s window mis-marginalizes the long LISA signal -> biased lnL + parser.add_argument("--d-max", type=float, default=100000.0) # LISA MBHBs reach cosmological distances + parser.add_argument("--d-min", type=float, default=1000.0) parser.add_argument("--event-time", type=float, default=0.0) parser.add_argument("--zero-likelihood", action="store_true") - parser.add_argument("--n-eff", type=int, default=2) - parser.add_argument("--n-max", type=int, default=20) - parser.add_argument("--n-chunk", type=int, default=10) + parser.add_argument("--no-adapt", action="store_true", + help="Disable adaptive extrinsic sampling (uniform). Loud signals need adaptation, so default is OFF.") + parser.add_argument("--ile-sampler-method", default="AV") + parser.add_argument("--n-eff", type=int, default=20) + parser.add_argument("--n-max", type=int, default=8000) + parser.add_argument("--n-chunk", type=int, default=500) parser.add_argument("--save-P", type=float, default=0.1) - parser.add_argument("--cip-fit-method", default="quadratic") + parser.add_argument("--cip-fit-method", default="rf") # paper uses random-forest; robust to non-quadratic lnL surfaces + parser.add_argument("--cip-sampler-method", default="AV") parser.add_argument("--cip-iterations", default="1") parser.add_argument("--cip-n-output-samples", type=int, default=100) - parser.add_argument("--cip-lnL-offset", type=float, default=100.0) + parser.add_argument("--cip-lnL-offset", type=float, default=2000.0) # keep all grid points for the fit (loud-signal lnL spread is large) + parser.add_argument("--cip-n-eff", type=int, default=100) + parser.add_argument("--cip-n-max", type=int, default=3000000) + parser.add_argument("--cip-m-max-cut", default="1e8", + help="CIP --M-max-cut (Msun). LISA MBHBs need a large value.") + parser.add_argument("--cip-sigma-cut", default="10.0", + help="CIP --sigma-cut. Relaxed for single-sample high-SNR demo integrals.") + parser.add_argument("--cip-mass-range-frac", type=float, default=0.0, + help="Explicit half-width (fractional) of the CIP mc/mtot range; if 0, auto = 3x the grid fractional width (brackets the grid).") parser.add_argument("--test-threshold", type=float, default=0.02) parser.add_argument("--cepp-exe", default="create_event_parameter_pipeline_BasicIteration") parser.add_argument("--ile-exe", default="integrate_likelihood_extrinsic_batchmode_lisa") @@ -158,11 +191,16 @@ def main(argv=None): workdir = os.path.abspath(opts.working_directory) os.makedirs(workdir, exist_ok=True) + if opts.channel_name is None: + opts.channel_name = ["A=fake_strain", "E=fake_strain", "T=fake_strain"] + if opts.psd_file is None: + opts.psd_file = ["A=A_psd.xml.gz", "E=E_psd.xml.gz", "T=T_psd.xml.gz"] input_grid = os.path.join(workdir, opts.input_grid) ile_args = os.path.join(workdir, opts.ile_args) cip_args_list = os.path.join(workdir, opts.cip_args_list) test_args = os.path.join(workdir, opts.test_args) + puff_args = os.path.join(workdir, opts.puff_args) transfer_file_list = os.path.join(workdir, opts.transfer_file_list) cepp_command_file = os.path.join(workdir, opts.cepp_command_file) @@ -182,11 +220,10 @@ def main(argv=None): ile_parts = [ "--LISA", "--h5-frame-FD", - "--lisa-fixed-sky", "1", - "--ecliptic-longitude", opts.ecliptic_longitude, - "--ecliptic-latitude", opts.ecliptic_latitude, + "--time-marginalization", "--lisa-reference-time", opts.lisa_reference_time, "--lisa-reference-frequency", opts.lisa_reference_frequency, + "--data-integration-window-half", opts.data_integration_window_half, "--modes", opts.modes, "--cache-file", opts.cache_file, "--event-time", opts.event_time, @@ -206,26 +243,116 @@ def main(argv=None): "--n-max", opts.n_max, "--n-chunk", opts.n_chunk, "--save-P", opts.save_P, - "--no-adapt", + "--sampler-method", opts.ile_sampler_method, "--internal-use-lnL", ] + # Adaptive sampling is REQUIRED for loud LISA signals: with uniform sampling + # (--no-adapt) the sharp extrinsic peak is single-sample-dominated (eff_samp + # collapses to 1). --force-adapt-all adapts every extrinsic dimension + # (distance, angles) so the sampler concentrates on the peak. + if opts.no_adapt: + ile_parts.append("--no-adapt") + else: + ile_parts.append("--force-adapt-all") + # ILE always evaluates at a FIXED sky per grid point. Known-sky: the single + # injected sky (hardcoded). Vary-sky: each grid point carries its own + # ecliptic_longitude/latitude, which --sim-grid feeds into the ILE -- so we + # still pass --lisa-fixed-sky 1 but let the per-row grid sky win (no hardcode). + if not opts.vary_sky: + ile_parts[3:3] = [ + "--lisa-fixed-sky", "1", + "--ecliptic-longitude", opts.ecliptic_longitude, + "--ecliptic-latitude", opts.ecliptic_latitude, + ] + else: + ile_parts[3:3] = ["--lisa-fixed-sky", "1"] if opts.zero_likelihood: ile_parts.append("--zero-likelihood") _write_arg_file(ile_args, ile_parts) + # CIP fits only the parameters that actually VARY across the grid. In + # known-sky mode the ecliptic sky location is fixed (constant columns), so + # fitting it is degenerate AND those coordinates are not understood by CIP's + # waveform-parameter machinery (-> "No attribute ecliptic_longitude"). Only + # add the sky as a fit parameter when it is varied (--vary-sky). + cip_params = ["--parameter", "mc", "--parameter", "eta"] + # LISA sky is fit as phi (ecliptic longitude) / theta (ecliptic latitude): + # CIP reads the ecliptic_longitude/latitude NAMED hyperpipeline columns into + # P.phi/P.theta (hyperpipeline_io alias) and fits them as ordinary + # coordinates -- no positional all.net special-casing. + sky_range_args = [] + if opts.vary_sky: + _skw = max(3.0 * opts.sky_grid_width, 0.06) + cip_params += ["--parameter", "phi", "--parameter", "theta"] + sky_range_args = [ + "--phi-range", "[{},{}]".format(opts.ecliptic_longitude - _skw, + opts.ecliptic_longitude + _skw), + "--theta-range", "[{},{}]".format(opts.ecliptic_latitude - _skw, + opts.ecliptic_latitude + _skw), + ] + # CIP's posterior MC sampler defaults to a STELLAR-mass chirp-mass range + # ([0.9, 250] Msun); for a LISA MBHB (mc ~ 1e4-1e7 Msun) the sampler would + # never place a point near the signal -> eff_samp=nan. Bracket mc and mtot + # around the injected masses (analogue of the paper's force-mc-range). + _mtot = opts.mass1 + opts.mass2 + _mc = (opts.mass1 * opts.mass2) ** 0.6 / _mtot ** 0.2 + # Bracket the GRID (plus margin): a range much wider than the grid samples + # mostly where the fit extrapolates -> eff_samp=nan; one matched to the grid + # keeps the CIP sampler where lnL is actually constrained. The injected + # mc/eta are measured to ~Fisher precision (<< grid), so grid-tied is also + # tight enough to bracket the posterior. + _w = max(opts.cip_mass_range_frac, 1.5 * opts.grid_fractional_width) + _eta = (opts.mass1 * opts.mass2) / _mtot ** 2 + # Bracket eta too: with only mc/mtot bounded, the CIP posterior drifts to the + # eta FLOOR (0.01) -> extreme q -> garbage m1 (mtot=mc/eta^0.6 blows up). The + # truth eta (~0.25 for near-equal MBHBs) is near the ceiling, so a tight + # grid-tied eta window is essential (analogue of the paper's force-eta-range). + cip_range_args = [ + "--mc-range", "[{},{}]".format(_mc * (1.0 - _w), _mc * (1.0 + _w)), + "--mtot-range", "[{},{}]".format(_mtot * (1.0 - _w), _mtot * (1.0 + _w)), + "--eta-range", "[{},{}]".format(_eta * (1.0 - _w), min(0.2499999, _eta * (1.0 + _w))), + "--n-eff", str(opts.cip_n_eff), "--n-max", str(opts.cip_n_max), + ] cip_line = _quote_join([ opts.cip_iterations, "--fit-method", opts.cip_fit_method, - "--parameter", "mc", - "--parameter", "eta", - "--parameter", "ecliptic_longitude", - "--parameter", "ecliptic_latitude", + # AV integrator works in lnL space automatically (no exp() overflow at + # loud-signal lnL) and avoids the default sampler's lsoda CDF inversion + # (which NaNs on the sharp high-SNR posterior). --internal-use-lnL too. + "--sampler-method", opts.cip_sampler_method, + "--internal-use-lnL", + *cip_params, + *cip_range_args, + *sky_range_args, "--n-output-samples", opts.cip_n_output_samples, "--lnL-offset", opts.cip_lnL_offset, + # LISA sources are massive black-hole binaries (M ~ 1e4-1e8 Msun), far + # above CIP's stellar-mass default (--M-max-cut 1e5) which would strip + # every grid point as "too massive". Likewise the synthetic high-SNR + # demo gives single-sample (n_eff=1) integrals, so relax CIP's own + # error cut (default 0.6) to keep those points. + "--M-max-cut", opts.cip_m_max_cut, + "--sigma-cut", opts.cip_sigma_cut, "--no-plots", ]) _write_cip_list(cip_args_list, [cip_line]) + # Puffball: between iterations, perturb the (very tight) CIP posterior so the + # next iteration's grid is not a near-degenerate cluster (which makes the CIP + # refit ill-conditioned and diverge). The CEPP wraps this with --inj-file / + # --inj-file-out; we supply the perturbation parameters + physical bounds. + puff_parts = [ + "--parameter", "mc", "--parameter", "eta", + "--puff-factor", opts.puff_factor, + "--mc-range", "[{},{}]".format(_mc * (1.0 - _w), _mc * (1.0 + _w)), + "--mtot-range", "[{},{}]".format(_mtot * (1.0 - _w), _mtot * (1.0 + _w)), + "--eta-range", "[{},{}]".format(_eta * (1.0 - _w), min(0.2499999, _eta * (1.0 + _w))), + ] + if opts.vary_sky: + # puff the sky too (phi/theta), else the next grid's sky collapses + puff_parts += ["--parameter", "phi", "--parameter", "theta"] + _write_arg_file(puff_args, puff_parts) + _write_arg_file(test_args, [ "--method", "lame", "--parameter", "mc", diff --git a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode_lisa b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode_lisa index 627819522..3754e5452 100755 --- a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode_lisa +++ b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode_lisa @@ -190,6 +190,7 @@ optp.add_option("-p", "--psd-file", action="append", help="instrument=psd-file, optp.add_option("-k", "--skymap-file", help="Use skymap stored in given FITS file.") optp.add_option("-x", "--coinc-xml", help="gstlal_inspiral XML file containing coincidence information.") optp.add_option("-I", "--sim-xml", help="XML file containing mass grid to be evaluated") +optp.add_option("--sim-grid", help="Hyperpipeline ASCII grid file (RIFT_HYPERPIPELINE_V1) of intrinsic points to be evaluated. Row index selected by --event. LISA-fork analogue of the main ILE --sim-grid; populates mass1/mass2/spin1z/spin2z and ecliptic_longitude/latitude from the selected row.") optp.add_option("-E", "--event", default=0,type=int, help="Event number used for this run") optp.add_option("--n-events-to-analyze", default=1,type=int, help="Number of events to analyze from this XML") optp.add_option("--soft-fail-event-range",action='store_true',help='Soft failure (exit 0) if event ID is out of range. This happens in pipelines, if we have pre-built a DAG attempting to analyze more points than we really have') @@ -521,6 +522,31 @@ print("######################################################################### # Struct to hold template parameters, P.deltaF comes from data (data.deltaF) P_list = None P=None # force allocation so I can use the preferred event later +if opts.sim_grid: + # Hyperpipeline intrinsic-grid handoff (the DAG/CEPP path passes one row per + # ILE worker via --sim-grid + --event). The LISA fork is single-point per + # job, so resolve the selected row here into the mass1/mass2/spin*/ecliptic + # opts and fall through to the standard single-point P construction below. + from RIFT.misc import hyperpipeline_io as _hpio + print("====Loading injection grid file:", opts.sim_grid, " event", opts.event, "=======") + _grid_arr, _grid_cols = _hpio.read_table(opts.sim_grid) + if opts.event is None or opts.event >= len(_grid_arr): + print(" Event index out of range for grid; soft exit") + sys.exit(0) + _row = _grid_arr[opts.event] + def _grid_get(name, default=0.0): + return float(_row[name]) if name in _grid_cols else default + opts.mass1 = _grid_get('m1') + opts.mass2 = _grid_get('m2') + opts.spin1z = _grid_get('a1z') + opts.spin2z = _grid_get('a2z') + if 'ecliptic_longitude' in _grid_cols: + opts.ecliptic_longitude = _grid_get('ecliptic_longitude') + if 'ecliptic_latitude' in _grid_cols: + opts.ecliptic_latitude = _grid_get('ecliptic_latitude') + print(" grid row -> m1={} m2={} s1z={} s2z={} ecl_long={} ecl_lat={}".format( + opts.mass1, opts.mass2, opts.spin1z, opts.spin2z, + opts.ecliptic_longitude, opts.ecliptic_latitude)) if opts.sim_xml: print(f"====Loading injection XML: {opts.sim_xml}, reading from event {opts.event} to event {opts.event+opts.n_events_to_analyze} =======") P_list = lalsimutils.xml_to_ChooseWaveformParams_array(str(opts.sim_xml)) @@ -1514,8 +1540,26 @@ def analyze_event_LISA(P_list, indx_event, data_dict, psd_dict, fmax, opts, inv_ event_id = -1 if opts.event == None: event_id = -1 - # Current response only applicable to quasicircular MBHB signals. Save sky location even if not varying. - numpy.savetxt(fname_output_txt, numpy.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y, P.s2z, lisa_sky_lamda, lisa_sky_beta, log_res+manual_avoid_overflow_logarithm, sqrt_var_over_res, sampler.ntotal, neff ]])) + # Hyperpipeline ASCII output (opt-in via RIFT_HYPERPIPELINE_FORMAT). The + # DAG/CEPP join (util_CleanILE_hyperpipeline.py) consumes self-describing + # header-bearing shards with lnL/sigma_lnL as columns 0/1. The legacy + # headerless savetxt below cannot be read by that path. + from RIFT.misc import hyperpipeline_io as _hpio + if _hpio.is_active(): + _cols = _hpio.build_column_list(use_sky=True) + _vals = { + "lnL": log_res + manual_avoid_overflow_logarithm, + "sigma_lnL": sqrt_var_over_res, + "m1": m1, "m2": m2, + "a1x": P.s1x, "a1y": P.s1y, "a1z": P.s1z, + "a2x": P.s2x, "a2y": P.s2y, "a2z": P.s2z, + "ecliptic_longitude": lisa_sky_lamda, + "ecliptic_latitude": lisa_sky_beta, + } + _hpio.write_row(fname_output_txt, _cols, [_vals.get(c, 0.0) for c in _cols]) + else: + # Current response only applicable to quasicircular MBHB signals. Save sky location even if not varying. + numpy.savetxt(fname_output_txt, numpy.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y, P.s2z, lisa_sky_lamda, lisa_sky_beta, log_res+manual_avoid_overflow_logarithm, sqrt_var_over_res, sampler.ntotal, neff ]])) # Comprehensive output (not yet provided) # Convert declination, inclination parameters in sampler if needed diff --git a/MonteCarloMarginalizeCode/Code/bin/util_CleanILE_hyperpipeline.py b/MonteCarloMarginalizeCode/Code/bin/util_CleanILE_hyperpipeline.py index 7f083dd8d..be580165a 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_CleanILE_hyperpipeline.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_CleanILE_hyperpipeline.py @@ -22,6 +22,7 @@ from __future__ import absolute_import, print_function import argparse +import os import sys # Import via an explicit file load so this script keeps working even when @@ -49,9 +50,11 @@ def main(argv=None): help="One or more hyperpipeline shard .dat files.") parser.add_argument("--output", "-o", default="-", help="Output filename, or '-' for stdout (default).") - parser.add_argument("--sigma-cut", type=float, default=0.9, + parser.add_argument("--sigma-cut", type=float, + default=float(os.environ.get("RIFT_ILE_SIGMA_CUT", 0.9)), help="Drop rows with sigma_lnL above this value " - "(default 0.9, mirrors util_CleanILE).") + "(default 0.9, mirrors util_CleanILE; override the " + "default via the RIFT_ILE_SIGMA_CUT env var).") parser.add_argument("--digits", type=int, default=5, help="Decimal-place precision used when grouping " "duplicate intrinsic rows (default 5, mirrors " diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ConstructEOSPosterior.py b/MonteCarloMarginalizeCode/Code/bin/util_ConstructEOSPosterior.py index e5513d23d..af9c1e0fa 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_ConstructEOSPosterior.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_ConstructEOSPosterior.py @@ -533,7 +533,10 @@ def fit_rf(x,y,y_errors=None,fname_export='nn_fit'): if y_errors is None: rf.fit(x,y) else: - rf.fit(x,y,sample_weight=1./y_errors**2) + # floor sigma so a zero error (placeholder rows can leak into the + # accumulated marg net with sigma=0) doesn't make sample_weight=1/sigma^2 + # infinite (sklearn rejects inf sample_weight). + rf.fit(x,y,sample_weight=1./np.maximum(np.asarray(y_errors,dtype=float),1e-3)**2) ### reject points with infinities : problems for inputs def fn_return(x_in,rf=rf): diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py index 0747c5e48..f2dcf7615 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py @@ -257,6 +257,8 @@ def extract_combination_from_LI(samples_LI, p): parser.add_argument("--mc-range",default=None,help="Chirp mass range [mc1,mc2]. Important if we have a low-mass object, to avoid wasting time sampling elsewhere.") parser.add_argument("--eta-range",default=None,help="Eta range. Important if we have a BNS or other item that has a strong constraint.") parser.add_argument("--mtot-range",default=None,help="Chirp mass range [mc1,mc2]. Important if we have a low-mass object, to avoid wasting time sampling elsewhere.") +parser.add_argument("--phi-range",default=None,help="LISA ecliptic-longitude (phi) range [lo,hi].") +parser.add_argument("--theta-range",default=None,help="LISA ecliptic-latitude (theta) range [lo,hi].") parser.add_argument("--trust-sample-parameter-box",action='store_true', help="If used, sets the prior range to the SAMPLE range for any parameters. NOT IMPLEMENTED. This should be automatically done for mc!") parser.add_argument("--plots-do-not-force-large-range",action='store_true', help = "If used, the plots do NOT automatically set the chieff range to [-1,1], the eta range to [0,1/4], etc") parser.add_argument("--downselect-parameter",action='append', help='Name of parameter to be used to eliminate grid points ') @@ -962,7 +964,10 @@ def normalized_zbar_prior(z): 'meanPerAno':meanPerAno_prior, 'chi_pavg':precession_prior, 'mu1': unnormalized_log_prior, - 'mu2': unnormalized_uniform_prior + 'mu2': unnormalized_uniform_prior, + # LISA ecliptic sky (phi=longitude, theta=latitude); uniform priors + 'phi': (lambda x: 1./(2*np.pi)), + 'theta': (lambda x: 1./np.pi), } prior_range_map = {"mtot": [1, 300], "q":[0.01,1], "s1z":[-0.999*chi_max,0.999*chi_max], "s2z":[-0.999*chi_small_max,0.999*chi_small_max], "mc":[0.9,250], "eta":[0.01,0.2499999],'delta_mc':[0,0.9], 'xi':[-chi_max,chi_max],'chi_eff':[-chi_max,chi_max],'delta':[-1,1], 's1x':[-chi_max,chi_max], @@ -998,12 +1003,18 @@ def normalized_zbar_prior(z): 'chi2_perp_u':[0,1], 's1z_bar':[-1,1], 's2z_bar':[-1,1], - 'mu1':[0.0001,1e3], # suboptimal, but something - 'mu2':[-300,1e3] + 'mu1':[0.0001,1e3], # suboptimal, but something + 'mu2':[-300,1e3], + 'phi':[0, 2*np.pi], # LISA ecliptic longitude (override via --phi-range) + 'theta':[-np.pi/2, np.pi/2] # LISA ecliptic latitude (override via --theta-range) } if not (opts.chiz_plus_range is None): print(" Warning: Overriding default chiz_plus range. USE WITH CARE", opts.chiz_plus_range) prior_range_map['chiz_plus']=eval(opts.chiz_plus_range) +if not (opts.phi_range is None): + prior_range_map['phi']=eval(opts.phi_range) +if not (opts.theta_range is None): + prior_range_map['theta']=eval(opts.theta_range) if not (opts.eta_range is None): print(f" Warning: Overriding default eta range to {eval(opts.eta_range)}. USE WITH CARE") @@ -1841,20 +1852,27 @@ def fit_gp_sparse(x): _use_tides = bool(opts.input_tides) or _has("lambda1") _use_eos = bool(opts.input_eos_index) or _has("eos_table_index") _use_dist = bool(opts.input_distance) or _has("distance") + # LISA sky: ecliptic_longitude/latitude are NAMED columns -> carry them + # through (aliased to P.phi/P.theta), so the sky is fit/sampled like any + # other coordinate. No positional all.net hacking. + _use_sky = _has("ecliptic_longitude") dat = _hpio.to_legacy_dat(_arr, use_eccentricity=_use_ecc, use_meanPerAno=_use_mpa, use_tides=_use_tides, use_eos_index=_use_eos, - use_distance=_use_dist) + use_distance=_use_dist, use_sky=_use_sky) _ix = _hpio.legacy_column_indices( use_eccentricity=_use_ecc, use_meanPerAno=_use_mpa, use_tides=_use_tides, use_eos_index=_use_eos, - use_distance=_use_dist) + use_distance=_use_dist, use_sky=_use_sky) col_lnL = _ix["lnL"] col_distance = _ix["distance"] col_lambda1 = _ix["lambda1"] col_eccentricity = _ix["eccentricity"] col_meanPerAno = _ix["meanPerAno"] + col_ecliptic_longitude = _ix["ecliptic_longitude"] + col_ecliptic_latitude = _ix["ecliptic_latitude"] else: + _use_sky = False dat = np.loadtxt(opts.fname) dat_orig = dat dat_orig = dat[dat[:,col_lnL].argsort()] # sort http://stackoverflow.com/questions/2828059/sorting-arrays-in-numpy-by-column @@ -1955,7 +1973,12 @@ def fit_gp_sparse(x): # P.meanPerAno = line[10] if opts.input_distance: P.dist = lal.PC_SI*1e6*line[col_distance] # 9. Previously incompatible with tides when hardcoded - + if _use_sky: + # LISA ecliptic sky stored as P.phi (longitude) / P.theta (latitude), + # matching hyperpipeline_io's column alias. Fit/sample as 'phi'/'theta'. + P.phi = line[col_ecliptic_longitude] + P.theta = line[col_ecliptic_latitude] + if opts.contingency_unevolved_neff == "quadpuff": P_copy = P.manual_copy() # prevent duplication P_list_in.append(P_copy) # store it, make sure distinct @@ -3026,7 +3049,7 @@ def parse_corr_params(my_str): _cols_out = _hpio.build_column_list( use_eccentricity=opts.use_eccentricity, use_meanPerAno=opts.use_meanPerAno, use_tides=opts.input_tides, use_eos_index=opts.input_eos_index, - use_distance=False) + use_distance=False, use_sky=_use_sky) _hpio.write_grid_from_P_list(opts.fname_output_samples, P_out_list[:n_output_size], _cols_out, @@ -3555,7 +3578,7 @@ def parse_corr_params(my_str): _cols_out = _hpio.build_column_list( use_eccentricity=opts.use_eccentricity, use_meanPerAno=opts.use_meanPerAno, use_tides=opts.input_tides, use_eos_index=opts.input_eos_index, - use_distance=False) + use_distance=False, use_sky=_use_sky) _hpio.write_grid_from_P_list(opts.fname_output_samples, P_list[:n_output_size], _cols_out, diff --git a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py index 835eafde1..281b99fd6 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py @@ -17,6 +17,8 @@ import numpy as np import argparse import os +import shlex +import subprocess import sys import lal import lalsimulation as lalsim @@ -154,6 +156,190 @@ def unsafe_parse_arg_string_dict(my_argstr): return dict_return +def _lisa_data_products_from_ini(opts): + """Fill the LISA data-product opts (channels / PSD files) from the conventional + production-ini sections so a LISA run can be driven by --use-ini like the + ground-based path. Scalars and lisa-* algorithm options are already populated + by the generic [rift-pseudo-pipe] parser (any CLI arg, by name); here we only + translate the per-channel *dict* products that don't map cleanly to a flat key: + + [data] channels = {'A': 'fake_strain', 'E': ..., 'T': ...} + [lalinference] psds = {'A': 'A_psd.xml.gz', ...} + + Values already set (e.g. via [rift-pseudo-pipe] lisa-channel-name) win, so the + flat CLI surface still overrides. Read-only; no LDG data-find is invoked. + """ + import configparser as _CfgP + cfg = _CfgP.ConfigParser() + cfg.optionxform = str + cfg.read(opts.use_ini) + + def _dict_to_assignments(section, key): + if not cfg.has_option(section, key): + return None + mapping = eval(cfg.get(section, key)) + return ["{}={}".format(ifo, val) for ifo, val in mapping.items()] + + if not opts.lisa_channel_name: + opts.lisa_channel_name = _dict_to_assignments("data", "channels") + if not opts.lisa_psd_file: + opts.lisa_psd_file = _dict_to_assignments("lalinference", "psds") + + +def run_lisa_known_sky_surface(opts): + if opts.approx is None: + print(" --lisa-known-sky requires --approx ") + sys.exit(1) + if opts.use_ini is not None: + # LISA production-ini path: scalars/algorithm options come from the + # generic [rift-pseudo-pipe] parser; fill the per-channel data products + # (channels, PSDs) from the conventional [data]/[lalinference] sections. + _lisa_data_products_from_ini(opts) + + bin_dir = os.path.dirname(os.path.abspath(__file__)) + helper = os.path.join(bin_dir, "helper_LISA_Events.py") + cepp = os.path.join(bin_dir, "create_event_parameter_pipeline_BasicIteration") + ile = os.path.join(bin_dir, "integrate_likelihood_extrinsic_batchmode_lisa") + + if opts.use_rundir: + workdir = os.path.abspath(opts.use_rundir) + else: + event_label = "manual_" + format_gps_time(opts.event_time) + sky_label = "variable_sky" if opts.lisa_vary_sky else "known_sky" + workdir = os.path.abspath( + event_label + "_LISA_" + opts.approx + "_" + sky_label + opts.manual_postfix + ) + os.makedirs(workdir, exist_ok=False) + + helper_cmd = [ + sys.executable, + helper, + "--working-directory", + workdir, + "--input-grid", + "proposed-grid.dat", + "--approximant", + opts.approx, + "--l-max", + str(opts.l_max), + "--event-time", + format_gps_time(opts.event_time), + "--cache-file", + opts.lisa_cache_file, + "--ecliptic-longitude", + str(opts.ecliptic_longitude), + "--ecliptic-latitude", + str(opts.ecliptic_latitude), + "--fmin-template", + str(opts.lisa_fmin_template), + "--fmax", + str(opts.lisa_fmax), + "--reference-freq", + str(opts.lisa_reference_freq), + "--srate", + str(opts.lisa_srate), + "--data-integration-window-half", + str(opts.lisa_data_integration_window_half), + "--grid-size", + str(opts.lisa_grid_size), + "--grid-fractional-width", + str(opts.lisa_grid_fractional_width), + "--sky-grid-width", + str(opts.lisa_sky_grid_width), + "--n-iterations", + str(opts.lisa_n_iterations), + "--n-samples-per-job", + str(opts.lisa_n_samples_per_job), + "--request-memory-ILE", + str(opts.internal_ile_request_memory), + "--request-memory-CIP", + str(opts.internal_cip_request_memory or 4096), + ] + if opts.lisa_vary_sky: + helper_cmd.append("--vary-sky") + if opts.lisa_zero_likelihood: + helper_cmd.append("--zero-likelihood") + for assignment in opts.lisa_channel_name or []: + helper_cmd.extend(["--channel-name", assignment]) + for assignment in opts.lisa_psd_file or []: + helper_cmd.extend(["--psd-file", assignment]) + if opts.extra_args_helper: + with open(opts.extra_args_helper) as extra: + helper_cmd.extend(shlex.split(extra.read())) + + print(" LISA known-sky helper command: ", " ".join(shlex.quote(x) for x in helper_cmd)) + subprocess.run(helper_cmd, check=True) + + if opts.lisa_skip_cepp_render: + print(" LISA helper bundle written in {}".format(workdir)) + return + + env = os.environ.copy() + env["RIFT_HYPERPIPELINE_FORMAT"] = "1" + env["PATH"] = bin_dir + os.pathsep + env.get("PATH", "") + env["PYTHONPATH"] = os.path.abspath(os.path.join(bin_dir, "..")) + os.pathsep + env.get("PYTHONPATH", "") + cepp_cmd = [ + sys.executable, + cepp, + "--ile-n-events-to-analyze", + "1", + "--input-grid", + os.path.join(workdir, "proposed-grid.dat"), + "--ile-exe", + ile, + "--ile-args", + os.path.join(workdir, "args_ile.txt"), + "--cip-args-list", + os.path.join(workdir, "args_cip_list.txt"), + "--test-args", + os.path.join(workdir, "args_test.txt"), + "--working-directory", + workdir, + "--n-iterations", + str(opts.lisa_n_iterations), + "--n-samples-per-job", + str(opts.lisa_n_samples_per_job), + "--n-copies", + str(opts.ile_copies), + "--request-memory-ILE", + str(opts.internal_ile_request_memory), + "--request-memory-CIP", + str(opts.internal_cip_request_memory or 4096), + "--transfer-file-list", + os.path.join(workdir, "helper_transfer_files.txt"), + ] + # Container: let write_ILE_sub_simple emit the singularity + file-transfer + # wiring (the LDG path's native mechanism) rather than any LISA-specific code. + # Needs SINGULARITY_RIFT_IMAGE (+ SINGULARITY_BASE_EXE_DIR) in the env. + if opts.lisa_use_singularity: + # write_ILE_sub_simple's singularity path requires the CEPP's --cache-file + # to be set (else "Need to specify frames_dir or cache_file to use + # singularity"); the LISA cache is otherwise only inside the ILE args. + cepp_cmd += ["--use-singularity", "--cache-file", opts.lisa_cache_file] + # Puffball between iterations: perturb the (very tight) CIP posterior so the + # next grid is not a near-degenerate cluster (else the CIP refit diverges). + if opts.lisa_n_iterations > 1 and not opts.lisa_no_puff: + cepp_cmd += [ + "--puff-exe", os.path.join(bin_dir, "util_ParameterPuffball.py"), + "--puff-args", os.path.join(workdir, "args_puff.txt"), + "--puff-cadence", "1", + "--puff-max-it", str(opts.lisa_n_iterations), + ] + # LISA reflected-sky-mode (vary-sky): reflect the grid at one iteration to + # explore the secondary sky mode (latitude bimodality). + if opts.lisa_search_reflected_sky_mode and opts.lisa_n_iterations > 1: + cepp_cmd += [ + "--search-reflected-sky-mode", + "--reflected-sky-mode-exe", os.path.join(bin_dir, "convert_primary_sky_mode_to_secondary"), + "--lisa-reference-time", str(opts.lisa_reference_time), + ] + if opts.lisa_search_reflected_sky_mode_iteration is not None: + cepp_cmd += ["--search-reflected-sky-mode-iteration", str(opts.lisa_search_reflected_sky_mode_iteration)] + print(" LISA known-sky CEPP command: ", " ".join(shlex.quote(x) for x in cepp_cmd)) + subprocess.run(cepp_cmd, check=True, cwd=workdir, env=env) + print(" LISA known-sky CEPP surface rendered in {}".format(workdir)) + + parser = argparse.ArgumentParser() parser.add_argument("--skip-reproducibility",action='store_true') @@ -200,6 +386,30 @@ def unsafe_parse_arg_string_dict(my_argstr): parser.add_argument("--use-legacy-gracedb",action='store_true') parser.add_argument("--internal-use-gracedb-bayestar",action='store_true',help="Retrieve BS skymap from gracedb (bayestar.fits), and use it internally in integration with --use-skymap bayestar.fits.") parser.add_argument("--event-time",default=None,type=float,help="Event time. Intended to override use of GracedbID. MUST provide --manual-initial-grid ") +parser.add_argument("--lisa-known-sky",action='store_true',help="Use the LISA helper to build a known-sky LISA CEPP surface and exit. Avoids the LDG event helper path.") +parser.add_argument("--lisa-vary-sky",action='store_true',help="With --lisa-known-sky, treat ecliptic sky location as intrinsic rather than pinning --lisa-fixed-sky.") +parser.add_argument("--lisa-skip-cepp-render",action='store_true',help="With --lisa-known-sky, only write the helper bundle; do not render the CEPP DAG.") +parser.add_argument("--lisa-cache-file",default="lisa.cache",help="With --lisa-known-sky, cache file passed to the LISA ILE.") +parser.add_argument("--lisa-channel-name",action="append",default=None,help="With --lisa-known-sky, channel assignment such as A=fake_strain. May be repeated.") +parser.add_argument("--lisa-psd-file",action="append",default=None,help="With --lisa-known-sky, PSD assignment such as A=A_psd.xml.gz. May be repeated.") +parser.add_argument("--ecliptic-longitude",default=1.0,type=float,help="With --lisa-known-sky, fixed ecliptic longitude.") +parser.add_argument("--ecliptic-latitude",default=0.3,type=float,help="With --lisa-known-sky, fixed ecliptic latitude.") +parser.add_argument("--lisa-fmin-template",default=1.0e-3,type=float,help="With --lisa-known-sky, template low-frequency cutoff.") +parser.add_argument("--lisa-fmax",default=0.125,type=float,help="With --lisa-known-sky, high-frequency cutoff.") +parser.add_argument("--lisa-reference-freq",default=5.0e-3,type=float,help="With --lisa-known-sky, waveform reference frequency.") +parser.add_argument("--lisa-srate",default=0.25,type=float,help="With --lisa-known-sky, sample rate. Kept as float for long-duration LISA data.") +parser.add_argument("--lisa-data-integration-window-half",default=300.0,type=float,help="With --lisa-known-sky, half-width of the ILE data integration window.") +parser.add_argument("--lisa-no-puff",action="store_true",help="Disable the inter-iteration puffball for the known-sky LISA path.") +parser.add_argument("--lisa-search-reflected-sky-mode",action="store_true",help="LISA vary-sky: at one iteration, reflect the grid to the secondary sky mode (handles the LISA latitude bimodality).") +parser.add_argument("--lisa-search-reflected-sky-mode-iteration",default=None,type=int,help="Iteration to reflect the sky (default n_iterations-2).") +parser.add_argument("--lisa-reference-time",default=0.0,type=float,help="LISA coalescence/reference time (for the reflected-sky transform).") +parser.add_argument("--lisa-use-singularity",action="store_true",help="Forward --use-singularity to the CEPP for the known-sky LISA path. The container wiring + transfer is then emitted by write_ILE_sub_simple, exactly as for the LDG path; set SINGULARITY_RIFT_IMAGE (osdf:// staged image preferred, so dag_utils file-transfers it) and SINGULARITY_BASE_EXE_DIR (dir of the LISA ILE *inside* the image).") +parser.add_argument("--lisa-grid-size",default=3,type=int,help="With --lisa-known-sky, number of synthetic initial-grid points.") +parser.add_argument("--lisa-grid-fractional-width",default=1.0e-3,type=float,help="With --lisa-known-sky, fractional mass width for the initial grid.") +parser.add_argument("--lisa-sky-grid-width",default=0.02,type=float,help="With --lisa-known-sky --lisa-vary-sky, ecliptic sky half-step scale for the initial grid.") +parser.add_argument("--lisa-n-iterations",default=1,type=int,help="With --lisa-known-sky, CEPP iteration count.") +parser.add_argument("--lisa-n-samples-per-job",default=1,type=int,help="With --lisa-known-sky, CEPP samples per job.") +parser.add_argument("--lisa-zero-likelihood",action='store_true',help="With --lisa-known-sky, pass --zero-likelihood through to the LISA ILE args.") parser.add_argument("--calibration",default="C00",type=str) parser.add_argument("--playground-data",action='store_true', help="Passed through to helper_LDG_events, and changes name prefix") parser.add_argument("--approx",default=None,type=str,help="Approximant. REQUIRED") @@ -421,6 +631,11 @@ def unsafe_parse_arg_string_dict(my_argstr): ile_condor_commands.append([item, val]) +if opts.lisa_known_sky: + run_lisa_known_sky_surface(opts) + sys.exit(0) + + if opts.use_osg: opts.condor_nogrid_nonworker = True # note we ALSO have to check this if we set use_osg in the ini file! Moved statement so flagged diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/lisa/BBH_lisa_demo.ini b/MonteCarloMarginalizeCode/Code/demo/rift/lisa/BBH_lisa_demo.ini new file mode 100644 index 000000000..5a2c76987 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/demo/rift/lisa/BBH_lisa_demo.ini @@ -0,0 +1,54 @@ +########################################################################################### +# LISA toy-demo config, production-ini form. +# +# Drives the SAME validated LISA known-sky workflow as the --lisa-* CLI flags, +# but through the production entry point: +# +# util_RIFT_pseudo_pipe.py --use-ini BBH_lisa_demo.ini --use-rundir +# +# Per-channel data products (channels, PSDs) are read from the conventional +# [data]/[lalinference] sections; every other option flows through the generic +# [rift-pseudo-pipe] parser by CLI-arg name (s/-/_/). Values are eval()'d, so +# strings/paths must be quoted and booleans are bare True/False. +# +# ConfigParser does no path interpolation: replace __BUNDLE_DIR__ below with the +# absolute path of the bundle written by make_synthetic_lisa_inputs.py (the dir +# holding lisa.cache and {A,E,T}_psd.xml.gz). run_lisa_ini_demo.sh does this for +# you. Stays IMRPhenomD -- toy demo. +########################################################################################### + +[data] +# strain channel name carried in the cache / h5 frames, per LISA channel A,E,T +channels = {'A': 'fake_strain', 'E': 'fake_strain', 'T': 'fake_strain'} + +[lalinference] +# per-channel PSD xml(.gz) +psds = {'A': '__BUNDLE_DIR__/A_psd.xml.gz', 'E': '__BUNDLE_DIR__/E_psd.xml.gz', 'T': '__BUNDLE_DIR__/T_psd.xml.gz'} + +[rift-pseudo-pipe] +### route to the LISA known-sky builder ### +lisa-known-sky=True +lisa-vary-sky=True +lisa-search-reflected-sky-mode=True +lisa-reference-time=0 + +### source / sky (injected ecliptic location) ### +approx="IMRPhenomD" +event-time=0 +ecliptic-longitude=1.0 +ecliptic-latitude=0.3 + +### data products / band ### +lisa-cache-file="__BUNDLE_DIR__/lisa.cache" +lisa-srate=0.25 +lisa-fmin-template=0.0001 +lisa-fmax=0.125 +lisa-reference-freq=0.005 +lisa-data-integration-window-half=300 + +### grid / iteration ### +lisa-grid-size=16 +lisa-n-iterations=4 +lisa-n-samples-per-job=16 +internal-ile-request-memory=2048 +internal-cip-request-memory=2048 diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/lisa/README.md b/MonteCarloMarginalizeCode/Code/demo/rift/lisa/README.md new file mode 100644 index 000000000..1254564f3 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/demo/rift/lisa/README.md @@ -0,0 +1,124 @@ +# LISA zero-likelihood CEPP demo + +This directory contains the first checked-in LISA/RIFT workflow scaffold. It is +intended to be small enough for CI contract testing while still looking like the +shape of a real LISA pipeline run. + +The demo uses `helper_LISA_Events.py` to write the files consumed by +`create_event_parameter_pipeline_BasicIteration`: + +- `proposed-grid.dat` +- `args_ile.txt` +- `args_cip_list.txt` +- `args_test.txt` +- `helper_transfer_files.txt` +- `command-cepp-lisa.sh` + +The default path uses `--zero-likelihood`, so it validates file formats, +hyperpipeline sky columns, executable handoffs, and DAG rendering without +requiring real Sangria frames or PSD products. + +## Run + +From a checkout with the normal RIFT runtime available: + +```bash +./MonteCarloMarginalizeCode/Code/demo/rift/lisa/run_lisa_zero_likelihood_cepp.sh +``` + +Useful environment overrides: + +- `RIFT_LISA_PYTHON`: Python executable to use. +- `RIFT_LISA_WORKDIR`: output directory; defaults under `/tmp`. +- `RIFT_LISA_RENDER_CEPP=0`: only write the helper bundle; do not render CEPP. + +Additional arguments are passed through to `helper_LISA_Events.py`, for example: + +```bash +RIFT_LISA_WORKDIR=/tmp/rift-lisa-demo \ + ./MonteCarloMarginalizeCode/Code/demo/rift/lisa/run_lisa_zero_likelihood_cepp.sh \ + --mass1 120000 --mass2 90000 --ecliptic-longitude 1.4 +``` + +Expected output is a helper bundle and, unless `RIFT_LISA_RENDER_CEPP=0`, CEPP +submit/DAG files in the work directory. This script does not submit the DAG. + +## Synthetic ILE surface + +The companion script builds tiny synthetic A/E/T inputs and runs the standalone +LISA ILE against them: + +```bash +./MonteCarloMarginalizeCode/Code/demo/rift/lisa/run_lisa_synthetic_ile.sh +``` + +It writes frequency-domain HDF5 frames, a `lisa.cache`, flat XML PSDs, helper +contract files, and `lisa_ile_0_.dat`. The run pins most extrinsic parameters +but leaves polarization open, so it exercises a nonzero LISA likelihood integral +without becoming a full PE run. + +This surface currently still uses XML PSDs because that is what the ILE path +loads today. The synthetic input builder keeps PSD generation local and +mechanical, making it a good target for a future ASCII-PSD path once the LISA +workflow no longer needs `lal.series` XML PSD documents. + +## Analytic PSD products + +For a closer analogue of the toy `generate_iligo_psd` examples, this directory +also provides: + +```bash +./MonteCarloMarginalizeCode/Code/demo/rift/lisa/make_lisa_psds.py \ + --output-directory /tmp/rift-lisa-psds --write-ascii +``` + +It writes analytic LISA A/E/T XML PSDs and, optionally, `LISA_psd.txt`. The +PP-style LISA surface in `MonteCarloMarginalizeCode/Code/test/pp_lisa` uses +this generator together with the synthetic frame builder and +`util_RIFT_pseudo_pipe.py --lisa-known-sky`. + +## Heavyweight end-to-end demo + +For a local end-to-end LISA exercise that uses synthetic data, analytic PSDs, +`pseudo_pipe`, CEPP rendering, a real direct ILE likelihood evaluation, and the +LISA run diagnostic: + +```bash +./MonteCarloMarginalizeCode/Code/demo/rift/lisa/run_lisa_end_to_end.sh +``` + +Useful environment overrides: + +- `RIFT_LISA_WORKDIR`: output directory; defaults under `/tmp`. +- `RIFT_LISA_RUN_ILE=0`: render inputs and CEPP files only. +- `RIFT_LISA_VARY_SKY=1`: render the `pseudo_pipe` surface with sky left as an + intrinsic hyperpipeline parameter. The direct ILE check remains fixed-sky so + the heavyweight demo has a stable, fast likelihood evaluation. + +Expected heavyweight outputs include `event_0/` synthetic data products, +`analysis_event_0/` CEPP files, `lisa_end_to_end_0_.dat`, and +`lisa_end_to_end_summary.json`. + +## Production-ini path + +The same known-sky workflow can be driven from an `.ini` file through the +production entry point, instead of `--lisa-*` CLI flags: + +```bash +util_RIFT_pseudo_pipe.py --use-ini BBH_lisa_demo.ini --use-rundir +``` + +`BBH_lisa_demo.ini` is a toy template (IMRPhenomD, vary-sky + reflected-sky). +Per-channel data products are read from the conventional `[data]` (`channels`) +and `[lalinference]` (`psds`) sections; every other option flows through the +generic `[rift-pseudo-pipe]` parser by CLI-arg name. ConfigParser does no path +interpolation, so substitute the `__BUNDLE_DIR__` placeholder for the absolute +bundle path. The convenience script does the build + substitute + render: + +```bash +./MonteCarloMarginalizeCode/Code/demo/rift/lisa/run_lisa_ini_demo.sh +``` + +The ini path is a thin front-end over the validated CLI machinery: it renders a +byte-identical workflow (see `test/test_lisa_ini_contract.py`), so it does not +submit and stays a drop-in alternative to the CLI form. diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/lisa/make_lisa_psds.py b/MonteCarloMarginalizeCode/Code/demo/rift/lisa/make_lisa_psds.py new file mode 100755 index 000000000..dae0f0b41 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/demo/rift/lisa/make_lisa_psds.py @@ -0,0 +1,75 @@ +#!/usr/bin/env python + +"""Generate lightweight analytic LISA PSD products for demo/test workflows.""" + +from argparse import ArgumentParser +import os + +from igwn_ligolw import utils as ligolw_utils +import lal +import lal.series +import numpy as np + +from RIFT.LISA.psd_generation import generate_LISA_psd +import RIFT.lalsimutils as lalsimutils + + +def build_parser(): + parser = ArgumentParser() + parser.add_argument("--output-directory", default=".") + parser.add_argument("--channels", default="A,E,T") + parser.add_argument("--fmax", type=float, default=0.125) + parser.add_argument("--npts", type=int, default=513) + parser.add_argument("--Tobs", type=float, default=0.5, help="Observation time in years.") + parser.add_argument("--NC", type=int, default=3, help="Number of LISA channels.") + parser.add_argument("--write-ascii", action="store_true") + return parser + + +def write_psd_xml(channel, path, fvals, psd_values): + deltaF = fvals[1] - fvals[0] + psd = lal.CreateREAL8FrequencySeries( + channel, + lal.LIGOTimeGPS(0), + fvals[0], + deltaF, + lalsimutils.lsu_HertzUnit, + len(fvals), + ) + psd.data.data[:] = psd_values + xmldoc = lal.series.make_psd_xmldoc({channel: psd}) + xmldoc.childNodes[0].attributes._attrs = {"Name": "psd"} + ligolw_utils.write_filename(xmldoc, path, compress="gz") + + +def main(argv=None): + opts = build_parser().parse_args(argv) + output_directory = os.path.abspath(opts.output_directory) + os.makedirs(output_directory, exist_ok=True) + + if opts.npts < 2: + raise ValueError("--npts must be at least 2") + fvals = np.linspace(0.0, opts.fmax, opts.npts) + positive_fvals, positive_psd = generate_LISA_psd.generate_psd( + fmin=fvals[1], + fmax=opts.fmax, + Tobs_years=opts.Tobs, + NC=opts.NC, + npts=opts.npts - 1, + ) + psd_values = np.zeros_like(fvals) + psd_values[1:] = positive_psd + + for channel in [item.strip() for item in opts.channels.split(",") if item.strip()]: + xml_path = os.path.join(output_directory, "{}_psd.xml.gz".format(channel)) + write_psd_xml(channel, xml_path, fvals, psd_values) + print("Wrote {}".format(xml_path)) + + if opts.write_ascii: + ascii_path = os.path.join(output_directory, "LISA_psd.txt") + np.savetxt(ascii_path, np.column_stack([positive_fvals, positive_psd])) + print("Wrote {}".format(ascii_path)) + + +if __name__ == "__main__": + main() diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/lisa/make_synthetic_lisa_inputs.py b/MonteCarloMarginalizeCode/Code/demo/rift/lisa/make_synthetic_lisa_inputs.py new file mode 100755 index 000000000..9905e1b81 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/demo/rift/lisa/make_synthetic_lisa_inputs.py @@ -0,0 +1,195 @@ +#!/usr/bin/env python + +"""Generate tiny synthetic LISA inputs for the demo analysis surface.""" + +from argparse import ArgumentParser +import os + +from igwn_ligolw import utils as ligolw_utils +import lal +import lal.series +import lalsimulation as lalsim +import numpy as np + +import RIFT.LISA.lalsimutils_compat as lisa_lalsimutils_compat +import RIFT.lalsimutils as lalsimutils +from RIFT.LISA.response import LISA_response +from RIFT.LISA.psd_generation import generate_LISA_psd + + +def build_parser(): + parser = ArgumentParser() + parser.add_argument("--output-directory", default=".") + parser.add_argument("--mass1", type=float, default=1.0e5) + parser.add_argument("--mass2", type=float, default=8.0e4) + parser.add_argument("--spin1z", type=float, default=0.1) + parser.add_argument("--spin2z", type=float, default=-0.05) + parser.add_argument("--distance-mpc", type=float, default=1.0e3) + parser.add_argument("--fmin", type=float, default=1.0e-3) + parser.add_argument("--fref", type=float, default=5.0e-3) + parser.add_argument("--fmax", type=float, default=0.125) + parser.add_argument("--deltaT", type=float, default=4.0) + parser.add_argument("--duration", type=float, default=4096.0) + parser.add_argument("--ecliptic-latitude", type=float, default=0.3) + parser.add_argument("--ecliptic-longitude", type=float, default=1.0) + parser.add_argument("--psi", type=float, default=0.2) + parser.add_argument("--inclination", type=float, default=0.4) + parser.add_argument("--phiref", type=float, default=0.1) + parser.add_argument("--psd-level", type=float, default=1.0e-40) + return parser + + +def synthetic_params(opts): + P = lalsimutils.ChooseWaveformParams() + P.m1 = opts.mass1 * lal.MSUN_SI + P.m2 = opts.mass2 * lal.MSUN_SI + P.s1z = opts.spin1z + P.s2z = opts.spin2z + P.dist = opts.distance_mpc * 1.0e6 * lal.PC_SI + P.fmin = opts.fmin + P.fref = opts.fref + P.fmax = opts.fmax + P.deltaT = float(opts.deltaT) + P.deltaF = 1.0 / float(opts.duration) + P.approx = lalsim.IMRPhenomD + P.theta = opts.ecliptic_latitude + P.phi = opts.ecliptic_longitude + P.psi = opts.psi + P.incl = opts.inclination + P.phiref = opts.phiref + return P + + +def write_cache(output_directory): + cache_path = os.path.join(output_directory, "lisa.cache") + rows = [ + ("A", "A-fake_strain-1000000-10000.h5"), + ("E", "E-fake_strain-1000000-10000.h5"), + ("T", "T-fake_strain-1000000-10000.h5"), + ] + with open(cache_path, "w") as out: + for channel, filename in rows: + path = os.path.abspath(os.path.join(output_directory, filename)) + out.write(f"{channel} {channel} 0 1 file://localhost{path}\n") + return cache_path + + +def _ldc_inference_psd(channel, fvals): + """The PSD people actually use for INFERENCE: the LISA Data Challenge noise + model (per A/E/T channel), via RIFT.LISA.utils.utils.get_ldc_psds. Requires + the optional `ldc` package. Returns None (so the caller falls back to the + analytic sensitivity curve) if ldc is unavailable -- the analytic curve is + sky-averaged and strictly a mismatch tool, but it matches the SciRDv1 model + to <~1% across the sensitive band, so it is an adequate CI proxy. (Guidance + from A. Jan; the production runs use the ldc PSD.) + """ + try: + import ldc.lisa.noise as _noise + except Exception: + return None + nm = _noise.get_noise_model("SciRDv1", fvals[1:]) + out = np.empty(len(fvals)) + out[0] = np.inf + out[1:] = nm.psd(fvals[1:], channel) + return out + + +def write_lisa_psd_xml(channel, output_directory, deltaF, length, Tobs_years=0.5, NC=3): + """Write a LISA PSD on the SAME frequency grid as the data (f0=0, this + deltaF, this length). + + Prefers the ldc INFERENCE PSD (per-channel SciRDv1); falls back to the + analytic sky-averaged sensitivity curve when `ldc` is not installed. + + A flat PSD is unphysical for LISA: at the SNRs of MBHBs the steeply-rising + low-frequency noise must be modelled or the low-f content dominates and the + posterior collapses to a delta (extrinsic eff_samp -> 1). Keeping the PSD on + the data's own grid also avoids ILE interpolation across a mismatched grid. + """ + fvals = np.arange(length) * deltaF + psd_values = _ldc_inference_psd(channel, fvals) + if psd_values is None: + R_exists, interp_func = generate_LISA_psd.response_interpolant(NC) + psd_values = np.empty(length) + # DC bin (f=0) is not physical for the sensitivity curve; set it huge so + # it carries zero weight in the likelihood. + psd_values[0] = np.inf + psd_values[1:] = generate_LISA_psd.Sn( + fvals[1:], Tobs_years * lal.YRSID_SI, NC, R_exists, interp_func) + psd = lal.CreateREAL8FrequencySeries( + channel, + lal.LIGOTimeGPS(0), + 0.0, + deltaF, + lalsimutils.lsu_HertzUnit, + length, + ) + psd.data.data[:] = psd_values + xmldoc = lal.series.make_psd_xmldoc({channel: psd}) + xmldoc.childNodes[0].attributes._attrs = {"Name": "psd"} + path = os.path.join(output_directory, f"{channel}_psd.xml.gz") + ligolw_utils.write_filename(xmldoc, path, compress="gz") + return path + + +def main(argv=None): + opts = build_parser().parse_args(argv) + output_directory = os.path.abspath(opts.output_directory) + os.makedirs(output_directory, exist_ok=True) + + P = synthetic_params(opts) + modes = [(2, 2)] + hlms = lisa_lalsimutils_compat.hlmoff_for_LISA(P, Lmax=2, modes=modes) + data = LISA_response.create_lisa_injections( + hlms, + P.fmax, + P.fref, + P.theta, + P.phi, + P.psi, + P.incl, + P.phiref, + tref=0.0, + ) + LISA_response.create_h5_files_from_data_dict(data, output_directory) + cache_path = write_cache(output_directory) + + psd_paths = {} + for channel, channel_data in data.items(): + psd_paths[channel] = write_lisa_psd_xml( + channel, + output_directory, + channel_data.deltaF, + channel_data.data.length, + ) + + summary_path = os.path.join(output_directory, "synthetic-params.env") + with open(summary_path, "w") as out: + out.write(f"MASS1={opts.mass1}\n") + out.write(f"MASS2={opts.mass2}\n") + out.write(f"SPIN1Z={opts.spin1z}\n") + out.write(f"SPIN2Z={opts.spin2z}\n") + out.write(f"DISTANCE_MPC={opts.distance_mpc}\n") + out.write(f"SRATE={1.0 / P.deltaT}\n") + out.write(f"DELTA_T={P.deltaT}\n") + out.write(f"DURATION={1.0 / P.deltaF}\n") + out.write(f"DELTA_F={P.deltaF}\n") + out.write(f"FMIN={P.fmin}\n") + out.write(f"FREF={P.fref}\n") + out.write(f"FMAX={P.fmax}\n") + out.write(f"ECLIPTIC_LATITUDE={P.theta}\n") + out.write(f"ECLIPTIC_LONGITUDE={P.phi}\n") + out.write(f"PSI={P.psi}\n") + out.write(f"INCLINATION={P.incl}\n") + out.write(f"PHIREF={P.phiref}\n") + out.write(f"CACHE_FILE={cache_path}\n") + for channel in ["A", "E", "T"]: + out.write(f"{channel}_PSD={psd_paths[channel]}\n") + + print(f"Wrote synthetic LISA inputs in {output_directory}") + print(f" cache: {cache_path}") + print(f" params: {summary_path}") + + +if __name__ == "__main__": + main() diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/lisa/run_lisa_end_to_end.sh b/MonteCarloMarginalizeCode/Code/demo/rift/lisa/run_lisa_end_to_end.sh new file mode 100755 index 000000000..c0947b935 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/demo/rift/lisa/run_lisa_end_to_end.sh @@ -0,0 +1,133 @@ +#!/usr/bin/env bash +set -euo pipefail + +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +CODE_DIR="$(cd "${SCRIPT_DIR}/../../.." && pwd)" +BIN_DIR="${CODE_DIR}/bin" +CHECK_SCRIPT="${CODE_DIR}/RIFT/LISA/run_checks/plot_RIFT.py" + +PYTHON_BIN="${RIFT_LISA_PYTHON:-}" +if [[ -z "${PYTHON_BIN}" ]]; then + PYTHON_BIN="$(command -v python3)" +fi + +WORKDIR="${RIFT_LISA_WORKDIR:-/tmp/rift-lisa-end-to-end-$(date +%s)}" +BUNDLE_DIR="${WORKDIR}/event_0" +RUNDIR="${WORKDIR}/analysis_event_0" +RUN_ILE="${RIFT_LISA_RUN_ILE:-1}" +VARY_SKY="${RIFT_LISA_VARY_SKY:-0}" + +SKY_ARGS=() +if [[ "${VARY_SKY}" == "1" ]]; then + SKY_ARGS=(--lisa-vary-sky --lisa-grid-size 3 --lisa-sky-grid-width 0.01) +else + SKY_ARGS=(--lisa-grid-size 1) +fi + +mkdir -p "${BUNDLE_DIR}" + +export PYTHONPATH="${CODE_DIR}${PYTHONPATH:+:${PYTHONPATH}}" +export PATH="${BIN_DIR}:${PATH}" +export MPLCONFIGDIR="${WORKDIR}/.matplotlib" +export RIFT_HYPERPIPELINE_FORMAT=1 + +"${PYTHON_BIN}" "${SCRIPT_DIR}/make_synthetic_lisa_inputs.py" \ + --output-directory "${BUNDLE_DIR}" \ + --duration 1024 \ + --deltaT 4 \ + "$@" + +"${PYTHON_BIN}" "${SCRIPT_DIR}/make_lisa_psds.py" \ + --output-directory "${BUNDLE_DIR}" \ + --fmax 0.125 \ + --npts 513 \ + --write-ascii + +set -a +source "${BUNDLE_DIR}/synthetic-params.env" +set +a + +"${PYTHON_BIN}" "${BIN_DIR}/util_RIFT_pseudo_pipe.py" \ + --lisa-known-sky \ + --use-rundir "${RUNDIR}" \ + --approx IMRPhenomD \ + --event-time 0 \ + --ecliptic-longitude "${ECLIPTIC_LONGITUDE}" \ + --ecliptic-latitude "${ECLIPTIC_LATITUDE}" \ + --lisa-cache-file "${CACHE_FILE}" \ + --lisa-channel-name A=fake_strain \ + --lisa-channel-name E=fake_strain \ + --lisa-channel-name T=fake_strain \ + --lisa-psd-file "A=${BUNDLE_DIR}/A_psd.xml.gz" \ + --lisa-psd-file "E=${BUNDLE_DIR}/E_psd.xml.gz" \ + --lisa-psd-file "T=${BUNDLE_DIR}/T_psd.xml.gz" \ + --lisa-srate "${SRATE}" \ + --lisa-fmin-template "${FMIN}" \ + --lisa-fmax "${FMAX}" \ + --lisa-reference-freq "${FREF}" \ + "${SKY_ARGS[@]}" \ + --lisa-n-iterations 1 \ + --lisa-n-samples-per-job 1 \ + --internal-ile-request-memory 1024 \ + --internal-cip-request-memory 1024 + +if [[ "${RUN_ILE}" == "0" ]]; then + echo "Rendered LISA end-to-end demo products in ${WORKDIR}" + exit 0 +fi + +"${PYTHON_BIN}" "${BIN_DIR}/integrate_likelihood_extrinsic_batchmode_lisa" \ + --LISA \ + --h5-frame-FD \ + --time-marginalization \ + --lisa-fixed-sky 1 \ + --ecliptic-longitude "${ECLIPTIC_LONGITUDE}" \ + --ecliptic-latitude "${ECLIPTIC_LATITUDE}" \ + --lisa-reference-time 0 \ + --lisa-reference-frequency "${FREF}" \ + --data-integration-window-half 8 \ + --modes "[(2,2)]" \ + --cache-file "${CACHE_FILE}" \ + --channel-name A=fake_strain \ + --channel-name E=fake_strain \ + --channel-name T=fake_strain \ + --psd-file "A=${BUNDLE_DIR}/A_psd.xml.gz" \ + --psd-file "E=${BUNDLE_DIR}/E_psd.xml.gz" \ + --psd-file "T=${BUNDLE_DIR}/T_psd.xml.gz" \ + --fmin-template "${FMIN}" \ + --fmin-ifo "A=${FMIN}" \ + --fmin-ifo "E=${FMIN}" \ + --fmin-ifo "T=${FMIN}" \ + --fmax "${FMAX}" \ + --reference-freq "${FREF}" \ + --srate "${SRATE}" \ + --l-max 2 \ + --approx IMRPhenomD \ + --mass1 "${MASS1}" \ + --mass2 "${MASS2}" \ + --spin1z "${SPIN1Z}" \ + --spin2z "${SPIN2Z}" \ + --d-max 5000 \ + --d-min 1 \ + --n-eff 2 \ + --n-max 40 \ + --n-chunk 20 \ + --save-P 1 \ + --no-adapt \ + --internal-use-lnL \ + --sampler-method AV \ + --inclination "${INCLINATION}" \ + --phi-orb "${PHIREF}" \ + --distance "${DISTANCE_MPC:-1000}" \ + --right-ascension 0 \ + --declination 0 \ + --internal-hard-fail-on-error \ + --output-file "${WORKDIR}/lisa_end_to_end" + +test -s "${WORKDIR}/lisa_end_to_end_0_.dat" +"${PYTHON_BIN}" "${CHECK_SCRIPT}" "${WORKDIR}/lisa_end_to_end_0_.dat" --json \ + > "${WORKDIR}/lisa_end_to_end_summary.json" + +echo "LISA end-to-end demo products: ${WORKDIR}" +echo "LISA end-to-end ILE output: ${WORKDIR}/lisa_end_to_end_0_.dat" +echo "LISA end-to-end summary: ${WORKDIR}/lisa_end_to_end_summary.json" diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/lisa/run_lisa_ini_demo.sh b/MonteCarloMarginalizeCode/Code/demo/rift/lisa/run_lisa_ini_demo.sh new file mode 100755 index 000000000..5c62c2b6f --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/demo/rift/lisa/run_lisa_ini_demo.sh @@ -0,0 +1,39 @@ +#!/usr/bin/env bash +# Render the LISA toy demo through the PRODUCTION-INI entry point +# (util_RIFT_pseudo_pipe.py --use-ini), to show the .ini path drives the same +# workflow as the --lisa-* CLI flags. Builds synthetic A/E/T inputs, instantiates +# BBH_lisa_demo.ini against that bundle, and renders the CEPP DAG. Does NOT submit. +set -euo pipefail + +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +CODE_DIR="$(cd "${SCRIPT_DIR}/../../.." && pwd)" +BIN_DIR="${CODE_DIR}/bin" + +PYTHON_BIN="${RIFT_LISA_PYTHON:-}" +if [[ -z "${PYTHON_BIN}" ]]; then PYTHON_BIN="$(command -v python3)"; fi + +WORKDIR="${RIFT_LISA_WORKDIR:-/tmp/rift-lisa-ini-$(date +%s)}" +BUNDLE_DIR="${WORKDIR}/event_0" +RUNDIR="${WORKDIR}/analysis_event_0" +mkdir -p "${BUNDLE_DIR}" + +export PYTHONPATH="${CODE_DIR}${PYTHONPATH:+:${PYTHONPATH}}" +export PATH="${BIN_DIR}:${PATH}" +export RIFT_HYPERPIPELINE_FORMAT=1 + +# Tiny synthetic A/E/T inputs (real grid-matched PSD + lisa.cache + frames). +"${PYTHON_BIN}" "${SCRIPT_DIR}/make_synthetic_lisa_inputs.py" \ + --output-directory "${BUNDLE_DIR}" --duration 16384 --deltaT 4 \ + --fmin 1e-4 --distance-mpc 20000 "$@" + +# Instantiate the template ini against this bundle (ConfigParser has no path +# interpolation, so substitute the placeholder). +INI="${WORKDIR}/BBH_lisa_demo.ini" +sed "s#__BUNDLE_DIR__#${BUNDLE_DIR}#g" "${SCRIPT_DIR}/BBH_lisa_demo.ini" > "${INI}" + +# Render the workflow from the ini. +"${PYTHON_BIN}" "${BIN_DIR}/util_RIFT_pseudo_pipe.py" \ + --use-ini "${INI}" --use-rundir "${RUNDIR}" + +echo "=== rendered LISA ini workflow in ${RUNDIR} ===" +ls "${RUNDIR}"/*.dag "${RUNDIR}"/args_ile.txt "${RUNDIR}"/args_cip_list.txt diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/lisa/run_lisa_synthetic_ile.sh b/MonteCarloMarginalizeCode/Code/demo/rift/lisa/run_lisa_synthetic_ile.sh new file mode 100755 index 000000000..275489124 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/demo/rift/lisa/run_lisa_synthetic_ile.sh @@ -0,0 +1,114 @@ +#!/usr/bin/env bash +set -euo pipefail + +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +CODE_DIR="$(cd "${SCRIPT_DIR}/../../.." && pwd)" +BIN_DIR="${CODE_DIR}/bin" + +PYTHON_BIN="${RIFT_LISA_PYTHON:-}" +if [[ -z "${PYTHON_BIN}" ]]; then + PYTHON_BIN="$(command -v python3)" +fi + +WORKDIR="${RIFT_LISA_WORKDIR:-/tmp/rift-lisa-synthetic-ile-$(date +%s)}" +RUN_ILE="${RIFT_LISA_RUN_ILE:-1}" + +mkdir -p "${WORKDIR}" + +export PYTHONPATH="${CODE_DIR}${PYTHONPATH:+:${PYTHONPATH}}" +export PATH="${BIN_DIR}:${PATH}" +export MPLCONFIGDIR="${WORKDIR}/.matplotlib" + +"${PYTHON_BIN}" "${SCRIPT_DIR}/make_synthetic_lisa_inputs.py" \ + --output-directory "${WORKDIR}" \ + "$@" + +set -a +source "${WORKDIR}/synthetic-params.env" +set +a + +"${PYTHON_BIN}" "${BIN_DIR}/helper_LISA_Events.py" \ + --working-directory "${WORKDIR}" \ + --cache-file "${CACHE_FILE}" \ + --psd-file "A=${A_PSD}" \ + --psd-file "E=${E_PSD}" \ + --psd-file "T=${T_PSD}" \ + --mass1 "${MASS1}" \ + --mass2 "${MASS2}" \ + --spin1z "${SPIN1Z}" \ + --spin2z "${SPIN2Z}" \ + --ecliptic-latitude "${ECLIPTIC_LATITUDE}" \ + --ecliptic-longitude "${ECLIPTIC_LONGITUDE}" \ + --fmin-template "${FMIN}" \ + --fmax "${FMAX}" \ + --reference-freq "${FREF}" \ + --lisa-reference-frequency "${FREF}" \ + --data-integration-window-half 8 \ + --srate "${SRATE}" \ + --d-min 1 \ + --d-max 5000 \ + --n-eff 2 \ + --n-max 40 \ + --n-chunk 20 \ + --save-P 1 \ + --grid-size 1 \ + --n-iterations 1 \ + --n-samples-per-job 1 \ + --request-memory-ILE 1024 \ + --request-memory-CIP 1024 + +if [[ "${RUN_ILE}" == "0" ]]; then + echo "Generated synthetic LISA analysis inputs in ${WORKDIR}" + exit 0 +fi + +"${PYTHON_BIN}" "${BIN_DIR}/integrate_likelihood_extrinsic_batchmode_lisa" \ + --LISA \ + --h5-frame-FD \ + --time-marginalization \ + --lisa-fixed-sky 1 \ + --ecliptic-longitude "${ECLIPTIC_LONGITUDE}" \ + --ecliptic-latitude "${ECLIPTIC_LATITUDE}" \ + --lisa-reference-time 0 \ + --lisa-reference-frequency "${FREF}" \ + --data-integration-window-half 8 \ + --modes "[(2,2)]" \ + --cache-file "${CACHE_FILE}" \ + --channel-name A=fake_strain \ + --channel-name E=fake_strain \ + --channel-name T=fake_strain \ + --psd-file "A=${A_PSD}" \ + --psd-file "E=${E_PSD}" \ + --psd-file "T=${T_PSD}" \ + --fmin-template "${FMIN}" \ + --fmin-ifo "A=${FMIN}" \ + --fmin-ifo "E=${FMIN}" \ + --fmin-ifo "T=${FMIN}" \ + --fmax "${FMAX}" \ + --reference-freq "${FREF}" \ + --srate "${SRATE}" \ + --l-max 2 \ + --approx IMRPhenomD \ + --mass1 "${MASS1}" \ + --mass2 "${MASS2}" \ + --spin1z "${SPIN1Z}" \ + --spin2z "${SPIN2Z}" \ + --d-max 5000 \ + --d-min 1 \ + --n-eff 2 \ + --n-max 40 \ + --n-chunk 20 \ + --save-P 1 \ + --no-adapt \ + --internal-use-lnL \ + --sampler-method AV \ + --inclination "${INCLINATION}" \ + --phi-orb "${PHIREF}" \ + --distance "${DISTANCE_MPC:-1000}" \ + --right-ascension 0 \ + --declination 0 \ + --internal-hard-fail-on-error \ + --output-file "${WORKDIR}/lisa_ile" + +test -s "${WORKDIR}/lisa_ile_0_.dat" +echo "Synthetic LISA ILE run wrote ${WORKDIR}/lisa_ile_0_.dat" diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/lisa/run_lisa_zero_likelihood_cepp.sh b/MonteCarloMarginalizeCode/Code/demo/rift/lisa/run_lisa_zero_likelihood_cepp.sh new file mode 100755 index 000000000..32d686d72 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/demo/rift/lisa/run_lisa_zero_likelihood_cepp.sh @@ -0,0 +1,53 @@ +#!/usr/bin/env bash +set -euo pipefail + +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +CODE_DIR="$(cd "${SCRIPT_DIR}/../../.." && pwd)" +BIN_DIR="${CODE_DIR}/bin" + +PYTHON_BIN="${RIFT_LISA_PYTHON:-}" +if [[ -z "${PYTHON_BIN}" ]]; then + PYTHON_BIN="$(command -v python3)" +fi + +WORKDIR="${RIFT_LISA_WORKDIR:-/tmp/rift-lisa-zero-likelihood-$(date +%s)}" +RENDER_CEPP="${RIFT_LISA_RENDER_CEPP:-1}" + +mkdir -p "${WORKDIR}" + +export RIFT_HYPERPIPELINE_FORMAT=1 +export PYTHONPATH="${CODE_DIR}${PYTHONPATH:+:${PYTHONPATH}}" +export PATH="${BIN_DIR}:${PATH}" + +"${PYTHON_BIN}" "${BIN_DIR}/helper_LISA_Events.py" \ + --working-directory "${WORKDIR}" \ + --zero-likelihood \ + --grid-size 1 \ + --n-iterations 1 \ + --n-samples-per-job 1 \ + --request-memory-ILE 1024 \ + --request-memory-CIP 1024 \ + "$@" + +echo "Generated LISA helper bundle in ${WORKDIR}" + +if [[ "${RENDER_CEPP}" == "0" ]]; then + exit 0 +fi + +"${PYTHON_BIN}" "${BIN_DIR}/create_event_parameter_pipeline_BasicIteration" \ + --ile-n-events-to-analyze 1 \ + --input-grid "${WORKDIR}/proposed-grid.dat" \ + --ile-exe "${BIN_DIR}/integrate_likelihood_extrinsic_batchmode_lisa" \ + --ile-args "${WORKDIR}/args_ile.txt" \ + --cip-args-list "${WORKDIR}/args_cip_list.txt" \ + --test-args "${WORKDIR}/args_test.txt" \ + --working-directory "${WORKDIR}" \ + --n-iterations 1 \ + --n-samples-per-job 1 \ + --n-copies 1 \ + --request-memory-ILE 1024 \ + --request-memory-CIP 1024 \ + --transfer-file-list "${WORKDIR}/helper_transfer_files.txt" + +echo "Rendered LISA zero-likelihood CEPP DAG in ${WORKDIR}" diff --git a/MonteCarloMarginalizeCode/Code/test/pp_lisa/README.md b/MonteCarloMarginalizeCode/Code/test/pp_lisa/README.md new file mode 100644 index 000000000..0f093f7ed --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/pp_lisa/README.md @@ -0,0 +1,29 @@ +# LISA PP Test Surface + +This directory is the LISA analogue of `test/pp`: it stages synthetic +injection-like data products and renders a small RIFT analysis surface from a +known truth. It is intentionally lightweight at first; full PP population +drivers can grow here without putting injection-generation workflows in the +main package path. + +The current smoke path builds one known-sky event: + +```bash +./run_pp_lisa_known_sky.sh +``` + +The driver writes, under `RIFT_PP_LISA_WORKDIR` or a temporary directory: + +- A/E/T frequency-domain HDF5 frame products +- `lisa.cache` +- analytic LISA A/E/T XML PSDs +- `synthetic-params.env` +- a `pseudo_pipe` known-sky run directory with hyperpipeline CEPP files + +Set `RIFT_PP_LISA_RUN_ILE=1` to also run the tiny direct ILE check after the +bundle and DAG are rendered. + +Set `RIFT_PP_LISA_VARY_SKY=1` to render the same surface with ecliptic sky +location left as an intrinsic hyperpipeline parameter. In that mode the helper +omits `--lisa-fixed-sky` from `args_ile.txt` and seeds a small three-point sky +grid. diff --git a/MonteCarloMarginalizeCode/Code/test/pp_lisa/run_pp_lisa_known_sky.sh b/MonteCarloMarginalizeCode/Code/test/pp_lisa/run_pp_lisa_known_sky.sh new file mode 100755 index 000000000..cd0726517 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/pp_lisa/run_pp_lisa_known_sky.sh @@ -0,0 +1,79 @@ +#!/usr/bin/env bash +set -euo pipefail + +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +CODE_DIR="$(cd "${SCRIPT_DIR}/../.." && pwd)" +DEMO_DIR="${CODE_DIR}/demo/rift/lisa" +BIN_DIR="${CODE_DIR}/bin" + +PYTHON_BIN="${RIFT_LISA_PYTHON:-}" +if [[ -z "${PYTHON_BIN}" ]]; then + PYTHON_BIN="$(command -v python3)" +fi + +WORKDIR="${RIFT_PP_LISA_WORKDIR:-/tmp/rift-pp-lisa-known-sky-$(date +%s)}" +BUNDLE_DIR="${WORKDIR}/event_0" +RUNDIR="${WORKDIR}/analysis_event_0" +RUN_ILE="${RIFT_PP_LISA_RUN_ILE:-0}" +VARY_SKY="${RIFT_PP_LISA_VARY_SKY:-0}" +SKY_ARGS=() +if [[ "${VARY_SKY}" == "1" ]]; then + SKY_ARGS=(--lisa-vary-sky --lisa-grid-size 3 --lisa-sky-grid-width 0.01) +else + SKY_ARGS=(--lisa-grid-size 1) +fi + +mkdir -p "${BUNDLE_DIR}" + +export PYTHONPATH="${CODE_DIR}${PYTHONPATH:+:${PYTHONPATH}}" +export PATH="${BIN_DIR}:${PATH}" +export MPLCONFIGDIR="${WORKDIR}/.matplotlib" + +"${PYTHON_BIN}" "${DEMO_DIR}/make_synthetic_lisa_inputs.py" \ + --output-directory "${BUNDLE_DIR}" \ + --duration 1024 \ + --deltaT 4 \ + "$@" + +"${PYTHON_BIN}" "${DEMO_DIR}/make_lisa_psds.py" \ + --output-directory "${BUNDLE_DIR}" \ + --fmax 0.125 \ + --npts 513 \ + --write-ascii + +set -a +source "${BUNDLE_DIR}/synthetic-params.env" +set +a + +"${PYTHON_BIN}" "${BIN_DIR}/util_RIFT_pseudo_pipe.py" \ + --lisa-known-sky \ + --use-rundir "${RUNDIR}" \ + --approx IMRPhenomD \ + --event-time 0 \ + --ecliptic-longitude "${ECLIPTIC_LONGITUDE}" \ + --ecliptic-latitude "${ECLIPTIC_LATITUDE}" \ + --lisa-cache-file "${CACHE_FILE}" \ + --lisa-channel-name A=fake_strain \ + --lisa-channel-name E=fake_strain \ + --lisa-channel-name T=fake_strain \ + --lisa-psd-file "A=${BUNDLE_DIR}/A_psd.xml.gz" \ + --lisa-psd-file "E=${BUNDLE_DIR}/E_psd.xml.gz" \ + --lisa-psd-file "T=${BUNDLE_DIR}/T_psd.xml.gz" \ + --lisa-srate "${SRATE}" \ + --lisa-fmin-template "${FMIN}" \ + --lisa-fmax "${FMAX}" \ + --lisa-reference-freq "${FREF}" \ + "${SKY_ARGS[@]}" \ + --lisa-n-iterations 1 \ + --lisa-n-samples-per-job 1 \ + --internal-ile-request-memory 1024 \ + --internal-cip-request-memory 1024 + +if [[ "${RUN_ILE}" == "1" ]]; then + RIFT_LISA_WORKDIR="${BUNDLE_DIR}" "${DEMO_DIR}/run_lisa_synthetic_ile.sh" \ + --duration 1024 \ + --deltaT 4 +fi + +echo "LISA PP known-sky bundle: ${BUNDLE_DIR}" +echo "LISA PP known-sky run: ${RUNDIR}" diff --git a/MonteCarloMarginalizeCode/Code/test/test_lisa_auxiliary_imports.py b/MonteCarloMarginalizeCode/Code/test/test_lisa_auxiliary_imports.py new file mode 100644 index 000000000..f72d90a07 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lisa_auxiliary_imports.py @@ -0,0 +1,87 @@ +"""Smoke and contract tests for imported LISA auxiliary modules.""" + +import os + +import numpy as np +import pytest + + +def test_lisa_auxiliary_modules_import_without_cli_side_effects(): + import RIFT.LISA.initial_grid.fisher_errors as fisher_errors + import RIFT.LISA.injections.LISA_injections as lisa_injections + import RIFT.LISA.injections.create_injections as create_injections + import RIFT.LISA.psd_generation.generate_LISA_psd as generate_LISA_psd + import RIFT.LISA.sangria_test.pycbc_to_rift as pycbc_to_rift + import RIFT.LISA.utils.utils as lisa_utils + + assert hasattr(fisher_errors, "get_error_bounds") + assert hasattr(lisa_injections, "generate_lisa_TDI_dict") + assert hasattr(create_injections, "parameter_dict_from_xml") + assert hasattr(generate_LISA_psd, "write_lisa_psd") + assert hasattr(pycbc_to_rift, "create_injection_from_time_series") + assert hasattr(pycbc_to_rift, "read_gwpy_channels") + assert hasattr(lisa_utils, "SSB_to_LISA") + + +def test_lisa_sky_frame_round_trip(): + from RIFT.LISA.utils import utils as lisa_utils + + t_ssb = np.array([1024.0]) + lam = np.array([1.2]) + beta = np.array([0.3]) + psi = np.array([0.4]) + + lisa_frame = lisa_utils.SSB_to_LISA(t_ssb, lam, beta, psi) + ssb_frame = lisa_utils.LISA_to_SSB( + lisa_frame[:, 0], + lisa_frame[:, 1], + lisa_frame[:, 2], + lisa_frame[:, 3], + ) + + assert np.allclose(ssb_frame[:, 0], t_ssb, rtol=0, atol=1e-4) + assert np.allclose(ssb_frame[:, 1], lam, rtol=0, atol=1e-10) + assert np.allclose(ssb_frame[:, 2], beta, rtol=0, atol=1e-10) + + +def test_lisa_psd_generator_writes_small_ascii_products(tmp_path): + from RIFT.LISA.psd_generation import generate_LISA_psd + + txt_path, png_path = generate_LISA_psd.write_lisa_psd( + os.fspath(tmp_path), + fmin=1.0e-4, + fmax=1.0e-3, + npts=16, + write_xml=False, + ) + + psd = np.loadtxt(txt_path) + assert psd.shape == (16, 2) + assert np.all(np.isfinite(psd)) + assert np.all(psd[:, 1] > 0) + assert os.path.exists(png_path) + + +def test_sangria_converter_prefers_gwpy_reader(monkeypatch): + from RIFT.LISA.sangria_test import pycbc_to_rift + + calls = [] + + def fake_read(frame_path, channel): + calls.append((frame_path, channel)) + return channel + + monkeypatch.setattr(pycbc_to_rift, "_read_gwpy_frame", fake_read) + assert pycbc_to_rift.read_gwpy_channels("frame.gwf", channels=("A", "E")) == {"A": "A", "E": "E"} + assert calls == [("frame.gwf", "A"), ("frame.gwf", "E")] + + +def test_sangria_converter_reports_missing_optional_gwpy(monkeypatch): + from RIFT.LISA.sangria_test import pycbc_to_rift + + def missing_gwpy(): + raise ImportError("gwpy is required to read Sangria frame files") + + monkeypatch.setattr(pycbc_to_rift, "_import_gwpy_timeseries", missing_gwpy) + with pytest.raises(ImportError, match="gwpy is required"): + pycbc_to_rift.read_gwpy_channels("missing.gwf") diff --git a/MonteCarloMarginalizeCode/Code/test/test_lisa_demo_contract.py b/MonteCarloMarginalizeCode/Code/test/test_lisa_demo_contract.py new file mode 100644 index 000000000..a8f956eee --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lisa_demo_contract.py @@ -0,0 +1,53 @@ +"""Contract test for the checked-in LISA zero-likelihood demo.""" + +import os +import subprocess + +from RIFT.misc import hyperpipeline_io + + +REPO_ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", "..")) +DEMO = os.path.join( + REPO_ROOT, + "MonteCarloMarginalizeCode", + "Code", + "demo", + "rift", + "lisa", + "run_lisa_zero_likelihood_cepp.sh", +) + + +def test_lisa_zero_likelihood_demo_renders_cepp_bundle(tmp_path): + env = os.environ.copy() + env["RIFT_LISA_WORKDIR"] = os.fspath(tmp_path) + subprocess.run([DEMO], check=True, env=env) + + expected = { + "proposed-grid.dat", + "args_ile.txt", + "args_cip_list.txt", + "args_test.txt", + "helper_transfer_files.txt", + "command-cepp-lisa.sh", + "ILE.sub", + "CIP.sub", + } + assert expected <= {path.name for path in tmp_path.iterdir()} + assert any(path.suffix == ".dag" for path in tmp_path.iterdir()) + + grid, columns = hyperpipeline_io.read_table(os.fspath(tmp_path / "proposed-grid.dat")) + assert grid.shape == (1,) + assert "ecliptic_longitude" in columns + assert "ecliptic_latitude" in columns + + ile_args = (tmp_path / "args_ile.txt").read_text() + assert "--zero-likelihood" in ile_args + assert "--LISA" in ile_args + assert "--cache-file lisa.cache" in ile_args + + cip_args = (tmp_path / "args_cip_list.txt").read_text() + assert "--parameter mc" in cip_args + assert "--parameter eta" in cip_args + assert "--parameter ecliptic_longitude" not in cip_args + assert "--parameter ecliptic_latitude" not in cip_args diff --git a/MonteCarloMarginalizeCode/Code/test/test_lisa_helper_contract.py b/MonteCarloMarginalizeCode/Code/test/test_lisa_helper_contract.py index 9e99ea097..40ca10cc3 100644 --- a/MonteCarloMarginalizeCode/Code/test/test_lisa_helper_contract.py +++ b/MonteCarloMarginalizeCode/Code/test/test_lisa_helper_contract.py @@ -49,7 +49,7 @@ def test_lisa_helper_writes_cepp_contract_files(tmp_path): assert expected <= {path.name for path in tmp_path.iterdir()} grid, columns = hyperpipeline_io.read_table(os.fspath(tmp_path / "proposed-grid.dat")) - assert grid.shape == (3,) + assert grid.shape[0] >= 3 assert "ecliptic_longitude" in columns assert "ecliptic_latitude" in columns @@ -57,7 +57,9 @@ def test_lisa_helper_writes_cepp_contract_files(tmp_path): assert ile_args.startswith("X ") assert "--LISA" in ile_args assert "--h5-frame-FD" in ile_args + assert "--time-marginalization" in ile_args assert "--zero-likelihood" in ile_args + assert "--data-integration-window-half 300.0" in ile_args assert "--cache-file lisa.cache" in ile_args assert "--channel-name A=fake_strain" in ile_args assert "--psd-file A=A_psd.xml.gz" in ile_args @@ -66,8 +68,10 @@ def test_lisa_helper_writes_cepp_contract_files(tmp_path): cip_args = (tmp_path / "args_cip_list.txt").read_text() assert cip_args.startswith("1 ") - assert "--parameter ecliptic_longitude" in cip_args - assert "--parameter ecliptic_latitude" in cip_args + assert "--parameter mc" in cip_args + assert "--parameter eta" in cip_args + assert "--parameter ecliptic_longitude" not in cip_args + assert "--parameter ecliptic_latitude" not in cip_args assert "--fname" not in cip_args test_args = (tmp_path / "args_test.txt").read_text() @@ -80,6 +84,70 @@ def test_lisa_helper_writes_cepp_contract_files(tmp_path): assert os.fspath(tmp_path / "proposed-grid.dat") in cepp_command +def test_lisa_helper_custom_data_products_replace_defaults(tmp_path): + _run_helper( + tmp_path, + "--cache-file", + os.fspath(tmp_path / "custom.cache"), + "--psd-file", + "A=/tmp/A.xml.gz", + "--psd-file", + "E=/tmp/E.xml.gz", + "--psd-file", + "T=/tmp/T.xml.gz", + "--channel-name", + "A=SYNTH", + "--channel-name", + "E=SYNTH", + "--channel-name", + "T=SYNTH", + ) + + ile_args = (tmp_path / "args_ile.txt").read_text() + assert "--cache-file {}/custom.cache".format(tmp_path) in ile_args + assert "--psd-file A=/tmp/A.xml.gz" in ile_args + assert "--psd-file A=A_psd.xml.gz" not in ile_args + assert "--channel-name A=SYNTH" in ile_args + assert "--channel-name A=fake_strain" not in ile_args + + transfer_files = (tmp_path / "helper_transfer_files.txt").read_text().splitlines() + assert os.fspath(tmp_path / "custom.cache") in transfer_files + assert "/tmp/A.xml.gz" in transfer_files + assert "A_psd.xml.gz" not in transfer_files + + +def test_lisa_helper_variable_sky_uses_grid_sky(tmp_path): + _run_helper( + tmp_path, + "--vary-sky", + "--grid-size", + "3", + "--sky-grid-width", + "0.01", + ) + + grid, _ = hyperpipeline_io.read_table(os.fspath(tmp_path / "proposed-grid.dat")) + assert grid.shape[0] >= 3 + assert len(set(grid["ecliptic_longitude"])) > 1 + assert len(set(grid["ecliptic_latitude"])) > 1 + # off-lattice sky jitter within +/- sky_grid_width (see pp_surface test for + # the collinearity rationale); assert it varies and stays in the width box. + assert all(abs(x - 1.0) <= 0.01 + 1e-9 for x in grid["ecliptic_longitude"]) + assert all(abs(x - 0.3) <= 0.01 + 1e-9 for x in grid["ecliptic_latitude"]) + + ile_args = (tmp_path / "args_ile.txt").read_text() + # vary-sky: ILE uses the per-row grid sky (--lisa-fixed-sky 1), no hardcode. + assert "--lisa-fixed-sky" in ile_args + assert "--ecliptic-longitude" not in ile_args + assert "--ecliptic-latitude" not in ile_args + + cip_args = (tmp_path / "args_cip_list.txt").read_text() + # CIP fits sky as phi/theta; hyperpipeline aliases the ecliptic_longitude/ + # latitude NAMED columns into P.phi/P.theta on read (no positional all.net). + assert "--parameter phi" in cip_args + assert "--parameter theta" in cip_args + + def test_lisa_helper_bundle_renders_basic_cepp_dag(tmp_path): _run_helper( tmp_path, diff --git a/MonteCarloMarginalizeCode/Code/test/test_lisa_ini_contract.py b/MonteCarloMarginalizeCode/Code/test/test_lisa_ini_contract.py new file mode 100644 index 000000000..7cd5dfe63 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lisa_ini_contract.py @@ -0,0 +1,142 @@ +#!/usr/bin/env python +"""Contract test for the LISA production-ini path (util_RIFT_pseudo_pipe --use-ini). + +The LISA known-sky workflow can be driven either from --lisa-* CLI flags or from +an .ini file (production form). This test renders BOTH ways with matched config +and asserts the generated workflow is identical -- so the ini path stays a thin +front-end over the validated CLI machinery and never silently diverges. + +The ini sources per-channel data products from the conventional sections +([data] channels, [lalinference] psds); everything else flows through the generic +[rift-pseudo-pipe] parser by CLI-arg name. +""" + +import os +import subprocess +import sys + +import pytest + +REPO_ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", "..")) +CODE_DIR = os.path.join(REPO_ROOT, "MonteCarloMarginalizeCode", "Code") +PSEUDO_PIPE = os.path.join(CODE_DIR, "bin", "util_RIFT_pseudo_pipe.py") + + +def _env(): + env = os.environ.copy() + env["PYTHONPATH"] = CODE_DIR + os.pathsep + env.get("PYTHONPATH", "") + env["PATH"] = os.path.join(CODE_DIR, "bin") + os.pathsep + env.get("PATH", "") + return env + + +# matched (CLI flag, ini key, value) for the scalar/algorithm options +_COMMON = dict( + approx="IMRPhenomD", + event_time="0", + ecliptic_longitude="1.0", + ecliptic_latitude="0.3", + srate="0.25", + fmin="0.0001", + fmax="0.125", + fref="0.005", + window="300", + grid_size="16", + n_iter="4", + n_samp="16", +) + + +def _write_ini(path, cache, psd): + path.write_text( + "[data]\n" + "channels = {{'A': 'fake_strain', 'E': 'fake_strain', 'T': 'fake_strain'}}\n" + "\n[lalinference]\n" + "psds = {{'A': '{A}', 'E': '{E}', 'T': '{T}'}}\n" + "\n[rift-pseudo-pipe]\n" + "lisa-known-sky=True\n" + "lisa-vary-sky=True\n" + "lisa-search-reflected-sky-mode=True\n" + "lisa-reference-time=0\n" + 'approx="{approx}"\n' + "event-time={event_time}\n" + "ecliptic-longitude={ecliptic_longitude}\n" + "ecliptic-latitude={ecliptic_latitude}\n" + 'lisa-cache-file="{cache}"\n' + "lisa-srate={srate}\n" + "lisa-fmin-template={fmin}\n" + "lisa-fmax={fmax}\n" + "lisa-reference-freq={fref}\n" + "lisa-data-integration-window-half={window}\n" + "lisa-grid-size={grid_size}\n" + "lisa-n-iterations={n_iter}\n" + "lisa-n-samples-per-job={n_samp}\n" + "internal-ile-request-memory=2048\n" + "internal-cip-request-memory=2048\n".format( + A=psd["A"], E=psd["E"], T=psd["T"], cache=cache, **_COMMON + ) + ) + + +def _render_cli(rundir, cache, psd): + cmd = [ + sys.executable, PSEUDO_PIPE, + "--lisa-known-sky", "--lisa-vary-sky", + "--lisa-search-reflected-sky-mode", "--lisa-reference-time", "0", + "--use-rundir", os.fspath(rundir), + "--approx", _COMMON["approx"], + "--event-time", _COMMON["event_time"], + "--ecliptic-longitude", _COMMON["ecliptic_longitude"], + "--ecliptic-latitude", _COMMON["ecliptic_latitude"], + "--lisa-cache-file", cache, + "--lisa-channel-name", "A=fake_strain", + "--lisa-channel-name", "E=fake_strain", + "--lisa-channel-name", "T=fake_strain", + "--lisa-psd-file", "A={}".format(psd["A"]), + "--lisa-psd-file", "E={}".format(psd["E"]), + "--lisa-psd-file", "T={}".format(psd["T"]), + "--lisa-srate", _COMMON["srate"], + "--lisa-fmin-template", _COMMON["fmin"], + "--lisa-fmax", _COMMON["fmax"], + "--lisa-reference-freq", _COMMON["fref"], + "--lisa-data-integration-window-half", _COMMON["window"], + "--lisa-grid-size", _COMMON["grid_size"], + "--lisa-n-iterations", _COMMON["n_iter"], + "--lisa-n-samples-per-job", _COMMON["n_samp"], + "--internal-ile-request-memory", "2048", + "--internal-cip-request-memory", "2048", + ] + subprocess.run(cmd, check=True, env=_env()) + + +def test_lisa_ini_path_matches_cli(tmp_path): + cache = os.fspath(tmp_path / "lisa.cache") + psd = {c: os.fspath(tmp_path / "{}_psd.xml.gz".format(c)) for c in ("A", "E", "T")} + + cli_dir = tmp_path / "cli" + ini_dir = tmp_path / "ini" + ini_file = tmp_path / "demo.ini" + _write_ini(ini_file, cache, psd) + + _render_cli(cli_dir, cache, psd) + subprocess.run( + [sys.executable, PSEUDO_PIPE, "--use-ini", os.fspath(ini_file), + "--use-rundir", os.fspath(ini_dir)], + check=True, env=_env(), + ) + + def norm(rundir, name): + return (rundir / name).read_text().replace(os.fspath(rundir), "RUNDIR") + + # the rendered ILE / CIP / test argument files must be byte-identical + for name in ("args_ile.txt", "args_cip_list.txt", "args_test.txt"): + assert norm(cli_dir, name) == norm(ini_dir, name), name + + # the ini path must produce the same workflow (incl. the reflect node) + cli_dag = next(cli_dir.glob("*.dag")).read_text() + ini_dag = next(ini_dir.glob("*.dag")).read_text() + assert cli_dag.count("\nJOB ") == ini_dag.count("\nJOB ") + assert "convert_primary_sky_mode_to_secondary" in ini_dag + # data products sourced from [data]/[lalinference] reach the ILE args + ile_args = (ini_dir / "args_ile.txt").read_text() + assert "--channel-name A=fake_strain" in ile_args + assert psd["A"] in ile_args diff --git a/MonteCarloMarginalizeCode/Code/test/test_lisa_pp_surface.py b/MonteCarloMarginalizeCode/Code/test/test_lisa_pp_surface.py new file mode 100644 index 000000000..170305e82 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lisa_pp_surface.py @@ -0,0 +1,120 @@ +"""Tests for the lightweight LISA PP-style surface.""" + +import os +import subprocess + +import pytest + +import RIFT.lalsimutils as lalsimutils + +from RIFT.misc import hyperpipeline_io + + +REPO_ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", "..")) +CODE_DIR = os.path.join(REPO_ROOT, "MonteCarloMarginalizeCode", "Code") +PP_LISA_DRIVER = os.path.join(CODE_DIR, "test", "pp_lisa", "run_pp_lisa_known_sky.sh") +MAKE_PSDS = os.path.join(CODE_DIR, "demo", "rift", "lisa", "make_lisa_psds.py") + + +def test_lisa_demo_psd_generator_writes_channel_xml(tmp_path): + subprocess.run( + [ + MAKE_PSDS, + "--output-directory", + os.fspath(tmp_path), + "--fmax", + "0.125", + "--npts", + "129", + "--write-ascii", + ], + check=True, + ) + + for channel in ["A", "E", "T"]: + psd_path = tmp_path / "{}_psd.xml.gz".format(channel) + assert psd_path.exists() + psd = lalsimutils.get_psd_series_from_xmldoc(os.fspath(psd_path), channel) + assert psd.data.length == 129 + assert psd.deltaF > 0 + assert psd.data.data[1] > 0 + + assert (tmp_path / "LISA_psd.txt").exists() + + +def test_lisa_pp_known_sky_surface_builds_bundle_and_dag(tmp_path): + env = os.environ.copy() + env["RIFT_PP_LISA_WORKDIR"] = os.fspath(tmp_path) + env["RIFT_PP_LISA_RUN_ILE"] = "0" + env["PYTHONPATH"] = CODE_DIR + os.pathsep + env.get("PYTHONPATH", "") + env["PATH"] = os.path.join(CODE_DIR, "bin") + os.pathsep + env.get("PATH", "") + + subprocess.run([PP_LISA_DRIVER], check=True, env=env) + + bundle_dir = tmp_path / "event_0" + rundir = tmp_path / "analysis_event_0" + + expected_bundle = { + "A-fake_strain-1000000-10000.h5", + "E-fake_strain-1000000-10000.h5", + "T-fake_strain-1000000-10000.h5", + "A_psd.xml.gz", + "E_psd.xml.gz", + "T_psd.xml.gz", + "LISA_psd.txt", + "lisa.cache", + "synthetic-params.env", + } + assert expected_bundle <= {path.name for path in bundle_dir.iterdir()} + + assert (rundir / "proposed-grid.dat").exists() + assert (rundir / "args_ile.txt").exists() + assert (rundir / "helper_transfer_files.txt").exists() + assert (rundir / "ILE.sub").exists() + assert (rundir / "CIP.sub").exists() + + _, columns = hyperpipeline_io.read_table(os.fspath(rundir / "proposed-grid.dat")) + assert "ecliptic_longitude" in columns + assert "ecliptic_latitude" in columns + + ile_args = (rundir / "args_ile.txt").read_text() + assert "--lisa-fixed-sky 1" in ile_args + assert "--cache-file {}".format(bundle_dir / "lisa.cache") in ile_args + assert "--psd-file A={}".format(bundle_dir / "A_psd.xml.gz") in ile_args + assert "--srate 0.25" in ile_args + + transfer_files = (rundir / "helper_transfer_files.txt").read_text().splitlines() + assert os.fspath(bundle_dir / "lisa.cache") in transfer_files + assert os.fspath(bundle_dir / "A_psd.xml.gz") in transfer_files + + +def test_lisa_pp_variable_sky_surface_builds_intrinsic_sky_grid(tmp_path): + env = os.environ.copy() + env["RIFT_PP_LISA_WORKDIR"] = os.fspath(tmp_path) + env["RIFT_PP_LISA_RUN_ILE"] = "0" + env["RIFT_PP_LISA_VARY_SKY"] = "1" + env["PYTHONPATH"] = CODE_DIR + os.pathsep + env.get("PYTHONPATH", "") + env["PATH"] = os.path.join(CODE_DIR, "bin") + os.pathsep + env.get("PATH", "") + + subprocess.run([PP_LISA_DRIVER], check=True, env=env) + + rundir = tmp_path / "analysis_event_0" + grid, _ = hyperpipeline_io.read_table(os.fspath(rundir / "proposed-grid.dat")) + assert grid.shape[0] >= 3 + assert len(set(grid["ecliptic_longitude"])) > 1 + assert len(set(grid["ecliptic_latitude"])) > 1 + # vary-sky jitters each point's sky OFF-LATTICE within +/- sky_grid_width so + # (mc, eta, phi, theta) span 4 independent dims -- a sky tied to the mass + # index would be collinear with mass1/mass2 -> singular CIP fit. Assert the + # sky varies and stays inside the +/- width box, not at exact lattice ends. + assert 0 < max(grid["ecliptic_longitude"]) - min(grid["ecliptic_longitude"]) <= 0.02 + 1e-9 + assert 0 < max(grid["ecliptic_latitude"]) - min(grid["ecliptic_latitude"]) <= 0.02 + 1e-9 + assert all(abs(x - 1.0) <= 0.01 + 1e-9 for x in grid["ecliptic_longitude"]) + assert all(abs(x - 0.3) <= 0.01 + 1e-9 for x in grid["ecliptic_latitude"]) + + ile_args = (rundir / "args_ile.txt").read_text() + # vary-sky: ILE evaluates each grid point at its OWN fixed sky (--lisa-fixed-sky + # 1, supplied per-row via --sim-grid); no single hardcoded ecliptic sky. + assert "--lisa-fixed-sky" in ile_args + assert "--ecliptic-longitude" not in ile_args + assert "--ecliptic-latitude" not in ile_args diff --git a/MonteCarloMarginalizeCode/Code/test/test_lisa_pseudo_pipe_contract.py b/MonteCarloMarginalizeCode/Code/test/test_lisa_pseudo_pipe_contract.py new file mode 100644 index 000000000..e07c73ff1 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lisa_pseudo_pipe_contract.py @@ -0,0 +1,166 @@ +#!/usr/bin/env python +"""Contract tests for the LISA known-sky pseudo_pipe surface.""" + +import os +import subprocess +import sys + +import pytest + +from RIFT.misc import hyperpipeline_io + + +REPO_ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", "..")) +CODE_DIR = os.path.join(REPO_ROOT, "MonteCarloMarginalizeCode", "Code") +PSEUDO_PIPE = os.path.join(CODE_DIR, "bin", "util_RIFT_pseudo_pipe.py") + + +def test_lisa_known_sky_pseudo_pipe_renders_cepp_surface(tmp_path): + rundir = tmp_path / "pseudo_lisa" + env = os.environ.copy() + env["PYTHONPATH"] = CODE_DIR + os.pathsep + env.get("PYTHONPATH", "") + env["PATH"] = os.path.join(CODE_DIR, "bin") + os.pathsep + env.get("PATH", "") + + cmd = [ + sys.executable, + PSEUDO_PIPE, + "--lisa-known-sky", + "--use-rundir", + os.fspath(rundir), + "--approx", + "IMRPhenomD", + "--event-time", + "1234.5", + "--ecliptic-longitude", + "1.25", + "--ecliptic-latitude", + "-0.4", + "--lisa-cache-file", + os.fspath(tmp_path / "lisa.cache"), + "--lisa-channel-name", + "A=SYNTH", + "--lisa-channel-name", + "E=SYNTH", + "--lisa-channel-name", + "T=SYNTH", + "--lisa-psd-file", + "A={}".format(tmp_path / "A_psd.xml.gz"), + "--lisa-psd-file", + "E={}".format(tmp_path / "E_psd.xml.gz"), + "--lisa-psd-file", + "T={}".format(tmp_path / "T_psd.xml.gz"), + "--lisa-srate", + "0.25", + "--lisa-fmin-template", + "0.001", + "--lisa-fmax", + "0.125", + "--lisa-grid-size", + "1", + "--lisa-n-iterations", + "1", + "--lisa-n-samples-per-job", + "1", + "--internal-ile-request-memory", + "1024", + "--internal-cip-request-memory", + "1024", + ] + subprocess.run(cmd, check=True, env=env) + + assert (rundir / "args_ile.txt").exists() + assert (rundir / "args_cip_list.txt").exists() + assert (rundir / "helper_transfer_files.txt").exists() + assert (rundir / "ILE.sub").exists() + assert (rundir / "CIP.sub").exists() + + grid, columns = hyperpipeline_io.read_table(os.fspath(rundir / "proposed-grid.dat")) + assert grid.shape == (1,) + assert "ecliptic_longitude" in columns + assert "ecliptic_latitude" in columns + assert grid["ecliptic_longitude"][0] == 1.25 + assert grid["ecliptic_latitude"][0] == -0.4 + + ile_args = (rundir / "args_ile.txt").read_text() + assert "--LISA" in ile_args + assert "--lisa-fixed-sky 1" in ile_args + assert "--ecliptic-longitude 1.25" in ile_args + assert "--ecliptic-latitude -0.4" in ile_args + assert "--srate 0.25" in ile_args + assert "--channel-name A=SYNTH" in ile_args + assert "A=fake_strain" not in ile_args + assert "--psd-file A={}".format(tmp_path / "A_psd.xml.gz") in ile_args + assert "A=A_psd.xml.gz" not in ile_args + + transfer_files = (rundir / "helper_transfer_files.txt").read_text().splitlines() + assert os.fspath(tmp_path / "lisa.cache") in transfer_files + assert os.fspath(tmp_path / "A_psd.xml.gz") in transfer_files + + +def test_lisa_variable_sky_pseudo_pipe_leaves_sky_intrinsic(tmp_path): + rundir = tmp_path / "pseudo_lisa_variable_sky" + env = os.environ.copy() + env["PYTHONPATH"] = CODE_DIR + os.pathsep + env.get("PYTHONPATH", "") + env["PATH"] = os.path.join(CODE_DIR, "bin") + os.pathsep + env.get("PATH", "") + + cmd = [ + sys.executable, + PSEUDO_PIPE, + "--lisa-known-sky", + "--lisa-vary-sky", + "--use-rundir", + os.fspath(rundir), + "--approx", + "IMRPhenomD", + "--event-time", + "1234.5", + "--ecliptic-longitude", + "1.25", + "--ecliptic-latitude", + "-0.4", + "--lisa-cache-file", + os.fspath(tmp_path / "lisa.cache"), + "--lisa-psd-file", + "A={}".format(tmp_path / "A_psd.xml.gz"), + "--lisa-psd-file", + "E={}".format(tmp_path / "E_psd.xml.gz"), + "--lisa-psd-file", + "T={}".format(tmp_path / "T_psd.xml.gz"), + "--lisa-srate", + "0.25", + "--lisa-grid-size", + "3", + "--lisa-sky-grid-width", + "0.01", + "--lisa-n-iterations", + "1", + "--lisa-n-samples-per-job", + "1", + "--internal-ile-request-memory", + "1024", + "--internal-cip-request-memory", + "1024", + ] + subprocess.run(cmd, check=True, env=env) + + grid, _ = hyperpipeline_io.read_table(os.fspath(rundir / "proposed-grid.dat")) + assert grid.shape[0] >= 3 + assert len(set(grid["ecliptic_longitude"])) > 1 + assert len(set(grid["ecliptic_latitude"])) > 1 + # off-lattice sky jitter within +/- sky_grid_width (see pp_surface test for + # the collinearity rationale); assert it varies and stays in the width box. + assert all(abs(x - 1.25) <= 0.01 + 1e-9 for x in grid["ecliptic_longitude"]) + assert all(abs(x + 0.40) <= 0.01 + 1e-9 for x in grid["ecliptic_latitude"]) + + ile_args = (rundir / "args_ile.txt").read_text() + assert "--LISA" in ile_args + # vary-sky: ILE uses the per-row grid sky (--lisa-fixed-sky 1), no hardcode. + assert "--lisa-fixed-sky" in ile_args + assert "--ecliptic-longitude" not in ile_args + assert "--ecliptic-latitude" not in ile_args + + cip_args = (rundir / "args_cip_list.txt").read_text() + # CIP fits sky as phi/theta; hyperpipeline aliases the ecliptic_longitude/ + # latitude NAMED columns into P.phi/P.theta on read (no positional all.net). + assert "--parameter phi" in cip_args + assert "--parameter theta" in cip_args diff --git a/MonteCarloMarginalizeCode/Code/test/test_lisa_run_checks.py b/MonteCarloMarginalizeCode/Code/test/test_lisa_run_checks.py new file mode 100644 index 000000000..44f518295 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lisa_run_checks.py @@ -0,0 +1,60 @@ +"""Tests for LISA run diagnostics.""" + +import json +import os +import subprocess +import sys + +import numpy as np + +from RIFT.LISA.run_checks import plot_RIFT + + +REPO_ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", "..")) +CHECK_SCRIPT = os.path.join( + REPO_ROOT, + "MonteCarloMarginalizeCode", + "Code", + "RIFT", + "LISA", + "run_checks", + "plot_RIFT.py", +) + + +def _write_lisa_table(path): + rows = np.array( + [ + [0, 1.0e5, 8.0e4, 0, 0, 0.1, 0, 0, -0.1, 1.0, 0.3, 10.0, 0.2, 40, 12], + [0, 1.1e5, 7.5e4, 0, 0, 0.0, 0, 0, 0.0, 1.1, 0.2, 14.0, 0.5, 50, 20], + [0, 0.9e5, 8.5e4, 0, 0, 0.2, 0, 0, 0.1, 0.9, 0.4, 13.0, 0.1, 60, 30], + ] + ) + np.savetxt(path, rows) + + +def test_lisa_run_summary_identifies_best_point(tmp_path): + table = tmp_path / "all.net" + _write_lisa_table(table) + + summary = plot_RIFT.summarize_lisa_ile(os.fspath(table), lnL_window=2.0, error_threshold=0.4) + assert summary["n_rows"] == 3 + assert summary["max_index"] == 1 + assert summary["max_lnL"] == 14.0 + assert summary["high_lnL_points"] == 2 + assert summary["high_lnL_low_error_points"] == 1 + assert summary["best"]["ecliptic_longitude"] == 1.1 + assert summary["best"]["ecliptic_latitude"] == 0.2 + + +def test_lisa_run_summary_cli_json(tmp_path): + table = tmp_path / "lisa_ile_0_.dat" + _write_lisa_table(table) + + output = subprocess.check_output( + [sys.executable, CHECK_SCRIPT, os.fspath(table), "--json"], + text=True, + ) + summary = json.loads(output) + assert summary["max_lnL"] == 14.0 + assert summary["best"]["n_eff"] == 20.0 diff --git a/MonteCarloMarginalizeCode/Code/test/test_lisa_synthetic_demo.py b/MonteCarloMarginalizeCode/Code/test/test_lisa_synthetic_demo.py new file mode 100644 index 000000000..61c8d53f1 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lisa_synthetic_demo.py @@ -0,0 +1,104 @@ +"""Tests for the LISA synthetic-data demo analysis surface.""" + +import os +import subprocess + +import numpy as np + +from RIFT.misc import hyperpipeline_io + + +REPO_ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", "..")) +DEMO_DIR = os.path.join(REPO_ROOT, "MonteCarloMarginalizeCode", "Code", "demo", "rift", "lisa") +MAKE_INPUTS = os.path.join(DEMO_DIR, "make_synthetic_lisa_inputs.py") +RUN_ILE = os.path.join(DEMO_DIR, "run_lisa_synthetic_ile.sh") +RUN_END_TO_END = os.path.join(DEMO_DIR, "run_lisa_end_to_end.sh") + + +def _read_env_file(path): + values = {} + with open(path) as inp: + for line in inp: + key, value = line.strip().split("=", 1) + values[key] = value + return values + + +def test_lisa_synthetic_input_builder_writes_analysis_products(tmp_path): + subprocess.run( + [ + MAKE_INPUTS, + "--output-directory", + os.fspath(tmp_path), + "--duration", + "1024", + "--deltaT", + "4", + ], + check=True, + ) + + expected = { + "A-fake_strain-1000000-10000.h5", + "E-fake_strain-1000000-10000.h5", + "T-fake_strain-1000000-10000.h5", + "A_psd.xml.gz", + "E_psd.xml.gz", + "T_psd.xml.gz", + "lisa.cache", + "synthetic-params.env", + } + assert expected <= {path.name for path in tmp_path.iterdir()} + + env = _read_env_file(tmp_path / "synthetic-params.env") + assert env["SRATE"] == "0.25" + assert env["DELTA_T"] == "4.0" + assert env["DURATION"] == "1024.0" + + +def test_lisa_synthetic_demo_wires_real_analysis_arguments(tmp_path): + env = os.environ.copy() + env["RIFT_LISA_WORKDIR"] = os.fspath(tmp_path) + env["RIFT_LISA_RUN_ILE"] = "0" + subprocess.run([RUN_ILE, "--duration", "1024"], check=True, env=env) + + _, columns = hyperpipeline_io.read_table(os.fspath(tmp_path / "proposed-grid.dat")) + assert "ecliptic_longitude" in columns + assert "ecliptic_latitude" in columns + + ile_args = (tmp_path / "args_ile.txt").read_text() + assert "--zero-likelihood" not in ile_args + assert "--time-marginalization" in ile_args + assert "--cache-file {}".format(tmp_path / "lisa.cache") in ile_args + assert "--psd-file A={}".format(tmp_path / "A_psd.xml.gz") in ile_args + assert "--psd-file A=A_psd.xml.gz" not in ile_args + assert "--srate 0.25" in ile_args + assert "--data-integration-window-half 8.0" in ile_args + + transfer_files = (tmp_path / "helper_transfer_files.txt").read_text().splitlines() + assert os.fspath(tmp_path / "lisa.cache") in transfer_files + assert os.fspath(tmp_path / "A_psd.xml.gz") in transfer_files + + +def test_lisa_synthetic_demo_runs_real_ile(tmp_path): + env = os.environ.copy() + env["RIFT_LISA_WORKDIR"] = os.fspath(tmp_path) + subprocess.run([RUN_ILE, "--duration", "1024"], check=True, env=env) + + output = np.loadtxt(tmp_path / "lisa_ile_0_.dat") + assert output.shape == (15,) + assert np.isfinite(output[11]) + assert output[13] > 0 + + +def test_lisa_end_to_end_demo_renders_heavyweight_surface(tmp_path): + env = os.environ.copy() + env["RIFT_LISA_WORKDIR"] = os.fspath(tmp_path) + env["RIFT_LISA_RUN_ILE"] = "0" + subprocess.run([RUN_END_TO_END], check=True, env=env) + + assert (tmp_path / "event_0" / "A-fake_strain-1000000-10000.h5").exists() + assert (tmp_path / "event_0" / "A_psd.xml.gz").exists() + assert (tmp_path / "event_0" / "LISA_psd.txt").exists() + assert (tmp_path / "analysis_event_0" / "args_ile.txt").exists() + assert (tmp_path / "analysis_event_0" / "ILE.sub").exists()