diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 677a7d428..646ae725b 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -87,6 +87,10 @@ import_check: stage: system tests script: - python .travis/test-all-mod.py + - python -m pytest -q MonteCarloMarginalizeCode/Code/test/test_lisa_response_import.py + - python -m pytest -q MonteCarloMarginalizeCode/Code/test/test_lisa_lalsimutils_compat.py + - python -m pytest -q MonteCarloMarginalizeCode/Code/test/test_lisa_operational_synthetic.py + - python -m pytest -q MonteCarloMarginalizeCode/Code/test/test_lisa_helper_contract.py sim_manager_check: stage: unit tests diff --git a/MonteCarloMarginalizeCode/Code/RIFT/LISA/__init__.py b/MonteCarloMarginalizeCode/Code/RIFT/LISA/__init__.py new file mode 100644 index 000000000..01d2971e3 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/LISA/__init__.py @@ -0,0 +1,2 @@ +"""LISA-specific response and analysis helpers for RIFT.""" + diff --git a/MonteCarloMarginalizeCode/Code/RIFT/LISA/lalsimutils_compat.py b/MonteCarloMarginalizeCode/Code/RIFT/LISA/lalsimutils_compat.py new file mode 100644 index 000000000..71441effb --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/LISA/lalsimutils_compat.py @@ -0,0 +1,283 @@ +"""Compatibility helpers for the standalone LISA ILE fork. + +These routines are adapted from Aasim Z. Jan's LISA-RIFT +``lisa_rift_paper`` branch. Keep LISA-only shims here while the normal +``RIFT.lalsimutils`` and ILE driver continue to evolve independently. +""" + +import os +import sys + +import h5py +import lal +import lalsimulation as lalsim +import numpy as np + +import RIFT.lalsimutils as lalsimutils + + +def _filter_modes(hlms, modes): + if modes is None: + return hlms + + requested_modes = {tuple(map(int, mode)) for mode in modes} + return {mode: hlm for mode, hlm in hlms.items() if mode in requested_modes} + + +def hlmoff_for_LISA( + P, + Lmax=4, + modes=None, + fd_standoff_factor=0.964, + fd_alignment_postevent_time=None, + path_to_NR_hdf5=None, + **kwargs +): + """Generate frequency-domain modes with the sizing expected by LISA ILE.""" + assert Lmax >= 2 + assert (not np.isnan(P.m1)) and (not np.isnan(P.m2)), "Masses are NaN." + + extra_waveform_args = {} + if "extra_waveform_args" in kwargs: + extra_waveform_args.update(kwargs["extra_waveform_args"]) + extra_params = P.to_lal_dict_extended(extra_args_dict=extra_waveform_args) + + fNyq = 0.5 / P.deltaT + TDlen = int(1.0 / (P.deltaT * P.deltaF)) + if fd_alignment_postevent_time: + if fd_alignment_postevent_time >= TDlen * P.deltaT / 2: + print( + " Warning: fd alignment postevent time requested incompatible with short duration ", + file=sys.stderr, + ) + + if path_to_NR_hdf5 is not None: + hlms_struct = lalsimutils.hlmoft_from_NRhdf5( + path_to_NR_hdf5, + P, + Lmax, + only_mode=modes, + taper_percent=10, + beta=8, + verbose=True, + ) + return {mode: lalsimutils.DataFourier(hlms_struct[mode]) for mode in hlms_struct} + + fd_mode_approximants = { + lalsimutils.lalIMRPhenomHM, + lalsimutils.lalIMRPhenomXPHM, + lalsimutils.lalIMRPhenomXHM, + } + if P.approx in fd_mode_approximants: + hlms_struct = lalsim.SimInspiralChooseFDModes( + P.m1, + P.m2, + P.s1x, + P.s1y, + P.s1z, + P.s2x, + P.s2y, + P.s2z, + P.deltaF, + P.fmin * fd_standoff_factor, + fNyq, + P.fref, + P.phiref, + P.dist, + P.incl, + extra_params, + P.approx, + ) + hlmsdict = lalsimutils.SphHarmFrequencySeries_to_dict(hlms_struct, Lmax) + hlmsdict = _filter_modes(hlmsdict, modes) + return { + mode: lal.ResizeCOMPLEX16FrequencySeries(hlm, 0, TDlen) + for mode, hlm in hlmsdict.items() + } + + if P.approx in {lalsimutils.lalNRHybSur3dq8, lalsimutils.lalIMRPhenomD}: + hlms_struct = lalsimutils.hlmoff(P, Lmax=Lmax) + hlmsdict = lalsimutils.SphHarmFrequencySeries_to_dict(hlms_struct, Lmax) + hlmsdict = _filter_modes(hlmsdict, modes) + for mode, hlm in list(hlmsdict.items()): + if not (1 / hlm.deltaF == P.deltaT * TDlen): + print( + "WARNING: RESIZING IN FD DOMAIN " + f"(1/deltaF = {1 / hlm.deltaF}, deltaT*TDlen = {P.deltaT * TDlen}), " + "THIS SHOULD NOT BE HAPPENING." + ) + hlmsdict[mode] = lal.ResizeCOMPLEX16FrequencySeries(hlm, 0, TDlen) + return hlmsdict + + raise ValueError(f"Unsupported LISA FD mode approximant: {P.approx}") + + +def _cache_path_for_channel(fname, channel): + cache_data = np.loadtxt(fname, dtype=str) + if cache_data.ndim == 1: + cache_data = cache_data.reshape(1, len(cache_data)) + + channel_prefix = channel[0] + for row in cache_data: + if row[0] == channel_prefix: + raw_path = row[-1] + if "localhost" in raw_path: + raw_path = raw_path.split("localhost", 1)[1] + elif raw_path.startswith("file://"): + raw_path = raw_path[len("file://"):] + return os.path.expanduser(raw_path) + + raise ValueError(f"Could not find channel {channel!r} in cache {fname!r}") + + +def frame_h5_to_hoff(fname, channel, start=None, stop=None, verbose=True): + """Read LISA frequency-domain HDF5 data from a cache entry.""" + if verbose: + print(" ++ Loading from cache ", fname, channel) + + path_to_h5 = _cache_path_for_channel(fname, channel) + if verbose: + print(f"Reading h5 file {path_to_h5}") + + with h5py.File(path_to_h5, "r") as data: + hoff = lal.CreateCOMPLEX16FrequencySeries( + "hoff", + data.attrs["epoch"], + data.attrs["f0"], + data.attrs["deltaF"], + lalsimutils.lsu_HertzUnit, + int(data.attrs["length"]), + ) + hoff.data.data = data["data"][()] + + return hoff + + +def frame_h5_to_hoft(fname, channel, start=None, stop=None, verbose=True): + """Read LISA time-domain HDF5 data from a cache entry.""" + if verbose: + print(" ++ Loading from cache ", fname, channel) + + path_to_h5 = _cache_path_for_channel(fname, channel) + if verbose: + print(f"Reading h5 file {path_to_h5}") + + with h5py.File(path_to_h5, "r") as data: + hoft = lal.CreateREAL8TimeSeries( + "hoft", + data.attrs["epoch"], + data.attrs["f0"], + data.attrs["deltaT"], + lal.DimensionlessUnit, + int(data.attrs["length"]), + ) + hoft.data.data = data["data"][()] + + return hoft + + +def frame_data_to_non_herm_hoff( + fname, + channel, + start=None, + stop=None, + TDlen=0, + window_shape=0.0, + verbose=True, + deltaT=None, + h5_frame=False, +): + """Read frame/HDF5 time-domain data and FFT to two-sided frequency data.""" + if not h5_frame: + return lalsimutils.frame_data_to_non_herm_hoff( + fname, + channel, + start=start, + stop=stop, + TDlen=TDlen, + window_shape=window_shape, + verbose=verbose, + deltaT=deltaT, + ) + + ht = frame_h5_to_hoft(fname, channel, start, stop, verbose) + + tmplen = ht.data.length + if TDlen == -1: + TDlen = tmplen + elif TDlen == 0: + TDlen = lalsimutils.nextPow2(tmplen) + else: + assert TDlen >= tmplen + + ht = lal.ResizeREAL8TimeSeries(ht, 0, TDlen) + hoftC = lal.CreateCOMPLEX16TimeSeries( + "h(t)", ht.epoch, ht.f0, ht.deltaT, ht.sampleUnits, TDlen + ) + hoftC.data.data = ht.data.data + 0j + fwdplan = lal.CreateForwardCOMPLEX16FFTPlan(TDlen, 0) + hf = lal.CreateCOMPLEX16FrequencySeries( + "Template h(f)", + ht.epoch, + ht.f0, + 1.0 / ht.deltaT / TDlen, + lalsimutils.lsu_HertzUnit, + TDlen, + ) + lal.COMPLEX16TimeFreqFFT(hf, hoftC, fwdplan) + if verbose: + print( + " ++ Loaded data h(f) of length n= ", + hf.data.length, + " (= ", + len(hf.data.data) * ht.deltaT, + "s) at sampling rate ", + 1.0 / ht.deltaT, + ) + return hf + + +def print_params_lisa(self, show_system_frame=False): + """LISA-flavored parameter printer used by the standalone ILE fork.""" + print("This ChooseWaveformParams has the following parameter values:") + print(f"m1 = {self.m1 / lalsimutils.lsu_MSUN / 1e3:0.3f} x 1e3 (Msun)") + print(f"m2 = {self.m2 / lalsimutils.lsu_MSUN / 1e3:0.3f} x 1e3 (Msun)") + print("s1x =", self.s1x) + print("s1y =", self.s1y) + print("s1z =", self.s1z) + print("s2x =", self.s2x) + print("s2y =", self.s2y) + print("s2z =", self.s2z) + if show_system_frame: + self.print_params(show_system_frame=show_system_frame) + print("lambda1 =", self.lambda1) + print("lambda2 =", self.lambda2) + print("inclination =", self.incl) + print("distance =", self.dist / 1.0e9 / lalsimutils.lsu_PC, "(Gpc)") + print("reference orbital phase =", self.phiref) + print("polarization angle =", self.psi) + print("eccentricity = ", self.eccentricity) + print("reference time = ", float(self.tref), "(s)") + print("detector is: LISA") + print("ecliptic latitude (beta):", self.theta, "(radians)") + print("ecliptic longitude (lambda):", self.phi, "(radians)") + print("starting frequency is =", self.fmin, "(Hz)") + print("reference frequency is =", self.fref, "(Hz)") + print("Max frequency is =", self.fmax, "(Hz)") + print("time step =", self.deltaT, "(s) <==>", 1.0 / self.deltaT, "(Hz) sample rate") + print("freq. bin size is =", self.deltaF, "(Hz)") + + +def install_choose_waveform_print_params_lisa(): + """Install the LISA printer on ChooseWaveformParams for the LISA fork.""" + lalsimutils.ChooseWaveformParams.print_params_lisa = print_params_lisa + + +__all__ = [ + "frame_data_to_non_herm_hoff", + "frame_h5_to_hoff", + "frame_h5_to_hoft", + "hlmoff_for_LISA", + "install_choose_waveform_print_params_lisa", + "print_params_lisa", +] diff --git a/MonteCarloMarginalizeCode/Code/RIFT/LISA/response/LISA_response.py b/MonteCarloMarginalizeCode/Code/RIFT/LISA/response/LISA_response.py new file mode 100644 index 000000000..a8855671f --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/LISA/response/LISA_response.py @@ -0,0 +1,727 @@ +"""This code is based on http://arxiv.org/abs/1806.10734 and 10.1103/PhysRevD.103.083011, with parts of code being taken from BBHx's C code. The response has been validated again BBHx with mismatch of around 10^(-10). This assumes fixed LISA arm length and is not the best formalism for Stellar mass binaries.""" + +# TO DO: +## Better documentation +## Impact of different fourier convention. Right now we follow theirs, not lal's. + +# Future work +## 1) Add a different response code for Stellar mass binaries. + +# Speed ups: +# https://stackoverflow.com/questions/49493482/numpy-np-multiply-vs-operator np.multiply vs * (no speedup) +# https://stackoverflow.com/questions/25870923/how-to-square-or-raise-to-a-power-elementwise-a-2d-numpy-array np.square vs **2 (no speedup for array of length 1e7). np.power is slow +# https://stackoverflow.com/questions/49459661/differences-between-numpy-divide-and-python-divide np.divide vs / (no speedup) + +import numpy as np +import lal +import RIFT.lalsimutils as lsu +import RIFT.LISA.lalsimutils_compat as lisa_lalsimutils_compat +import sys +import h5py + +__author__ = "A. Jan" +########################################################################################### +# Constants +########################################################################################### +e = 0.004824185218078991 +omega0 = 1.9909865927683788 * 10**(-7) #1/seconds +a = 149597870700. #meters +C_SI = 299792458. +L = 2*np.sqrt(3)*a*e +YRSID_SI = 31558149.763545603 + +########################################################################################### +# FUNCTIONS +########################################################################################### +def create_lal_frequency_series(frequency_values, frequency_series, deltaF, epoch = 950000000, f0 = 0.0): + """A helper function to create lal COMPLEX16FrequencySeries. + Args: + frequency_values (numpy.array): Frequency values at which the series is defined. + frequency_series (numpy.array): Corresponding strain in frequency domain. + epoch (float): Needed to create COMPLEX16FrequencySeries, by default it is 950000000. + f0 (float): Needed to create COMPLEX16FrequencySeries, by default it is 0.0 . + Output: + lal.COMPLEX16FrequencySeries object""" + assert len(frequency_values) == len(frequency_series), "frequency_values and frequency_series don't have the same length." + hf_lal = lal.CreateCOMPLEX16FrequencySeries("hf", epoch, f0, deltaF, lal.HertzUnit, len(frequency_values)) + hf_lal.data.data = frequency_series + return hf_lal + +def create_lal_COMPLEX16TimeSeries(deltaT, time_series, epoch = 950000000, f0 = 0.0, data_is_real = True): + """A helper function to create lal COMPLEX16TimeSeries. + Args: + deltaT (float): time step (1/sampling rate). + time_series (numpy.array): strain in time domain. + epoch (float): Needed to create COMPLEX16TimeSeries, by default it is 950000000. + f0 (float): Needed to create COMPLEX16TimeSeries, by default it is 0.0 . + Output: + lal.COMPLEX16TimeSeries object""" + ht_lal = lal.CreateCOMPLEX16TimeSeries("ht_lal", epoch, f0, deltaT, lal.DimensionlessUnit, len(time_series)) + if data_is_real: + ht_lal.data.data = time_series + 0j + else: + ht_lal.data.data = time_series + return ht_lal + + +def convert_double_sided_to_single_sided(frequency_values, frequency_series, data_defined="negative"): + assert len(frequency_values) == len(frequency_series), "frequency_values and frequency_series don't have the same length." + if data_defined == "negative": + print("Negative") + index = np.argwhere(frequency_values<=0).flatten() + hf_onesided = create_lal_frequency_series(frequency_values[index], np.conj(frequency_series[index][::-1])) # do I need to conjugate? + assert len(frequency_series)//2 + 1 == hf_onesided.data.length + return hf_onesided + elif data_defined == "positive": + print("Positve") + index = np.argwhere(frequency_values>=0).flatten() + hf_onesided = create_lal_frequency_series(frequency_values[index], frequency_series[index]) + assert len(frequency_series)//2 + 1 == hf_onesided.data.length + return hf_onesided + else: + raise ValueError("data_defined must be either 'negative' or 'positive'") + +def frequency_series_double_sided(frequency_values, frequency_series, data_defined = "negative"): + hf_small = convert_double_sided_to_single_sided(frequency_values, frequency_series, data_defined) + tmp = np.zeros(len(frequency_series), dtype=complex) + tmp[:len(hf_small.data.data)] = np.conj(hf_small.data.data[::-1]) + tmp[len(hf_small.data.data)-1:] = (hf_small.data.data[:-1]) + hf = create_lal_frequency_series(frequency_values, tmp) + return hf + +def get_fvals(frequency_series): + """A function to evaulate frequency values of a COMPLEX16FrequencySeries. Goes from [-fNyq, fNyq - deltaF]. + Args: + frequency_series (COMPLEX16FrequencySeries): + Output: + frequency array (numpy.array)""" + fvals = -frequency_series.deltaF*np.arange(frequency_series.data.length//2, -frequency_series.data.length//2, -1) + return fvals + + + +def get_Ylm(inclination, phiref, l ,m, s = -2): + """A function to call spherical harmonics, should change to calling the GPU version in future but for now, I am only using to create injections so it should be fine. + Args: + inclination (float): inclination in SSB frame, + phiref (float): phase at coalescence, + l, m: modes, + s: spin weight, -2 by default. + Output: + Spin weighted spherical harmonics (complex)""" + return lal.SpinWeightedSphericalHarmonic(inclination, phiref, s, l, m) + +def get_closest_index(array, value): + """A function that gives you the index at which the array has a value closest to the one you want. + Args: + array (np.array): + value (float): + Output: + index (float)""" + return np.argmin(np.abs(array - value)) + +def transformed_Hplus_Hcross(beta, lamda, psi, theta, phiref, l, m): + """This function transforms the plus and cross polarization from wave frame to SSB frame. + Args: + beta (float) = ecliptic latitude, + lamda (float) = ecliptic longitude, + psi (float) = polarization + Returns: + transformed Plm (used to create the transfer function, equation 21 (http://arxiv.org/abs/2003.00357) + """ + + # defining cos and sin for ease (checked) + cl, sl = np.cos(lamda), np.sin(lamda) + cb, sb = np.cos(beta), np.sin(beta) + cp, sp = np.cos(psi), np.sin(psi) + + # polarization in waveframe (checked) + Hp = [[1,0,0], [0,-1,0], [0,0,0]] + Hc = [[0,1,0], [1,0,0] , [0,0,0]] + + # Rotation matrix (checked) + O1 = np.zeros((3,3)) + + O1[0,0] = cp*sl - cl*sb*sp + O1[0,1] = -cl*cp*sb - sl*sp + O1[0,2] = -cb*cl + + O1[1,0] = -cl*cp - sb*sl*sp + O1[1,1] = -cp*sb*sl + cl*sp + O1[1,2] = -cb*sl + + O1[2,0] = cb*sp + O1[2,1] = cb*cp + O1[2,2] = -sb + + # Transpose of Rotation matrix (checked) + TO1 = np.zeros((3,3)) + + TO1[0,0] = cp*sl - cl*sb*sp + TO1[0,1] = -cp*cl - sl*sb*sp + TO1[0,2] = cb*sp + + TO1[1,0] = -cl*cp*sb - sl*sp + TO1[1,1] = -sl*cp*sb + cl*sp + TO1[1,2] = cb*cp + + TO1[2,0] = -cb*cl + TO1[2,1] = -cb*sl + TO1[2,2] = -sb + + # Ylm factors (checked) + # For injection, we shouldn't take -phiref as we do in marginalization. Confirm with ROS. + Ylm = get_Ylm(theta, phiref, l, m, -2) + Y_lm = (-1)**(l) * np.conj(get_Ylm(theta, phiref, l, -m, -2)) + + Yfactorplus = 0.5*(Ylm + Y_lm) + Yfactorcross = 0.5*1j*(Ylm - Y_lm) + + Plm = Yfactorplus*np.matmul(O1, np.matmul(Hp,TO1)) + Yfactorcross*np.matmul(O1, np.matmul(Hc,TO1)) + return Plm + + +def get_tf_from_phase(hlm, fmax, debug = False):#tested + """This function differentiates phase to get tf. Similar to pycbc's time_from_frequencyseries (waveforms/utils.py) function but does not include their discontinuity check. (Now it does have those checks). + Args: + hlm (COMPLEX16FrequencySeries): The mode for which you are calculating tf , + fmax (float): maximum frequency (fNyq for RIFT). + Returns: + tf (numpy.array): Time, + frequency (numpy.array): Frequency (numpy.array). + """ + # send in hlm not shifted in time + + # get frequency and mode data + freq = np.arange(-fmax, fmax+hlm.deltaF, hlm.deltaF) + if debug: + print(f"len(freq) = {len(freq)}, freq[0] = {freq[0]} Hz, freq[-1] = {freq[-1]} Hz, len(hlm) = {hlm.data.length}") + + # get amplitude and phase + phase = np.unwrap(np.angle(hlm.data.data)) + # compute tf = -1/(2pi) * d(phase)/df + phase = phase - phase[0] #Pycbc does this, doesn't change answer as expected. + dphi = np.unwrap(np.diff(phase)) + time = -dphi / (2.*np.pi*np.diff(freq)) + # diff reduces len by 1 so artifically increasing it by adding an extra zero at the end + tmp = np.zeros(len(time)+1) + tmp[:-1] = time + time = tmp + + # only focusing on f bins where data exists + nzidx = np.nonzero(abs(hlm.data.data))[0] + if debug: + print(nzidx) + print(f"tf[0](after stripping zeros) = {time[kmin]}, tf[-1](after stripping zeros) {time[kmax]}") + kmin, kmax = nzidx[0], nzidx[-2] + time[:kmin] = time[kmin] + time[kmax:] = time[kmax] + if debug: + print(f"len(time) = {len(time)}, len(freq) = {len(freq)}") + + # saving data + return time, freq[::-1] #inverting frequency array to match http://arxiv.org/abs/1806.10734i fourier convention, will change it as we validate these codes. + +def Evaluate_Gslr(tf, f, H, beta, lamda): + """This function takes in tf, f, Plm (from transformed_Hplus_Hcross), beta and lamda to generate transfer function for a given mode. yslr = Sum_l^r Gslr * hlm. + Args: + tf (numpy.array)= -1/2pi d(phase)/df, + f (numpy.array)= frequency array, + H (numpy.array)= The plus and cross polarization matrices transformed from wave to SSB frame , + beta (float) = ecliptic latitude, + lamda (float) = ecliptic longitude + Returns: + Transfer function L1 (numpy.array), Transfer function L2 (numpy.array), Transfer function L3 (numpy.array) + """ + alpha = omega0*tf + c, s = np.cos(alpha), np.sin(alpha) + k = np.array([-np.cos(beta)*np.cos(lamda), -np.cos(beta)*np.sin(lamda), -np.sin(beta)]) + p0 = np.array([a*c, a*s, np.zeros(len(tf))]) # (3, N) + kR = np.dot(k, p0) # (N,) + phaseRdelay = 2.*np.pi/C_SI *f*kR #(N,) + + p1L = np.array([-a*e*(1 + s*s), a*e*c*s, -a*e*np.sqrt(3)*c]) # (3, N) + p2L = np.array([a*e/2*(np.sqrt(3)*c*s + (1 + s*s)), a*e/2*(-c*s - np.sqrt(3)*(1 + c*c)), -a*e*np.sqrt(3)/2*(np.sqrt(3)*s - c)]) # (3, N) + p3L = np.array([a*e/2*(-np.sqrt(3)*c*s + (1 + s*s)), a*e/2*(-c*s + np.sqrt(3)*(1 + c*c)), -a*e*np.sqrt(3)/2*(-np.sqrt(3)*s - c)]) # (3, N) + + + n1 = np.array([-1./2*c*s, 1./2*(1 + c*c), np.sqrt(3)/2*s]) # (3, N) + kn1= np.dot(k, n1) #(N,) + n1Hn1 = np.einsum("ij,ji->i",n1.T, np.einsum("ij,jk", H, n1)) + + n2 = 1./4. * np.array([c*s - np.sqrt(3)*(1 + s*s), np.sqrt(3)*c*s - (1 + c*c), -np.sqrt(3)*s - 3*c]) + kn2= np.dot(k, n2) + n2Hn2 = np.einsum("ij,ji->i",n2.T, np.einsum("ij,jk", H, n2)) + + n3 = 1./4*np.array([c*s + np.sqrt(3)*(1 + s*s), -np.sqrt(3)*c*s - (1 + c*c), -np.sqrt(3)*s + 3*c]) + kn3= np.dot(k, n3) + n3Hn3 = np.einsum("ij,ji->i",n3.T, np.einsum("ij,jk", H, n3)) + + + + kp1Lp2L = np.dot(k, (p1L+p2L)) + kp2Lp3L = np.dot(k, (p2L+p3L)) + kp3Lp1L = np.dot(k, (p3L+p1L)) + kp0 = np.dot(k, p0) + + factorcexp0 = np.exp(1j*2.*np.pi*f/C_SI * kp0) + prefactor = np.pi*f*L/C_SI + + factorcexp12 = np.exp(1j*prefactor * (1.+kp1Lp2L/L)) + factorcexp23 = np.exp(1j*prefactor * (1.+kp2Lp3L/L)) + factorcexp31 = np.exp(1j*prefactor * (1.+kp3Lp1L/L)) + + factorsinc12 = np.sinc( (prefactor * (1.-kn3))/np.pi) + factorsinc21 = np.sinc( (prefactor * (1.+kn3))/np.pi) + factorsinc23 = np.sinc( (prefactor * (1.-kn1))/np.pi) + factorsinc32 = np.sinc( (prefactor * (1.+kn1))/np.pi) + factorsinc31 = np.sinc( (prefactor * (1.-kn2))/np.pi) + factorsinc13 = np.sinc( (prefactor * (1.+kn2))/np.pi) + + commonfac = 1j*prefactor*factorcexp0 + G12 = commonfac * n3Hn3 * factorsinc12 * factorcexp12 * np.exp(-1j*phaseRdelay) + G21 = commonfac * n3Hn3 * factorsinc21 * factorcexp12 * np.exp(-1j*phaseRdelay) + G23 = commonfac * n1Hn1 * factorsinc23 * factorcexp23 * np.exp(-1j*phaseRdelay) + G32 = commonfac * n1Hn1 * factorsinc32 * factorcexp23 * np.exp(-1j*phaseRdelay) + G31 = commonfac * n2Hn2 * factorsinc31 * factorcexp31 * np.exp(-1j*phaseRdelay) + G13 = commonfac * n2Hn2 * factorsinc13 * factorcexp31 * np.exp(-1j*phaseRdelay) + + x = np.pi*f*L/C_SI + z = np.exp(1j*2.*x) + + factor_convention = 2 + factorAE = 1j*np.sqrt(2)*np.sin(2.*x)*z + factorT = 2.*np.sqrt(2)*np.sin(2.*x)*np.sin(x)*np.exp(1j*3.*x) + + Araw = 0.5 * ( (1.+z)*(G31 + G13) - G23 - z*G32 - G21 - z*G12 ) + Eraw = 0.5*1/np.sqrt(3) * ( (1.-z)*(G13 - G31) + (2.+z)*(G12 - G32) + (1.+2.*z)*(G21 - G23) ) + Traw = 1/np.sqrt(6) * (G21 - G12 + G32 - G23 + G13 - G31) + transferL1 = factor_convention * factorAE * Araw + transferL2 = factor_convention * factorAE * Eraw + transferL3 = factor_convention * factorT * Traw + + return transferL1, transferL2, transferL3 + + +########################################################################################### +# FUNCTIONS USED IN RECOVERY +########################################################################################### +def get_amplitude_phase(hf): #tested + # NOTE: Send in hlm not shifted in time, I add time shift in precomputation. + """This function splits h(f) into Amplitude A(f) and phase phase(f). We then express h_lm as A_lm(f) * exp(i*phase(f)). Note the sign in the exponent. + Args: + hf (CreateCOMPLEX16FrequencySeries): Frequency domain waveform. + Returns: + Amplitude as a function of frequency (numpy.array), phase as a function of frequency (numpy.array) + """ + return np.abs(hf.data.data), np.unwrap(np.angle(hf.data.data)) + + + +def get_tf_from_phase_dict(hlm, fmax, fref=None, debug=True, shift=True):#tested + """This function differentiates phase for each mode to get tf. Similar to pycbc's time_from_frequencyseries (waveforms/utils.py) function. + Args: + hlm (dict): mode dict generated by std_and_conj_hlmoff or any other hlm(f) function in RIFT, + fmax (float): maximum frequency (fNyq for RIFT). + Returns: + tf (dict): tf for each mode in hlm, + frequency (dict): corresponding frequency for each mode in hlm, + amplitude (dict): frequency domain amplitude for each mode, + phase (dict): frequency domain phase for each mode. + """ + # send in hlm not shifted in time + tf_dict = {} + freq_dict = {} + amp_dict = {} + phase_dict = {} + print("Computing time frequency correspondence for mode") + modes = np.array(list(hlm.keys())) + #freq = -lsu.evaluate_fvals(hlm[tuple(modes[0])]) # THIS CONSUMES MOST TIME + freq = -hlm[tuple(modes[0])].deltaF*np.arange(hlm[tuple(modes[0])].data.length//2, -hlm[tuple(modes[0])].data.length//2, -1) # this matches evaluate_fvals without using for loops + for mode in modes: + print(f"\n\tMode = {mode}") + mode = tuple(mode) + # get amplitude and phase + amp, phase = get_amplitude_phase(hlm[mode]) + # compute tf = -1/(2pi) * d(phase)/df + dphi = np.unwrap(np.diff(phase)) + time = np.divide(-dphi, (2.*np.pi*hlm[mode].deltaF)) + # diff reduces len by 1 so artifically increasing it by adding an extra zero at the end + tmp = np.zeros(len(time)+1) + tmp[:-1] = time + time = tmp + + # only focusing on f bins where data exists + # I had to introduce this statement since sometimes a mode doesn't have data (odd m modes are not excited for q=1, so the mode content is all zero.) + try: + nzidx = np.nonzero(abs(hlm[mode].data.data))[0] + kmin, kmax = nzidx[0], nzidx[-2] + time[:kmin] = time[kmin] + time[kmax:] = time[kmax] + except: + print(f"No data for {mode}") + pass + # saving data + tf_dict[mode] = time + freq_dict[mode] = freq[::-1] + amp_dict[mode] = amp + phase_dict[mode] = phase + + if shift: + modes = list(hlm.keys()) + print(f"Shifting of time and phase with fref = {fref}.") + assert (2,2) in modes, "(2,2) mode needs to be present." + # phase and tf shifts + if not fref: + # if fref not provided, set it to frequency at max (f^2 * A_{2,2}(f)) (BBHx) + fref = freq_dict[2,2][np.argmax(freq_dict[2,2]**2 * amp_dict[2,2])] # frequency at max (f^2 * A_{2,2}(f)) + + # find tf at fref + index_at_fref = get_closest_index(freq_dict[2,2], fref) + tf_22_current = tf_dict[2,2][index_at_fref] + phase_22_current = phase_dict[2,2][index_at_fref] + + time_shift = tf_22_current + reference_phase = 0.0 + + # for loop needs to start with (2,2) mode + modes.remove((2,2)) + modes.insert(0, (2,2)) + if debug: + print(f"tf[2,2] at fref ({freq_dict[2,2][index_at_fref]} Hz) before shift is {tf_22_current}s (phase[2,2] = {phase_22_current}).") + + # subtract that from all modes. tf for (2,2) needs to be zero at fref, I will add t_ref to all modes later (create_lisa_injections for injections and precompute for recovery), making tf=t_ref at fref. + for mode in modes: + if debug: + print(f"\tShifting {mode}") + tf_dict[mode] = tf_dict[mode] - time_shift # confirmed that I don't need to set all modes tf as 0. Conceptually, for the same time the other modes will be at a different frequency. + phase_dict[mode] = phase_dict[mode] - 2*np.pi*time_shift*freq_dict[mode] + if mode == (2,2): + phase_22_current = phase_dict[2,2][index_at_fref] + difference = reference_phase - phase_22_current + phase_dict[mode] = phase_dict[mode] + mode[1]/2 * difference + print(f"{mode}, phase = {phase_dict[mode][index_at_fref]}") + if debug: + print(f"tf[2,2] at fref ({fref} Hz) after shift is {tf_dict[2,2][index_at_fref]} (phase[2,2] = {phase_dict[2,2][index_at_fref]}).") + + return tf_dict, freq_dict, amp_dict, phase_dict + + +def get_beta_lamda_psi_terms_Hp(beta, lamda, psi): + """This function gives beta lamda psi terms for each term when we split up n_l * P_lm * n_l in equation 21 of http://arxiv.org/abs/2003.00357. We need this to bring out the psi dependence to marginalize over it. This gives those terms after we transform the plus polarization's frame to SSB frame. + Args: + beta (float): ecliptic latitude, + lamda (float): ecliptic longitude, + psi (numpy.array of shape (n,1): polarization angle array. + Return: + xx, xy, xz, yy, yz, zz terms (each is a numpy array of shape (n,1) + + """ + xx_term = (np.cos(psi)*np.sin(lamda) - np.cos(lamda)*np.sin(beta)*np.sin(psi))**2 + \ + (-np.cos(lamda)*np.cos(psi)*np.sin(beta) - np.sin(lamda)*np.sin(psi))*(np.cos(lamda)*np.cos(psi)*np.sin(beta) + np.sin(lamda)*np.sin(psi)) + + xy_term = ((np.cos(psi)*np.sin(beta)*np.sin(lamda)-np.cos(lamda)*np.sin(psi)) * (-np.cos(lamda)*np.cos(psi)*np.sin(beta) - np.sin(lamda)*np.sin(psi)) + + (np.cos(psi)*np.sin(lamda)-np.cos(lamda)*np.sin(beta)*np.sin(psi))*(-np.cos(lamda)*np.cos(psi) - np.sin(beta)*np.sin(lamda)* np.sin(psi))) + \ + ((-np.cos(psi)*np.sin(beta)*np.sin(lamda)+ np.cos(lamda)*np.sin(psi))*(np.cos(lamda)*np.cos(psi)*np.sin(beta) + np.sin(lamda)*np.sin(psi)) + + (np.cos(psi)*np.sin(lamda)-np.cos(lamda)*np.sin(beta)*np.sin(psi))*(-np.cos(lamda)*np.cos(psi) - np.sin(beta)*np.sin(lamda)*np.sin(psi))) + + xz_term = (np.cos(beta)*np.sin(psi)*(np.cos(psi)*np.sin(lamda) - np.cos(lamda)*np.sin(beta)*np.sin(psi)) + + np.cos(beta)*np.cos(psi)*(np.cos(lamda)*np.cos(psi)*np.sin(beta) + np.sin(lamda)*np.sin(psi))) + \ + (np.cos(beta)*np.sin(psi)*(np.cos(psi)*np.sin(lamda) - np.cos(lamda)*np.sin(beta)*np.sin(psi)) - + np.cos(beta)*np.cos(psi)*(-np.cos(lamda)*np.cos(psi)*np.sin(beta) - np.sin(lamda)*np.sin(psi))) + + yy_term = (np.cos(psi)*np.sin(beta)*np.sin(lamda) - np.cos(lamda)*np.sin(psi))*(-np.cos(psi)*np.sin(beta)*np.sin(lamda)+ np.cos(lamda)*np.sin(psi)) + \ + (-np.cos(lamda)*np.cos(psi) - np.sin(beta)*np.sin(lamda)*np.sin(psi))**2 + + yz_term = (np.cos(beta)*np.cos(psi)*(np.cos(psi)*np.sin(beta)*np.sin(lamda) - np.cos(lamda)*np.sin(psi)) + + np.cos(beta)*np.sin(psi)*(-np.cos(lamda)*np.cos(psi) - np.sin(beta)*np.sin(lamda)*np.sin(psi))) + \ + (-np.cos(beta)*np.cos(psi)*(-np.cos(psi)*np.sin(beta)*np.sin(lamda) + np.cos(lamda)*np.sin(psi)) + + np.cos(beta)*np.sin(psi)*(-np.cos(lamda)*np.cos(psi) - np.sin(beta)*np.sin(lamda)*np.sin(psi))) + + zz_term = (-np.cos(beta)**2 * np.cos(psi)**2 + np.cos(beta)**2 * np.sin(psi)**2) + + combined = np.vstack([xx_term, xy_term, xz_term, yy_term, yz_term, zz_term]) + return combined + +def get_beta_lamda_psi_terms_Hc(beta, lamda, psi): + """This function gives beta lamda psi terms for each term when we split up n_l * P_lm * n_l in equation 21 of http://arxiv.org/abs/2003.00357. We need this to bring out the psi dependence to marginalize over +it. This gives those terms after we transform the cross polarization's frame to SSB frame. + Args: + beta (float): ecliptic latitude, + lamda (float): ecliptic longitude, + psi (numpy.array of shape (n,1): polarization angle array. + Return: + xx, xy, xz, yy, yz, zz terms (each is a numpy array of shape (n,1) + + """ + xx_term = 2*(np.cos(psi)*np.sin(lamda) - np.cos(lamda)*np.sin(beta)*np.sin(psi))*(-np.cos(lamda)*np.cos(psi)*np.sin(beta) - np.sin(lamda)*np.sin(psi)) + + xy_term = ((-np.cos(psi)*np.sin(beta)*np.sin(lamda) + np.cos(lamda)*np.sin(psi))*(np.cos(psi)*np.sin(lamda) - np.cos(lamda)*np.sin(beta)*np.sin(psi)) + + (-np.cos(lamda)*np.cos(psi)*np.sin(beta) - np.sin(lamda)*np.sin(psi))*(-np.cos(lamda)*np.cos(psi) - np.sin(beta)*np.sin(lamda)*np.sin(psi))) + \ + ((-np.cos(psi)*np.sin(beta)*np.sin(lamda) + np.cos(lamda)*np.sin(psi))*(np.cos(psi)*np.sin(lamda) - np.cos(lamda)*np.sin(beta)*np.sin(psi)) + + (-np.cos(lamda)*np.cos(psi)*np.sin(beta) - np.sin(lamda)*np.sin(psi))*(-np.cos(lamda)*np.cos(psi) - np.sin(beta)*np.sin(lamda)*np.sin(psi))) + + xz_term = (np.cos(beta)*np.cos(psi)*(np.cos(psi)*np.sin(lamda) - np.cos(lamda)*np.sin(beta)*np.sin(psi)) + + np.cos(beta)*np.sin(psi)*(-np.cos(lamda)*np.cos(psi)*np.sin(beta) - np.sin(lamda)*np.sin(psi))) + \ + (np.cos(beta)*np.cos(psi)*(np.cos(psi)*np.sin(lamda) - np.cos(lamda)*np.sin(beta)*np.sin(psi)) + + np.cos(beta)*np.sin(psi)*(-np.cos(lamda)*np.cos(psi)*np.sin(beta) - np.sin(lamda)*np.sin(psi))) + + yy_term = 2*(-np.cos(psi)*np.sin(beta)*np.sin(lamda) + np.cos(lamda)*np.sin(psi))*(-np.cos(lamda)*np.cos(psi) - np.sin(beta)*np.sin(lamda)*np.sin(psi)) + + yz_term = (np.cos(beta)*np.sin(psi)*(-np.cos(psi)*np.sin(beta)*np.sin(lamda) + np.cos(lamda)*np.sin(psi)) + + np.cos(beta)*np.cos(psi)*(-np.cos(lamda)*np.cos(psi) - np.sin(beta)*np.sin(lamda)*np.sin(psi))) + \ + (np.cos(beta)*np.sin(psi)*(-np.cos(psi)*np.sin(beta)*np.sin(lamda) + np.cos(lamda)*np.sin(psi)) + + np.cos(beta)*np.cos(psi)*(-np.cos(lamda)*np.cos(psi) - np.sin(beta)*np.sin(lamda)*np.sin(psi))) + + zz_term = 2 * np.cos(beta)**2 * np.cos(psi) * np.sin(psi) + + combined = np.vstack([xx_term, xy_term, xz_term, yy_term, yz_term, zz_term]) + return combined + + +def Evaluate_Gslr_test_2(tf, f, beta, lamda): + """This is the main function, takes in tf, f, beta and lamda to generate transfer function for a given mode for each xx, xy, xz, yy, yz and zz term. (need to explain this in paper) + Args: + tf (numpy.array)= -1/2pi d(phase)/df, + f (numpy.array)= frequency array, + beta (float) = ecliptic latitude, + lamda (float) = ecliptic longitude + Returns: + Transfer function L1 (numpy.array with xx, xy, xz, yy, yz, zz), Transfer function L2 (numpy.array with xx, xy, xz, yy, yz, zz), Transfer function L3 (numpy.array with xx, xy, xz, yy, yz, zz) + """ + alpha = omega0*tf + c, s = np.cos(alpha), np.sin(alpha) + k = np.array([-np.cos(beta)*np.cos(lamda), -np.cos(beta)*np.sin(lamda), -np.sin(beta)]) + p0 = np.array([a*c, a*s, np.zeros(len(tf))]) # (3, N) + kR = np.dot(k, p0) # (N,) + phaseRdelay = 2.*np.pi/C_SI *f*kR #(N,) + + p1L =np.array([-a*e*(1 + s*s), a*e*c*s, -a*e*np.sqrt(3)*c]) # (3, N) + p2L =np.array([a*e/2*(np.sqrt(3)*c*s + (1 + s*s)), a*e/2*(-c*s - np.sqrt(3)*(1 + c*c)), -a*e*np.sqrt(3)/2*(np.sqrt(3)*s - c)]) # (3, N) + p3L =np.array([a*e/2*(-np.sqrt(3)*c*s + (1 + s*s)), a*e/2*(-c*s + np.sqrt(3)*(1 + c*c)), -a*e*np.sqrt(3)/2*(-np.sqrt(3)*s - c)]) # (3, N) + + + n1 = np.array([-1./2*c*s, 1./2*(1 + c*c), np.sqrt(3)/2*s]) # (3, N) + kn1= np.dot(k, n1) #(N,) + + + n2 = 1./4. * np.array([c*s - np.sqrt(3)*(1 + s*s), np.sqrt(3)*c*s - (1 + c*c), -np.sqrt(3)*s - 3*c]) + kn2= np.dot(k, n2) + + + n3 = 1./4*np.array([c*s + np.sqrt(3)*(1 + s*s), -np.sqrt(3)*c*s - (1 + c*c), -np.sqrt(3)*s + 3*c]) + kn3= np.dot(k, n3) + + + kp1Lp2L = np.dot(k, (p1L+p2L)) + kp2Lp3L = np.dot(k, (p2L+p3L)) + kp3Lp1L = np.dot(k, (p3L+p1L)) + kp0 = np.dot(k, p0) + + factorcexp0 = np.exp(1j*2.*np.pi*f/C_SI * kp0) + prefactor = np.pi*f*L/C_SI + + factorcexp12 = np.exp(1j*prefactor * (1.+kp1Lp2L/L)) + factorcexp23 = np.exp(1j*prefactor * (1.+kp2Lp3L/L)) + factorcexp31 = np.exp(1j*prefactor * (1.+kp3Lp1L/L)) + + factorsinc12 = np.sinc( (prefactor * (1.-kn3))/np.pi) + factorsinc21 = np.sinc( (prefactor * (1.+kn3))/np.pi) + factorsinc23 = np.sinc( (prefactor * (1.-kn1))/np.pi) + factorsinc32 = np.sinc( (prefactor * (1.+kn1))/np.pi) + factorsinc31 = np.sinc( (prefactor * (1.-kn2))/np.pi) + factorsinc13 = np.sinc( (prefactor * (1.+kn2))/np.pi) + + commonfac = 1j*prefactor*factorcexp0 + G12_term = commonfac * factorcexp12 * np.exp(-1j*phaseRdelay) + G23_term = commonfac * factorcexp23 * np.exp(-1j*phaseRdelay) + G31_term = commonfac * factorcexp31 * np.exp(-1j*phaseRdelay) + + G12xx = G12_term * n3[0,:]*n3[0,:] * factorsinc12 + G21xx = G12_term * n3[0,:]*n3[0,:] * factorsinc21 + G23xx = G23_term * n1[0,:]*n1[0,:] * factorsinc23 + G32xx = G23_term * n1[0,:]*n1[0,:] * factorsinc32 + G31xx = G31_term * n2[0,:]*n2[0,:] * factorsinc31 + G13xx = G31_term * n2[0,:]*n2[0,:] * factorsinc13 + + G12xy = G12_term * n3[0,:]*n3[1,:] * factorsinc12 + G21xy = G12_term * n3[0,:]*n3[1,:] * factorsinc21 + G23xy = G23_term * n1[0,:]*n1[1,:] * factorsinc23 + G32xy = G23_term * n1[0,:]*n1[1,:] * factorsinc32 + G31xy = G31_term * n2[0,:]*n2[1,:] * factorsinc31 + G13xy = G31_term * n2[0,:]*n2[1,:] * factorsinc13 + + G12xz = G12_term * n3[0,:]*n3[2,:] * factorsinc12 + G21xz = G12_term * n3[0,:]*n3[2,:] * factorsinc21 + G23xz = G23_term * n1[0,:]*n1[2,:] * factorsinc23 + G32xz = G23_term * n1[0,:]*n1[2,:] * factorsinc32 + G31xz = G31_term * n2[0,:]*n2[2,:] * factorsinc31 + G13xz = G31_term * n2[0,:]*n2[2,:] * factorsinc13 + + G12yy = G12_term * n3[1,:]*n3[1,:] * factorsinc12 + G21yy = G12_term * n3[1,:]*n3[1,:] * factorsinc21 + G23yy = G23_term * n1[1,:]*n1[1,:] * factorsinc23 + G32yy = G23_term * n1[1,:]*n1[1,:] * factorsinc32 + G31yy = G31_term * n2[1,:]*n2[1,:] * factorsinc31 + G13yy = G31_term * n2[1,:]*n2[1,:] * factorsinc13 + + G12yz = G12_term * n3[1,:]*n3[2,:] * factorsinc12 + G21yz = G12_term * n3[1,:]*n3[2,:] * factorsinc21 + G23yz = G23_term * n1[1,:]*n1[2,:] * factorsinc23 + G32yz = G23_term * n1[1,:]*n1[2,:] * factorsinc32 + G31yz = G31_term * n2[1,:]*n2[2,:] * factorsinc31 + G13yz = G31_term * n2[1,:]*n2[2,:] * factorsinc13 + + G12zz = G12_term * n3[2,:]*n3[2,:] * factorsinc12 + G21zz = G12_term * n3[2,:]*n3[2,:] * factorsinc21 + G23zz = G23_term * n1[2,:]*n1[2,:] * factorsinc23 + G32zz = G23_term * n1[2,:]*n1[2,:] * factorsinc32 + G31zz = G31_term * n2[2,:]*n2[2,:] * factorsinc31 + G13zz = G31_term * n2[2,:]*n2[2,:] * factorsinc13 + + x = np.pi*f*L/C_SI + z = np.exp(1j*2.*x) + + factor_convention = 2 + #factor_convention = 2/np.sqrt(2) # for radler dataset + + factorAE = 1j*np.sqrt(2)*np.sin(2.*x)*z + factorT = 2.*np.sqrt(2)*np.sin(2.*x)*np.sin(x)*np.exp(1j*3.*x) + + Araw_xx = 0.5 * ( (1.+z)*(G31xx + G13xx) - G23xx - z*G32xx - G21xx - z*G12xx ) + Araw_xy = 0.5 * ( (1.+z)*(G31xy + G13xy) - G23xy - z*G32xy - G21xy - z*G12xy ) + Araw_xz = 0.5 * ( (1.+z)*(G31xz + G13xz) - G23xz - z*G32xz - G21xz - z*G12xz ) + Araw_yy = 0.5 * ( (1.+z)*(G31yy + G13yy) - G23yy - z*G32yy - G21yy - z*G12yy ) + Araw_yz = 0.5 * ( (1.+z)*(G31yz + G13yz) - G23yz - z*G32yz - G21yz - z*G12yz ) + Araw_zz = 0.5 * ( (1.+z)*(G31zz + G13zz) - G23zz - z*G32zz - G21zz - z*G12zz ) + + Eraw_xx = 0.5*1/np.sqrt(3) * ( (1.-z)*(G13xx - G31xx) + (2.+z)*(G12xx - G32xx) + (1.+2.*z)*(G21xx - G23xx) ) + Eraw_xy = 0.5*1/np.sqrt(3) * ( (1.-z)*(G13xy - G31xy) + (2.+z)*(G12xy - G32xy) + (1.+2.*z)*(G21xy - G23xy) ) + Eraw_xz = 0.5*1/np.sqrt(3) * ( (1.-z)*(G13xz - G31xz) + (2.+z)*(G12xz - G32xz) + (1.+2.*z)*(G21xz - G23xz) ) + Eraw_yy = 0.5*1/np.sqrt(3) * ( (1.-z)*(G13yy - G31yy) + (2.+z)*(G12yy - G32yy) + (1.+2.*z)*(G21yy - G23yy) ) + Eraw_yz = 0.5*1/np.sqrt(3) * ( (1.-z)*(G13yz - G31yz) + (2.+z)*(G12yz - G32yz) + (1.+2.*z)*(G21yz - G23yz) ) + Eraw_zz = 0.5*1/np.sqrt(3) * ( (1.-z)*(G13zz - G31zz) + (2.+z)*(G12zz - G32zz) + (1.+2.*z)*(G21zz - G23zz) ) + + Traw_xx = 1/np.sqrt(6) * (G21xx - G12xx + G32xx - G23xx + G13xx - G31xx) + Traw_xy = 1/np.sqrt(6) * (G21xy - G12xy + G32xy - G23xy + G13xy - G31xy) + Traw_xz = 1/np.sqrt(6) * (G21xz - G12xz + G32xz - G23xz + G13xz - G31xz) + Traw_yy = 1/np.sqrt(6) * (G21yy - G12yy + G32yy - G23yy + G13yy - G31yy) + Traw_yz = 1/np.sqrt(6) * (G21yz - G12yz + G32yz - G23yz + G13yz - G31yz) + Traw_zz = 1/np.sqrt(6) * (G21zz - G12zz + G32zz - G23zz + G13zz - G31zz) + + AE_term = factor_convention * factorAE + transferL1_xx = AE_term * Araw_xx + transferL1_xy = AE_term * Araw_xy + transferL1_xz = AE_term * Araw_xz + transferL1_yy = AE_term * Araw_yy + transferL1_yz = AE_term * Araw_yz + transferL1_zz = AE_term * Araw_zz + + transferL2_xx = AE_term * Eraw_xx + transferL2_xy = AE_term * Eraw_xy + transferL2_xz = AE_term * Eraw_xz + transferL2_yy = AE_term * Eraw_yy + transferL2_yz = AE_term * Eraw_yz + transferL2_zz = AE_term * Eraw_zz + + T_term = factor_convention * factorT + transferL3_xx = T_term * Traw_xx + transferL3_xy = T_term * Traw_xy + transferL3_xz = T_term * Traw_xz + transferL3_yy = T_term * Traw_yy + transferL3_yz = T_term * Traw_yz + transferL3_zz = T_term * Traw_zz + + return np.array([transferL1_xx, transferL1_xy, transferL1_xz, transferL1_yy, transferL1_yz, transferL1_zz]), np.array([transferL2_xx, transferL2_xy, transferL2_xz, transferL2_yy, transferL2_yz, transferL2_zz]), np.array([transferL3_xx, transferL3_xy, transferL3_xz, transferL3_yy, transferL3_yz, transferL3_zz]) + +########################################################################################### +# FOR INJECTIONS +########################################################################################### +def create_lisa_injections(hlmf, fmax, fref, beta, lamda, psi, inclination, phi_ref, tref, return_response = False): + print(f"create_lisa_injections function has been called with following arguments:\n{locals()}") + tf_dict, f_dict, amp_dict, phase_dict = get_tf_from_phase_dict(hlmf, fmax, fref) + A = 0.0 + E = 0.0 + T = 0.0 + modes = list(hlmf.keys()) + response = {} + mode_TDI = {} + for mode in modes: + H_0 = transformed_Hplus_Hcross(beta, lamda, psi, inclination, -phi_ref, mode[0], mode[1]) + L1, L2, L3 = Evaluate_Gslr(tf_dict[mode] + tref, f_dict[mode], H_0, beta, lamda) + time_shifted_phase = phase_dict[mode] + 2*np.pi*tref*f_dict[mode] + tmp_data = amp_dict[mode] * np.exp(1j*time_shifted_phase) + # I belive BBHx conjugates because the formalism is define for A*exp(-1jphase), but I need to check with ROS and Mike Katz. + A += np.conj(tmp_data * L1) + E += np.conj(tmp_data * L2) + T += np.conj(tmp_data * L3) + response[mode], mode_TDI[mode] = {}, {} + response[mode]["L1"], response[mode]["L2"], response[mode]["L3"] = np.conj(L1), np.conj(L2), np.conj(L3) + mode_TDI[mode]["L1"], mode_TDI[mode]["L2"], mode_TDI[mode]["L3"] = np.conj(tmp_data*L1), np.conj(tmp_data*L2), np.conj(tmp_data*L3) + A_lal, E_lal, T_lal = create_lal_frequency_series(f_dict[modes[0]], A, hlmf[modes[0]].deltaF), create_lal_frequency_series(f_dict[modes[0]], E, hlmf[modes[0]].deltaF), create_lal_frequency_series(f_dict[modes[0]], T, hlmf[modes[0]].deltaF) + data_dict = {} + data_dict["A"], data_dict["E"], data_dict["T"] = A_lal, E_lal, T_lal + if return_response: + return data_dict, response, mode_TDI + else: + return data_dict + +def generate_lisa_TDI(P_inj, lmax=4, modes=None, tref=0.0, fref=None, return_response=False, path_to_NR_hdf5=None): + print(f"generate_lisa_TDI function has been called with following arguments:\n{locals()}") + P = lsu.ChooseWaveformParams() + + P.m1 = P_inj.m1 + P.m2 = P_inj.m2 + P.s1z = P_inj.s1z + P.s2z = P_inj.s2z + P.dist = P_inj.dist + P.fmin = P_inj.fmin + P.fmax = 0.5/P_inj.deltaT + P.deltaF = P_inj.deltaF + P.deltaT = P_inj.deltaT + + + P.phiref = 0.0 + P.inclination = 0.0 + P.psi = 0.0 + P.fref = P_inj.fref + P.tref = 0.0 + + P.approx = P_inj.approx + hlmf = lisa_lalsimutils_compat.hlmoff_for_LISA( + P, Lmax=lmax, modes=modes, path_to_NR_hdf5=path_to_NR_hdf5 + ) + modes = list(hlmf.keys()) + + # create TDI + output = create_lisa_injections(hlmf, P.fmax, fref, P_inj.theta, P_inj.phi, P_inj.psi, P_inj.incl, P_inj.phiref, tref, return_response) + + if return_response: + return output[0], output[1], output[2] + else: + return output + + +def create_h5_files_from_data_dict(data_dict, save_path): + """This function takes in data dictionary and creates h5 files from them. Assumes the data is stores as COMPLEX16FrequencySeries. + Args: + data_dict (dictonary): contains data for A, E, T channels, + save_path (string): path to where you want to save the h5 files. + Output: + None""" + A_h5_file = h5py.File(f'{save_path}/A-fake_strain-1000000-10000.h5', 'w') + A_h5_file.create_dataset('data', data=data_dict["A"].data.data) + A_h5_file.attrs["deltaF"], A_h5_file.attrs["epoch"], A_h5_file.attrs["length"], A_h5_file.attrs["f0"] = data_dict["A"].deltaF, float(data_dict["A"].epoch), data_dict["A"].data.length, data_dict["A"].f0 + A_h5_file.close() + + E_h5_file = h5py.File(f'{save_path}/E-fake_strain-1000000-10000.h5', 'w') + E_h5_file.create_dataset('data', data=data_dict["E"].data.data) + E_h5_file.attrs["deltaF"], E_h5_file.attrs["epoch"], E_h5_file.attrs["length"], E_h5_file.attrs["f0"] = data_dict["E"].deltaF, float(data_dict["E"].epoch), data_dict["E"].data.length, data_dict["E"].f0 + E_h5_file.close() + + T_h5_file = h5py.File(f'{save_path}/T-fake_strain-1000000-10000.h5', 'w') + T_h5_file.create_dataset('data', data=data_dict["T"].data.data) + T_h5_file.attrs["deltaF"], T_h5_file.attrs["epoch"], T_h5_file.attrs["length"], T_h5_file.attrs["f0"] = data_dict["T"].deltaF, float(data_dict["T"].epoch), data_dict["T"].data.length, data_dict["T"].f0 + T_h5_file.close() + + return None diff --git a/MonteCarloMarginalizeCode/Code/RIFT/LISA/response/__init__.py b/MonteCarloMarginalizeCode/Code/RIFT/LISA/response/__init__.py new file mode 100644 index 000000000..5a7346780 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/LISA/response/__init__.py @@ -0,0 +1,2 @@ +"""Time-dependent LISA response utilities.""" + diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood_LISA.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood_LISA.py new file mode 100644 index 000000000..57e3a880b --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood_LISA.py @@ -0,0 +1,431 @@ +from __future__ import print_function +import sys +from RIFT.LISA.response.LISA_response import * +import lal +import lalsimulation as lalsim +import RIFT.lalsimutils as lsu # problem of relative comprehensive import - dangerous due to package name +import numpy as np + +__author__ = "A. Jan" + +try: + import cupy + from . import optimized_gpu_tools + from . import Q_inner_product + xpy_default=cupy + junk_to_check_installed = cupy.array(5) # this will fail if GPU not installed correctly +except: + print(' no cupy (factored)') + cupy=np #import numpy as cupy # make sure pointer is identical + optimized_gpu_tools=None + Q_inner_product=None + xpy_default=np + +from .SphericalHarmonics_gpu import SphericalHarmonicsVectorized + + +from scipy import integrate + +import os +if 'PROFILE' not in os.environ: + def profile(fn): + return fn + +has_GWS=False # make sure defined in top-level scope +try: + import RIFT.physics.GWSignal as rgws + has_GWS=True +except: + has_GWS=False + +if not( 'RIFT_LOWLATENCY' in os.environ): + # Dont support external packages in low latency + try: + import NRWaveformCatalogManager3 as nrwf + useNR =True + print(" factored_likelihood.py : NRWaveformCatalogManager3 available ") + except ImportError: + useNR=False + + + try: + import RIFT.physics.ROMWaveformManager as romwf + print(" factored_likelihood.py: ROMWaveformManager as romwf") + useROM=True + rom_basis_scale = 1.0*1e-21 # Fundamental problem: Inner products with ROM basis vectors/Sh are tiny. Need to rescale to avoid overflow/underflow and simplify comparisons + except ImportError: + useROM=False + print(" factored_likelihood.py: - no ROM - ") + rom_basis_scale =1 + + try: + hasEOB=True + import RIFT.physics.EOBTidalExternalC as eobwf + # import EOBTidalExternal as eobwf + except: + hasEOB=False + print(" factored_likelihood: no EOB ") +else: + hasEOB=False + useROM=False; rom_basis_scale=1 + useNR=False + +distGpcRef = 200 # a fiducial distance for the template source. +tWindowExplore = [-0.15, 0.15] # Not used in main code. Provided for backward compatibility for ROS. Should be consistent with t_ref_wind in ILE. +rosDebugMessages = True +rosDebugMessagesDictionary = {} # Mutable after import (passed by reference). Not clear if it can be used by caling routines + # BUT if every module has a `PopulateMessagesDictionary' module, I can set their internal copies +rosDebugMessagesDictionary["DebugMessages"] = False +rosDebugMessagesDictionary["DebugMessagesLong"] = False + + + +########################################################################################### +# Main functions +########################################################################################### +def ComputeIPTimeSeries(IP, hf, data, N_shift, N_window, analyticPSD_Q=False, + inv_spec_trunc_Q=False, T_spec=0., debug=True): + r""" + Compute the complex-valued overlap between + each member of a SphHarmFrequencySeries 'hlms' + and the interferometer data COMPLEX16FrequencySeries 'data', + weighted the power spectral density REAL8FrequencySeries 'psd'. + + The integrand is non-zero in the range: [-fNyq, -fmin] union [fmin, fNyq]. + This integrand is then inverse-FFT'd to get the inner product + at a discrete series of time shifts. + + Returns a SphHarmTimeSeries object containing the complex inner product + for discrete values of the reference time tref. The epoch of the + SphHarmTimeSeries object is set to account for the transformation + """ + assert data.deltaF == hf.deltaF, f"DeltaF of template is{data.deltaF}, deltaF of data is{hf.delta}. They should be the same, this could potentionally be caused by numerical precision." + assert data.data.length == hf.data.length + + rho, rhoTS, rhoIdx, rhoPhase = IP.ip(hf, data) + rhoTS.epoch = data.epoch - hf.epoch + tmp= lsu.DataRollBins(rhoTS, N_shift) # restore functionality for bidirectional shifts: waveform need not start at t=0 + rho_time_series =lal.CutCOMPLEX16TimeSeries(rhoTS, 0, N_window) + if debug: + print(f"Max in the original series = {np.max(rhoTS.data.data)}, max in the truncated series = {np.max(rho_time_series.data.data)}, max index in the original series = {np.argmax(rhoTS.data.data) + N_shift}, max occurs at time shift of {rhoIdx * rhoTS.deltaT}s") + return rho_time_series + +def PrecomputeAlignedSpinLISA(tref, fref, t_window, hlms, hlms_conj, data_dict, psd_dict, flow, fNyq, fhigh, deltaT, beta, lamda, analyticPSD_Q=False, inv_spec_trunc_Q=False, T_spec=0.): + print(f"PrecomputeAlignedSpinLISA has been called with the following arguments: \n{locals()}") + # GENERATE DETECTOR RESPONSE + # Compute time truncation (this assumes no existing time shifts, so don't inlcude them) + # Compute Qlm (order of entires in the IP, conj, t_ref) + # Compute Ulm (sensitive to t_ref being used) + + # TO precompute innerproducts, we need tref for time shift, t_window to select what time values we need + # hlms, hlms_conj, data_dict, psd_dict, flow, fhigh, deltaT (integration range) + # beta, lamdba for detector response + + # Let's work in indexes for time shifts + # LISA response is time dependent, large shifts in time won't be good so your t_ref needs to be right. Assuming that, the peak in time series should be near the first index + + N_shift = int(t_window/deltaT) + N_window = int(2 * t_window/deltaT) + + # first get 6 terms per mode, and multiply with detector response + tf_dict, f_dict, amp_dict, phase_dict = get_tf_from_phase_dict(hlms, fNyq, fref) #here we need fmax, but RIFT has fhigh + index_at_fref = get_closest_index(f_dict[2,2], fref) + + collect_mode_terms = {} + modes = list(hlms.keys()) + + reference_phase = 0.0 + modes.remove((2,2)) + modes.insert(0, (2,2)) + + collect_mode_terms["A"] = {} + collect_mode_terms["E"] = {} + collect_mode_terms["T"] = {} + for mode in modes: + A_terms, E_terms, T_terms = Evaluate_Gslr_test_2(tf_dict[mode]+tref, f_dict[mode], beta, lamda) # NOTE: I added t_ref + + amp, phase = amp_dict[mode], phase_dict[mode] + shifted_phase = (phase + 2*np.pi*f_dict[mode]*tref) #take care of convention + #if mode == (2,2): + # phase_22_current = shifted_phase[index_at_fref] + # difference = reference_phase - phase_22_current + #shifted_phase = shifted_phase + mode[1]/2 * difference + #print(f"Precompute: {mode}, phase = {shifted_phase[index_at_fref]}, time = {tf_dict[mode][index_at_fref]+tref}") + + tmp_mode_data = (amp * np.exp(1j*shifted_phase)).reshape(1, -1) #take care of convention + + # tmp_mode_data = hlms[mode].data.data + collect_mode_terms["A"][mode] = {} + tmp_mode_here = np.conj(A_terms * tmp_mode_data) + collect_mode_terms["A"][mode][0] = create_lal_frequency_series(f_dict[mode], tmp_mode_here[0], hlms[mode].deltaF) # A_term (6,n), hlm (n,1) + collect_mode_terms["A"][mode][1] = create_lal_frequency_series(f_dict[mode], tmp_mode_here[1], hlms[mode].deltaF) + collect_mode_terms["A"][mode][2] = create_lal_frequency_series(f_dict[mode], tmp_mode_here[2], hlms[mode].deltaF) + collect_mode_terms["A"][mode][3] = create_lal_frequency_series(f_dict[mode], tmp_mode_here[3], hlms[mode].deltaF) + collect_mode_terms["A"][mode][4] = create_lal_frequency_series(f_dict[mode], tmp_mode_here[4], hlms[mode].deltaF) + collect_mode_terms["A"][mode][5] = create_lal_frequency_series(f_dict[mode], tmp_mode_here[5], hlms[mode].deltaF) + + + collect_mode_terms["E"][mode] = {} + tmp_mode_here = np.conj(E_terms * tmp_mode_data) + collect_mode_terms["E"][mode][0] = create_lal_frequency_series(f_dict[mode], tmp_mode_here[0], hlms[mode].deltaF) # A_term (6,n), hlm (n,1) + collect_mode_terms["E"][mode][1] = create_lal_frequency_series(f_dict[mode], tmp_mode_here[1], hlms[mode].deltaF) + collect_mode_terms["E"][mode][2] = create_lal_frequency_series(f_dict[mode], tmp_mode_here[2], hlms[mode].deltaF) + collect_mode_terms["E"][mode][3] = create_lal_frequency_series(f_dict[mode], tmp_mode_here[3], hlms[mode].deltaF) + collect_mode_terms["E"][mode][4] = create_lal_frequency_series(f_dict[mode], tmp_mode_here[4], hlms[mode].deltaF) + collect_mode_terms["E"][mode][5] = create_lal_frequency_series(f_dict[mode], tmp_mode_here[5], hlms[mode].deltaF) + + + collect_mode_terms["T"][mode] = {} + tmp_mode_here = np.conj(T_terms * tmp_mode_data) + collect_mode_terms["T"][mode][0] = create_lal_frequency_series(f_dict[mode], tmp_mode_here[0], hlms[mode].deltaF) # A_term (6,n), hlm (n,1) + collect_mode_terms["T"][mode][1] = create_lal_frequency_series(f_dict[mode], tmp_mode_here[1], hlms[mode].deltaF) + collect_mode_terms["T"][mode][2] = create_lal_frequency_series(f_dict[mode], tmp_mode_here[2], hlms[mode].deltaF) + collect_mode_terms["T"][mode][3] = create_lal_frequency_series(f_dict[mode], tmp_mode_here[3], hlms[mode].deltaF) + collect_mode_terms["T"][mode][4] = create_lal_frequency_series(f_dict[mode], tmp_mode_here[4], hlms[mode].deltaF) + collect_mode_terms["T"][mode][5] = create_lal_frequency_series(f_dict[mode], tmp_mode_here[5], hlms[mode].deltaF) + # calculate + IP_time_series_A= lsu.ComplexOverlap(flow, fhigh, fNyq, hlms[2,2].deltaF, psd_dict["A"], analyticPSD_Q, inv_spec_trunc_Q, T_spec, full_output =True) # Incase the three arms have different PSDs. Assume all arms have same PSD for now and 2,2 mode is present + IP_time_series_E= lsu.ComplexOverlap(flow, fhigh, fNyq, hlms[2,2].deltaF, psd_dict["E"], analyticPSD_Q, inv_spec_trunc_Q, T_spec, full_output =True) + IP_time_series_T= lsu.ComplexOverlap(flow, fhigh, fNyq, hlms[2,2].deltaF, psd_dict["T"], analyticPSD_Q, inv_spec_trunc_Q, T_spec, full_output =True) + + Q_lm = {} + for channel in ["A", "E", "T"]: + if channel == "A": + inner_product= IP_time_series_A + if channel == "E": + inner_product= IP_time_series_E + if channel == "T": + inner_product= IP_time_series_T + for mode in hlms.keys(): + l, m = mode[0], mode[1] + + Q_lm[f"{channel}_{l}_{m}_xx"] = ((ComputeIPTimeSeries(inner_product, data_dict[channel], collect_mode_terms[channel][mode][0], -N_shift, N_window))) + Q_lm[f"{channel}_{l}_{m}_xy"] = ((ComputeIPTimeSeries(inner_product, data_dict[channel], collect_mode_terms[channel][mode][1], -N_shift, N_window))) + Q_lm[f"{channel}_{l}_{m}_xz"] = ((ComputeIPTimeSeries(inner_product, data_dict[channel], collect_mode_terms[channel][mode][2], -N_shift, N_window))) + Q_lm[f"{channel}_{l}_{m}_yy"] = ((ComputeIPTimeSeries(inner_product, data_dict[channel], collect_mode_terms[channel][mode][3], -N_shift, N_window))) + Q_lm[f"{channel}_{l}_{m}_yz"] = ((ComputeIPTimeSeries(inner_product, data_dict[channel], collect_mode_terms[channel][mode][4], -N_shift, N_window))) + Q_lm[f"{channel}_{l}_{m}_zz"] = ((ComputeIPTimeSeries(inner_product, data_dict[channel], collect_mode_terms[channel][mode][5], -N_shift, N_window))) + + # calculate + U_lm_pq = {} + IP_A = lsu.ComplexIP(flow, fhigh, fNyq, hlms[2,2].deltaF, psd_dict["A"], analyticPSD_Q, inv_spec_trunc_Q, T_spec) # Incase the three arms have different PSDs. Assume all arms have same PSD for now and 2,2 mode is present + IP_E = lsu.ComplexIP(flow, fhigh, fNyq, hlms[2,2].deltaF, psd_dict["E"], analyticPSD_Q, inv_spec_trunc_Q, T_spec) + IP_T = lsu.ComplexIP(flow, fhigh, fNyq, hlms[2,2].deltaF, psd_dict["T"], analyticPSD_Q, inv_spec_trunc_Q, T_spec) + for i in np.arange(len(modes)): + for j in np.arange(len(modes))[i:]: + # print(modes[i], modes[j]) + l, m, p, q = modes[i][0], modes[i][1], modes[j][0], modes[j][1] + for channel in ["A", "E", "T"]: + if channel == "A": + inner_product= IP_A + if channel == "E": + inner_product= IP_E + if channel == "T": + inner_product= IP_T + U_lm_pq[f"{channel}_{l}_{m}_xx_{p}_{q}_xx"] = inner_product.ip(collect_mode_terms[channel][modes[i]][0], collect_mode_terms[channel][modes[j]][0]) + U_lm_pq[f"{channel}_{l}_{m}_xx_{p}_{q}_xy"] = inner_product.ip(collect_mode_terms[channel][modes[i]][0], collect_mode_terms[channel][modes[j]][1]) + U_lm_pq[f"{channel}_{l}_{m}_xx_{p}_{q}_xz"] = inner_product.ip(collect_mode_terms[channel][modes[i]][0], collect_mode_terms[channel][modes[j]][2]) + U_lm_pq[f"{channel}_{l}_{m}_xx_{p}_{q}_yy"] = inner_product.ip(collect_mode_terms[channel][modes[i]][0], collect_mode_terms[channel][modes[j]][3]) + U_lm_pq[f"{channel}_{l}_{m}_xx_{p}_{q}_yz"] = inner_product.ip(collect_mode_terms[channel][modes[i]][0], collect_mode_terms[channel][modes[j]][4]) + U_lm_pq[f"{channel}_{l}_{m}_xx_{p}_{q}_zz"] = inner_product.ip(collect_mode_terms[channel][modes[i]][0], collect_mode_terms[channel][modes[j]][5]) + + U_lm_pq[f"{channel}_{l}_{m}_xy_{p}_{q}_xx"] = inner_product.ip(collect_mode_terms[channel][modes[i]][1], collect_mode_terms[channel][modes[j]][0]) + U_lm_pq[f"{channel}_{l}_{m}_xy_{p}_{q}_xy"] = inner_product.ip(collect_mode_terms[channel][modes[i]][1], collect_mode_terms[channel][modes[j]][1]) + U_lm_pq[f"{channel}_{l}_{m}_xy_{p}_{q}_xz"] = inner_product.ip(collect_mode_terms[channel][modes[i]][1], collect_mode_terms[channel][modes[j]][2]) + U_lm_pq[f"{channel}_{l}_{m}_xy_{p}_{q}_yy"] = inner_product.ip(collect_mode_terms[channel][modes[i]][1], collect_mode_terms[channel][modes[j]][3]) + U_lm_pq[f"{channel}_{l}_{m}_xy_{p}_{q}_yz"] = inner_product.ip(collect_mode_terms[channel][modes[i]][1], collect_mode_terms[channel][modes[j]][4]) + U_lm_pq[f"{channel}_{l}_{m}_xy_{p}_{q}_zz"] = inner_product.ip(collect_mode_terms[channel][modes[i]][1], collect_mode_terms[channel][modes[j]][5]) + + U_lm_pq[f"{channel}_{l}_{m}_xz_{p}_{q}_xx"] = inner_product.ip(collect_mode_terms[channel][modes[i]][2], collect_mode_terms[channel][modes[j]][0]) + U_lm_pq[f"{channel}_{l}_{m}_xz_{p}_{q}_xy"] = inner_product.ip(collect_mode_terms[channel][modes[i]][2], collect_mode_terms[channel][modes[j]][1]) + U_lm_pq[f"{channel}_{l}_{m}_xz_{p}_{q}_xz"] = inner_product.ip(collect_mode_terms[channel][modes[i]][2], collect_mode_terms[channel][modes[j]][2]) + U_lm_pq[f"{channel}_{l}_{m}_xz_{p}_{q}_yy"] = inner_product.ip(collect_mode_terms[channel][modes[i]][2], collect_mode_terms[channel][modes[j]][3]) + U_lm_pq[f"{channel}_{l}_{m}_xz_{p}_{q}_yz"] = inner_product.ip(collect_mode_terms[channel][modes[i]][2], collect_mode_terms[channel][modes[j]][4]) + U_lm_pq[f"{channel}_{l}_{m}_xz_{p}_{q}_zz"] = inner_product.ip(collect_mode_terms[channel][modes[i]][2], collect_mode_terms[channel][modes[j]][5]) + + U_lm_pq[f"{channel}_{l}_{m}_yy_{p}_{q}_xx"] = inner_product.ip(collect_mode_terms[channel][modes[i]][3], collect_mode_terms[channel][modes[j]][0]) + U_lm_pq[f"{channel}_{l}_{m}_yy_{p}_{q}_xy"] = inner_product.ip(collect_mode_terms[channel][modes[i]][3], collect_mode_terms[channel][modes[j]][1]) + U_lm_pq[f"{channel}_{l}_{m}_yy_{p}_{q}_xz"] = inner_product.ip(collect_mode_terms[channel][modes[i]][3], collect_mode_terms[channel][modes[j]][2]) + U_lm_pq[f"{channel}_{l}_{m}_yy_{p}_{q}_yy"] = inner_product.ip(collect_mode_terms[channel][modes[i]][3], collect_mode_terms[channel][modes[j]][3]) + U_lm_pq[f"{channel}_{l}_{m}_yy_{p}_{q}_yz"] = inner_product.ip(collect_mode_terms[channel][modes[i]][3], collect_mode_terms[channel][modes[j]][4]) + U_lm_pq[f"{channel}_{l}_{m}_yy_{p}_{q}_zz"] = inner_product.ip(collect_mode_terms[channel][modes[i]][3], collect_mode_terms[channel][modes[j]][5]) + + U_lm_pq[f"{channel}_{l}_{m}_yz_{p}_{q}_xx"] = inner_product.ip(collect_mode_terms[channel][modes[i]][4], collect_mode_terms[channel][modes[j]][0]) + U_lm_pq[f"{channel}_{l}_{m}_yz_{p}_{q}_xy"] = inner_product.ip(collect_mode_terms[channel][modes[i]][4], collect_mode_terms[channel][modes[j]][1]) + U_lm_pq[f"{channel}_{l}_{m}_yz_{p}_{q}_xz"] = inner_product.ip(collect_mode_terms[channel][modes[i]][4], collect_mode_terms[channel][modes[j]][2]) + U_lm_pq[f"{channel}_{l}_{m}_yz_{p}_{q}_yy"] = inner_product.ip(collect_mode_terms[channel][modes[i]][4], collect_mode_terms[channel][modes[j]][3]) + U_lm_pq[f"{channel}_{l}_{m}_yz_{p}_{q}_yz"] = inner_product.ip(collect_mode_terms[channel][modes[i]][4], collect_mode_terms[channel][modes[j]][4]) + U_lm_pq[f"{channel}_{l}_{m}_yz_{p}_{q}_zz"] = inner_product.ip(collect_mode_terms[channel][modes[i]][4], collect_mode_terms[channel][modes[j]][5]) + + U_lm_pq[f"{channel}_{l}_{m}_zz_{p}_{q}_xx"] = inner_product.ip(collect_mode_terms[channel][modes[i]][5], collect_mode_terms[channel][modes[j]][0]) + U_lm_pq[f"{channel}_{l}_{m}_zz_{p}_{q}_xy"] = inner_product.ip(collect_mode_terms[channel][modes[i]][5], collect_mode_terms[channel][modes[j]][1]) + U_lm_pq[f"{channel}_{l}_{m}_zz_{p}_{q}_xz"] = inner_product.ip(collect_mode_terms[channel][modes[i]][5], collect_mode_terms[channel][modes[j]][2]) + U_lm_pq[f"{channel}_{l}_{m}_zz_{p}_{q}_yy"] = inner_product.ip(collect_mode_terms[channel][modes[i]][5], collect_mode_terms[channel][modes[j]][3]) + U_lm_pq[f"{channel}_{l}_{m}_zz_{p}_{q}_yz"] = inner_product.ip(collect_mode_terms[channel][modes[i]][5], collect_mode_terms[channel][modes[j]][4]) + U_lm_pq[f"{channel}_{l}_{m}_zz_{p}_{q}_zz"] = inner_product.ip(collect_mode_terms[channel][modes[i]][5], collect_mode_terms[channel][modes[j]][5]) + + if (l,m) != (p,q): + U_lm_pq[f"{channel}_{p}_{q}_xx_{l}_{m}_xx"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_xx_{p}_{q}_xx"]) + U_lm_pq[f"{channel}_{p}_{q}_xx_{l}_{m}_xy"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_xy_{p}_{q}_xx"]) + U_lm_pq[f"{channel}_{p}_{q}_xx_{l}_{m}_xz"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_xz_{p}_{q}_xx"]) + U_lm_pq[f"{channel}_{p}_{q}_xx_{l}_{m}_yy"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_yy_{p}_{q}_xx"]) + U_lm_pq[f"{channel}_{p}_{q}_xx_{l}_{m}_yz"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_yz_{p}_{q}_xx"]) + U_lm_pq[f"{channel}_{p}_{q}_xx_{l}_{m}_zz"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_zz_{p}_{q}_xx"]) + + U_lm_pq[f"{channel}_{p}_{q}_xy_{l}_{m}_xx"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_xx_{p}_{q}_xy"]) + U_lm_pq[f"{channel}_{p}_{q}_xy_{l}_{m}_xy"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_xy_{p}_{q}_xy"]) + U_lm_pq[f"{channel}_{p}_{q}_xy_{l}_{m}_xz"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_xz_{p}_{q}_xy"]) + U_lm_pq[f"{channel}_{p}_{q}_xy_{l}_{m}_yy"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_yy_{p}_{q}_xy"]) + U_lm_pq[f"{channel}_{p}_{q}_xy_{l}_{m}_yz"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_yz_{p}_{q}_xy"]) + U_lm_pq[f"{channel}_{p}_{q}_xy_{l}_{m}_zz"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_zz_{p}_{q}_xy"]) + + U_lm_pq[f"{channel}_{p}_{q}_xz_{l}_{m}_xx"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_xx_{p}_{q}_xz"]) + U_lm_pq[f"{channel}_{p}_{q}_xz_{l}_{m}_xy"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_xy_{p}_{q}_xz"]) + U_lm_pq[f"{channel}_{p}_{q}_xz_{l}_{m}_xz"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_xz_{p}_{q}_xz"]) + U_lm_pq[f"{channel}_{p}_{q}_xz_{l}_{m}_yy"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_yy_{p}_{q}_xz"]) + U_lm_pq[f"{channel}_{p}_{q}_xz_{l}_{m}_yz"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_yz_{p}_{q}_xz"]) + U_lm_pq[f"{channel}_{p}_{q}_xz_{l}_{m}_zz"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_zz_{p}_{q}_xz"]) + + U_lm_pq[f"{channel}_{p}_{q}_yy_{l}_{m}_xx"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_xx_{p}_{q}_yy"]) + U_lm_pq[f"{channel}_{p}_{q}_yy_{l}_{m}_xy"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_xy_{p}_{q}_yy"]) + U_lm_pq[f"{channel}_{p}_{q}_yy_{l}_{m}_xz"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_xz_{p}_{q}_yy"]) + U_lm_pq[f"{channel}_{p}_{q}_yy_{l}_{m}_yy"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_yy_{p}_{q}_yy"]) + U_lm_pq[f"{channel}_{p}_{q}_yy_{l}_{m}_yz"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_yz_{p}_{q}_yy"]) + U_lm_pq[f"{channel}_{p}_{q}_yy_{l}_{m}_zz"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_zz_{p}_{q}_yy"]) + + U_lm_pq[f"{channel}_{p}_{q}_yz_{l}_{m}_xx"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_xx_{p}_{q}_yz"]) + U_lm_pq[f"{channel}_{p}_{q}_yz_{l}_{m}_xy"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_xy_{p}_{q}_yz"]) + U_lm_pq[f"{channel}_{p}_{q}_yz_{l}_{m}_xz"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_xz_{p}_{q}_yz"]) + U_lm_pq[f"{channel}_{p}_{q}_yz_{l}_{m}_yy"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_yy_{p}_{q}_yz"]) + U_lm_pq[f"{channel}_{p}_{q}_yz_{l}_{m}_yz"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_yz_{p}_{q}_yz"]) + U_lm_pq[f"{channel}_{p}_{q}_yz_{l}_{m}_zz"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_zz_{p}_{q}_yz"]) + + U_lm_pq[f"{channel}_{p}_{q}_zz_{l}_{m}_xx"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_xx_{p}_{q}_zz"]) + U_lm_pq[f"{channel}_{p}_{q}_zz_{l}_{m}_xy"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_xy_{p}_{q}_zz"]) + U_lm_pq[f"{channel}_{p}_{q}_zz_{l}_{m}_xz"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_xz_{p}_{q}_zz"]) + U_lm_pq[f"{channel}_{p}_{q}_zz_{l}_{m}_yy"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_yy_{p}_{q}_zz"]) + U_lm_pq[f"{channel}_{p}_{q}_zz_{l}_{m}_yz"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_yz_{p}_{q}_zz"]) + U_lm_pq[f"{channel}_{p}_{q}_zz_{l}_{m}_zz"] = np.conj(U_lm_pq[f"{channel}_{l}_{m}_zz_{p}_{q}_zz"]) + + # Match RIFT's output + guess_snr = None # right now, it is painful to evaluate this + rholms_intp = None + V_lm_pq = None + rest = None + # current RIFT output for precompute is: rholms_intp, cross_terms, cross_terms_V, rholms, guess_snr, rest + return rholms_intp, U_lm_pq, V_lm_pq, Q_lm, guess_snr, rest + + +def FactoredLogLikelihoodAlignedSpinLISA(Q_lm, U_lm_pq, beta, lam, psi, inclination, phi_ref, distance, modes, reference_distance, return_lnLt=False, only_positive_modes=True): + # Calculated marginalized likelihood, marginalized over time, psi, inlincation, phiref and distance. + # call psi terms + plus_terms = get_beta_lamda_psi_terms_Hp(beta, lam, psi) + cross_terms = get_beta_lamda_psi_terms_Hc(beta, lam, psi) + terms = 0.5*(plus_terms - 1j*cross_terms) + conj_terms = 0.5*(plus_terms + 1j*cross_terms) + + # call spherical harmonics + factor = np.ones(modes.shape) + factor[:,1] = -factor[:,1] + negative_m_modes = modes * factor + spherical_harmonics = SphericalHarmonicsVectorized(modes, inclination, -phi_ref) + negative_m_harmonics = SphericalHarmonicsVectorized(negative_m_modes, inclination, -phi_ref) + + term_lm_conj_conjterm_lm__ = {} + conjterm_lm_term_lm__conj = {} + + Qlm = {} + Qlm["A"], Qlm["E"], Qlm["T"] = 0, 0, 0 + for i, mode in enumerate(modes): + l, m = mode[0], mode[1] + + # This is used in calculating Ulmpq too, so we will save it now instead of calculating it again + Sp_harm_lm, Sp_harm_lm_ = spherical_harmonics[:,i], negative_m_harmonics[:,i] + Sp_harm_lm_conjugate, Sp_harm_lm__conjugate = np.conj((Sp_harm_lm)), np.conj(Sp_harm_lm_) + + term_lm_conj_conjterm_lm__[f"{l}_{m}"] = (terms)*Sp_harm_lm_conjugate + (-1)**(l)*(conj_terms)*(Sp_harm_lm_) + conjterm_lm_term_lm__conj[f"{l}_{m}"] = conj_terms*(Sp_harm_lm) + (-1)**(l)*(terms)*Sp_harm_lm__conjugate + + for channel in ["A", "E", "T"]: + Qlm[channel] += Q_lm[f"{channel}_{l}_{m}_xx"].data.data.reshape(-1,1) * (term_lm_conj_conjterm_lm__[f"{l}_{m}"][0]) + Qlm[channel] += Q_lm[f"{channel}_{l}_{m}_xy"].data.data.reshape(-1,1) * (term_lm_conj_conjterm_lm__[f"{l}_{m}"][1]) + Qlm[channel] += Q_lm[f"{channel}_{l}_{m}_xz"].data.data.reshape(-1,1) * (term_lm_conj_conjterm_lm__[f"{l}_{m}"][2]) + Qlm[channel] += Q_lm[f"{channel}_{l}_{m}_yy"].data.data.reshape(-1,1) * (term_lm_conj_conjterm_lm__[f"{l}_{m}"][3]) + Qlm[channel] += Q_lm[f"{channel}_{l}_{m}_yz"].data.data.reshape(-1,1) * (term_lm_conj_conjterm_lm__[f"{l}_{m}"][4]) + Qlm[channel] += Q_lm[f"{channel}_{l}_{m}_zz"].data.data.reshape(-1,1) * (term_lm_conj_conjterm_lm__[f"{l}_{m}"][5]) + # Qlm[channel] += Q_lm[f"{channel}_{l}_{m}_xx"].data.data.reshape(-1,1) * (0.5*(plus_terms[0] - 1j*cross_terms[0])*Sp_harm_lm_conjugate + 0.5*(-1)**(l)*(plus_terms[0] + 1j*cross_terms[0])*(Sp_harm_lm_)) + # Qlm[channel] += Q_lm[f"{channel}_{l}_{m}_xy"].data.data.reshape(-1,1) * (0.5*(plus_terms[1] - 1j*cross_terms[1])*Sp_harm_lm_conjugate + 0.5*(-1)**(l)*(plus_terms[1] + 1j*cross_terms[1])*(Sp_harm_lm_)) + # Qlm[channel] += Q_lm[f"{channel}_{l}_{m}_xz"].data.data.reshape(-1,1) * (0.5*(plus_terms[2] - 1j*cross_terms[2])*Sp_harm_lm_conjugate + 0.5*(-1)**(l)*(plus_terms[2] + 1j*cross_terms[2])*(Sp_harm_lm_)) + # Qlm[channel] += Q_lm[f"{channel}_{l}_{m}_yy"].data.data.reshape(-1,1) * (0.5*(plus_terms[3] - 1j*cross_terms[3])*Sp_harm_lm_conjugate + 0.5*(-1)**(l)*(plus_terms[3] + 1j*cross_terms[3])*(Sp_harm_lm_)) + # Qlm[channel] += Q_lm[f"{channel}_{l}_{m}_yz"].data.data.reshape(-1,1) * (0.5*(plus_terms[4] - 1j*cross_terms[4])*Sp_harm_lm_conjugate + 0.5*(-1)**(l)*(plus_terms[4] + 1j*cross_terms[4])*(Sp_harm_lm_)) + # Qlm[channel] += Q_lm[f"{channel}_{l}_{m}_zz"].data.data.reshape(-1,1) * (0.5*(plus_terms[5] - 1j*cross_terms[5])*Sp_harm_lm_conjugate + 0.5*(-1)**(l)*(plus_terms[5] + 1j*cross_terms[5])*(Sp_harm_lm_)) + + Ulmpq = {} + Ulmpq["A"], Ulmpq["E"], Ulmpq["T"] = 0, 0, 0 + for i in np.arange(len(modes)): + l, m = modes[i][0], modes[i][1] + for j in np.arange(len(modes))[i:]: + p, q = modes[j][0], modes[j][1] + for channel in ["A", "E", "T"]: + temp_value = 0.0 # so we can exploit conjugate symmetry! (tested) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_xx_{p}_{q}_xx"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][0]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][0]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_xx_{p}_{q}_xy"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][0]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][1]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_xx_{p}_{q}_xz"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][0]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][2]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_xx_{p}_{q}_yy"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][0]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][3]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_xx_{p}_{q}_yz"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][0]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][4]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_xx_{p}_{q}_zz"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][0]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][5]) + + temp_value += U_lm_pq[f"{channel}_{l}_{m}_xy_{p}_{q}_xx"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][1]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][0]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_xy_{p}_{q}_xy"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][1]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][1]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_xy_{p}_{q}_xz"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][1]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][2]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_xy_{p}_{q}_yy"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][1]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][3]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_xy_{p}_{q}_yz"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][1]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][4]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_xy_{p}_{q}_zz"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][1]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][5]) + + temp_value += U_lm_pq[f"{channel}_{l}_{m}_xz_{p}_{q}_xx"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][2]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][0]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_xz_{p}_{q}_xy"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][2]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][1]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_xz_{p}_{q}_xz"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][2]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][2]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_xz_{p}_{q}_yy"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][2]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][3]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_xz_{p}_{q}_yz"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][2]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][4]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_xz_{p}_{q}_zz"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][2]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][5]) + + temp_value += U_lm_pq[f"{channel}_{l}_{m}_yy_{p}_{q}_xx"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][3]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][0]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_yy_{p}_{q}_xy"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][3]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][1]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_yy_{p}_{q}_xz"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][3]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][2]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_yy_{p}_{q}_yy"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][3]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][3]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_yy_{p}_{q}_yz"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][3]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][4]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_yy_{p}_{q}_zz"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][3]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][5]) + + temp_value += U_lm_pq[f"{channel}_{l}_{m}_yz_{p}_{q}_xx"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][4]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][0]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_yz_{p}_{q}_xy"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][4]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][1]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_yz_{p}_{q}_xz"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][4]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][2]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_yz_{p}_{q}_yy"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][4]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][3]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_yz_{p}_{q}_yz"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][4]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][4]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_yz_{p}_{q}_zz"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][4]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][5]) + + temp_value += U_lm_pq[f"{channel}_{l}_{m}_zz_{p}_{q}_xx"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][5]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][0]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_zz_{p}_{q}_xy"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][5]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][1]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_zz_{p}_{q}_xz"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][5]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][2]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_zz_{p}_{q}_yy"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][5]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][3]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_zz_{p}_{q}_yz"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][5]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][4]) + temp_value += U_lm_pq[f"{channel}_{l}_{m}_zz_{p}_{q}_zz"]* (conjterm_lm_term_lm__conj[f"{l}_{m}"][5]) * (term_lm_conj_conjterm_lm__[f"{p}_{q}"][5]) + if (l,m) != (p,q): + temp_value += np.conj(temp_value) + Ulmpq[channel] += temp_value + total_lnL = 0 + for channel in ["A", "E", "T"]: + total_lnL += np.real(reference_distance/distance * Qlm[channel] - ((reference_distance/distance)**2)*0.5*Ulmpq[channel]) + # if not including -m modes, then the innerproduct is 4*integral, but all innerproducts evaluated so far are 2*integral + if only_positive_modes: + total_lnL *= 2 + # for time sampling, return likelihood time series. + if return_lnLt: + return total_lnL + # shape (time terms, extrinsic_params), integrating in time --> axis 0 + Q_lm_term = list(Q_lm.keys())[0] + L_t = np.exp(total_lnL - np.max(total_lnL, axis=0)) + L = integrate.simpson(L_t, dx = Q_lm[Q_lm_term].deltaT, axis=0) #P.deltaT + lnL = np.max(total_lnL, axis=0) + np.log(L) + return lnL + + diff --git a/MonteCarloMarginalizeCode/Code/RIFT/misc/hyperpipeline_io.py b/MonteCarloMarginalizeCode/Code/RIFT/misc/hyperpipeline_io.py index 681b40e9e..2461f8029 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/misc/hyperpipeline_io.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/misc/hyperpipeline_io.py @@ -37,6 +37,7 @@ * ``eccentricity`` (and optionally ``meanPerAno``) * ``lambda1 lambda2`` (and optionally ``eos_table_index``) * ``distance`` (in Mpc) +* ``ecliptic_longitude ecliptic_latitude`` (sky position treated as intrinsic) The reader returns a structured ``numpy.ndarray`` so callers can address columns by name regardless of order, and a small adaptor function is provided @@ -81,6 +82,8 @@ "lambda2", "eos_table_index", "distance", + "ecliptic_longitude", + "ecliptic_latitude", ) @@ -100,14 +103,14 @@ def is_active(env=None): def build_column_list(use_eccentricity=False, use_meanPerAno=False, use_tides=False, use_eos_index=False, - use_distance=False): + use_distance=False, use_sky=False): """Compose the column-name tuple for a given physics configuration. The resulting tuple always begins with :data:`DEFAULT_BASE_COLUMNS` and appends the requested optional groups in a fixed canonical order:: (eccentricity, [meanPerAno,] [lambda1, lambda2,] - [eos_table_index,] [distance]) + [eos_table_index,] [distance,] [ecliptic_longitude, ecliptic_latitude]) """ cols = list(DEFAULT_BASE_COLUMNS) if use_eccentricity: @@ -120,6 +123,8 @@ def build_column_list(use_eccentricity=False, use_meanPerAno=False, cols.append("eos_table_index") if use_distance: cols.append("distance") + if use_sky: + cols.extend(["ecliptic_longitude", "ecliptic_latitude"]) return tuple(cols) @@ -267,7 +272,8 @@ def read_table(fname): # --------------------------------------------------------------------------- def to_legacy_dat(arr, use_eccentricity=False, use_meanPerAno=False, - use_tides=False, use_eos_index=False, use_distance=False): + use_tides=False, use_eos_index=False, use_distance=False, + use_sky=False): """Reshape a hyperpipeline structured array into the legacy CIP layout. The legacy CIP reader expects a plain ``(N, K)`` ``ndarray`` whose @@ -276,6 +282,7 @@ def to_legacy_dat(arr, use_eccentricity=False, use_meanPerAno=False, [event_id, m1, m2, s1x, s1y, s1z, s2x, s2y, s2z, (distance?), (lambda1, lambda2, (eos_index)?)?, (eccentricity, (meanPerAno)?)?, + (ecliptic_longitude, ecliptic_latitude)?, lnL, sigma_lnL] The ordering of optional groups before ``lnL`` mirrors what the legacy @@ -297,6 +304,8 @@ def to_legacy_dat(arr, use_eccentricity=False, use_meanPerAno=False, cols.append("eccentricity") if use_meanPerAno: cols.append("meanPerAno") + if use_sky: + cols.extend(["ecliptic_longitude", "ecliptic_latitude"]) cols.extend(["lnL", "sigma_lnL"]) out = np.zeros((n, len(cols)), dtype=float) @@ -425,6 +434,8 @@ def consolidate(arr, columns, sigma_cut=0.9, digits=5): COLUMN_ALIAS_DISK_TO_ATTR = { "a1x": "s1x", "a1y": "s1y", "a1z": "s1z", "a2x": "s2x", "a2y": "s2y", "a2z": "s2z", + "ecliptic_longitude": "phi", + "ecliptic_latitude": "theta", } COLUMN_ALIAS_ATTR_TO_DISK = {v: k for k, v in COLUMN_ALIAS_DISK_TO_ATTR.items()} @@ -614,7 +625,7 @@ def read_grid_to_P_list(fname, P_factory, lal_module=None, def legacy_column_indices(use_eccentricity=False, use_meanPerAno=False, use_tides=False, use_eos_index=False, - use_distance=False): + use_distance=False, use_sky=False): """Return the positional column indices the legacy CIP loop expects. Mirrors the layout produced by :func:`to_legacy_dat` so callers can @@ -623,13 +634,15 @@ def legacy_column_indices(use_eccentricity=False, use_meanPerAno=False, CIP indexing logic without having to recompute them. Returns a dict keyed by ``'lnL'``, ``'sigma_lnL'``, ``'distance'``, - ``'lambda1'``, ``'eccentricity'``, ``'meanPerAno'``. Any column not + ``'lambda1'``, ``'eccentricity'``, ``'meanPerAno'``, + ``'ecliptic_longitude'``, ``'ecliptic_latitude'``. Any column not present in the configuration maps to ``None``. """ # event_id m1 m2 a1x a1y a1z a2x a2y a2z = 9 leading columns. idx = 9 out = {"lnL": None, "sigma_lnL": None, "distance": None, - "lambda1": None, "eccentricity": None, "meanPerAno": None} + "lambda1": None, "eccentricity": None, "meanPerAno": None, + "ecliptic_longitude": None, "ecliptic_latitude": None} if use_distance: out["distance"] = idx idx += 1 @@ -644,6 +657,10 @@ def legacy_column_indices(use_eccentricity=False, use_meanPerAno=False, if use_meanPerAno: out["meanPerAno"] = idx idx += 1 + if use_sky: + out["ecliptic_longitude"] = idx + out["ecliptic_latitude"] = idx + 1 + idx += 2 out["lnL"] = idx out["sigma_lnL"] = idx + 1 return out diff --git a/MonteCarloMarginalizeCode/Code/bin/helper_LISA_Events.py b/MonteCarloMarginalizeCode/Code/bin/helper_LISA_Events.py new file mode 100755 index 000000000..44afeda45 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/bin/helper_LISA_Events.py @@ -0,0 +1,269 @@ +#!/usr/bin/env python +"""Emit a minimal CEPP argument bundle for LISA-RIFT workflows. + +This helper is intentionally separate from ``helper_LDG_Events.py``. It does +not select observing-run data, skymaps, or astrophysical fit strategies; it +only writes the files consumed by ``create_event_parameter_pipeline_*`` for a +LISA-specific ILE run. +""" + +from __future__ import print_function + +import argparse +import os +import shlex + +import numpy as np + +from RIFT.misc import hyperpipeline_io + + +def _split_assignment(text, option_name): + if "=" not in text: + raise argparse.ArgumentTypeError( + "{} entries must have the form NAME=value".format(option_name) + ) + name, value = text.split("=", 1) + if not name or not value: + raise argparse.ArgumentTypeError( + "{} entries must have the form NAME=value".format(option_name) + ) + return name, value + + +def _quote_join(parts): + return " ".join(shlex.quote(str(part)) for part in parts if str(part) != "") + + +def _write_arg_file(path, parts): + with open(path, "w") as out: + out.write("X " + _quote_join(parts) + "\n") + + +def _write_cip_list(path, lines): + with open(path, "w") as out: + out.write("\n".join(lines) + "\n") + + +def _write_transfer_file(path, names): + with open(path, "w") as out: + for name in names: + out.write(str(name) + "\n") + + +def _write_initial_grid(path, opts): + columns = hyperpipeline_io.build_column_list(use_sky=True) + base = np.array([ + 0.0, + 1.0, + opts.mass1, + opts.mass2, + opts.spin1x, + opts.spin1y, + opts.spin1z, + opts.spin2x, + opts.spin2y, + opts.spin2z, + opts.ecliptic_longitude, + opts.ecliptic_latitude, + ]) + + 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) + hyperpipeline_io.write_table(path, columns, np.array(rows)) + + +def build_parser(): + parser = argparse.ArgumentParser( + description="Write CEPP contract files for a standalone LISA-RIFT run." + ) + parser.add_argument("--working-directory", default=".") + parser.add_argument("--input-grid", default="proposed-grid.dat") + 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("--transfer-file-list", default="helper_transfer_files.txt") + parser.add_argument("--cepp-command-file", default="command-cepp-lisa.sh") + + parser.add_argument("--cache-file", default="lisa.cache") + parser.add_argument( + "--channel-name", + action="append", + default=["A=fake_strain", "E=fake_strain", "T=fake_strain"], + 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"], + help="PSD assignment, e.g. A=A_psd.xml.gz.", + ) + + parser.add_argument("--mass1", type=float, default=1.0e5) + parser.add_argument("--mass2", type=float, default=8.0e4) + parser.add_argument("--spin1x", type=float, default=0.0) + parser.add_argument("--spin1y", type=float, default=0.0) + parser.add_argument("--spin1z", type=float, default=0.0) + parser.add_argument("--spin2x", type=float, default=0.0) + parser.add_argument("--spin2y", type=float, default=0.0) + 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("--grid-size", type=int, default=3) + parser.add_argument("--grid-fractional-width", type=float, default=1.0e-3) + + parser.add_argument("--approximant", default="IMRPhenomD") + parser.add_argument("--fmin-template", type=float, default=1.0e-3) + 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) + parser.add_argument("--l-max", type=int, default=2) + 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("--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("--save-P", type=float, default=0.1) + + parser.add_argument("--cip-fit-method", default="quadratic") + 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("--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") + parser.add_argument("--n-samples-per-job", type=int, default=1) + parser.add_argument("--n-iterations", type=int, default=1) + parser.add_argument("--request-memory-ILE", type=int, default=4096) + parser.add_argument("--request-memory-CIP", type=int, default=4096) + return parser + + +def main(argv=None): + opts = build_parser().parse_args(argv) + if opts.grid_size < 1: + raise ValueError("--grid-size must be positive") + + workdir = os.path.abspath(opts.working_directory) + os.makedirs(workdir, exist_ok=True) + + 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) + transfer_file_list = os.path.join(workdir, opts.transfer_file_list) + cepp_command_file = os.path.join(workdir, opts.cepp_command_file) + + _write_initial_grid(input_grid, opts) + + channel_args = [] + transfer_files = [opts.cache_file] + for assignment in opts.channel_name: + name, value = _split_assignment(assignment, "--channel-name") + channel_args.extend(["--channel-name", "{}={}".format(name, value)]) + psd_args = [] + for assignment in opts.psd_file: + name, value = _split_assignment(assignment, "--psd-file") + psd_args.extend(["--psd-file", "{}={}".format(name, value)]) + transfer_files.append(value) + + ile_parts = [ + "--LISA", + "--h5-frame-FD", + "--lisa-fixed-sky", "1", + "--ecliptic-longitude", opts.ecliptic_longitude, + "--ecliptic-latitude", opts.ecliptic_latitude, + "--lisa-reference-time", opts.lisa_reference_time, + "--lisa-reference-frequency", opts.lisa_reference_frequency, + "--modes", opts.modes, + "--cache-file", opts.cache_file, + "--event-time", opts.event_time, + ] + channel_args + psd_args + [ + "--fmin-template", opts.fmin_template, + "--fmin-ifo", "A={}".format(opts.fmin_template), + "--fmin-ifo", "E={}".format(opts.fmin_template), + "--fmin-ifo", "T={}".format(opts.fmin_template), + "--fmax", opts.fmax, + "--reference-freq", opts.reference_freq, + "--srate", opts.srate, + "--l-max", opts.l_max, + "--approx", opts.approximant, + "--d-max", opts.d_max, + "--d-min", opts.d_min, + "--n-eff", opts.n_eff, + "--n-max", opts.n_max, + "--n-chunk", opts.n_chunk, + "--save-P", opts.save_P, + "--no-adapt", + "--internal-use-lnL", + ] + if opts.zero_likelihood: + ile_parts.append("--zero-likelihood") + _write_arg_file(ile_args, ile_parts) + + cip_line = _quote_join([ + opts.cip_iterations, + "--fit-method", opts.cip_fit_method, + "--parameter", "mc", + "--parameter", "eta", + "--parameter", "ecliptic_longitude", + "--parameter", "ecliptic_latitude", + "--n-output-samples", opts.cip_n_output_samples, + "--lnL-offset", opts.cip_lnL_offset, + "--no-plots", + ]) + _write_cip_list(cip_args_list, [cip_line]) + + _write_arg_file(test_args, [ + "--method", "lame", + "--parameter", "mc", + "--parameter", "eta", + "--iteration", "$(macroiteration)", + "--threshold", opts.test_threshold, + "--always-succeed", + ]) + + _write_transfer_file(transfer_file_list, transfer_files) + + cepp_parts = [ + "RIFT_HYPERPIPELINE_FORMAT=1", + opts.cepp_exe, + "--ile-n-events-to-analyze", 1, + "--input-grid", input_grid, + "--ile-exe", opts.ile_exe, + "--ile-args", ile_args, + "--cip-args-list", cip_args_list, + "--test-args", test_args, + "--working-directory", workdir, + "--n-iterations", opts.n_iterations, + "--n-samples-per-job", opts.n_samples_per_job, + "--n-copies", 1, + "--request-memory-ILE", opts.request_memory_ILE, + "--request-memory-CIP", opts.request_memory_CIP, + "--transfer-file-list", transfer_file_list, + ] + with open(cepp_command_file, "w") as out: + out.write(_quote_join(cepp_parts) + "\n") + + print("Wrote LISA CEPP helper bundle in {}".format(workdir)) + print(" input grid: {}".format(input_grid)) + print(" ILE args: {}".format(ile_args)) + print(" CIP args list: {}".format(cip_args_list)) + print(" test args: {}".format(test_args)) + print(" CEPP command: {}".format(cepp_command_file)) + + +if __name__ == "__main__": + main() diff --git a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode_lisa b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode_lisa new file mode 100755 index 000000000..627819522 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode_lisa @@ -0,0 +1,2482 @@ +#!/usr/bin/env python + +# LISA-specific ILE fork. +# +# This executable is intentionally separate from +# integrate_likelihood_extrinsic_batchmode to avoid conflicts with the normal +# ILE and the incoming calmarg-in-loop work. It was seeded from Aasim Jan's +# LISA-RIFT lisa_rift_paper branch and should be treated as an experimental +# LISA test surface until the newer driver architecture has a stable hook. + +# This code marginalizes over extrinsic parameters by precomputing terms, setting up a sampler (based on adaptive importance sampling) and creating a likelihood function which\ +# is passed downstream to the sampler. + +# FOLLOWUP TEST (if known injection) +# util_PrintInjection.py --inj mdc.xml.gz --event 0 --verbose +# util_PrintInjection.py --inj maxpt_zero_noise.xml.gz.xml.gz --event 0 +""" +Integrate the extrinsic parameters of the prefactored likelihood function. +""" +import functools +from optparse import OptionParser, OptionGroup +import sys +import numpy +import numpy as np +import os +print("###########################################################################################") +print("# Loading packages and samplers") +print("###########################################################################################") +try: + import cupy + xpy_default=cupy + identity_convert = cupy.asnumpy + identity_convert_togpu = cupy.asarray + identity_convert_lnL= identity_convert # see later, will be assigned something different if cupy active AND AC AND gpu + junk_to_check_installed = cupy.array(5) # this will fail if GPU not installed correctly + if not('RIFT_LOWLATENCY' in os.environ): + print(cupy.show_config()) # print provenance/debugging information + + # Check memory allocation + mem_info = cupy.cuda.Device().mem_info # memory available in bytes + n_mb_total = mem_info[1]/1024./1024. + n_mb_free = mem_info[0]/1024/1024. + print( " cupy memory [total, available] {} {} Mb".format(n_mb_total,n_mb_free)) + # Add failure mode test if memory insufficient! + if n_mb_free < 3000: + print( " - low cupy memory - ") + cupy_success=True +except: + print( ' no cupy') +# import numpy as cupy # will automatically replace cupy calls with numpy! + xpy_default=numpy # just in case, to make replacement clear and to enable override + identity_convert = lambda x: x # trivial return itself + identity_convert_togpu = lambda x: x + identity_convert_lnL= lambda x:x # see later + cupy_success=False + +import lal +from igwn_ligolw import utils, lsctables, table, ligolw +from igwn_ligolw.utils import process +import glue.lal + +import RIFT.lalsimutils as lalsimutils +import RIFT.LISA.lalsimutils_compat as lisa_lalsimutils_compat +import RIFT.likelihood.factored_likelihood as factored_likelihood +import RIFT.likelihood.factored_likelihood_LISA as factored_likelihood_LISA +import RIFT.integrators.mcsampler as mcsampler +import RIFT.misc.sky_rotations as sky_rotations +try: + import RIFT.integrators.mcsamplerEnsemble as mcsamplerEnsemble + mcsampler_gmm_ok = True +except: + print(" No mcsamplerEnsemble ") + mcsampler_gmm_ok = False +try: + import RIFT.integrators.mcsamplerGPU as mcsamplerGPU + mcsampler_gpu_ok = True +except: + print( " No mcsamplerGPU ") + mcsampler_gpu_ok = False +try: + import RIFT.integrators.mcsamplerAdaptiveVolume as mcsamplerAdaptiveVolume + mcsampler_AV_ok = True +except: + print(" No mcsamplerAV ") + mcsampler_AV_ok = False +mcsampler_Portfolio_ok=False +try: + import RIFT.integrators.mcsamplerPortfolio as mcsamplerPortfolio + mcsampler_Portfolio_ok = True +except: + print(" No mcsamplerPortolfio ") +lisa_lalsimutils_compat.install_choose_waveform_print_params_lisa() + +import RIFT.likelihood.priors_utils as priors_utils +import RIFT.misc.xmlutils as xmlutils + +########################################################################################### +class EvenBivariateLinearInterpolator: + def __init__(self, x0, dx, y0, dy, f): + self._x0 = x0 + self._dx = dx + self._y0 = y0 + self._dy = dy + self._fgrid = xpy_default.asarray(f) + + self._dx_inv = 1.0 / self._dx + self._dy_inv = 1.0 / self._dy + + self._N, self._M = self._fgrid.shape + + def __call__(self, x, y): + # Compute the fractional indices into the lookup table where the free + # parameters lie. + i_mid = self._dx_inv * (x - self._x0) + j_mid = self._dy_inv * (y - self._y0) + + # Compute the floor and ceiling of the fractional indices to get the + # indices of the boundaries of the bin `x` and `y` lie in. + # NOTE: In the case where `x` or `y` lie directly on a boundary, the + # floor and ceiling will be equal, but the output is written in such a + # way that we'd still get the correct result. + i_lo = xpy_default.floor(i_mid).astype(int) + j_lo = xpy_default.floor(j_mid).astype(int) + i_hi = xpy_default.ceil(i_mid).astype(int) + j_hi = xpy_default.ceil(j_mid).astype(int) + + # Compute just the fractional part of each index from the low and high + # points. + p = i_mid - i_lo + q = j_mid - j_lo + p_ = 1-p + q_ = 1-q + + # Compute the interpolated result. + f_approx = p_*q_ * self._fgrid[i_lo,j_lo] + f_approx += p *q_ * self._fgrid[i_hi,j_lo] + f_approx += p_*q * self._fgrid[i_lo,j_hi] + f_approx += p *q * self._fgrid[i_hi,j_hi] + + return f_approx + + +__author__ = "Evan Ochsner , Chris Pankow , R. O'Shaughnessy" + + + +# +# Pinnable parameters -- for command line processing +# +LIKELIHOOD_PINNABLE_PARAMS = ["right_ascension", "declination", "psi", "distance", "phi_orb", "t_ref", "inclination"] + +def get_pinned_params(opts): + """ + Retrieve a dictionary of user pinned parameters and their pin values. + """ + return dict([(p,v) for p, v in opts.__dict__.items() if p in LIKELIHOOD_PINNABLE_PARAMS and v is not None]) + +def get_unpinned_params(opts, params): + """ + Retrieve a set of unpinned parameters. + """ + return params - set([p for p, v in opts.__dict__.items() if p in LIKELIHOOD_PINNABLE_PARAMS and v is not None]) + +def zero_like(*args,**kwargs): + if len(kwargs)>0: + arg0 = kwargs['psi'] + else: + arg0 = args[0] + return xpy_default.zeros(len(arg0)) + +def unit_like(*args,**kwargs): + if len(kwargs)>0: + arg0 = kwargs['psi'] + else: + arg0 = args[0] + return xpy_default.ones(len(arg0)) + + +def sampler_param_tuple(sampler,args): + param_list = sampler.params_ordered + return tuple( [param_list.index(x) for x in args]) + + +####################### Option parsing ################################### +optp = OptionParser() +optp.add_option( "--zero-likelihood", action='store_true', help="Run with exactly zero likelihood. Prior test (e.g., alternative distance priors)") +optp.add_option("-c", "--cache-file", default=None, help="LIGO cache file containing all data needed.") +optp.add_option("-C", "--channel-name", action="append", help="instrument=channel-name, e.g. H1=FAKE-STRAIN. Can be given multiple times for different instruments.") +optp.add_option("-p", "--psd-file", action="append", help="instrument=psd-file, e.g. H1=H1_PSD.xml.gz. Can be given multiple times for different instruments.") +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("-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') +optp.add_option("-f", "--reference-freq", type=float, default=100.0, help="Waveform reference frequency. Required, default is 100 Hz.") +optp.add_option("--fmin-template", dest='fmin_template', type=float, default=40, help="Waveform starting frequency. Default is 40 Hz. Also equal to starting frequency for integration") +optp.add_option("--fmin-template-correct-for-lmax",action='store_true',help="Modify amount of data selected, waveform starting frequency to account for l-max, to better insure all requested modes start within the targeted band") +optp.add_option("--fmin-ifo", action='append' , help="Minimum frequency for each IFO. Implemented by setting the PSD=0 below this cutoff. Use with care.") +#optp.add_option("--nr-params",default=None, help="List of specific NR parameters and groups (and masses?) to use for the grid.") +#optp.add_option("--nr-index",type=int,default=-1,help="Index of specific NR simulation to use [integer]. Mass used: mtot= m1+m2") +optp.add_option('--nr-group', default=None,help="If using a *ssingle specific simulation* specified on the command line, provide it here") +optp.add_option('--nr-param', default=None,help="If using a *ssingle specific simulation* specified on the command line, provide it here") +optp.add_option("--nr-lookup",action='store_true', help=" Look up parameters from an NR catalog, instead of using the approximant specified") +optp.add_option("--nr-lookup-group",action='append', help="Restriction on 'group' for NR lookup") +optp.add_option("--nr-hybrid-use",action='store_true',help="Enable use of NR (or ROM!) hybrid, using --approx as the default approximant and with a frequency fmin") +optp.add_option("--nr-hybrid-method",default="taper_add",help="Hybridization method for NR (or ROM!). Passed through to LALHybrid. pseudo_aligned_from22 will provide ad-hoc higher modes, if the early-time hybridization model only includes the 22 mode") +optp.add_option("--rom-group",default=None) +optp.add_option("--rom-param",default=None) +optp.add_option("--rom-use-basis",default=False,action='store_true',help="Use the ROM basis for inner products.") +optp.add_option("--rom-limit-basis-size-to",default=None,type=int) +optp.add_option("--rom-integrate-intrinsic",default=False,action='store_true',help='Integrate over intrinsic variables. REQUIRES rom_use_basis at present. ONLY integrates in mass ratio as present') +optp.add_option("--nr-perturbative-extraction",default=False,action='store_true') +optp.add_option("--nr-perturbative-extraction-full",default=False,action='store_true') +optp.add_option("--nr-use-provided-strain",default=False,action='store_true') +optp.add_option("--no-memory",default=False,action='store_true', help="At present, turns off m=0 modes. Use with EXTREME caution only if requested by model developer") +optp.add_option("--restricted-mode-list-file",default=None,help="A list of ALL modes to use in likelihood. Incredibly dangerous. Only use when comparing with models which provide restricted mode sets, or otherwise to isolate the effect of subsets of modes on the whole") +optp.add_option("--use-gwsignal",default=False,action='store_true',help='Use gwsignal. In this case the approx name is passed as a string to the lalsimulation.gwsignal interface') +optp.add_option("--use-gwsignal-lmax-nyquist",default=None,type=int,help='Passes lmax_nyquist integer to the gwsignal waveform interface') +optp.add_option("--use-external-EOB",default=False,action='store_true') +optp.add_option("--maximize-only",default=False, action='store_true',help="After integrating, attempts to find the single best fitting point") +optp.add_option("--dump-lnL-time-series",default=False, action='store_true',help="(requires --sim-xml) Dump lnL(t) at the injected parameters") +optp.add_option("-a", "--approximant", default="TaylorT4", help="Waveform family to use for templates. Any approximant implemented in LALSimulation is valid.") +optp.add_option("-A", "--amp-order", type=int, default=0, help="Include amplitude corrections in template waveforms up to this e.g. (e.g. 5 <==> 2.5PN), default is Newtonian order.") +optp.add_option("--l-max", type=int, default=2, help="Include all (l,m) modes with l less than or equal to this value.") +optp.add_option("-s", "--data-start-time", type=float, default=None, help="GPS start time of data segment. If given, must also give --data-end-time. If not given, sane start and end time will automatically be chosen.") +optp.add_option("-e", "--data-end-time", type=float, default=None, help="GPS end time of data segment. If given, must also give --data-start-time. If not given, sane start and end time will automatically be chosen.") +optp.add_option("--data-integration-window-half",default=75*1e-3,type=float,help="Only change this window size if you are an expert. The window for time integration is -/+ this quantity around the event time") +optp.add_option("-F", "--fmax", type=float, help="Upper frequency of signal integration. Default is use PSD's maximum frequency.") +optp.add_option("--srate",default=16384,type=float,help="Sampling rate. Change ONLY IF YOU ARE ABSOLUTELY SURE YOU KNOW WHAT YOU ARE DOING.") +optp.add_option("-t", "--event-time", type=float, help="GPS time of the event --- probably the end time. Required if --coinc-xml not given.") +optp.add_option("-i", "--inv-spec-trunc-time", type=float, default=8., help="Timescale of inverse spectrum truncation in seconds (Default is 8 - give 0 for no truncation)") +optp.add_option("-w", "--window-shape", type=float, default=0, help="Shape of Tukey window to apply to data (default is no windowing)") +optp.add_option("--psd-window-shape", type=float, default=0, help="Shape of Tukey window that *was* applied to the PSD being passed. If nonzero, we will rescale the PSD by the ratio of window shape results") +optp.add_option("-m", "--time-marginalization", action="store_true", help="Perform marginalization over time via direct numerical integration. Default is false.") +#optp.add_option("--n-fairdraw-extrinsic-samples",default=None,type=int,help="Extracts a concrete number of fair draw extrinsic samples, bounded above by n_eff") +optp.add_option("--resample-time-marginalization",action='store_true', help="If time-marginalizaiton is true (and should almost always be true), at the end export step use resampling. REQUIRES using fairdraw-extrinsic-output") +optp.add_option("-d", "--distance-marginalization", action="store_true", help="Perform marginalization over distance via a look-up table. Default is false.") +optp.add_option("-l", "--distance-marginalization-lookup-table", default=None, help="Look-up table for distance marginalization.") +optp.add_option("--vectorized", action="store_true", help="Perform manipulations of lm and timeseries using numpy arrays, not LAL data structures. (Combine with --gpu to enable GPU use, where available)") +optp.add_option("--gpu", action="store_true", help="Perform manipulations of lm and timeseries using numpy arrays, CONVERTING TO GPU when available. You MUST use this option with --vectorized (otherwise it is a no-op). You MUST have a suitable version of cupy installed, your cuda operational, etc") +optp.add_option("--force-gpu-only", action="store_true", help="Hard fail if no GPU present (assessed by cupy not loading)") +optp.add_option("--force-xpy", action="store_true", help="Use the xpy code path. Use with --vectorized --gpu to use the fallback CPU-based code path. Useful for debugging.") +optp.add_option("-o", "--output-file", help="Save result to this file.") +optp.add_option("-O", "--output-format", default='xml', help="[xml|hdf5]") +optp.add_option("-S", "--save-samples", action="store_true", help="Save sample points to output-file. Requires --output-file to be defined.") +optp.add_option("-L", "--save-deltalnL", type=float, default=float("Inf"), help="Threshold on deltalnL for points preserved in output file. Requires --output-file to be defined") +optp.add_option("-P", "--save-P", type=float,default=0.1, help="Threshold on cumulative probability for points preserved in output file. Requires --output-file to be defined") +optp.add_option("--internal-hard-fail-on-error",action='store_true',help='If true, fails with exit code 1 if any point is unsuccessful') +optp.add_option("--internal-soft-fail-on-cuda-error",action='store_true',help='If true, returns with exit code 0 on any CUDA error. Use with care (e.g., if many jobs failing)') +optp.add_option("--internal-make-empty-file-on-error",action='store_true',help='If true, failed points generate empty output file. Protects against OSG workflow problems') +optp.add_option("--internal-waveform-fd-L-frame",action='store_true',help='If true, passes extra_waveform_kwargs = {fd_L_frame=True} to lalsimutils hlmoft. Impacts outputs of ChooseFDWaveform calls only.') +optp.add_option("--internal-waveform-fd-no-condition",action='store_true',help='If true, adds extra_waveform_kwargs = {no_condition=True} to lalsimutils hlmoft. Impacts outputs of ChooseFDWaveform calls only. Provided to enable controlled tests of conditioning impact on PE') +optp.add_option("--internal-waveform-extra-lalsuite-args",type=str,default=None) +optp.add_option("--verbose",action='store_true') +optp.add_option("--save-eccentricity", action="store_true") + +# +# Add LISA options +# +lisa_params = OptionGroup(optp, "LISA arguments") +lisa_params.add_option("--LISA", action="store_true", help = "Let the code know that you are analysing LISA data. For now makes the code use analyze_event_LISA instead of analyze_event.") +lisa_params.add_option("--h5-frame", action="store_true", help = "LISA injection frames are in h5 format. Data is assumed to be in time domain when this option is used.") +lisa_params.add_option("--h5-frame-FD", action="store_true", help = "LISA injection frames are in h5 format. Data is assumed to be in frequency domain when this option is used.") +lisa_params.add_option("--lisa-reference-time", default=0.0, type=float, help = "Time in seconds at reference frequency, if reference frequency is not provided then reference frequency is the frequency at max(f^2 * A_22(f)).") +lisa_params.add_option("--lisa-reference-frequency", default=None, type=float, help = "Reference frequency in Hz, if reference frequency is not provided then reference frequency is the frequency at max(f^2 * A_22(f)).") +lisa_params.add_option("--modes", default=None, help = "If you need specific modes, set modes to an array eg --modes '[(2,2),(3,3)]") +lisa_params.add_option("--lisa-fixed-sky", default=False, help="Not varying skylocation") +lisa_params.add_option("--ecliptic-latitude", default=0, help="Value of ecliptic latitude (beta) if sky location is fixed") +lisa_params.add_option("--ecliptic-longitude", default=0, help="Value of ecliptic longitude (lambda) if sky location is fixed") +optp.add_option_group(lisa_params) + +# +# Add the integration options +# +integration_params = OptionGroup(optp, "Integration Parameters", "Control the integration with these options.") +# Default is actually None, but that tells the integrator to go forever or until n_eff is hit. +integration_params.add_option("--n-max", type=int, help="Total number of samples points to draw. If this number is hit before n_eff, then the integration will terminate. Default is 'infinite'.",default=1e7) +integration_params.add_option("--n-eff", type=int, default=100, help="Total number of effective samples points to calculate before the integration will terminate. Default is 100") +integration_params.add_option("--fairdraw-extrinsic-output", action='store_true' , help="Output is fair draw, rather than being comprehensive") +integration_params.add_option("--n-chunk", type=int, help="Chunk'.",default=10000) +integration_params.add_option("--convergence-tests-on",default=False,action='store_true') +integration_params.add_option("--seed", type=int, help="Random seed to use. Default is to not seed the RNG.") +integration_params.add_option("--no-adapt", action="store_true", help="Turn off adaptive sampling. Adaptive sampling is on by default.") +integration_params.add_option("--force-adapt-all", action="store_true", help="Force adaptive sampling for all parameters.") +integration_params.add_option("--force-reset-all", action="store_true", help="Force reset of sampling every iteration. (Recommended if AC and not using no-adapt-after-first)") +integration_params.add_option("--no-adapt-distance", action="store_true", help="Turn off adaptive sampling, just for distance. Adaptive sampling is on by default.") +integration_params.add_option("--no-adapt-after-first",action='store_true',help="Disables adaptation after first iteration with significant lnL") +integration_params.add_option("--adapt-weight-exponent", type=float, default=1.0, help="Exponent to use with weights (likelihood integrand) when doing adaptive sampling. Used in tandem with --adapt-floor-level to prevent overconvergence. Default is 1.0.") +integration_params.add_option("--adapt-floor-level", type=float, default=0.1, help="Floor to use with weights (likelihood integrand) when doing adaptive sampling. This is necessary to ensure the *sampling* prior is non zero during adaptive sampling and to prevent overconvergence. Default is 0.1 (no floor)") +integration_params.add_option("--adapt-adapt",action='store_true',help="Adapt the tempering exponent") +integration_params.add_option("--adapt-log",action='store_true',help="Use a logarithmic tempering exponent") +integration_params.add_option("--interpolate-time", default=False,help="If using time marginalization, compute using a continuously-interpolated array. (Default=false)") +integration_params.add_option("--d-prior",default='Euclidean' ,type=str,help="Distance prior for dL. Options are dL^2 (Euclidean), uniform and 'pseudo_cosmo' .") +integration_params.add_option("--d-max", default=10000,type=float,help="Maximum distance in volume integral. Used to SET THE PRIOR; changing this value changes the numerical answer.") +integration_params.add_option("--d-min", default=1,type=float,help="Minimum distance in volume integral. Used to SET THE PRIOR; changing this value changes the numerical answer.") +integration_params.add_option("--declination-cosine-sampler",action='store_true',help="If specified, the parameter used for declination is cos(dec), not dec") +integration_params.add_option("--inclination-cosine-sampler",action='store_true',help="If specified, the parameter used for inclination is cos(dec), not dec") +integration_params.add_option("--internal-rotate-phase", action='store_true',help="If specified, the integration sampler uses phase_p ==phi+psi and phase_m == phi-psi as sampling coordinates, both ranging from 0 to 4 pi. The prior is twice as large.") +integration_params.add_option("--internal-sky-network-coordinates",action='store_true',help="If specified, perform integration in sky coordinates aligned with the first two IFOs provided") +integration_params.add_option("--internal-sky-network-coordinates-raw",action='store_true',help="If specified, does not attempt to organize IFO network sensibly, uses them AS PROVIDED IN ORDER.") +integration_params.add_option("--manual-logarithm-offset",type=float,default=0,help="Target value of logarithm lnL. Integrand is reduced by exp(-manual_logarithm_offset). Important for high-SNR sources! Should be set dynamically") +integration_params.add_option("--auto-logarithm-offset",action='store_true',help="Use the 'guess_snr' field returned in the precompute stage to change --manual-logarithm-offset for each event.") +integration_params.add_option("--internal-use-lnL",action='store_true',help="likelihood returns lnL, and integrator integrates lnL") +integration_params.add_option("--sampler-method",default="adaptive_cartesian_gpu",help="adaptive_cartesian|GMM|adaptive_cartesian_gpu") +integration_params.add_option("--sampler-portfolio",default=None,action='append',type=str,help="comma-separated strings, matching sampler methods other than portfolio") +integration_params.add_option("--sampler-portfolio-args",default=None, action='append', type=str, help='eval-able dictionaryo to be passed to that sampler') +integration_params.add_option("--sampler-xpy",default=None,help="numpy|cupy if the adaptive_cartesian_gpu sampler is active, use that.") +integration_params.add_option("--supplementary-likelihood-factor-code", default=None,type=str,help="Import a module (in your pythonpath!) containing a supplementary factor for the likelihood. Used to impose supplementary external priors of arbitrary complexity and external dependence (e.g., EM observations). EXPERTS-ONLY") +integration_params.add_option("--supplementary-likelihood-factor-function", default=None,type=str,help="With above option, specifies the specific function used as an external prior. EXPERTS ONLY") +integration_params.add_option("--supplementary-likelihood-factor-ini", default=None,type=str,help="With above option, specifies an ini file that is parsed (here) and passed to the preparation code, called when the module is first loaded, to configure the module. EXPERTS ONLY") +optp.add_option_group(integration_params) + +# +# Add the intrinsic parameters +# +intrinsic_params = OptionGroup(optp, "Intrinsic Parameters", "Intrinsic parameters (e.g component mass) to use.") +intrinsic_params.add_option("--pin-distance-to-sim",action='store_true', help="Pin *distance* value to sim entry. Used to enable source frame reconstruction with NR.") +intrinsic_params.add_option("--mass1", type=float, help="Value of first component mass, in solar masses. Required if not providing coinc tables.") +intrinsic_params.add_option("--mass2", type=float, help="Value of second component mass, in solar masses. Required if not providing coinc tables.") +intrinsic_params.add_option("--spin1z", type=float, help="Value of first component spin (aligned with angular momentum), dimensionless.") +intrinsic_params.add_option("--spin2z", type=float, help="Value of second component spin (aligned with angular momentum), dimensionless.") +intrinsic_params.add_option("--eff-lambda", type=float, help="Value of effective tidal parameter. Optional, ignored if not given.") +intrinsic_params.add_option("--deff-lambda", type=float, help="Value of second effective tidal parameter. Optional, ignored if not given") +intrinsic_params.add_option("--export-eos-index",action='store_true') +optp.add_option_group(intrinsic_params) + + +# +# Add options to integrate over intrinsic parameters. Same conventions as util_ManualOverlapGrid.py. +# Parameters have special names, and we adopt priors that use those names. +# NOTE: Only 'q' implemented +# +intrinsic_int_params = OptionGroup(optp, "Intrinsic integrated parameters", "Intrinsic parameters to integrate over. ONLY currently used with ROM version") +intrinsic_int_params.add_option("--parameter",action='append') +intrinsic_int_params.add_option("--parameter-range",action='append',type=str) +intrinsic_int_params.add_option("--adapt-intrinsic",action='store_true') +optp.add_option_group(intrinsic_int_params) + +# +# Add the pinnable parameters +pinnable = OptionGroup(optp, "Pinnable Parameters", "Specifying these command line options will pin the value of that parameter to the specified value with a probability of unity.") +for pin_param in LIKELIHOOD_PINNABLE_PARAMS: + option = "--" + pin_param.replace("_", "-") + pinnable.add_option(option, type=float, help="Pin the value of %s." % pin_param) +optp.add_option_group(pinnable) + +opts, args = optp.parse_args() + +# +# Failure modes +# +ok_lnL_methods = ['GMM', 'adaptive_cartesian', 'adaptive_cartesian_gpu','AV','portfolio'] +if opts.internal_use_lnL and not(opts.sampler_method in ok_lnL_methods ): + print(" OPTION MISMATCH : --internal-use-lnL not compatible with", opts.sampler_method, " can only use ", ok_lnL_methods) + sys.exit(99) +# if we are on a GPU and using a GPU-accelerated likelihood, don't send the likelihood data back from the GPU needlessly +if cupy_success and opts.gpu and opts.sampler_method == 'adaptive_cartesian_gpu': + identity_convert_lnL = lambda x:x + + +if opts.resample_time_marginalization: + import scipy.special +if opts.resample_time_marginalization and not(opts.fairdraw_extrinsic_output): + raise Exception(" Resampled time output requires --fairdraw-extrinsic-output ") + +# Fairdraw is NOT YET IMPLEMENTED for these other integrators! +#if opts.fairdraw_extrinsic_output: +# if opts.sampler_method == 'adaptive_cartesian_gpu': +# raise Exception(" Fairdraw not available for this sampler") +# if opts.sampler_method == 'GMM': +# raise Exception(" Fairdraw not available for this sampler") + + +supplemental_ln_likelihood= None +supplemental_ln_likelhood_prep=None +supplemental_ln_likelhood_parsed_ini=None +# Supplemental likelihood factor. Must have identical call sequence to 'likelihood_function'. Called with identical raw inputs (including cosines/etc) +if opts.supplementary_likelihood_factor_code and opts.supplementary_likelihood_factor_function: + print(" EXTERNAL SUPPLEMENTARY LIKELIHOOD FACTOR : {}.{} ".format(opts.supplementary_likelihood_factor_code,opts.supplementary_likelihood_factor_function)) + __import__(opts.supplementary_likelihood_factor_code) + external_likelihood_module = sys.modules[opts.supplementary_likelihood_factor_code] + supplemental_ln_likelihood = getattr(external_likelihood_module,opts.supplementary_likelihood_factor_function) + name_prep = "prepare_"+opts.supplementary_likelihood_factor_function + if hasattr(external_likelihood_module,name_prep): + supplemental_ln_likelhood_prep=getattr(external_likelihood_module,name_prep) + # Check for and load in ini file associated with external library + if opts.supplementary_likelihood_factor_ini: + import configparser as ConfigParser + config = ConfigParser.ConfigParser() + config.optionxform=str # force preserve case! + config.read(opts.supplementary_likelihood_factor_ini) + supplemental_ln_likelhood_parsed_ini=config + +if opts.distance_marginalization: + lookup_table = np.load(opts.distance_marginalization_lookup_table) + opts.maximize_only = False + +if opts.gpu is None: + opts.gpu = False # None can be treated differently from false? + +if opts.gpu and opts.force_gpu_only and not cupy_success: + print("GPU requested but no GPU/cupy found: Hard fail by request") + sys.exit(35) # unique code to make resuming due to this error easier to identify + +if opts.gpu and xpy_default is numpy: + print( " Override --gpu (not available); use --force-xpy to require the identical code path is used (with xpy =[np|cupy]") + opts.gpu=False + if opts.force_xpy: + opts.gpu=True + +manual_avoid_overflow_logarithm=opts.manual_logarithm_offset +manual_avoid_overflow_logarithm_default = manual_avoid_overflow_logarithm + +# LISA check (Sampling rate) +deltaT = None +fSample= opts.srate # change sampling rate +if not(fSample is None): + deltaT =1./fSample + + +# Load in restricted mode set, if available +restricted_mode_list=None +if not(opts.restricted_mode_list_file is None): + modes =numpy.loadtxt(opts.restricted_mode_list_file,dtype=int) # columns are l m. Must contain all. Only integers obviously + restricted_mode_list = [ (l,m) for l,m in modes] + print( " RESTRICTED MODE LIST target :", restricted_mode_list) + +intrinsic_param_names = opts.parameter +valid_intrinsic_param_names = ['q'] +if intrinsic_param_names: + for param in intrinsic_param_names: + + # Check if in the valid list + if not(param in valid_intrinsic_param_names): + print( ' Invalid param ', param, ' not in ', valid_intrinsic_param_names) + sys.exit(1) + param_ranges = [] + if len(intrinsic_param_names) == len(opts.parameter_range): + param_ranges = numpy.array(map(eval, opts.parameter_range)) + # Rescale mass-dependent ranges to SI units + for p in ['mc', 'm1', 'm2', 'mtot']: + if p in intrinsic_param_names: + indx = intrinsic_param_names.index(p) + param_ranges[indx]= numpy.array(param_ranges[indx])* lal.MSUN_SI + + +# LISA check (what if I am working in FD?) +# Check both or neither of --data-start/end-time given +if opts.data_start_time is None and opts.data_end_time is not None: + raise ValueError("You must provide both or neither of --data-start-time and --data-end-time.") +if opts.data_end_time is None and opts.data_start_time is not None: + raise ValueError("You must provide both or neither of --data-start-time and --data-end-time.") + +# +# Import NR grid +# +NR_template_group=None +NR_template_param=None +if opts.nr_group and opts.nr_param: + import NRWaveformCatalogManager3 as nrwf + NR_template_group = opts.nr_group + if nrwf.internal_ParametersAreExpressions[NR_template_group]: + NR_template_param = eval(opts.nr_param) + else: + NR_template_param = opts.nr_param + + +# +# Hardcoded variables +# +template_min_freq = opts.fmin_template # minimum frequency of template +#t_ref_wind = 50e-3 # Interpolate in a window +/- this width about event time. +t_ref_wind = opts.data_integration_window_half +T_safety = 2. # Safety buffer (in sec) for wraparound corruption + +# +# Inverse spectrum truncation control +# +T_spec = opts.inv_spec_trunc_time +if T_spec == 0.: # Do not do inverse spectrum truncation + inv_spec_trunc_Q = False + T_safety += 8. # Add a bit more safety buffer in this case +else: + inv_spec_trunc_Q = True + +# +# Integrator options +# +n_max = opts.n_max # Max number of extrinsic points to evaluate at +n_eff = opts.n_eff # Effective number of points evaluated + +# +# Initialize the RNG, if needed +# +# TODO: Do we seed a given instance of the integrator, or set it for all +# or both? +if opts.seed is not None: + numpy.random.seed(opts.seed) + +# LISA check, reference time instead of event time +if not(opts.LISA): + if opts.event_time is not None: + event_time = glue.lal.LIGOTimeGPS(opts.event_time) + print( "Event time from command line: %s" % str(event_time)) + else: + print( " Error! ") + sys.exit(1) + fiducial_epoch = lal.LIGOTimeGPS() + fiducial_epoch = event_time.seconds + 1e-9*event_time.nanoseconds # no more direct access to gpsSeconds +else: + fiducial_epoch = opts.lisa_reference_time +# +# Template descriptors +# + + +print("###########################################################################################") +print("# Reading grid, data and psd") +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_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)) + if len(P_list) < opts.event: +#+opts.n_events_to_analyze: + print( " Event list of range; soft exit") + sys.exit(0) + n_event_max= np.min([len(P_list), opts.event+opts.n_events_to_analyze]) + P_list = P_list[opts.event:n_event_max] + if len(P_list) ==0: + print(" No events to analyze, terminating with success ") + sys.exit(0) + for P in P_list: + P.radec =False # do NOT propagate the epoch later + # LISA + if opts.LISA: + #P.fref = opts.lisa_reference_frequency + P.fref = opts.reference_freq # we are distinguishing between tempalte fref and lisa_fref + P.tref = 0.0 + P.dist = factored_likelihood_LISA.distGpcRef * 1.e9 * lal.PC_SI + else: + P.fref = opts.reference_freq + P.tref = fiducial_epoch # the XML table + P.dist = factored_likelihood.distMpcRef * 1.e6 * lal.PC_SI # use *nonstandard* distance + P.fmin = template_min_freq + m1 = P.m1/lal.MSUN_SI + m2 =P.m2/lal.MSUN_SI + lambda1, lambda2 = P.lambda1,P.lambda2 + P.phiref=0.0 + P.psi=0.0 + P.incl = 0.0 # only works for aligned spins. Be careful. + #P.fref = opts.reference_freq # twice? + if opts.approximant != "TaylorT4": # not default setting + P.approx = lalsimutils.lalsim.GetApproximantFromString(opts.approximant) # allow user to override the approx setting. Important for NR followup, where no approx set in sim_xml! + if opts.approximant == "EccentricTD": + P.phaseO = 3 + + P = P_list[0] # Load in the physical parameters of the injection. +else: + if opts.mass1 is not None and opts.mass2 is not None: + m1, m2 = opts.mass1, opts.mass2 + else: + raise RuntimeError('Missing --mass1 or --mass2') + s1z=s2z=0. + if opts.spin1z is not None: + s1z = opts.spin1z + if opts.spin2z is not None: + s2z = opts.spin2z + +# LISA check (should I change things here? Seems like it is for Neutron stars, so maybe not) + lambda1, lambda2 = 0, 0 + if opts.eff_lambda is not None: + lambda1, lambda2 = lalsimutils.tidal_lambda_from_tilde(m1, m2, opts.eff_lambda, opts.deff_lambda or 0) + P = lalsimutils.ChooseWaveformParams( + approx = lalsimutils.lalsim.GetApproximantFromString(opts.approximant), + fmin = template_min_freq, + radec = False, # do NOT propagate the epoch later + incl = 0.0, # only works for aligned spins. Be careful. + phiref = 0.0, + theta = 0.0, + phi = 0.0, + psi = 0.0, + m1 = m1 * lal.MSUN_SI, + m2 = m2 * lal.MSUN_SI, + lambda1 = lambda1, + lambda2 = lambda2, + s1z = s1z, + s2z = s2z, + ampO = opts.amp_order, + fref = opts.reference_freq, + tref = fiducial_epoch, + dist = factored_likelihood.distMpcRef * 1.e6 * lal.PC_SI + ) + P_list= [P] + +# Time information not needed for FD data +if not(opts.h5_frame_FD): + # User requested bounds for data segment + if not (opts.data_start_time == None) and not (opts.data_end_time == None): + start_time = opts.data_start_time + end_time = opts.data_end_time + print( "Fetching data segment with start=", start_time) + print( " end=", end_time) + + # Automatically choose data segment bounds so region of interest isn't corrupted + # FIXME: Use estimate, instead of painful full waveform generation call here. + else: + approxTmp = P.approx + LmaxEff = 2 + if opts.fmin_template_correct_for_lmax: + LmaxEff=opts.l_max + T_tmplt = lalsimutils.estimateWaveformDuration(P,LmaxEff) + 4 # Much more robust than the previous (slow, prone-to-crash approach) + # if (m1+m2)<30: + # P.approx = lalsimutils.lalsim.TaylorT4 # should not impact length much. IMPORTANT because EOB calls will fail at the default sampling rate + # htmplt = lalsimutils.hoft(P) # Horribly wasteful waveform gneeration solely to estimate duration. Will also crash if spins are used. + # P.approx = approxTmp + # T_tmplt = - float(htmplt.epoch) # ASSUMES time returned by hlm is a RELATIVE time that does not add tref. P.radec=False ! + # print fiducial_epoch, event_time, htmplt.epoch + # else: + # print " Using estimate for waveform length in a high mass regime: beware!" + # T_tmplt = lalsimutils.estimateWaveformDuration(P) + 4 + T_seg = T_tmplt + T_spec + T_safety # Amount before and after event time + # if opts.use_external_EOB: + # T_seg *=2 # extra safety factor + start_time = float(event_time) - T_seg + end_time = float(event_time) + T_seg + print( "Fetching data segment with start=", start_time) + print( " end=", end_time) + print( "\t\tEvent time is: ", float(event_time)) + print( "\t\tT_seg is: ", T_seg) + +# +# Load in data and PSDs +# +data_dict, psd_dict = {}, {} + +# this if statement will handle LISA FD data. +if opts.h5_frame_FD: + for inst, chan in map(lambda c: c.split("="), opts.channel_name): + print("Reading channel %s from cache %s" % (inst+":"+chan, opts.cache_file)) + data_dict[inst] = lisa_lalsimutils_compat.frame_h5_to_hoff(opts.cache_file,inst) + print("Frequency binning: %e (1/deltaF = %f), length %d" % (data_dict[inst].deltaF, 1/data_dict[inst].deltaF, data_dict[inst].data.length) ) + +# this can handle both LIGO and LISA TD data +else: + for inst, chan in map(lambda c: c.split("="), opts.channel_name): + print( "Reading channel %s from cache %s" % (inst+":"+chan, opts.cache_file)) + data_dict[inst] = lisa_lalsimutils_compat.frame_data_to_non_herm_hoff(opts.cache_file, + inst+":"+chan, start=start_time, stop=end_time, + window_shape=opts.window_shape, deltaT=deltaT, h5_frame=opts.h5_frame) + print( "Frequency binning: %f, length %d" % (data_dict[inst].deltaF, + data_dict[inst].data.length) ) + +flow_ifo_dict = {} +if opts.fmin_ifo: + for inst, freq_str in map(lambda c: c.split("="), opts.fmin_ifo): + freq_low_here = float(freq_str) + print( "Reading low frequency cutoff for instrument %s from %s" % (inst, freq_str), freq_low_here) + flow_ifo_dict[inst] = freq_low_here + +for inst, psdf in map(lambda c: c.split("="), opts.psd_file): + print( "Reading PSD for instrument %s from %s" % (inst, psdf)) + psd_dict[inst] = lalsimutils.get_psd_series_from_xmldoc(psdf, inst) + + deltaF = data_dict[inst].deltaF + psd_dict[inst] = lalsimutils.resample_psd_series(psd_dict[inst], deltaF) + print( "PSD deltaF after interpolation %e" % psd_dict[inst].deltaF) + + # Implement PSD window rescaling: see T1900249 + if opts.psd_window_shape > 0 or opts.window_shape > 0: + # Note this is just a constant scale factor in the likelihood, but it *can* have some effect on parameters if one windowing size is small (as it stretches/compresses the likleihood) + window_fac_psd = lalsimutils.psd_windowing_factor(opts.psd_window_shape, len(psd_dict[inst].data.data)) # assume the windowing factor IS accounted for, and we have to undo it + window_fac_data = lalsimutils.psd_windowing_factor(opts.window_shape, len(data_dict[inst].data.data)) + psd_dict[inst].data.data *= window_fac_data/window_fac_psd # scale the windowing factor to be appropriate for our actual data window: the PSD was measured using a windowing different than the one we use + + # implement cutoff. + if inst in flow_ifo_dict.keys(): + if isinstance(psd_dict[inst], lal.REAL8FrequencySeries): + psd_fvals = psd_dict[inst].f0 + deltaF*numpy.arange(psd_dict[inst].data.length) + psd_dict[inst].data.data[ psd_fvals < flow_ifo_dict[inst]] = 0 # + else: + print('FAIL on PSD import') + sys.exit(1) +# elif isinstance(psd_dict[inst], pylal.xlal.datatypes.real8frequencyseries.REAL8FrequencySeries): # for backward compatibility +# psd_fvals = psd_dict[inst].f0 + deltaF*numpy.arange(len(psd_dict[inst].data)) +# psd_dict[inst].data[psd_fvals < ifo_dict[inst]] =0 + + + assert psd_dict[inst].deltaF == deltaF + + # Highest freq. at which PSD is defined + # if isinstance(psd_dict[inst], + # pylal.xlal.datatypes.real8frequencyseries.REAL8FrequencySeries): + # fmax = psd_dict[inst].f0 + deltaF * (len(psd_dict[inst].data) - 1) + if isinstance(psd_dict[inst], lal.REAL8FrequencySeries): + fmax = psd_dict[inst].f0 + deltaF * (psd_dict[inst].data.length - 1) + + # Assert upper limit of IP integral does not go past where PSD defined + assert opts.fmax is None or opts.fmax<= fmax + # Allow us to target a smaller upper limit than provided by the PSD. Important for numerical PSDs that turn over at high frequency + if opts.fmax and opts.fmax < fmax: + fmax = opts.fmax # fmax is now the upper freq. of IP integral + +# Ensure data and PSDs keyed to same detectors +# For LISA we may only provde one PSD, so this check is not necessary. +if sorted(psd_dict.keys()) != sorted(data_dict.keys()) and not(opts.LISA): + print >>sys.stderr, "Got a different set of instruments based on data and PSDs provided." + +# Ensure waveform has same sample rate, padded length as data +# +# N.B. This assumes all detector data has same sample rate, length +# +# data_dict holds 2-sided FrequencySeries, so their length is the same as +# that of the original TimeSeries that was FFT'd = Nsamples +# Also, deltaF = 1/T, with T = the duration (in sec) of the original TimeSeries +# Therefore 1/(data.length*deltaF) = T/Nsamples = deltaT +P.deltaT = 1./ (data_dict[list(data_dict.keys())[0]].data.length * deltaF) +P.deltaF = deltaF +for Psig in P_list: + Psig.deltaT = P.deltaT + Psig.deltaf = P.deltaF + + +# +# Set up parameters and bounds +# + +# PROBLEM: if too large, you can MISS a source. Does NOT need to be fixed for all masses *IF* the problem really has strong support +dmin = opts.d_min # min distance +dmax = opts.d_max # max distance FOR ANY SOURCE EVER. EUCLIDEAN + + + +dmax_sampling_guess = dmax +distBoundGuess = dmax + +print( "Recommended distance for sampling ", dmax_sampling_guess, " and probably near ", distBoundGuess, " smaller than ", dmax) +print( " (recommendation not yet used) ") +param_limits = { "psi": (0, 2*numpy.pi), + "phi_orb": (0, 2*numpy.pi), + "distance": (dmin, dmax), # CAN LEAD TO CATASTROPHIC FAILURE if dmax is too large (adaptive fails - too few bins) + "right_ascension": (0, 2*numpy.pi), + "declination": (-numpy.pi/2, numpy.pi/2), + "t_ref": (-t_ref_wind, t_ref_wind), + "inclination": (0, numpy.pi) +} +if opts.internal_rotate_phase: + param_limits['psi'] = (0, 4*numpy.pi) + param_limits['phi_orb'] = (0, 4*numpy.pi) + +# +# Parameter integral sampling strategy +# +print("\n###########################################################################################\nInitiating integrator\n###########################################################################################\n") +# Oracle for portfolio if needed + + +# Portfolio +use_portfolio=False +params = {} +sampler = mcsampler.MCSampler() +xpy_asarray_already = functools.partial(xpy_default.asarray,dtype=np.float64) + +if opts.sampler_method == "adaptive_cartesian_gpu": + print(" ILE: {}".format(opts.sampler_method)) + sampler = mcsamplerGPU.MCSampler() + sampler.xpy = xpy_default + sampler.identity_convert=identity_convert + mcsampler = mcsamplerGPU # force use of routines in that file, for properly configured GPU-accelerated code as needed + + xpy_asarray_already = lambda x: x # do nothing because we are already on the board for GPU-generated + + if opts.sampler_xpy == "numpy": + mcsampler.set_xpy_to_numpy() + sampler.xpy= numpy + sampler.identity_convert= lambda x: x + +elif opts.sampler_method == "GMM": + print(" ILE: {}".format(opts.sampler_method)) + sampler = mcsamplerEnsemble.MCSampler() + +elif opts.sampler_method == 'AV': + print(" ILE: {}".format(opts.sampler_method)) + opts.internal_use_lnL=True + + sampler = mcsamplerAdaptiveVolume.MCSampler(n_chunk=opts.n_chunk) # note larger is better, but keep GPU mem limits in mind + sampler.xpy = xpy_default + sampler.identity_convert=identity_convert +# mcsampler = mcsampler # force use of routines in that file, for properly configured GPU-accelerated code as needed + + xpy_asarray_already = lambda x: x # do nothing because we are already on the board for GPU-generated + + if opts.sampler_xpy == "numpy": + mcsampler.set_xpy_to_numpy() + sampler.xpy= numpy + sampler.identity_convert= lambda x: x + +elif opts.sampler_method == "portfolio": + use_portfolio=True + opts.internal_use_lnL=True # required, we only implement those scenarios right now + sampler_list = [] + sampler_types = opts.sampler_portfolio + + # prep xpy, etc + my_xpy = xpy_default + my_identity_convert=identity_convert + my_identity_convert_togpu=identity_convert_togpu + print(" PORTFOLIO ", opts.sampler_portfolio) + xpy_asarray_already = lambda x: x # do nothing because we are already on the board for GPU-generated + if opts.sampler_xpy == "numpy": + mcsampler.set_xpy_to_numpy() + my_xpy= numpy + my_identity_convert= lambda x: x + my_identity_convert_togpu= lambda x: x + for name in sampler_types: + if name =='AV': + sampler = mcsamplerAdaptiveVolume.MCSampler(n_chunk=opts.n_chunk) # enforce now, so provided for setup phase + if name =='GMM': + sampler = mcsamplerEnsemble.MCSampler() + # following override means sampler_method is CHANGED, so THIS MUST BE LAST, and can't condition on portfolio + opts.sampler_method = 'GMM' # this will force the creation/parsing of GMM-specific arguments below, so they are properly passed + if name == "adaptive_cartesian_gpu": + sampler = mcsamplerGPU.MCSampler() + mcsampler = mcsamplerGPU # force use of routines in that file, for properly configured GPU-accelerated code as needed + print('PORTFOLIO: adding {} '.format(name)) + # enable xpy for low level sampler as needed + if hasattr(sampler, 'xpy'): + sampler.xpy = my_xpy + sampler.identity_convert= my_identity_convert + sampler.identity_convert_togpu=identity_convert_togpu + sampler_list.append(sampler) + sampler = mcsamplerPortfolio.MCSampler(portfolio=sampler_list) + sampler.xpy = my_xpy + sampler.identity_convert= my_identity_convert + sampler.identity_convert_togpu= my_identity_convert_togpu + # sampler weights will be CPU-typed, so don't change them + +else: + print(" ILE: original sampler") + print(" ILE: {}".format(opts.sampler_method)) + +# +# Psi -- polarization angle +# sampler: uniform in [0, pi) +# +psi_sampler = mcsampler.ret_uniform_samp_vector_alt( + param_limits["psi"][0], param_limits["psi"][1]) +psi_sampler_cdf_inv = functools.partial(mcsampler.uniform_samp_cdf_inv_vector, + param_limits["psi"][0], param_limits["psi"][1]) +sampler.add_parameter("psi", + pdf = psi_sampler, + cdf_inv = psi_sampler_cdf_inv, + left_limit = param_limits["psi"][0], + right_limit = param_limits["psi"][1], + prior_pdf = mcsampler.uniform_samp_psi, + adaptive_sampling=opts.internal_rotate_phase or opts.force_adapt_all) + +# +# Phi - orbital phase +# sampler: uniform in [0, 2*pi) +# +if not (opts.distance_marginalization and lookup_table["phase_marginalization"]): + phi_sampler = mcsampler.ret_uniform_samp_vector_alt( + param_limits["phi_orb"][0], param_limits["phi_orb"][1]) + phi_sampler_cdf_inv = functools.partial(mcsampler.uniform_samp_cdf_inv_vector, + param_limits["phi_orb"][0], param_limits["phi_orb"][1]) + sampler.add_parameter("phi_orb", + pdf = phi_sampler, + cdf_inv = phi_sampler_cdf_inv, + left_limit = param_limits["phi_orb"][0], + right_limit = param_limits["phi_orb"][1], + prior_pdf = mcsampler.uniform_samp_phase, + adaptive_sampling=opts.internal_rotate_phase or opts.force_adapt_all) + + +# +# inclination - angle of system angular momentum with line of sight +# sampler: cos(incl) uniform in [-1, 1) +# + +adapt_extra_extrinsic=False +if (opts.sampler_method == "adaptive_cartesian_gpu" or opts.sampler_method == 'AV' ) or use_portfolio or opts.force_adapt_all: # this is a better/more stable/faster adaptive code, trust it to adapt in more extrinsic dimensions + adapt_extra_extrinsic=True + +if not opts.inclination_cosine_sampler: + incl_sampler = mcsampler.cos_samp_vector # this is NOT dec_samp_vector, because the angular zero point is different! + incl_sampler_cdf_inv = mcsampler.cos_samp_cdf_inv_vector + sampler.add_parameter("inclination", + pdf = incl_sampler, + cdf_inv = incl_sampler_cdf_inv, + left_limit = param_limits["inclination"][0], + right_limit = param_limits["inclination"][1], + prior_pdf = mcsampler.uniform_samp_theta) # do not adapt in parameter going to zero at edge +else: + incl_sampler = mcsampler.ret_uniform_samp_vector_alt(-1.0,1.0) + incl_sampler_cdf_inv = lambda x: x*2.0-1. #functools.partial(mcsampler.uniform_samp_cdf_inv_vector,-1,1) + sampler.add_parameter("inclination", + pdf = incl_sampler, + cdf_inv = incl_sampler_cdf_inv, + left_limit = -1, + right_limit = 1, + prior_pdf = incl_sampler, + adaptive_sampling=adapt_extra_extrinsic) + +# +# Distance - luminosity distance to source in parsecs +# sampler: uniform distance over [dmin, dmax), adaptive sampling +# +print(f" Using {opts.d_prior} prior in distance.") +if not opts.distance_marginalization: + dist_sampler = mcsampler.ret_uniform_samp_vector_alt( param_limits["distance"][0], param_limits["distance"][1]) + dist_sampler_cdf_inv = functools.partial(mcsampler.uniform_samp_cdf_inv_vector, param_limits["distance"][0], param_limits["distance"][1]) + #dist_sampler=functools.partial( mcsampler.uniform_samp_withfloor_vector, numpy.min([distBoundGuess,param_limits["distance"][1]]), param_limits["distance"][1], 0.001) + dist_prior_pdf = lambda x: x**2/(param_limits["distance"][1]**3/3. - param_limits["distance"][0]**3/3.) + if opts.d_prior == 'pseudo_cosmo': + nm = priors_utils.dist_prior_pseudo_cosmo_eval_norm(param_limits["distance"][0],param_limits["distance"][1]) + dist_prior_pdf =functools.partial( priors_utils.dist_prior_pseudo_cosmo, nm=nm,xpy=xpy_default) + elif opts.d_prior == 'uniform': + dist_prior_pdf = dist_sampler + elif opts.d_prior != 'Euclidean': + print(" ==== WARNING UNKNOWN DISTANCE PRIOR === ") + raise Exception('distance prior') + #dist_sampler_cdf_inv=None + sampler.add_parameter("distance", + pdf = dist_sampler, + cdf_inv = dist_sampler_cdf_inv, + left_limit = param_limits["distance"][0], + right_limit = param_limits["distance"][1], + prior_pdf = dist_prior_pdf, + adaptive_sampling = (adapt_extra_extrinsic and not (opts.no_adapt or opts.no_adapt_distance)) or opts.force_adapt_all) + +# +# Rotate sky coordinates +# +if opts.internal_sky_network_coordinates: + ifo_list = list(psd_dict) + if not(opts.internal_sky_network_coordinates_raw): + # remove V and K : the sky ring is almost always only HL + if 'K1' in ifo_list: + ifo_list.remove('K1') + if 'V1' in ifo_list: + ifo_list.remove('V1') + # problem: ordering of psd_dict is rarely rational; almost always we want an HL network, rarely V + if len(ifo_list) <2 or opts.LISA: + opts.internal_sky_network_coordinates = False # stop using this code + else: + sky_rotations.assign_sky_frame(ifo_list[0], ifo_list[1], fiducial_epoch) + frm = identity_convert_togpu(sky_rotations.frm) + my_rotation = functools.partial(lalsimutils.polar_angles_in_frame_alt,frm,xpy=xpy_default) + my_rotation_cpu = functools.partial(lalsimutils.polar_angles_in_frame_alt,sky_rotations.frm,xpy=np) + +# +# Intrinsic parameters +# +sampler_lookup = {} +sampler_inv_lookup = {} +sampler_lookup['q'] = mcsampler.q_samp_vector # only one intrinsic parameter possible +sampler_lookup['M'] = mcsampler.M_samp_vector +sampler_inv_lookup['q'] = mcsampler.q_cdf_inv_vector +sampler_inv_lookup['M'] = None # mcsampler.M_cdf_inv_vector + +if opts.rom_use_basis and opts.rom_integrate_intrinsic: + for p in intrinsic_param_names: + indx = intrinsic_param_names.index(p) + qmin,qmax = param_ranges[indx] + q_pdf =mcsampler.ret_uniform_samp_vector_alt( qmin, qmax) # sample uniformly by default +# q_pdf =functools.partial(sampler_lookup[p], qmin, qmax) # sample uniformly by default + q_cdf_inv = functools.partial(mcsampler.uniform_samp_cdf_inv_vector, qmin,qmax) # sample uniformly by default +# q_cdf_inv= None + q_pdf_prior = functools.partial(sampler_lookup[param], qmin, qmax) # true prior, from lookup + sampler.add_parameter(p, + pdf=q_pdf, + cdf_inv=q_cdf_inv, + left_limit=qmin, right_limit=qmax,prior_pdf=q_pdf_prior, + adaptive_sampling = opts.adapt_intrinsic) + + + +if False: #opts.skymap_file is not None: + from ligo.skymap.io import fits as bfits + # + # Right ascension and declination -- use a provided skymap + # + smap, _ = bfits.read_sky_map(opts.skymap_file) + # FIXME: Uncomment for 'mixed' map + #smap = 0.9*smap + 0.1*numpy.ones(len(smap))/len(smap) + ss_sampler = mcsampler.HealPixSampler(smap) + #isotropic_2d_sampler = numpy.vectorize(lambda dec, ra: mcsampler.dec_samp_vector(dec)/2/numpy.pi) + isotropic_bstar_sampler = numpy.vectorize(lambda dec, ra: 1.0/len(smap)) + + # FIXME: Should the left and right limits be modified? + sampler.add_parameter(("declination", "right_ascension"), + pdf = ss_sampler.pseudo_pdf, + cdf_inv = ss_sampler.pseudo_cdf_inverse, + left_limit = (param_limits["declination"][0], param_limits["right_ascension"][0]), + right_limit = (param_limits["declination"][1], param_limits["right_ascension"][1]), + prior_pdf = isotropic_bstar_sampler) + +else: + # + # Right ascension - angle in radians from prime meridian plus hour angle + # sampler: uniform in [0, 2pi), adaptive sampling + # + ra_sampler = mcsampler.ret_uniform_samp_vector_alt( + param_limits["right_ascension"][0], param_limits["right_ascension"][1]) + ra_sampler_cdf_inv = functools.partial(mcsampler.uniform_samp_cdf_inv_vector, + param_limits["right_ascension"][0], param_limits["right_ascension"][1]) + sampler.add_parameter("right_ascension", + pdf = ra_sampler, + cdf_inv = ra_sampler_cdf_inv, + left_limit = param_limits["right_ascension"][0], + right_limit = param_limits["right_ascension"][1], + prior_pdf = mcsampler.uniform_samp_phase, + adaptive_sampling = opts.force_adapt_all or ((not opts.no_adapt) and (not opts.internal_sky_network_coordinates))) # TOO DANGEROUS to double-adapt in sky, ends up overconverging too easily. Just leave the whole sky ring in place if we want to do this, it's usually there and we don't lose that much. If we have full localization, we should just remove the sky network coordinates argument! + + # + # declination - angle in radians from the north pole piercing the celestial + # sky sampler: cos(dec) uniform in [-1, 1), adaptive sampling + # + if not opts.declination_cosine_sampler: + dec_sampler = mcsampler.dec_samp_vector + dec_sampler_cdf_inv = mcsampler.dec_samp_cdf_inv_vector + sampler.add_parameter("declination", + pdf = dec_sampler, + cdf_inv = dec_sampler_cdf_inv, + left_limit = param_limits["declination"][0], + right_limit = param_limits["declination"][1], + prior_pdf = mcsampler.uniform_samp_dec, + adaptive_sampling = opts.force_adapt_all or (not opts.no_adapt)) + else: + # Sample uniformly in cos(polar_theta), =1 for north pole, -1 for south pole. + # Propagate carefully in conversions: time of flight libraries use RA,DEC + dec_sampler = mcsampler.ret_uniform_samp_vector_alt(-1.0,1.0) + dec_sampler_cdf_inv = lambda x: x*2.0-1. # functools.partial(mcsampler.uniform_samp_cdf_inv_vector,-1,1) + sampler.add_parameter("declination", + pdf = dec_sampler, + cdf_inv = dec_sampler_cdf_inv, + left_limit = -1, + right_limit = 1, + prior_pdf = dec_sampler, + adaptive_sampling = opts.force_adapt_all or (not opts.no_adapt)) + +# LISA: What will change here if we don't marginalize in time for LISA? +if not opts.time_marginalization: + # + # tref - GPS time of geocentric end time + # sampler: uniform in +/-2 ms window around estimated end time + # + tref_sampler = mcsampler.ret_uniform_samp_vector_alt( + param_limits["t_ref"][0], param_limits["t_ref"][1]) + + tref_sampler_cdf_inv = functools.partial(mcsampler.uniform_samp_cdf_inv_vector, + param_limits["t_ref"][0], param_limits["t_ref"][1]) + sampler.add_parameter("t_ref", + pdf = tref_sampler, + cdf_inv = None, + left_limit = param_limits["t_ref"][0], + right_limit = param_limits["t_ref"][1], + prior_pdf = functools.partial(mcsampler.uniform_samp_vector, param_limits["t_ref"][0], param_limits["t_ref"][1])) + + +# skymap oracle +oracleRS=None +if opts.skymap_file: + # Read file + from RIFT.misc.reference_samples import ReferenceSamples + rs_object = ReferenceSamples() + if opts.skymap_file.endswith('fits'): + rs_object.from_skymap_fits(fname=opts.skymap_file,cos_dec=opts.declination_cosine_sampler) + else: + rs_object.from_ascii(fname=opts.skymap_file, reference_params=['ra', 'dec']) + if opts.declination_cosine_sampler: + rs_object.reference_params[:,1] = np.arccos(rs.reference_params[:,1]) + # relabel + rs_object.reference_params = ['right_ascension', 'declination'] + + # Pass to resampling oracle + from RIFT.integrators.unreliable_oracle.resampling import ResamplingOracle + oracleRS = ResamplingOracle() + for p in sampler.params_ordered: + oracleRS.add_parameter(p,pdf=None, left_limit =sampler.llim[p], right_limit = sampler.rlim[p]) + oracleRS.setup(reference_samples=rs_object.reference_samples, reference_params=rs_object.reference_params) + +# +# Determine pinned and non-pinned parameters +# + +# LISA +# I am calling ecliptic latitude (beta) as declination and ecliptic longitude (lambda) as right-ascension. + +pinned_params = get_pinned_params(opts) +unpinned_params = get_unpinned_params(opts, sampler.params) + +if opts.LISA and opts.lisa_fixed_sky: + lisa_sky_lamda = float(opts.ecliptic_longitude) + lisa_sky_beta = float(opts.ecliptic_latitude) + pinned_params["right_ascension"] = 0.0 + pinned_params["declination"] = 0.0 +if opts.LISA and not(opts.lisa_fixed_sky): # read the grid for sky location co-ordinate, but fixing it to prevent sampler from wasting time sampling skylocation. + pinned_params["right_ascension"] = 0.0 + pinned_params["declination"] = 0.0 + + +print( "{0:<25s} {1:>5s} {2:>5s} {3:>20s} {4:<10s}".format("parameter", "lower limit", "upper limit", "pinned?", "pin value")) +plen = len(sorted(sampler.params, key=lambda p: len(p))[-1]) +for p in sampler.params: + if p in pinned_params: + pinned, value = True, "%1.3g" % pinned_params[p] + else: + pinned, value = False, "" + + if isinstance(p, tuple): + for subp, subl, subr in zip(p, sampler.llim[p], sampler.rlim[p]): + subp = subp + " "*min(0, plen-len(subp)) + print( "|{0:<25s} {1:>1.3g} {2:>1.3g} {3:>20s} {4:<10s}".format(subp, subl, subr, str(False), "")) + else: + p = p + " "*min(0, plen-len(p)) + print("{0:<25s} {1:>5g} {2:>5g} {3:>20s} {4:<10s}".format(p, sampler.llim[p], sampler.rlim[p], str(pinned), value)) + +# Special case: t_ref is assumed to be relative to the epoch +if "t_ref" in pinned_params: + pinned_params["t_ref"] -= float(fiducial_epoch) + +# +# Provide convergence tests +# FIXME: Currently using hardcoded thresholds, poorly hand-tuned +# +test_converged = {} + +# +# Merge options into one big ol' kwargs dict +# + +pinned_params.update({ + # Iteration settings and termination conditions + "n": min(opts.n_chunk, n_max), # Number of samples in a chunk + "nmax": n_max, # Total number of samples to draw before termination + "neff": n_eff, # Total number of effective samples to collect before termination + + "convergence_tests" : test_converged, # Dictionary of convergence tests + + # Adaptive sampling settings + "tempering_exp": opts.adapt_weight_exponent if not opts.no_adapt else 0.0, # Weights will be raised to this power to prevent overconvergence + "tempering_log": opts.adapt_log, + "tempering_adapt": opts.adapt_adapt, + + "floor_level": opts.adapt_floor_level if not opts.no_adapt else 0.0, # The new sampling distribution at the end of each chunk will be floor_level-weighted average of a uniform distribution and the (L^tempering_exp p/p_s)-weighted histogram of sampled points. + "history_mult": 10, # Multiplier on 'n' - number of samples to estimate marginalized 1-D histograms + "n_adapt": 100 if not opts.no_adapt else 0, # Number of chunks to allow adaption over + + # Verbosity settings + "verbose": True, #not opts.rom_integrate_intrinsic, + "extremely_verbose": False, + + # Sample caching + "save_intg": opts.save_samples, # Cache the samples (and integrand values)? + "igrand_threshold_deltalnL": opts.save_deltalnL, # Threshold on distance from max L to save sample + "igrand_threshold_p": opts.save_P, # Threshold on cumulative probability contribution to cache sample + "igrand_fairdraw_samples": opts.fairdraw_extrinsic_output, + "igrand_fairdraw_samples_max": opts.n_eff +}) +if opts.sampler_method == "adaptive_cartesian_gpu": + pinned_params.update({"save_no_samples":True}) # do not exhaust GPU memory with MC samples! +return_lnL=False +if opts.sampler_method=="GMM" and opts.internal_use_lnL: + return_lnL=True + pinned_params.update({"use_lnL":True,"return_lnI":True}) +if opts.sampler_method =="adaptive_cartesian_gpu" and opts.internal_use_lnL: + return_lnL=True + pinned_params.update({"use_lnL":True}) +if opts.sampler_method =="portfolio": + return_lnL=True + pinned_params.update({"use_lnL":True}) +if opts.sampler_method == "GMM": + n_step =pinned_params["n"] + n_max_blocks = ((1.0*int(opts.n_max))/n_step) + # pairing coordinates for adaptive integration: see definition of order below + # (distance,inclination) + # (ra, dec) + # (psi) (orb_phase) # only do 1d adaptivity there, +# gmm_dict = {tuple([2]):None,tuple([4]):None,(3,5):None,(0,1):None} +# comp_dict = {tuple([2]):1,tuple([4]):1,(3,5):3,(0,1):4} + default_phipsi = mcsamplerEnsemble.create_wide_single_component_prior( [ param_limits['psi'], param_limits['phi_orb']]) +# pair_sky = sampler_param_tuple(sampler, ['right_ascension','declination']) + pair_ra_dec = sampler_param_tuple(sampler, ["right_ascension", "declination"]) + n_d = 2 # 2 distance-inclination lobes by default + n_sky = 4 + n_phase =4 + adapt_phase = False + if 'distance' in sampler.params: + pair_d_incl = sampler_param_tuple(sampler, ['distance','inclination']) + else: + pair_d_incl = sampler_param_tuple(sampler, ['inclination']) + n_d = 1 # one peak in distance by default + if 'phi_orb' in sampler.params: + pair_phi_psi = sampler_param_tuple(sampler, ['psi', "phi_orb"]) + else: + pair_phi_psi = sampler_param_tuple(sampler, ['psi']) + adapt_phase = True + n_phase = 2 + gmm_dict = {pair_ra_dec:None,pair_d_incl:None,pair_phi_psi:default_phipsi} + gmm_adapt = {pair_ra_dec: True, pair_d_incl: True, pair_phi_psi: False} + if opts.internal_rotate_phase: + gmm_adapt[pair_phi_psi] = True +# gmm_dict = {tuple([2]):None,tuple([4]):None,(3,5):None,tuple([0]):None,tuple([1]):None} +# if len(psd_dict.keys()) > 2: +# n_sky=2 # presumably we have localized better, avoid failure modes +# comp_dict = {tuple([2]):1,tuple([4]):1,(3,5):3,tuple([0]):n_sky,tuple([1]):n_sky} + comp_dict = {pair_ra_dec:n_sky,pair_d_incl:n_d,pair_phi_psi:n_phase} + extra_args = {'n_comp':comp_dict,'max_iter':n_max_blocks,'gmm_dict':gmm_dict, 'gmm_adapt':gmm_adapt} # made up for now, should adjust + print("GMM:",extra_args) + print("GMM:",sampler.params_ordered) + # if opts.distance_marginalization: + # if lookup_table["phase_marginalization"]: + # gmm_dict = {pair_ra_dec:None,(3,4):None,(0,):None} + # gmm_adapt = {pair_ra_dec: True, (3,4): True, (0,): True} + # comp_dict = {pair_ra_dec:n_sky,tuple([4]):1,(0,):2} + # else: + # gmm_dict = {pair_ra_dec:None,tuple([4]):None,(0,1):None} + # gmm_adapt = {pair_ra_dec:True,tuple([4]):True,(0,1):True} + # comp_dict = {pair_ra_dec:n_sky,tuple([4]):1,(0,1):4} +# extra_args = {'n_comp':comp_dict,'max_iter':n_max_blocks,'gmm_dict':gmm_dict} # made up for now, should adjust + + # force adapt in all parameters, if we request it. Requires pass-by-reference in above + if opts.force_adapt_all: + for param in gmm_adapt: + gmm_adapt[param] = True + pinned_params.update(extra_args) + + # cannot pin params at present + # unpinned params MUST be full call signature of likelihood, IN ORDER + if not(opts.time_marginalization): + unpinned_params =['right_ascension','declination', 't_ref','phi_orb','inclination','psi','distance'] + else: + unpinned_params =['right_ascension','declination', 'phi_orb','inclination','psi','distance'] + if opts.distance_marginalization: + unpinned_params.remove('distance') + if lookup_table["phase_marginalization"]: + unpinned_params.remove('phi_orb') +if opts.sampler_method == "AV": + return_lnL=True + pinned_params.update( { 'enforce_bounds':True}) # don't go out of range : choose integer bin sizes + + +# set up sampler, as needed. Mainly for portfolio integrator +if use_portfolio: + print(" PORTFOLIO : setup") + if opts.sampler_portfolio_args: + print(" PRE_EVAL", opts.sampler_portfolio_args) + #opts.sampler_portfolio_args = list(map(lambda x: eval(' "{}" '.format(x)), opts.sampler_portfolio_args)) + opts.sampler_portfolio_args = list(map(eval, opts.sampler_portfolio_args)) + # confirm all are dict + for indx in range(len(opts.sampler_portfolio_args)): + if not(isinstance(opts.sampler_portfolio_args[indx], dict)): + print(indx,opts.sampler_portfolio_args[indx]) + print(" ARGS ", opts.sampler_portfolio_args) + sampler.setup(portfolio_args=opts.sampler_portfolio_args, **pinned_params) # directly pass all parameters set above to low-level portfolios. In particular, GMM setup + +# initialize sampler, before we call integrate, so we can seed it +if opts.sampler_method == 'adaptive_cartesian_gpu' and opts.skymap_file: + sampler.setup() + + +def resample_samples_LISA(my_samples, rholms, cross_terms, right_ascension, declination, P, modes, reference_distance): + """This function takes in extrinsic samples and for each sample samples a time shift. This is done by generating a likelihood time series at an extrinsic sample and then weighted sampling in time.""" + # How many time samples? Same as the extrinsic samples being passed + n_samples = len(my_samples['longitude']) + print(" Time resampling size : {} ".format(n_samples)) + + # Hardcoded time sampling limits, should change it in future + low_t_lim, high_t_lim, delta_t = -5, 5, 0.001 + + # t_ref_wind is defined as data_window_integration_half + + # This is the time series over which the lnL series is defined. + tvals = xpy_default.linspace(-t_ref_wind, t_ref_wind - P.deltaT,int(t_ref_wind*2/P.deltaT)) # the array should have P.deltaT as deltaT between entries and for that you need to subtract P.deltaT from the final point. + + # Generating likelihood time series for these samples + P.phi = identity_convert_togpu(my_samples['right_ascension']) # cast to float + P.theta = identity_convert_togpu(my_samples['declination']) + P.tref = float(fiducial_epoch) # defined earlier as lisa-reference-time + P.phiref = identity_convert_togpu(my_samples['coa_phase']) + P.incl = identity_convert_togpu(my_samples['inclination']) + P.psi = identity_convert_togpu(my_samples['psi']) + P.dist = identity_convert_togpu(my_samples['distance']* 1.e6 * lalsimutils.lsu_PC) # luminosity distance + + if not(cupy_success): + # convert objects + for name in ['phi','theta', 'phiref','incl', 'psi','dist']: + setattr(P,name, getattr(P,name).astype(float) ) + + # plugging the selected samples + lnLt = factored_likelihood_LISA.FactoredLogLikelihoodAlignedSpinLISA(rholms, cross_terms, declination, right_ascension, P.psi, P.incl, P.phiref, P.dist, modes, reference_distance, return_lnLt = True) + + # axis 0 is time axis 1 is each point, transposing it + lnLt = lnLt.T + lnLt = identity_convert(lnLt) # back to CPU. Note we have removed offsets + if opts.zero_likelihood: + lnLt =np.zeros(lnLt.shape) # zero likelihood + lnLt_norm = scipy.special.logsumexp(lnLt,axis=-1) + tvals = identity_convert(tvals) # back to CPU + + # Loop over and resample in time, picking index according to weights + t_out = np.zeros(n_samples) + lnL_out = np.zeros(n_samples) + # identifiy max lnL(t) point, and then interpolate around that point to avoid Nan weights. + index_max_lnLt = np.argmax(lnLt[0]).flatten() + tval_at_max = tvals[index_max_lnLt] + print(f"Identified max lnL(t) as {tval_at_max}s. Interpolating around that point. Sampling for time from {tval_at_max + low_t_lim}s to {tval_at_max + high_t_lim}s with deltaT {delta_t}s.") + + # tvals over which we will resample + interp_tvals = xpy_default.arange(low_t_lim + tval_at_max, high_t_lim + tval_at_max, delta_t) + indx_list = np.arange(len(interp_tvals)) + # Loop over each extrinsic sample + for indx in np.arange(n_samples): + interp = scipy.interpolate.interp1d(tvals, lnLt[indx] - lnLt_norm[indx], kind='cubic') + # interpolating + interp_lnLt = interp(interp_tvals) + # weighted resampling + indx_choose = np.random.choice(indx_list, p=np.exp(interp_lnLt)/np.sum(np.exp(interp_lnLt))) + t_out[indx] = interp_tvals[indx_choose] + lnL_out[indx] = interp_lnLt[indx_choose] + print(t_out) + # saving + my_samples['t_ref'] = fiducial_epoch + t_out # add sample time jitter from reweighting to samples + my_samples["lnL_raw"] = lnL_out # export likelihoods (equivalent to SNR). Note needs downstream code filters to catch this and put it somewhere. + return my_samples + + +def resample_samples(my_samples, + lookupNKDict=None, rholmArrayDict=None, ctUArrayDict=None, ctVArrayDict=None,epochDict=None): # uses LOTS of global variables, don't pass them all + # will look a LOT like the likelihood function definitions, unfortunately + if not opts.vectorized: + raise Exception( ' resample_samples currently only for vectorized ') + if opts.distance_marginalization: + raise Exception( ' resample_samples currently only final extrinsic samples; you should NOT have distance marginalization on ') + + n_samples = len(my_samples['longitude']) + print(" Time resampling size : {} ".format(n_samples)) + + + # if we are using GPU-based generation this is ok; if itis mcsampler, we will have a lot of 'object' casts to fix, arg + + tvals = xpy_default.linspace(-t_ref_wind,t_ref_wind,int((t_ref_wind)*2/P.deltaT)) # choose an array at the target sampling rate. P is inherited globally + P.phi = identity_convert_togpu(my_samples['right_ascension']) # cast to float + P.theta = identity_convert_togpu(my_samples['declination']) + P.tref = float(fiducial_epoch) + P.phiref = identity_convert_togpu(my_samples['coa_phase']) + P.incl = identity_convert_togpu(my_samples['inclination']) + P.psi = identity_convert_togpu(my_samples['psi']) + P.dist = identity_convert_togpu(my_samples['distance']* 1.e6 * lalsimutils.lsu_PC) # luminosity distance + if not(cupy_success): + # convert objects + for name in ['phi','theta', 'phiref','incl', 'psi','dist']: + setattr(P,name, getattr(P,name).astype(float) ) + + lnLt = factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop(tvals, + P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max,xpy=xpy_default,return_lnLt=True) + + lnLt = identity_convert(lnLt) # back to CPU. Note we have removed offsets + if opts.zero_likelihood: + lnLt =np.zeros(lnLt.shape) # zero likelihood + lnLt_norm = scipy.special.logsumexp(lnLt,axis=-1) + tvals = identity_convert(tvals) # back to CPU +# print(lnLt.shape, lnLt_norm.shape,tvals.shape) + # Loop over and resample in time, picking index according to weights + t_out = np.zeros(n_samples) + lnL_out = np.zeros(n_samples) + indx_list =np.arange(len(tvals)) + for indx in np.arange(n_samples): + indx_choose = np.random.choice(indx_list, p=np.exp(lnLt[indx] - lnLt_norm[indx])) + t_out[indx] = tvals[indx_choose] + lnL_out[indx] = lnLt[indx][indx_choose] +# print(' Resampled time offset {} '.format(t_out[indx])) #, lnLt[indx]-lnLt_norm[indx]) + + my_samples['t_ref'] = fiducial_epoch+t_out # add sample time jitter from reweighting to samples + my_samples["lnL_raw"] = lnL_out # export likelihoods (equivalent to SNR). Note needs downstream code filters to catch this and put it somewhere. + + return my_samples + +# +# Main analysis functions for LIGO (analyze_event) and LISA (analyze_event_LISA). This precomputes terms (eq 23 10.1103/PhysRevD.92.023002) and then based on your settings\ +# creates a likelihood function which is then passed to a sampler. +# + +def analyze_event_LISA(P_list, indx_event, data_dict, psd_dict, fmax, opts, inv_spec_trunc_Q=inv_spec_trunc_Q, T_spec=T_spec): + print("\n###########################################################################################\nPrecomputing\n###########################################################################################") + nEvals=0 + P = P_list[indx_event] + # if pin-distance-to-sim, change the distance prior accordingly + if opts.pin_distance_to_sim: + pinned_params['distance']=P.dist/(1.e6 * lal.PC_SI) + + # Call external likleihood preparation if any + if supplemental_ln_likelhood_prep: + supplemental_ln_likelhood_prep(P=P,config=supplemental_ln_likelhood_parsed_ini) + + print("Generating modes for parameters:") + P.print_params_lisa() + # Generate modes + if opts.modes: + modes=np.array(eval(opts.modes)) + print(f"modes = \n{list(modes)}") + P.deltaF=deltaF # why do we need this? why is it becoming None? + hlms_FD = lisa_lalsimutils_compat.hlmoff_for_LISA(P, opts.l_max, modes) + fNyq = 0.5/P.deltaT + modes = np.array(list(hlms_FD.keys())) + reference_distance = P.dist + print(f"Reference distance = {reference_distance/lal.PC_SI/1e9} Gpc, 1/deltaF = {1/P.deltaF} s, deltaT = {P.deltaT} s") + + # Are we varying skylocation? If we are, read it from the grid and pass it to precompute. If not, take in the value as is passed as an argument to ILE. + if not(opts.lisa_fixed_sky): + lisa_sky_beta = P.theta + lisa_sky_lamda = P.phi + if opts.lisa_fixed_sky: + lisa_sky_lamda = float(opts.ecliptic_longitude) + lisa_sky_beta = float(opts.ecliptic_latitude) + + # Precompute + # fNyq is the max resolvable frequency for a waveform. It is 0.5/deltaT. RIFT needs deltaT, deltaF for waveform generation (information present in P) and for integration it needs fmax (fmax <= fNyq) + print(f"Sky location lambda = {lisa_sky_lamda}, sky location beta = {lisa_sky_beta}") + + rholms_intp, cross_terms, cross_terms_V, rholms, guess_snr, rest = factored_likelihood_LISA.PrecomputeAlignedSpinLISA(opts.lisa_reference_time, opts.lisa_reference_frequency, opts.data_integration_window_half, hlms_FD, None, data_dict, psd_dict, flow_ifo_dict["A"], fNyq, fmax, P.deltaT, lisa_sky_beta, lisa_sky_lamda, analyticPSD_Q=False, inv_spec_trunc_Q=False, T_spec=0.0) + + # reset to default. Should not be needed, but weird python scoping error + manual_avoid_overflow_logarithm = manual_avoid_overflow_logarithm_default + + def likelihood_function(right_ascension, declination, phi_orb, inclination, psi, distance): + P.phi = xpy_asarray_already(right_ascension) + P.theta = xpy_asarray_already(declination) + P.phiref = xpy_asarray_already(phi_orb) + if opts.inclination_cosine_sampler: + P.incl = xpy_default.arccos(xpy_asarray_already(inclination)) + else: + P.incl = xpy_asarray_already(inclination) + P.psi = xpy_asarray_already(psi) + P.dist = xpy_asarray_already(distance* 1.e6 * lalsimutils.lsu_PC) # luminosity distance + + # rotate phase if needed + # make copies of arrays + # Sampling assumes P.phiref == phi+ \psi , P.psi == phi - psi + if opts.internal_rotate_phase: + phi_orb_true = (P.phiref + P.psi)/2. + psi_true = (P.phiref - P.psi)/2. + P.psi= psi_true + P.phiref = phi_orb_true + #print(f"Sky location lambda = {lisa_sky_lamda}, sky location beta = {lisa_sky_beta}") + + + lnL = factored_likelihood_LISA.FactoredLogLikelihoodAlignedSpinLISA(rholms, cross_terms, lisa_sky_beta, lisa_sky_lamda, P.psi, P.incl, P.phiref, P.dist, modes, reference_distance) + return lnL + #return identity_convert_lnL(xpy_default.exp(lnL-manual_avoid_overflow_logarithm)) + + if opts.sampler_method == "adaptive_cartesian_gpu": + # reset sampling parameter as needed (distance, inclination but not sky location) + # distance (and inclination) can conceivably correlate with mass, and we don't want to truncate in distance/inclination prematurely from history effects + # method ONLY implemented for adaptcart! + if 'distance' in sampler.params: + sampler.reset_sampling('distance') + sampler.reset_sampling('inclination') + elif opts.sampler_method == "GMM": + if 'distance' in sampler.params: + pair_d_incl = sampler_param_tuple(sampler, ['distance','inclination']) + if pair_d_incl in gmm_dict: + gmm_dict[pair_d_incl] = None # reset inclination/distance sampling + single_incl = sampler_param_tuple(sampler, ['inclination']) + if (opts.distance_marginalization) and single_incl in gmm_dict: + gmm_dict[single_incl] = None # reset inclination sampling + print("\n###########################################################################################\nIntegrating\n###########################################################################################") + # Integrate + args = likelihood_function.__code__.co_varnames[:likelihood_function.__code__.co_argcount] + print( " --------> Arguments ", args) + # if exactly zero likelihood function + like_to_integrate = likelihood_function + if opts.zero_likelihood: + if opts.internal_use_lnL: + like_to_integrate = zero_like + else: + like_to_integrate = unit_like + if oracleRS: + print(" ORACLE - seeding sampling with" , oracleRS.params_ordered) + if hasattr(sampler, 'update_sampling_prior'): + rvs_train = {} + _, _, rv_oracle = oracleRS.draw_simplified(opts.n_chunk) + for indx,p in enumerate(sampler.params_ordered): + rvs_train[p] = rv_oracle[:,indx] + # train with equal weight - no likelihood information + lnL_oracles = np.zeros(opts.n_chunk) + sampler.update_sampling_prior(lnL_oracles, opts.n_chunk, external_rvs=rvs_train,log_scale_weights=True,floor_integrated_probability=opts.adapt_floor_level) + + res, var, neff, dict_return = sampler.integrate(like_to_integrate, *unpinned_params, **pinned_params) + + if not(res): # no resut + raise ValueError(" No integral result returned") + + if not(opts.internal_use_lnL): + log_res = numpy.log(res) + sqrt_var_over_res = numpy.sqrt(var)/res + else: + log_res = res + sqrt_var_over_res = numpy.exp(var/2 - log_res) + + # Report results + if opts.output_file: + fname_output_txt = opts.output_file +"_"+str(indx_event)+"_" + ".dat" + m1 =P.m1/lal.MSUN_SI + m2 =P.m2/lal.MSUN_SI + if opts.sim_xml: + event_id = opts.event + else: + 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 ]])) + + # Comprehensive output (not yet provided) + # Convert declination, inclination parameters in sampler if needed + if opts.save_samples and opts.output_file: + import copy + samples = copy.deepcopy(sampler._rvs) # deep copy: avoid modifying structures and having side effect on integrator, which loops over keys Expensive! + # Insert reference distance if it was marginalized over + if "distance" not in samples: + # Not distance output is the same as internal calculations: in *Mpc* + samples["distance"] = np.full_like( + samples["psi"], + factored_likelihood.distMpcRef, #*1e6*lal.PC_SI, + ) + # Insert reference phase if it was marginalized over + if "phi_orb" not in samples: + samples["phi_orb"] = np.full_like(samples["psi"], 0.) + if opts.inclination_cosine_sampler: + samples["inclination"] = numpy.arccos(samples["inclination"].astype(numpy.float64)) + if opts.declination_cosine_sampler: + samples["declination"] = numpy.pi/2 - numpy.arccos(samples["declination"].astype(numpy.float64)) + xmldoc = ligolw.Document() + xmldoc.appendChild(ligolw.LIGO_LW()) + process.register_to_xmldoc(xmldoc, sys.argv[0], opts.__dict__) + if not(opts.resample_time_marginalization): + if not opts.time_marginalization: + samples["t_ref"] += float(fiducial_epoch) + else: + samples["t_ref"] = float(fiducial_epoch)*numpy.ones(len(samples["psi"])) + # rotate sky to recover physical coordinates. Note on CPU + if opts.internal_sky_network_coordinates: + tmp_th, tmp_ph = my_rotation_cpu(np.pi/2 - samples['declination'],samples['right_ascension']) + samples['declination'] = np.pi/2 - tmp_th + samples['right_ascension'] = np.mod(tmp_ph, 2*np.pi) + + # recover true phase, polarization samples + if opts.internal_rotate_phase: + phi_orb_true = np.mod((samples['phi_orb'] + samples['psi'])/2., 2*np.pi) + psi_true = np.mod((samples['phi_orb'] - samples['psi'])/2., 2*np.pi) # keep as 2 pi range, to be consistent with past work + samples['psi']= psi_true + samples['phi_orb'] = phi_orb_true + samples["polarization"] = samples["psi"] + samples["coa_phase"] = samples["phi_orb"] + if ("declination", "right_ascension") in sampler.params: + samples["latitude"], samples["longitude"] = samples[("declination", "right_ascension")] + else: + samples["latitude"] = samples["declination"] + samples["longitude"] = samples["right_ascension"] + if "log_integrand" in samples: + samples["loglikelihood"] = samples["log_integrand"] + manual_avoid_overflow_logarithm + else: + samples["loglikelihood"] = numpy.log(samples["integrand"]) + manual_avoid_overflow_logarithm # export with consistent offset + if not opts.rom_integrate_intrinsic: + # ILE mode: insert fixed model parameters + samples["mass1"] = numpy.ones(samples["psi"].shape)*m1 # opts.mass1 + samples["mass2"] = numpy.ones(samples["psi"].shape)*m2 # opts.mass2 + samples["spin1x"] =numpy.ones(samples["psi"].shape)*P.s1x + samples["spin1y"] =numpy.ones(samples["psi"].shape)*P.s1y + samples["spin1z"] =numpy.ones(samples["psi"].shape)*P.s1z + samples["spin2x"] =numpy.ones(samples["psi"].shape)*P.s2x + samples["spin2y"] =numpy.ones(samples["psi"].shape)*P.s2y + samples["spin2z"] =numpy.ones(samples["psi"].shape)*P.s2z + samples["alpha4"] =numpy.ones(samples["psi"].shape)*P.eccentricity + samples["alpha5"] =numpy.ones(samples["psi"].shape)*P.lambda1 + samples["alpha6"] =numpy.ones(samples["psi"].shape)*P.lambda2 + # Below exist solely to placate XML export; new issue as of latest lalsuite say 7.15+ or so + samples["alpha2"] = numpy.zeros(samples["psi"].shape) + samples["alpha3"] = numpy.zeros(samples["psi"].shape) + + if opts.resample_time_marginalization: + samples = resample_samples_LISA(samples, rholms, cross_terms, lisa_sky_lamda, lisa_sky_beta, P, modes, reference_distance) + samples["loglikelihood" ] = samples["lnL_raw"] # export the non-time-marginalized likelihood, if we are in the final stages +# print(samples['t_ref'] - fiducial_epoch, len(samples['t_ref'])) + xmlutils.append_samples_to_xmldoc(xmldoc, samples) + # Extra metadata + dict_out={"mass1": m1, "mass2": m2, "spin1z": P.s1z, "spin2z": P.s2z, "alpha4": P.eccentricity, "alpha5":P.lambda1, "alpha6":P.lambda2, "event_duration": sqrt_var_over_res, "ttotal": sampler.ntotal} +# if 'distance' in pinned_params: +# dict_out['distance'] = pinned_params["distance"] + converged_result = False + if "convergence_test_results" in dict_return: + converged_result = dict_return["convergence_test_results"]["normal_integral"] + xmlutils.append_likelihood_result_to_xmldoc(xmldoc, log_res+manual_avoid_overflow_logarithm, neff=neff, converged=converged_result, **dict_out) + fname_output_xml = opts.output_file +"_"+str(indx_event)+"_" + ".xml.gz" + utils.write_filename(xmldoc, fname_output_xml, compress="gz") + + + + + if opts.maximize_only and opts.output_file: + # Pick the best extrinsic parameters, except for time (assumed not set: time marginalization) + if "integrand" in sampler._rvs: + indx_guess = numpy.argmax(sampler._rvs["integrand"]) # start search near maximum-likelihood point. (WARNING: can be very close by) + else: + indx_guess = numpy.argmax(sampler._rvs["log_integrand"]) + P.radec=True + P.phi = sampler._rvs["right_ascension"][indx_guess] + P.theta = sampler._rvs["declination"][indx_guess] + P.phiref = sampler._rvs["phi_orb"][indx_guess] + P.incl = sampler._rvs["inclination"][indx_guess] + P.psi = sampler._rvs["psi"][indx_guess] + P.dist = sampler._rvs["distance"][indx_guess]*1e6*lal.PC_SI + print( " ---- Best extrinsic paramers in MC ---- ") + P.print_params() + lalsimutils.ChooseWaveformParams_array_to_xml([P],"notime_raw_maxpt_"+opts.output_file) # best point, not including time + + import scipy.optimize + def fn_scaled(x): + P.phi = float(x[0]*2*numpy.pi) # right ascension + P.theta = float((x[1])*numpy.pi) # declination, really polar angle. NOT zero on equator + P.tref = fiducial_epoch + float((x[2]-0.5)*2*t_ref_wind) # ref. time (rel to epoch for data taking) + P.phiref = float(x[3]*numpy.pi*2) # ref. orbital phase + P.incl = float(x[4]*numpy.pi) # inclination + P.psi = float(x[5]*numpy.pi) # polarization angle + P.dist = x[6]*dmax* 1.e6 * lalsimutils.lsu_PC # luminosity distance + + return -1.0* factored_likelihood.FactoredLogLikelihood( + P, rholms,rholms_intp, cross_terms, cross_terms_V, opts.l_max) + + # Pick the best extrinsic parameters, except for time (assumed not set: time marginalization) + x0 =numpy.array( [ \ + sampler._rvs["right_ascension"][indx_guess]/(2*numpy.pi) , \ + (sampler._rvs["declination"][indx_guess]/numpy.pi), \ + 0.5, \ + (sampler._rvs["phi_orb"][indx_guess]/(2*numpy.pi)), \ + sampler._rvs["inclination"][indx_guess]/(numpy.pi),\ + sampler._rvs["psi"][indx_guess]/numpy.pi,\ + sampler._rvs["distance"][indx_guess]/dmax\ + ],dtype=numpy.float128) + x0 = numpy.fmod(x0,numpy.ones(len(x0))) # had BETTER be defined on this range! + # Pick the best starting time. BRUTE FORCE METHOD: use grid + def fn_scaled_t(t,x0): + return fn_scaled( [x0[0], x0[1], t, x0[3], x0[4], x0[5], x0[6]]) + npts_guess = int(t_ref_wind*2/(0.5*1e-4)) # Need to have enough points to fully explore the peak, timing to sub-ms accuracy + print( " Using ", npts_guess, " time points to select the best time, fixing the remaining extrinsic parameters ") + tvals_scaled_guess = numpy.linspace(0,1,npts_guess) + lnLvals = numpy.array([-1*fn_scaled_t(t,x0) for t in tvals_scaled_guess]) # note the two -1's cancel + # from matplotlib import pyplot as plt; plt.plot(tvals_scaled_guess,lnLvals,'o'); plt.show(); numpy.savetxt("dump-lnL.dat", numpy.array([tvals_scaled_guess*2*t_ref_wind,lnLvals]).T) + tbest = tvals_scaled_guess[numpy.argmax(lnLvals)] + print( " ---- Best extrinsic paramers in MC, after time offset ---- ") + P.tref = fiducial_epoch + tbest + P.print_params() + lalsimutils.ChooseWaveformParams_array_to_xml([P],"withtime_raw_maxpt_"+opts.output_file) # best point, not including time + x0[2] = tbest + x0p = x0 + # Refine best starting time with a search + # t_scaled_est = scipy.optimize.brent(fn_scaled_t, brack=(tbest-0.01,tbest,tbest+0.01),args=(x0),maxiter=500) + # print "Scaled starting time estimate after, before ", t_scaled_est, tbest + # x0p[2] = t_scaled_est + t_best_now = lal.LIGOTimeGPS() # create correct data type + t_best_now = fiducial_epoch + float((x0p[2]-0.5)*2*t_ref_wind) + t_best_now_s = int(numpy.floor(t_best_now)) + t_best_now_ns = int((t_best_now - t_best_now_s)*1e9) + print( " Fitting best (geocentric) time : ", str(t_best_now_s)+ '.'+ str(t_best_now_ns) , " relative time ", (x0p[2]-0.5)*2*t_ref_wind) + x0 = x0p + # Full multi-d search for best point + print( " Starting point [dimensionless] ", x0) + print( " Starting point [physical] ", [x0[0]*2*numpy.pi, numpy.pi*(x0[1]), t_ref_wind*2*(x0[2]-0.5), 2*numpy.pi*x0[3], numpy.pi*x0[4], numpy.pi*x0[5], x0[6]*dmax]) + lnLstart = -1*fn_scaled(x0) # note the two -1's cancel. Compare to 'rho^2/2' value reported by code + print( " lnL at start :", lnLstart, " [note can be lower than peak because of time offset]: best reported by ILE (including weights) is ", numpy.log(numpy.max(sampler._rvs["integrand"]))) + x = scipy.optimize.fmin(fn_scaled,x0, xtol=1e-5,ftol=1e-3,maxiter=opts.n_max,maxfun=opts.n_max) + lnLmax = -1.0*fn_scaled(x) + if lnLmax < lnLstart: + print( " Maximization failed to improve our initial point ! ") + x = x0p + lnLmax = lnLstart + t_best_now = fiducial_epoch + float((x[2]-0.5)*2*t_ref_wind) + print( "Best lnL = ", lnLmax) + print( "Best time =", t_best_now) + print( "Best point", [x[0]*2*numpy.pi, numpy.pi*(x[1]), t_ref_wind*2*(x[2]-0.5), 2*numpy.pi*x[3], numpy.pi*x[4], numpy.pi*x[5], x[6]*dmax]) + # Output best fit. Note STANDARD output will occur as usual + P_max = P.manual_copy() + P_max.tref = t_best_now + P_max.phi = float(x[0]*2*numpy.pi) + P_max.theta = float(numpy.pi*(x[1])) + P_max.phiref =float(x[3]*numpy.pi*2) + P_max.incl = float(x[4]*numpy.pi) + P_max.psi = float(x[5]*numpy.pi) + P_max.dist = x[6]*dmax*1e6*lal.PC_SI + P_max.m1 = lal.MSUN_SI*m1 + P_max.m2 = lal.MSUN_SI*m2 + print( " Sanity check: log likelihood for this set is ", factored_likelihood.FactoredLogLikelihood( + P_max, rholms, rholms_intp, cross_terms, cross_terms_V, opts.l_max)) + print( " ---- Best extrinsic paramers after polishing ---- ") + P_max.print_params() + P_list = [P_max] + lalsimutils.ChooseWaveformParams_array_to_xml(P_list,"maxpt_"+opts.output_file) + + # Clear sampler _rvs, to avoid side effects when called again + for key in sampler._rvs.keys(): + sampler._rvs[key] = [] + + return res + +def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts, inv_spec_trunc_Q=inv_spec_trunc_Q, T_spec=T_spec): + nEvals=0 + P = P_list[indx_event] + # if pin-distance-to-sim, change the distance prior accordingly + if opts.pin_distance_to_sim: + pinned_params['distance']=P.dist/(1.e6 * lal.PC_SI) + + # Call external likleihood preparation if any + if supplemental_ln_likelhood_prep: + supplemental_ln_likelhood_prep(P=P,config=supplemental_ln_likelhood_parsed_ini) + + extra_waveform_kwargs = {} + extra_waveform_kwargs['fd_alignment_postevent_time'] = 2 # 2 seconds after merger for ChooseFDModes + if opts.internal_waveform_fd_L_frame: + extra_waveform_kwargs['fd_L_frame'] = True + if opts.internal_waveform_fd_no_condition: + extra_waveform_kwargs['no_condition'] = True + if opts.rom_group: + extra_waveform_kwargs['rom_taper_start'] = True # really for NRHyb3dq8 given discussion with Aasim for high-SNR sources, but valuable generally + if opts.use_gwsignal_lmax_nyquist: + extra_waveform_kwargs['lmax_nyquist'] = int(opts.use_gwsignal_lmax_nyquist) + if opts.internal_waveform_extra_lalsuite_args: + extra_args_dict = eval(opts.internal_waveform_extra_lalsuite_args) # should only do this once and for all, not in loop! + print(" Waveform interface: extra args passed ", extra_args_dict) + extra_waveform_kwargs['extra_waveform_args'] = extra_args_dict + # Precompute + t_window = 0.15 + rholms_intp, cross_terms, cross_terms_V, rholms, guess_snr, rest=factored_likelihood.PrecomputeLikelihoodTerms( + fiducial_epoch, t_window, P, data_dict, psd_dict, opts.l_max, fmax, + False, inv_spec_trunc_Q, T_spec, + NR_group=NR_template_group,NR_param=NR_template_param, + use_gwsignal=opts.use_gwsignal, + use_gwsignal_approx=opts.approximant, + use_external_EOB=opts.use_external_EOB,nr_lookup=opts.nr_lookup,nr_lookup_valid_groups=opts.nr_lookup_group,perturbative_extraction=opts.nr_perturbative_extraction,perturbative_extraction_full=opts.nr_perturbative_extraction_full,use_provided_strain=opts.nr_use_provided_strain,hybrid_use=opts.nr_hybrid_use,hybrid_method=opts.nr_hybrid_method,ROM_group=opts.rom_group,ROM_param=opts.rom_param,ROM_use_basis=opts.rom_use_basis,verbose=opts.verbose,quiet=not opts.verbose,ROM_limit_basis_size=opts.rom_limit_basis_size_to,no_memory=opts.no_memory,skip_interpolation=opts.vectorized, extra_waveform_kwargs=extra_waveform_kwargs) + + # skip nan ! Something horrible has happened + if np.isnan(guess_snr): + print(" --- NAN SNR GUESS, ABORTING THIS EVENT {} --".format(indx_event)) + raise Exception(" NAN SNR ") + + if opts.auto_logarithm_offset and guess_snr: + # important: this only impacts *this* analysis + # important: if we have a very loud signal, it is important to change this dynamically to avoid *underflow*, which sometimes happens otherwise if we fix the scale early on. + print(" : naive overflow protection: updating lnL overflow based on SNR guess of {} ".format(guess_snr)) + print(" : reminder, diagnostics for lnLmax below are offset by this amount! ") + # Note by changing the offset, we change how adapt-weight-exponent is being applied (because we've offset the range before squashing it!). Use with care: pipeline auto-sets the weight exponent + manual_avoid_overflow_logarithm = guess_snr**2/2 - 100 # more conservative than helper, if I were to auto-set it from --hint-snr. So the integral peak is more likely positive + # Don't scale so much that peak guess likelihood is negative: we do fine for low-ampltiude sources + if manual_avoid_overflow_logarithm < 0: + manual_avoid_overflow_logarithm = 0 + elif opts.auto_logarithm_offset: + print(" PROBLEM: guess_snr not being returned, but auto-tuning requested ") + # reset to default. Should not be needed, but weird python scoping error + manual_avoid_overflow_logarithm = manual_avoid_overflow_logarithm_default + else: + # reset to default. Should not be needed, but weird python scoping error + manual_avoid_overflow_logarithm = manual_avoid_overflow_logarithm_default + + if opts.vectorized: + lookupNKDict = {} + lookupKNDict={} + lookupKNconjDict={} + ctUArrayDict = {} + ctVArrayDict={} + rholmArrayDict={} + rholms_intpArrayDict={} + epochDict={} + for det in rholms_intp.keys(): + print( " Packing ", det) + lookupNKDict[det],lookupKNDict[det], lookupKNconjDict[det], ctUArrayDict[det], ctVArrayDict[det], rholmArrayDict[det], rholms_intpArrayDict[det], epochDict[det] = factored_likelihood.PackLikelihoodDataStructuresAsArrays( rholms[det].keys(), rholms_intp[det], rholms[det], cross_terms[det],cross_terms_V[det]) + if opts.gpu and (not xpy_default is np): + lookupNKDict[det] = cupy.asarray(lookupNKDict[det]) + rholmArrayDict[det] = cupy.asarray(rholmArrayDict[det]) + ctUArrayDict[det] = cupy.asarray(ctUArrayDict[det]) + ctVArrayDict[det] = cupy.asarray(ctVArrayDict[det]) + epochDict[det] = cupy.asarray(epochDict[det]) + + + + # Likelihood + if not opts.time_marginalization: + + def likelihood_function(right_ascension, declination, t_ref, phi_orb, + inclination, psi, distance): + + dec = numpy.copy(declination).astype(numpy.float64) + if opts.declination_cosine_sampler: + dec = numpy.pi/2 - numpy.arccos(dec) + incl = numpy.copy(inclination).astype(numpy.float64) + if opts.inclination_cosine_sampler: + incl = numpy.arccos(incl) + + # use EXTREMELY many bits + lnL = numpy.zeros(right_ascension.shape,dtype=numpy.float128) + i = 0 + for ph, th, tr, phr, ic, ps, di in zip(right_ascension, dec, + t_ref, phi_orb, incl, psi, distance): + P.phi = ph # right ascension + P.theta = th # declination + P.tref = fiducial_epoch + tr # ref. time (rel to epoch for data taking) + P.phiref = phr # ref. orbital phase + P.incl = ic # inclination + P.psi = ps # polarization angle + P.dist = di* 1.e6 * lalsimutils.lsu_PC # luminosity distance + + lnL[i] = factored_likelihood.FactoredLogLikelihood( + P, rholms_intp, cross_terms, cross_terms_V, opts.l_max) + i+=1 + if return_lnL: + return lnL - manual_avoid_overflow_logarithm + return numpy.exp(lnL - manual_avoid_overflow_logarithm) + + else: # Sum over time for every point in other extrinsic params + if not (opts.rom_integrate_intrinsic or opts.vectorized): + def likelihood_function(right_ascension, declination, phi_orb, inclination, + psi, distance): + dec = numpy.copy(declination).astype(numpy.float64) # get rid of 'object', and allocate space + if opts.declination_cosine_sampler: + dec = numpy.pi/2 - numpy.arccos(dec) + incl = numpy.copy(inclination).astype(numpy.float64) + if opts.inclination_cosine_sampler: + incl = numpy.arccos(incl) + + # use EXTREMELY many bits + lnL = numpy.zeros(right_ascension.shape,dtype=numpy.float128) + i = 0 + tvals = numpy.linspace(-t_ref_wind,t_ref_wind,int((t_ref_wind)*2/P.deltaT)) # choose an array at the target sampling rate. P is inherited globally + + for ph, th, phr, ic, ps, di in zip(right_ascension, dec, + phi_orb, inclination, psi, distance): + P.phi = ph # right ascension + P.theta = th # declination + P.tref = fiducial_epoch # see 'tvals', above + P.phiref = phr # ref. orbital phase + P.incl = ic # inclination + P.psi = ps # polarization angle + P.dist = di* 1.e6 * lalsimutils.lsu_PC # luminosity distance + + + lnL[i] = factored_likelihood.FactoredLogLikelihoodTimeMarginalized(tvals, + P, rholms_intp, rholms, cross_terms, cross_terms_V, + opts.l_max,interpolate=opts.interpolate_time) + i+=1 + if supplemental_ln_likelihood: + lnL += supplemental_ln_likelihood(right_ascension, declination, phi_orb,inclination, psi, distance) + if return_lnL: + return lnL -manual_avoid_overflow_logarithm + return numpy.exp(lnL -manual_avoid_overflow_logarithm) + + elif opts.vectorized: # use array-based multiplications, fewer for loops + if (not opts.gpu): + if opts.distance_marginalization: + print(" **Warning**: distance marginalization not being used, this is the old vectorized code path not the xpy path. Check your code path; you may want --force-xpy ") + def likelihood_function(right_ascension, declination, phi_orb, inclination, + psi, distance): +# global nEvals + tvals = numpy.linspace(-t_ref_wind,t_ref_wind,int((t_ref_wind)*2/P.deltaT)) # choose an array at the target sampling rate. P is inherited globally + dec = numpy.copy(declination).astype(numpy.float64) + if opts.declination_cosine_sampler: + dec = numpy.pi/2 - numpy.arccos(dec) + incl = numpy.copy(inclination).astype(numpy.float64) + if opts.inclination_cosine_sampler: + incl = numpy.arccos(incl) + + # use EXTREMELY many bits + lnL = numpy.zeros(right_ascension.shape,dtype=numpy.float128) + P.phi = right_ascension.astype(float) # cast to float + P.theta = dec #declination.astype(float) + P.tref = float(fiducial_epoch) + P.phiref = phi_orb.astype(float) + P.incl = incl #inclination.astype(float) + P.psi = psi.astype(float) + P.dist = (distance* 1.e6 * lalsimutils.lsu_PC).astype(float) # luminosity distance + + # rotate sky if needed + if opts.internal_sky_network_coordinates: + P.theta,P.phi = my_rotation(np.pi/2 - P.theta,P.phi) + P.theta = np.pi/2 - P.theta + P.phi = xpy_default.mod(P.phi, 2*np.pi) + + # rotate phase if needed + # make copies of arrays + # Sampling assumes P.phiref == phi+ \psi , P.psi == phi - psi + if opts.internal_rotate_phase: + phi_orb_true = (P.phiref + P.psi)/2. + psi_true = (P.phiref - P.psi)/2. + P.psi= psi_true + P.phiref = phi_orb_true + + lnL = factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVector(tvals, + P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max) +# nEvals +=len(right_ascension) + if supplemental_ln_likelihood: + lnL += supplemental_ln_likelihood(P.phi, P.theta, P.phiref ,P.incl, P.psi, P.dist) # use these variables so they are already float-type + if return_lnL: + return lnL -manual_avoid_overflow_logarithm + return numpy.exp(lnL-manual_avoid_overflow_logarithm) + else: # vectorized and gpu either available or being forced to use xpy code path + print( " Using CUDA GPU likelihood, if cupy available ") + if not opts.distance_marginalization: + def likelihood_function(right_ascension, declination, phi_orb, inclination, + psi, distance): +# global nEvals + tvals = xpy_default.linspace(-t_ref_wind,t_ref_wind,int((t_ref_wind)*2/P.deltaT)) # choose an array at the target sampling rate. P is inherited globally + P.phi = xpy_asarray_already(right_ascension) # cast to float + if opts.declination_cosine_sampler: + P.theta = numpy.pi/2 - xpy_default.arccos(xpy_asarray_already(declination)) + else: + P.theta = xpy_asarray_already(declination) + P.tref = float(fiducial_epoch) + P.phiref = xpy_asarray_already(phi_orb) + if opts.inclination_cosine_sampler: + P.incl = xpy_default.arccos(xpy_asarray_already(inclination)) + else: + P.incl = xpy_asarray_already(inclination) + P.psi = xpy_asarray_already(psi) + P.dist = xpy_asarray_already(distance* 1.e6 * lalsimutils.lsu_PC) # luminosity distance + + # rotate sky if needed + if opts.internal_sky_network_coordinates: + P.theta,P.phi = my_rotation(np.pi/2 - P.theta,P.phi) + P.theta = np.pi/2 - P.theta + P.phi = xpy_default.mod(P.phi, 2*np.pi) + + # rotate phase if needed + # make copies of arrays + # Sampling assumes P.phiref == phi+ \psi , P.psi == phi - psi + if opts.internal_rotate_phase: + phi_orb_true = (P.phiref + P.psi)/2. + psi_true = (P.phiref - P.psi)/2. + P.psi= psi_true + P.phiref = phi_orb_true + + + lnL = factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop(tvals, + P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max,xpy=xpy_default) +# nEvals +=len(right_ascension) + if supplemental_ln_likelihood: + lnL += supplemental_ln_likelihood(P.phi, P.theta, P.phiref ,P.incl, P.psi, P.dist,xpy=xpy_default) # use these variables so they are already float-type + if return_lnL: + return identity_convert_lnL(lnL -manual_avoid_overflow_logarithm) + return identity_convert_lnL(xpy_default.exp(lnL-manual_avoid_overflow_logarithm)) + else: + print( " Using direct distance marginalization ") + xmin = factored_likelihood.distMpcRef / dmax + xmax = factored_likelihood.distMpcRef / dmin + bmax = xpy_default.asarray(lookup_table["bmax"]) + sqrt_bmax = xpy_default.sqrt(bmax) + bref = xpy_default.asarray(lookup_table["bref"]) + s_array = xpy_default.asarray(lookup_table["s_array"]) + smin = s_array[0] + smax = s_array[-1] + t_array = xpy_default.asarray(lookup_table["t_array"]) + tmax = t_array[-1] + lnI_array = xpy_default.asarray(lookup_table["lnI_array"]) + + intp = EvenBivariateLinearInterpolator(s_array[0], s_array[1] - s_array[0], t_array[0], t_array[1] - t_array[0], lnI_array) + + def exponent_max(x0, b): + x0_expmax = xpy_default.clip(x0, a_min=xmin, a_max=xmax) + return b * x0_expmax * (x0 - 0.5*x0_expmax) + + def b_to_t(b): + # TODO: this function is duplicate with what is in util_InitMargTable +# return np.arcsinh(b / bref) + b_by_bref = b / bref + return xpy_default.arcsinh(b_by_bref, out=b_by_bref) + + + def x0_to_s(x0): + # TODO: this function is duplicate with what is in util_InitMargTable +# return np.arcsinh(np.sqrt(bmax) * (x0 - xmin)) - np.arcsinh(np.sqrt(bmax) * (xmax - x0)) + A = x0 - xmin + A *= sqrt_bmax + xpy_default.arcsinh(A, out=A) + + B = xmax - x0 + B *= sqrt_bmax + xpy_default.arcsinh(B, out=B) + + A -= B + return A + + def distmarg_loglikelihood(kappa_sq, rho_sq): + x0 = kappa_sq / rho_sq +# lnI = np.ones(shape=x0.shape) * -np.inf + lnI = xpy_default.full_like(x0, xpy_default.NINF) + s = x0_to_s(x0) + t = b_to_t(rho_sq) +# in_bounds = (s > smin) * (s < smax) * (t < tmax) + in_bounds = (s > smin) & (s < smax) & (t < tmax) + lnI[in_bounds] = intp(s[in_bounds], t[in_bounds]) + return exponent_max(x0, rho_sq) + lnI + + if lookup_table["phase_marginalization"]: + + print( " Using direct phase marginalization ") + for det in lookupNKDict: + if set((lm[0], lm[1]) for lm in lookupNKDict[det]) != {(2, 2), (2, -2)}: + raise Exception( + " Phase marginalization is implemented only for 2-2 modes, " + f"while the modes consired here are {lookupNKDict[det]}." + ) + + def likelihood_function(right_ascension, declination, inclination, psi): +# global nEvals + tvals = xpy_default.linspace(-t_ref_wind,t_ref_wind,int((t_ref_wind)*2/P.deltaT)) # choose an array at the target sampling rate. P is inherited globally + P.phi = xpy_default.asarray(right_ascension, dtype=np.float64) # cast to float + if opts.declination_cosine_sampler: + P.theta = numpy.pi/2 - xpy_default.arccos(xpy_default.asarray(declination,dtype=np.float64)) + else: + P.theta = xpy_default.asarray(declination,dtype=np.float64) + P.tref = float(fiducial_epoch) + P.phiref = xpy_default.full_like(inclination, 0., dtype=np.float64) + if opts.inclination_cosine_sampler: + P.incl = xpy_default.arccos(xpy_default.asarray(inclination, dtype=np.float64)) + else: + P.incl = xpy_default.asarray(inclination, dtype=np.float64) + P.psi = xpy_default.asarray(psi, dtype=np.float64) + P.dist = xpy_default.asarray(factored_likelihood.distMpcRef * 1.e6 * lalsimutils.lsu_PC,dtype=np.float64) # luminosity distance + + # rotate sky if needed + if opts.internal_sky_network_coordinates: + P.theta,P.phi = my_rotation(np.pi/2 - P.theta,P.phi) + P.theta = np.pi/2 - P.theta + P.phi = xpy_default.mod(P.phi, 2*np.pi) + + # rotate phase if needed + # make copies of arrays + # Sampling assumes P.phiref == phi+ \psi , P.psi == phi - psi + if opts.internal_rotate_phase: + phi_orb_true = (P.phiref + P.psi)/2. + psi_true = (P.phiref - P.psi)/2. + P.psi= psi_true + P.phiref = phi_orb_true + + lnL = factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop(tvals, + P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max,xpy=xpy_default, loglikelihood=distmarg_loglikelihood, phase_marginalization=True) +# nEvals +=len(right_ascension) + if supplemental_ln_likelihood: + lnL += supplemental_ln_likelihood(P.phi, P.theta, P.phiref ,P.incl, P.psi, 0,xpy=xpy_default) # Same API + if return_lnL: + return identity_convert_lnL(lnL -manual_avoid_overflow_logarithm) + return identity_convert_lnL(xpy_default.exp(lnL-manual_avoid_overflow_logarithm)) + + else: + + def likelihood_function(right_ascension, declination, phi_orb, inclination, psi): +# global nEvals + tvals = xpy_default.linspace(-t_ref_wind,t_ref_wind,int((t_ref_wind)*2/P.deltaT)) # choose an array at the target sampling rate. P is inherited globally + P.phi = xpy_default.asarray(right_ascension, dtype=np.float64) # cast to float + if opts.declination_cosine_sampler: + P.theta = numpy.pi/2 - xpy_default.arccos(xpy_default.asarray(declination,dtype=np.float64)) + else: + P.theta = xpy_default.asarray(declination,dtype=np.float64) + P.tref = float(fiducial_epoch) + P.phiref = xpy_default.asarray(phi_orb, dtype=np.float64) + if opts.inclination_cosine_sampler: + P.incl = xpy_default.arccos(xpy_default.asarray(inclination, dtype=np.float64)) + else: + P.incl = xpy_default.asarray(inclination, dtype=np.float64) + P.psi = xpy_default.asarray(psi, dtype=np.float64) + P.dist = xpy_default.asarray(factored_likelihood.distMpcRef * 1.e6 * lalsimutils.lsu_PC,dtype=np.float64) # luminosity distance + + # rotate sky if needed + if opts.internal_sky_network_coordinates: + P.theta,P.phi = my_rotation(np.pi/2 - P.theta,P.phi) + P.theta = np.pi/2 - P.theta + P.phi = xpy_default.mod(P.phi, 2*np.pi) + + # rotate phase if needed + # make copies of arrays + # Sampling assumes P.phiref == phi+ \psi , P.psi == phi - psi + if opts.internal_rotate_phase: + phi_orb_true = (P.phiref + P.psi)/2. + psi_true = (P.phiref - P.psi)/2. + P.psi= psi_true + P.phiref = phi_orb_true + + lnL = factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop(tvals, + P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max,xpy=xpy_default, loglikelihood=distmarg_loglikelihood) +# nEvals +=len(right_ascension) + if supplemental_ln_likelihood: + lnL += supplemental_ln_likelihood(P.phi, P.theta, P.phiref ,P.incl, P.psi, 0,xpy=xpy_default) # Same API + if return_lnL: + return identity_convert_lnL(lnL -manual_avoid_overflow_logarithm) + return identity_convert_lnL(xpy_default.exp(lnL-manual_avoid_overflow_logarithm)) + + else: # integrate over intrinsic variables. Right now those variables ahave HARDCODED NAMES, alas + def likelihood_function(right_ascension, declination, phi_orb, inclination, + psi, distance,q): + dec = numpy.copy(declination).astype(numpy.float64) + if opts.declination_cosine_sampler: + dec = numpy.pi/2 - numpy.arccos(dec) + incl = numpy.copy(inclination).astype(numpy.float64) + if opts.inclination_cosine_sampler: + incl = numpy.arccos(incl) + + lnL = numpy.zeros(len(right_ascension),dtype=numpy.float128) +# i = 0 + tvals = numpy.linspace(-t_ref_wind,t_ref_wind,int((t_ref_wind)*2/P.deltaT)) # choose an array at the target sampling rate. P is inherited globally + + +# t_start =lal.GPSTimeNow() + for ph, th, phr, ic, ps, di,qi in zip(right_ascension, dec, + phi_orb, incl, psi, distance,q): + # Reconstruct U,V using ROM fits. PROBABLY should do this once for every q, rather than deep on the loop + P.assign_param('q',qi) # mass ratio + rholms_intp_A, cross_terms_A, cross_terms_V_A, rholms_A, rest_A = factored_likelihood.ReconstructPrecomputedLikelihoodTermsROM(P, rest, rholms_intp, cross_terms, cross_terms_V, rholms,verbose=False) + # proceed for rest + P.phi = ph # right ascension + P.theta = th # declination + P.tref = fiducial_epoch # see 'tvals', above + P.phiref = phr # ref. orbital phase + P.incl = ic # inclination + P.psi = ps # polarization angle + P.dist = di* 1.e6 * lalsimutils.lsu_PC # luminosity distance + + + lnL[i] = factored_likelihood.FactoredLogLikelihoodTimeMarginalized(tvals, + P, rholms_intp_A, rholms_A, cross_terms_A, cross_terms_V_A, + opts.l_max,interpolate=opts.interpolate_time) + if numpy.isnan(lnL[i]) or lnL[i]<-200: + lnL[i] = -200 # regularize : a hack, for now, to deal with rare ROM problems. Only on the ROM logic fork + i+=1 +# t_end =lal.GPSTimeNow() +# print " Cost per evaluation ", (t_end - t_start)/len(q) +# print " Max lnL for this iteration ", numpy.max(lnL) + if return_lnL: + return lnL -manual_avoid_overflow_logarithm + return numpy.exp(lnL - manual_avoid_overflow_logarithm) + + if opts.sampler_method == "adaptive_cartesian_gpu": + # reset sampling parameter as needed (distance, inclination but not sky location) + # distance (and inclination) can conceivably correlate with mass, and we don't want to truncate in distance/inclination prematurely from history effects + # method ONLY implemented for adaptcart! + if 'distance' in sampler.params: + sampler.reset_sampling('distance') + sampler.reset_sampling('inclination') + elif opts.sampler_method == "GMM": + if 'distance' in sampler.params: + pair_d_incl = sampler_param_tuple(sampler, ['distance','inclination']) + if pair_d_incl in gmm_dict: + gmm_dict[pair_d_incl] = None # reset inclination/distance sampling + single_incl = sampler_param_tuple(sampler, ['inclination']) + if (opts.distance_marginalization) and single_incl in gmm_dict: + gmm_dict[single_incl] = None # reset inclination sampling + + + # Integrate + args = likelihood_function.__code__.co_varnames[:likelihood_function.__code__.co_argcount] + print( " --------> Arguments ", args) + # if exactly zero likelihood function + like_to_integrate = likelihood_function + if opts.zero_likelihood: + if opts.internal_use_lnL: + like_to_integrate = zero_like + else: + like_to_integrate = unit_like + if oracleRS: + print(" ORACLE - seeding sampling with" , oracleRS.params_ordered) + if hasattr(sampler, 'update_sampling_prior'): + rvs_train = {} + _, _, rv_oracle = oracleRS.draw_simplified(opts.n_chunk) + for indx,p in enumerate(sampler.params_ordered): + rvs_train[p] = rv_oracle[:,indx] + # train with equal weight - no likelihood information + lnL_oracles = np.zeros(opts.n_chunk) + sampler.update_sampling_prior(lnL_oracles, opts.n_chunk, external_rvs=rvs_train,log_scale_weights=True,floor_integrated_probability=opts.adapt_floor_level) + + res, var, neff, dict_return = sampler.integrate(like_to_integrate, *unpinned_params, **pinned_params) + + if not(res): # no resut + raise ValueError(" No integral result returned") + + if not(opts.internal_use_lnL): + log_res = numpy.log(res) + sqrt_var_over_res = numpy.sqrt(var)/res + else: + log_res = res + sqrt_var_over_res = numpy.exp(var/2 - log_res) + + # Report results + if opts.output_file: + fname_output_txt = opts.output_file +"_"+str(indx_event)+"_" + ".dat" + m1 =P.m1/lal.MSUN_SI + m2 =P.m2/lal.MSUN_SI + if opts.sim_xml: + event_id = opts.event + else: + event_id = -1 + if opts.event == None: + event_id = -1 + if opts.save_eccentricity: + # output format when eccentricity is being used + numpy.savetxt(fname_output_txt, numpy.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y, P.s2z, P.eccentricity, log_res+manual_avoid_overflow_logarithm, sqrt_var_over_res,sampler.ntotal, neff ]])) #dict_return["convergence_test_results"]["normal_integral]" + elif not (P.lambda1>0 or P.lambda2>0): + # output format when lambda is NOT used + if not opts.pin_distance_to_sim: + numpy.savetxt(fname_output_txt, numpy.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y, P.s2z, log_res+manual_avoid_overflow_logarithm, sqrt_var_over_res,sampler.ntotal, neff ]])) #dict_return["convergence_test_results"]["normal_integral]" + else: + numpy.savetxt(fname_output_txt, numpy.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y, P.s2z, pinned_params["distance"], log_res+manual_avoid_overflow_logarithm, sqrt_var_over_res,sampler.ntotal, neff ]])) #dict_return["convergence_test_results"]["normal_integral]" + else: + if not(opts.export_eos_index): + # Alternative output format if lambda is active + numpy.savetxt(fname_output_txt, numpy.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y, P.s2z, P.lambda1, P.lambda2, log_res+manual_avoid_overflow_logarithm, sqrt_var_over_res,sampler.ntotal, neff ]])) #dict_return["convergence_test_results"]["normal_integral]" + else: + numpy.savetxt(fname_output_txt, numpy.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y, P.s2z, P.lambda1, P.lambda2, P.eos_table_index, log_res+manual_avoid_overflow_logarithm, sqrt_var_over_res,sampler.ntotal, neff ]])) #dict_return["convergence_test_results"]["normal_integral]" + + + # Comprehensive output (not yet provided) + # Convert declination, inclination parameters in sampler if needed + if opts.save_samples and opts.output_file: + import copy + samples = copy.deepcopy(sampler._rvs) # deep copy: avoid modifying structures and having side effect on integrator, which loops over keys Expensive! + # Insert reference distance if it was marginalized over + if "distance" not in samples: + # Not distance output is the same as internal calculations: in *Mpc* + samples["distance"] = np.full_like( + samples["psi"], + factored_likelihood.distMpcRef, #*1e6*lal.PC_SI, + ) + # Insert reference phase if it was marginalized over + if "phi_orb" not in samples: + samples["phi_orb"] = np.full_like(samples["psi"], 0.) + if opts.inclination_cosine_sampler: + samples["inclination"] = numpy.arccos(samples["inclination"].astype(numpy.float64)) + if opts.declination_cosine_sampler: + samples["declination"] = numpy.pi/2 - numpy.arccos(samples["declination"].astype(numpy.float64)) + xmldoc = ligolw.Document() + xmldoc.appendChild(ligolw.LIGO_LW()) + process.register_to_xmldoc(xmldoc, sys.argv[0], opts.__dict__) + if not(opts.resample_time_marginalization): + if not opts.time_marginalization: + samples["t_ref"] += float(fiducial_epoch) + else: + samples["t_ref"] = float(fiducial_epoch)*numpy.ones(len(samples["psi"])) + # rotate sky to recover physical coordinates. Note on CPU + if opts.internal_sky_network_coordinates: + tmp_th, tmp_ph = my_rotation_cpu(np.pi/2 - samples['declination'],samples['right_ascension']) + samples['declination'] = np.pi/2 - tmp_th + samples['right_ascension'] = np.mod(tmp_ph, 2*np.pi) + + # recover true phase, polarization samples + if opts.internal_rotate_phase: + phi_orb_true = np.mod((samples['phi_orb'] + samples['psi'])/2., 2*np.pi) + psi_true = np.mod((samples['phi_orb'] - samples['psi'])/2., 2*np.pi) # keep as 2 pi range, to be consistent with past work + samples['psi']= psi_true + samples['phi_orb'] = phi_orb_true + samples["polarization"] = samples["psi"] + samples["coa_phase"] = samples["phi_orb"] + if ("declination", "right_ascension") in sampler.params: + samples["latitude"], samples["longitude"] = samples[("declination", "right_ascension")] + else: + samples["latitude"] = samples["declination"] + samples["longitude"] = samples["right_ascension"] + if "log_integrand" in samples: + samples["loglikelihood"] = samples["log_integrand"] + manual_avoid_overflow_logarithm + else: + samples["loglikelihood"] = numpy.log(samples["integrand"]) + manual_avoid_overflow_logarithm # export with consistent offset + if not opts.rom_integrate_intrinsic: + # ILE mode: insert fixed model parameters + samples["mass1"] = numpy.ones(samples["psi"].shape)*m1 # opts.mass1 + samples["mass2"] = numpy.ones(samples["psi"].shape)*m2 # opts.mass2 + samples["spin1x"] =numpy.ones(samples["psi"].shape)*P.s1x + samples["spin1y"] =numpy.ones(samples["psi"].shape)*P.s1y + samples["spin1z"] =numpy.ones(samples["psi"].shape)*P.s1z + samples["spin2x"] =numpy.ones(samples["psi"].shape)*P.s2x + samples["spin2y"] =numpy.ones(samples["psi"].shape)*P.s2y + samples["spin2z"] =numpy.ones(samples["psi"].shape)*P.s2z + samples["alpha4"] =numpy.ones(samples["psi"].shape)*P.eccentricity + samples["alpha5"] =numpy.ones(samples["psi"].shape)*P.lambda1 + samples["alpha6"] =numpy.ones(samples["psi"].shape)*P.lambda2 + # Below exist solely to placate XML export; new issue as of latest lalsuite say 7.15+ or so + samples["alpha2"] = numpy.zeros(samples["psi"].shape) + samples["alpha3"] = numpy.zeros(samples["psi"].shape) + + if opts.resample_time_marginalization: + samples = resample_samples(samples,lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict) + samples["loglikelihood" ] = samples["lnL_raw"] # export the non-time-marginalized likelihood, if we are in the final stages +# print(samples['t_ref'] - fiducial_epoch, len(samples['t_ref'])) + xmlutils.append_samples_to_xmldoc(xmldoc, samples) + # Extra metadata + dict_out={"mass1": m1, "mass2": m2, "spin1z": P.s1z, "spin2z": P.s2z, "alpha4": P.eccentricity, "alpha5":P.lambda1, "alpha6":P.lambda2, "event_duration": sqrt_var_over_res, "ttotal": sampler.ntotal} +# if 'distance' in pinned_params: +# dict_out['distance'] = pinned_params["distance"] + converged_result = False + if "convergence_test_results" in dict_return: + converged_result = dict_return["convergence_test_results"]["normal_integral"] + xmlutils.append_likelihood_result_to_xmldoc(xmldoc, log_res+manual_avoid_overflow_logarithm, neff=neff, converged=converged_result, **dict_out) + fname_output_xml = opts.output_file +"_"+str(indx_event)+"_" + ".xml.gz" + utils.write_filename(xmldoc, fname_output_xml, compress="gz") + + + + + if opts.maximize_only and opts.output_file: + # Pick the best extrinsic parameters, except for time (assumed not set: time marginalization) + if "integrand" in sampler._rvs: + indx_guess = numpy.argmax(sampler._rvs["integrand"]) # start search near maximum-likelihood point. (WARNING: can be very close by) + else: + indx_guess = numpy.argmax(sampler._rvs["log_integrand"]) + P.radec=True + P.phi = sampler._rvs["right_ascension"][indx_guess] + P.theta = sampler._rvs["declination"][indx_guess] + P.phiref = sampler._rvs["phi_orb"][indx_guess] + P.incl = sampler._rvs["inclination"][indx_guess] + P.psi = sampler._rvs["psi"][indx_guess] + P.dist = sampler._rvs["distance"][indx_guess]*1e6*lal.PC_SI + print( " ---- Best extrinsic paramers in MC ---- ") + P.print_params() + lalsimutils.ChooseWaveformParams_array_to_xml([P],"notime_raw_maxpt_"+opts.output_file) # best point, not including time + + import scipy.optimize + def fn_scaled(x): + P.phi = float(x[0]*2*numpy.pi) # right ascension + P.theta = float((x[1])*numpy.pi) # declination, really polar angle. NOT zero on equator + P.tref = fiducial_epoch + float((x[2]-0.5)*2*t_ref_wind) # ref. time (rel to epoch for data taking) + P.phiref = float(x[3]*numpy.pi*2) # ref. orbital phase + P.incl = float(x[4]*numpy.pi) # inclination + P.psi = float(x[5]*numpy.pi) # polarization angle + P.dist = x[6]*dmax* 1.e6 * lalsimutils.lsu_PC # luminosity distance + + return -1.0* factored_likelihood.FactoredLogLikelihood( + P, rholms,rholms_intp, cross_terms, cross_terms_V, opts.l_max) + + # Pick the best extrinsic parameters, except for time (assumed not set: time marginalization) + x0 =numpy.array( [ \ + sampler._rvs["right_ascension"][indx_guess]/(2*numpy.pi) , \ + (sampler._rvs["declination"][indx_guess]/numpy.pi), \ + 0.5, \ + (sampler._rvs["phi_orb"][indx_guess]/(2*numpy.pi)), \ + sampler._rvs["inclination"][indx_guess]/(numpy.pi),\ + sampler._rvs["psi"][indx_guess]/numpy.pi,\ + sampler._rvs["distance"][indx_guess]/dmax\ + ],dtype=numpy.float128) + x0 = numpy.fmod(x0,numpy.ones(len(x0))) # had BETTER be defined on this range! + # Pick the best starting time. BRUTE FORCE METHOD: use grid + def fn_scaled_t(t,x0): + return fn_scaled( [x0[0], x0[1], t, x0[3], x0[4], x0[5], x0[6]]) + npts_guess = int(t_ref_wind*2/(0.5*1e-4)) # Need to have enough points to fully explore the peak, timing to sub-ms accuracy + print( " Using ", npts_guess, " time points to select the best time, fixing the remaining extrinsic parameters ") + tvals_scaled_guess = numpy.linspace(0,1,npts_guess) + lnLvals = numpy.array([-1*fn_scaled_t(t,x0) for t in tvals_scaled_guess]) # note the two -1's cancel + # from matplotlib import pyplot as plt; plt.plot(tvals_scaled_guess,lnLvals,'o'); plt.show(); numpy.savetxt("dump-lnL.dat", numpy.array([tvals_scaled_guess*2*t_ref_wind,lnLvals]).T) + tbest = tvals_scaled_guess[numpy.argmax(lnLvals)] + print( " ---- Best extrinsic paramers in MC, after time offset ---- ") + P.tref = fiducial_epoch + tbest + P.print_params() + lalsimutils.ChooseWaveformParams_array_to_xml([P],"withtime_raw_maxpt_"+opts.output_file) # best point, not including time + x0[2] = tbest + x0p = x0 + # Refine best starting time with a search + # t_scaled_est = scipy.optimize.brent(fn_scaled_t, brack=(tbest-0.01,tbest,tbest+0.01),args=(x0),maxiter=500) + # print "Scaled starting time estimate after, before ", t_scaled_est, tbest + # x0p[2] = t_scaled_est + t_best_now = lal.LIGOTimeGPS() # create correct data type + t_best_now = fiducial_epoch + float((x0p[2]-0.5)*2*t_ref_wind) + t_best_now_s = int(numpy.floor(t_best_now)) + t_best_now_ns = int((t_best_now - t_best_now_s)*1e9) + print( " Fitting best (geocentric) time : ", str(t_best_now_s)+ '.'+ str(t_best_now_ns) , " relative time ", (x0p[2]-0.5)*2*t_ref_wind) + x0 = x0p + # Full multi-d search for best point + print( " Starting point [dimensionless] ", x0) + print( " Starting point [physical] ", [x0[0]*2*numpy.pi, numpy.pi*(x0[1]), t_ref_wind*2*(x0[2]-0.5), 2*numpy.pi*x0[3], numpy.pi*x0[4], numpy.pi*x0[5], x0[6]*dmax]) + lnLstart = -1*fn_scaled(x0) # note the two -1's cancel. Compare to 'rho^2/2' value reported by code + print( " lnL at start :", lnLstart, " [note can be lower than peak because of time offset]: best reported by ILE (including weights) is ", numpy.log(numpy.max(sampler._rvs["integrand"]))) + x = scipy.optimize.fmin(fn_scaled,x0, xtol=1e-5,ftol=1e-3,maxiter=opts.n_max,maxfun=opts.n_max) + lnLmax = -1.0*fn_scaled(x) + if lnLmax < lnLstart: + print( " Maximization failed to improve our initial point ! ") + x = x0p + lnLmax = lnLstart + t_best_now = fiducial_epoch + float((x[2]-0.5)*2*t_ref_wind) + print( "Best lnL = ", lnLmax) + print( "Best time =", t_best_now) + print( "Best point", [x[0]*2*numpy.pi, numpy.pi*(x[1]), t_ref_wind*2*(x[2]-0.5), 2*numpy.pi*x[3], numpy.pi*x[4], numpy.pi*x[5], x[6]*dmax]) + # Output best fit. Note STANDARD output will occur as usual + P_max = P.manual_copy() + P_max.tref = t_best_now + P_max.phi = float(x[0]*2*numpy.pi) + P_max.theta = float(numpy.pi*(x[1])) + P_max.phiref =float(x[3]*numpy.pi*2) + P_max.incl = float(x[4]*numpy.pi) + P_max.psi = float(x[5]*numpy.pi) + P_max.dist = x[6]*dmax*1e6*lal.PC_SI + P_max.m1 = lal.MSUN_SI*m1 + P_max.m2 = lal.MSUN_SI*m2 + print( " Sanity check: log likelihood for this set is ", factored_likelihood.FactoredLogLikelihood( + P_max, rholms, rholms_intp, cross_terms, cross_terms_V, opts.l_max)) + print( " ---- Best extrinsic paramers after polishing ---- ") + P_max.print_params() + P_list = [P_max] + lalsimutils.ChooseWaveformParams_array_to_xml(P_list,"maxpt_"+opts.output_file) + + # Clear sampler _rvs, to avoid side effects when called again + for key in sampler._rvs.keys(): + sampler._rvs[key] = [] + + return res + + +lnL_sofar = -np.inf +no_adapt_sky = False +for indx in numpy.arange(len(P_list)): + try: + if opts.LISA: + res = analyze_event_LISA(P_list, indx, data_dict, psd_dict, fmax, opts) + else: + res = analyze_event(P_list, indx, data_dict, psd_dict, fmax, opts) + # abort if horrible (nan event) - done with 'raise' + lnL_sofar = np.max([lnL_sofar,res]) + if opts.force_reset_all: # depends on integrator! May not always be availble + if opts.sampler_method == "adaptive_cartesian_gpu": + for name in sampler.params: + sampler.reset_sampling(name) + elif opts.sampler_method == "GMM": + # reset the GMM dictionary + for component in gmm_dict: + gmm_dict[component] = None + elif opts.sampler_method == 'AV': + print(" AV always resets every iteration ! ") + else: + print(" force-reset-all not defined for this integrator ") + sys.exit(1) + if opts.no_adapt_after_first and (not no_adapt_sky): + if lnL_sofar > 20: # Use absolute threshold. Expect this will give modest sky localization. + # remove right_ascension, declination from adaptive parameters + params_adapt = sampler.adaptive + params_adapt = list(set(params_adapt) - set(['right_ascension','declination'])) + sampler.adaptive = params_adapt + # Disable saving of the integrand -- saves on memory and will reduce speed. + # Note this needs to be done in a few places + pinned_params.update({"force_no_adapt":True,"save_intg":False, "igrand_threshold_deltalnL":20}) # massively reduce memory usage in logic branch, don't save all. Highly redundant sequence to self-document all related logic branches + except Exception as exception_failure: + print(" ===> FAILED ANALYSIS <==== ") + print(exception_failure) + if opts.internal_make_empty_file_on_error: + fname_output_txt = opts.output_file +"_"+str(indx)+"_" + ".dat" + open(fname_output_txt,'a').close() # create empty file + if opts.internal_hard_fail_on_error: + sys.exit(1) +# if len(exception_failure) >0: +# if "CUBLAS" in exception_failure[0]: # Hard fail if a cuda error! +# sys.exit(1) + if ("CUDA" in str(exception_failure)) or ('CUBLAS' in str(exception_failure)) or ('cuda' in str(exception_failure)) or ('compilation' in (str(exception_failure)) ): # Hard fail if a cuda error with a catchable error code + if (opts.internal_soft_fail_on_cuda_error): + sys.exit(0) + sys.exit(62) + if 'Out of memory' in str(exception_failure): # should never happen, most likely a crappy node or failure to set correct memory limit for ILE. + sys.exit(63) + str_err = " {} ".format(exception_failure) + if ('Zero prior failure' in str_err) or ('effective samples' in str_err): + # common failures that require reset of sampler + # - zero prior : very rare case where the prior is exactly zero for some reason. User error most likely (floors, etc). Should never happen + # - effective samples = nan : error with AC where the integrator gets confused. Reset + # reset all variables sampling, so we don't contaminate subsequent versions + # - again, this should only be a problem if we are using multiply adaptive sampling, so it should NEVER happen on the first time through + for name in sampler.params: + sampler.reset_sampling(param) + + #a sketch of how I would do custom failure modes, but for sake of speed I'm just setting the above right now + #for i,failure_mode in enumerate(opts.custom_fails): + # if failure_mode in str(exception_failure): + # sys.exit(opts.custom_fail_codes[i]) + + print( " Probable reasons: SEOB nyquist or starting frequency limit or signal duration ") + print( " Skipping the following binary! ") + # Zero out extrinsic parameters -- these are CUDA-populated / meaningless, but could cause errors if populated + P_list[indx].incl = P_list[indx].tref = P_list[indx].dist = P_list[indx].phiref = P_list[indx].psi =P_list[indx].theta = P_list[indx].phi =0 + P_list[indx].print_params() + diff --git a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py index 7afac2f6c..835eafde1 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py @@ -28,6 +28,28 @@ else: assume_lowlatency=False +# ---------------------------------------------------------------------- +# Hyperpipeline ASCII grid format (opt-in via env var). When set, every +# grid file pseudo_pipe touches (target_params, proposed-grid, the +# --input-grid handed to BasicIteration, and the --sim-xml command-single +# sanity-check invocation) is .dat instead of .xml.gz, and the +# downstream pipeline runs in hyperpipeline mode end-to-end. When unset, +# behaviour is identical to the legacy XML pipeline. +# +# By design, pseudo_pipe does NOT convert formats internally -- the +# entire process operates cohesively in one mode or the other. In +# hyperpipeline mode the user is responsible for staging any external +# inputs (e.g. --manual-initial-grid) as hyperpipeline .dat; the +# auto-generated AMR / template-bank seed-grid paths still emit XML and +# will fail downstream unless --manual-initial-grid is supplied. +# ---------------------------------------------------------------------- +_use_hpip_pp = str(os.environ.get("RIFT_HYPERPIPELINE_FORMAT", "")).strip().lower() in ("1", "true", "yes", "on") +grid_suffix_pp = "dat" if _use_hpip_pp else "xml.gz" +sim_grid_flag_pp = "--sim-grid" if _use_hpip_pp else "--sim-xml" +if _use_hpip_pp: + print(" === pseudo_pipe: hyperpipeline ASCII grid format active (RIFT_HYPERPIPELINE_FORMAT) ===") + print(" Inter-stage grids will be .{}, command-single will use {}".format(grid_suffix_pp, sim_grid_flag_pp)) + # Backward compatibility from RIFT.misc.dag_utils_generic import which ligolw_prefix = 'igwn_' @@ -320,7 +342,7 @@ def unsafe_parse_arg_string_dict(my_argstr): parser.add_argument("--export-distance-slices-skip-threshold",default=None,type=float,help="Passthrough: --distance-slice-skip-threshold for the .dslice export (absolute peak-lnL detectability cut).") parser.add_argument("--ile-additional-files-to-transfer",default=None,help="Comma-separated list of filenames. To append to the transfer file list for ILE jobs (only). Intended for surrogates in LAL_DATA_PATH for wide-ranging use") parser.add_argument("--internal-cip-use-lnL",action='store_true') -parser.add_argument("--manual-initial-grid",default=None,type=str,help="Filename (full path) to initial grid. Copied into proposed-grid.xml.gz, overwriting any grid assignment done here") +parser.add_argument("--manual-initial-grid",default=None,type=str,help="Filename (full path) to initial grid. Copied into proposed-grid., overwriting any grid assignment done here. Suffix is .xml.gz by default and .dat when RIFT_HYPERPIPELINE_FORMAT is set; the source file's format must match the active mode.") parser.add_argument("--manual-initial-grid-supplements",action='store_true', help="Manual inital grid used to SUPPLEMENT output of the default helper grid.") parser.add_argument("--manual-extra-ile-args",default=None,type=str,help="Avenue to adjoin extra ILE arguments. Needed for unusual configurations (e.g., if channel names are not being selected, etc)") parser.add_argument("--internal-ile-force-adapt-all",action='store_true', help="Syntactic sugar to prevent need to add manual-extra-ile-args for this: easier on user") @@ -660,8 +682,17 @@ def unsafe_parse_arg_string_dict(my_argstr): P.eccentricity = event_dict["eccentricity"] if not(event_dict['meanPerAno'] is None): P.meanPerAno = event_dict["meanPerAno"] - # Write 'target_params.xml.gz' file - lalsimutils.ChooseWaveformParams_array_to_xml([P], "target_params") + # Write 'target_params' file -- hyperpipeline .dat or legacy XML. + if _use_hpip_pp: + from RIFT.misc import hyperpipeline_io as _hpio + _cols = _hpio.build_column_list( + use_eccentricity=(P.eccentricity != 0), + use_meanPerAno=(P.meanPerAno != 0)) + _hpio.write_grid_from_P_list("target_params", [P], _cols, + lal_module=lal, + lalsimutils_module=lalsimutils) + else: + lalsimutils.ChooseWaveformParams_array_to_xml([P], "target_params") if opts.use_production_defaults: # use more workers for high-q triggers @@ -852,7 +883,7 @@ def unsafe_parse_arg_string_dict(my_argstr): cmd += " --use-cvmfs-frames " # only run with CVMFS data, otherwise very very painful if opts.use_ini: cmd += " --use-ini " + opts.use_ini - cmd += " --sim-xml {}/target_params.xml.gz --event 0 ".format(base_dir + "/"+ dirname_run) # full path to target_params.xml.gz + cmd += " {} {}/target_params.{} --event 0 ".format(sim_grid_flag_pp, base_dir + "/"+ dirname_run, grid_suffix_pp) # full path to target_params (xml.gz or .dat) if (opts.event_time is None): cmd += " --event-time " + str(event_dict["tref"]) # @@ -1522,11 +1553,23 @@ def unsafe_parse_arg_string_dict(my_argstr): # Overwrite grid if needed if not (opts.manual_initial_grid is None): if opts.manual_initial_grid_supplements: + if _use_hpip_pp: + # ligolw_add is XML-only -- the equivalent for hyperpipeline is + # the head-line + cat | sort | uniq | shuf shell pattern from + # the BasicIteration join_grids.sh. We refuse rather than + # silently produce a broken proposed-grid. + raise SystemExit( + "pseudo_pipe: --manual-initial-grid-supplements is XML-only " + "(uses ligolw_add); incompatible with RIFT_HYPERPIPELINE_FORMAT. " + "Pre-merge your supplements into a single .dat grid and pass it " + "via --manual-initial-grid instead.") cmd_add = '{}ligolw_add {} proposed-grid.xml.gz --output tmp.xml.gz'.format(ligolw_prefix,opts.manual_initial_grid) os.system(cmd_add) shutil.copyfile('tmp.xml.gz', "proposed-grid.xml.gz") else: - shutil.copyfile(opts.manual_initial_grid, "proposed-grid.xml.gz") + # shutil.copyfile is format-agnostic: works for either .xml.gz or + # .dat as long as the source matches the active mode. + shutil.copyfile(opts.manual_initial_grid, "proposed-grid." + grid_suffix_pp) # override npts_it if needed if opts.internal_n_evaluations_per_iteration: @@ -1551,7 +1594,7 @@ def unsafe_parse_arg_string_dict(my_argstr): print(" WARNING: --pipeline-builder {} overrides --use-subdags routing; AMR/subdag runs require AlternateIteration ".format(opts.pipeline_builder)) cepp = "create_event_parameter_pipeline_" + opts.pipeline_builder print(" Pipeline builder (create_event_parameter_pipeline_*): ", cepp) -cmd =cepp+ " --ile-n-events-to-analyze {} --input-grid proposed-grid.xml.gz --ile-exe `which integrate_likelihood_extrinsic_batchmode` --ile-args `pwd`/args_ile.txt --cip-args-list args_cip_list.txt --test-args args_test.txt --request-memory-CIP {} --request-memory-ILE {} --n-samples-per-job ".format(n_jobs_per_worker,cip_mem,ile_mem) + str(npts_it) + " --working-directory `pwd` --n-iterations " + str(n_iterations) + " --n-iterations-subdag-max {} ".format(opts.internal_n_iterations_subdag_max) + " --n-copies {} ".format(opts.ile_copies) + " --ile-retries "+ str(opts.ile_retries) + " --general-retries " + str(opts.general_retries) +cmd =cepp+ " --ile-n-events-to-analyze {} --input-grid proposed-grid.{} --ile-exe `which integrate_likelihood_extrinsic_batchmode` --ile-args `pwd`/args_ile.txt --cip-args-list args_cip_list.txt --test-args args_test.txt --request-memory-CIP {} --request-memory-ILE {} --n-samples-per-job ".format(n_jobs_per_worker,grid_suffix_pp,cip_mem,ile_mem) + str(npts_it) + " --working-directory `pwd` --n-iterations " + str(n_iterations) + " --n-iterations-subdag-max {} ".format(opts.internal_n_iterations_subdag_max) + " --n-copies {} ".format(opts.ile_copies) + " --ile-retries "+ str(opts.ile_retries) + " --general-retries " + str(opts.general_retries) if opts.ile_jobs_per_worker_first: cmd += " --ile-n-events-to-analyze-first {} ".format(opts.ile_jobs_per_worker_first) if opts.assume_matter or opts.assume_eccentric: @@ -1603,6 +1646,18 @@ def unsafe_parse_arg_string_dict(my_argstr): cmd += " --cip-exe-G `which util_ConstructIntrinsicPosterior_GaussianResampling.py ` " if opts.internal_use_amr: print(" AMR prototype: Using hardcoded aligned-spin settings, assembling grid, requires coinc!") + if _use_hpip_pp and opts.manual_initial_grid is None: + # The AMR seed-grid generators (util_AMRGrid.py / + # util_GridSubsetOfTemplateBank.py) emit XML and have not been + # converted to hyperpipeline. Per the project policy of + # operating cohesively in one mode or the other, we refuse + # rather than producing a proposed-grid.xml.gz the rest of the + # hyperpipeline workflow can't consume. + raise SystemExit( + "pseudo_pipe: --internal-use-amr seed-grid auto-generation is " + "XML-only and incompatible with RIFT_HYPERPIPELINE_FORMAT. " + "Stage your initial grid as a hyperpipeline .dat and pass it " + "via --manual-initial-grid, or run without RIFT_HYPERPIPELINE_FORMAT.") cmd += " --cip-exe `which util_AMRGrid.py ` " coinc_file = "coinc.xml" if not(os.path.exists("coinc.xml")) and not(opts.use_coinc): diff --git a/MonteCarloMarginalizeCode/Code/test/test_hyperpipeline_io.py b/MonteCarloMarginalizeCode/Code/test/test_hyperpipeline_io.py new file mode 100644 index 000000000..844de6eb9 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_hyperpipeline_io.py @@ -0,0 +1,513 @@ +""" +Tests for RIFT/misc/hyperpipeline_io.py. + +These intentionally avoid importing lalsuite or any of the heavy ILE/CIP +modules so the round-trip can be exercised in a minimal environment:: + + python test/test_hyperpipeline_io.py + pytest test/test_hyperpipeline_io.py +""" + +import os +import sys +import tempfile +import numpy as np + +# Load hyperpipeline_io directly from its file so this test does not +# trigger the RIFT.__init__ import chain (which pulls in lalsimutils + +# scipy + lalsuite -- not needed for testing pure-numpy I/O code). +HERE = os.path.dirname(os.path.abspath(__file__)) +_HPIO_PATH = os.path.normpath(os.path.join( + HERE, "..", "RIFT", "misc", "hyperpipeline_io.py")) +import importlib.util as _ilu +_spec = _ilu.spec_from_file_location("hyperpipeline_io", _HPIO_PATH) +hpio = _ilu.module_from_spec(_spec) +_spec.loader.exec_module(hpio) + + +def _roundtrip(columns, values): + """Write *values* under *columns*, read back, return structured array.""" + with tempfile.NamedTemporaryFile("w", suffix=".dat", delete=False) as fp: + fname = fp.name + try: + hpio.write_table(fname, columns, np.asarray(values, dtype=float)) + arr, hdr = hpio.read_table(fname) + assert hdr == tuple(columns), (hdr, tuple(columns)) + return arr + finally: + os.unlink(fname) + + +def test_default_roundtrip(): + cols = hpio.build_column_list() # default: 10-col bbh + rows = [ + [-12.34, 0.05, 30.0, 25.0, 0.1, 0.2, 0.3, -0.1, -0.2, -0.3], + [-13.50, 0.07, 28.0, 27.0, 0.0, 0.0, 0.4, 0.0, 0.0, -0.4], + ] + arr = _roundtrip(cols, rows) + assert len(arr) == 2 + np.testing.assert_allclose(arr["lnL"], [-12.34, -13.50]) + np.testing.assert_allclose(arr["sigma_lnL"], [0.05, 0.07]) + np.testing.assert_allclose(arr["m1"], [30.0, 28.0]) + np.testing.assert_allclose(arr["a2z"], [-0.3, -0.4]) + print("test_default_roundtrip: OK") + + +def test_eccentricity_columns(): + cols = hpio.build_column_list(use_eccentricity=True, use_meanPerAno=True) + assert cols[-2:] == ("eccentricity", "meanPerAno") + rows = [[-1.0, 0.01, 1.4, 1.4, 0,0,0, 0,0,0, 0.3, 1.7]] + arr = _roundtrip(cols, rows) + np.testing.assert_allclose(arr["eccentricity"], [0.3]) + np.testing.assert_allclose(arr["meanPerAno"], [1.7]) + print("test_eccentricity_columns: OK") + + +def test_tides_with_eos_index(): + cols = hpio.build_column_list(use_tides=True, use_eos_index=True) + assert "lambda1" in cols and "eos_table_index" in cols + rows = [[-0.5, 0.02, 1.4, 1.35, 0,0,0, 0,0,0, 500.0, 480.0, 7.0]] + arr = _roundtrip(cols, rows) + np.testing.assert_allclose(arr["lambda1"], [500.0]) + np.testing.assert_allclose(arr["eos_table_index"], [7.0]) + print("test_tides_with_eos_index: OK") + + +def test_sky_columns(): + cols = hpio.build_column_list(use_sky=True) + assert cols[-2:] == ("ecliptic_longitude", "ecliptic_latitude") + rows = [[-0.5, 0.02, 1.4, 1.35, 0,0,0, 0,0,0, 1.25, -0.4]] + arr = _roundtrip(cols, rows) + np.testing.assert_allclose(arr["ecliptic_longitude"], [1.25]) + np.testing.assert_allclose(arr["ecliptic_latitude"], [-0.4]) + + legacy = hpio.to_legacy_dat(arr, use_sky=True) + ix = hpio.legacy_column_indices(use_sky=True) + np.testing.assert_allclose(legacy[:, ix["ecliptic_longitude"]], [1.25]) + np.testing.assert_allclose(legacy[:, ix["ecliptic_latitude"]], [-0.4]) + assert ix["lnL"] == 11 + print("test_sky_columns: OK") + + +def test_to_legacy_dat_default(): + """Hyperpipeline -> legacy positional matrix, default 10-col case.""" + cols = hpio.build_column_list() + rows = [[-12.34, 0.05, 30.0, 25.0, 0.1, 0.2, 0.3, -0.1, -0.2, -0.3]] + arr = _roundtrip(cols, rows) + legacy = hpio.to_legacy_dat(arr) + # Layout: event_id, m1, m2, a1xyz, a2xyz, lnL, sigma_lnL = 11 columns + assert legacy.shape == (1, 11), legacy.shape + # event_id is synthesized to -1 + assert legacy[0, 0] == -1.0 + # m1, m2 + np.testing.assert_allclose(legacy[0, 1:3], [30.0, 25.0]) + # spins + np.testing.assert_allclose(legacy[0, 3:9], [0.1, 0.2, 0.3, -0.1, -0.2, -0.3]) + # lnL, sigma_lnL at the end + np.testing.assert_allclose(legacy[0, 9:11], [-12.34, 0.05]) + print("test_to_legacy_dat_default: OK") + + +def test_legacy_column_indices_consistency(): + """The col_lnL etc. indices must match the actual layout produced by + to_legacy_dat across all combinations of optional groups.""" + combos = [ + dict(), + dict(use_distance=True), + dict(use_tides=True), + dict(use_tides=True, use_eos_index=True), + dict(use_eccentricity=True), + dict(use_eccentricity=True, use_meanPerAno=True), + dict(use_distance=True, use_tides=True), + ] + for kw in combos: + cols = hpio.build_column_list(**kw) + # synthesize a single row of zeros + rows = [[0.0] * len(cols)] + arr = _roundtrip(cols, rows) + legacy = hpio.to_legacy_dat(arr, **kw) + ix = hpio.legacy_column_indices(**kw) + assert ix["lnL"] < legacy.shape[1], (kw, ix, legacy.shape) + assert ix["sigma_lnL"] < legacy.shape[1], (kw, ix, legacy.shape) + if kw.get("use_distance"): + assert ix["distance"] is not None + assert ix["distance"] < legacy.shape[1] + if kw.get("use_tides"): + assert ix["lambda1"] is not None + if kw.get("use_eccentricity"): + assert ix["eccentricity"] is not None + print("test_legacy_column_indices_consistency: OK") + + +def test_sniff_distinguishes_legacy(): + """sniff() must return False on a legacy-style file (no header).""" + with tempfile.NamedTemporaryFile("w", suffix=".dat", delete=False) as fp: + # Plain numeric line -- looks like a legacy ILE shard. + fp.write("0 30.0 25.0 0 0 0.5 0 0 -0.3 -12.34 0.05 1000 200\n") + fname = fp.name + try: + assert not hpio.sniff(fname), "legacy file misidentified as hyperpipeline" + finally: + os.unlink(fname) + print("test_sniff_distinguishes_legacy: OK") + + +def test_sniff_recognizes_new_format(): + cols = hpio.build_column_list() + rows = [[0.0] * len(cols)] + with tempfile.NamedTemporaryFile("w", suffix=".dat", delete=False) as fp: + fname = fp.name + try: + hpio.write_table(fname, cols, rows) + assert hpio.sniff(fname), "hyperpipeline file not recognized" + finally: + os.unlink(fname) + print("test_sniff_recognizes_new_format: OK") + + +def test_env_flag(): + saved = os.environ.pop(hpio.ENV_FLAG, None) + try: + assert not hpio.is_active() + for v in ("1", "true", "YES", "On"): + os.environ[hpio.ENV_FLAG] = v + assert hpio.is_active(), v + for v in ("0", "false", "no", "", "wat"): + os.environ[hpio.ENV_FLAG] = v + assert not hpio.is_active(), v + finally: + os.environ.pop(hpio.ENV_FLAG, None) + if saved is not None: + os.environ[hpio.ENV_FLAG] = saved + print("test_env_flag: OK") + + +def test_concatenated_shards(): + """Concatenate two single-row hyperpipeline files (as ILE produces) and + confirm the reader still parses them. Each shard carries its own + magic+header lines, both of which must be treated as comments.""" + cols = hpio.build_column_list() + with tempfile.NamedTemporaryFile("w", suffix=".dat", delete=False) as fp: + a = fp.name + with tempfile.NamedTemporaryFile("w", suffix=".dat", delete=False) as fp: + b = fp.name + with tempfile.NamedTemporaryFile("w", suffix=".composite", delete=False) as fp: + c = fp.name + try: + hpio.write_row(a, cols, [-1.0, 0.1, 30, 25, 0,0,0, 0,0,0]) + hpio.write_row(b, cols, [-2.0, 0.2, 31, 26, 0,0,0, 0,0,0]) + with open(c, "w") as out: + for src in (a, b): + with open(src) as inp: + out.write(inp.read()) + arr, hdr = hpio.read_table(c) + assert len(arr) == 2, len(arr) + np.testing.assert_allclose(sorted(arr["lnL"]), [-2.0, -1.0]) + finally: + for f in (a, b, c): + os.unlink(f) + print("test_concatenated_shards: OK") + + +def test_read_many_skips_empties_and_mismatches(): + cols_a = hpio.build_column_list() + cols_b = hpio.build_column_list(use_eccentricity=True) + good1 = tempfile.NamedTemporaryFile("w", suffix=".dat", delete=False).name + good2 = tempfile.NamedTemporaryFile("w", suffix=".dat", delete=False).name + bad_cols = tempfile.NamedTemporaryFile("w", suffix=".dat", delete=False).name + empty = tempfile.NamedTemporaryFile("w", suffix=".dat", delete=False).name + missing = "/tmp/this_file_does_not_exist_12345.dat" + try: + hpio.write_row(good1, cols_a, [-1, 0.1, 30, 25, 0,0,0, 0,0,0]) + hpio.write_row(good2, cols_a, [-2, 0.2, 31, 26, 0,0,0, 0,0,0]) + hpio.write_row(bad_cols, cols_b, [-3, 0.3, 30, 25, 0,0,0, 0,0,0, 0.1]) + # `empty` already exists with zero bytes + arr, hdr = hpio.read_many([good1, empty, good2, bad_cols, missing]) + assert hdr == cols_a + assert len(arr) == 2 # bad/empty/missing all skipped + np.testing.assert_allclose(sorted(arr["lnL"]), [-2.0, -1.0]) + finally: + for f in (good1, good2, bad_cols, empty): + if os.path.exists(f): + os.unlink(f) + print("test_read_many_skips_empties_and_mismatches: OK") + + +def test_consolidate_weighted_average(): + """Two rows with identical intrinsic params should weighted-average.""" + cols = hpio.build_column_list() + f1 = tempfile.NamedTemporaryFile("w", suffix=".dat", delete=False).name + f2 = tempfile.NamedTemporaryFile("w", suffix=".dat", delete=False).name + f3 = tempfile.NamedTemporaryFile("w", suffix=".dat", delete=False).name + try: + # Two rows at the SAME intrinsic point with very different sigma -- + # the consolidated lnL should be pulled toward the lower-sigma one. + hpio.write_row(f1, cols, [-10.0, 0.01, 30, 25, 0,0,0.5, 0,0,-0.3]) + hpio.write_row(f2, cols, [-12.0, 0.50, 30, 25, 0,0,0.5, 0,0,-0.3]) + # And one row at a DIFFERENT intrinsic point; should pass through. + hpio.write_row(f3, cols, [-5.0, 0.05, 28, 27, 0,0,0.4, 0,0,-0.4]) + arr, hdr = hpio.read_many([f1, f2, f3]) + assert len(arr) == 3 + out, out_cols = hpio.consolidate(arr, hdr) + # Two unique intrinsic points -> 2 consolidated rows. + assert len(out) == 2, len(out) + # Sorted by lnL desc; check the order. + assert out["lnL"][0] >= out["lnL"][1] + # The weighted-average for the duplicated point: with sigma=0.01 vs + # 0.50, the weight is overwhelmingly on the first row (~2500x). + # So the consolidated lnL should be very close to -10.0, not the + # arithmetic mean of -11.0. + cons = out[(out["m1"] == 30.0) & (out["m2"] == 25.0)][0] + assert abs(cons["lnL"] - (-10.0)) < 0.05, cons["lnL"] + # The pass-through row should retain its original values. + passthru = out[(out["m1"] == 28.0) & (out["m2"] == 27.0)][0] + np.testing.assert_allclose(passthru["lnL"], -5.0) + np.testing.assert_allclose(passthru["sigma_lnL"], 0.05) + finally: + for f in (f1, f2, f3): + os.unlink(f) + print("test_consolidate_weighted_average: OK") + + +class _FakeP(object): + """Stand-in for lalsimutils.ChooseWaveformParams. + + Carries m1, m2 in kg and dist in metres (matching the real class's + SI-internal convention). Spin / lambda / eccentricity components are + plain dimensionless floats. Only the attribute / setattr surface used + by hyperpipeline_io is implemented. + """ + def __init__(self): + self.m1 = 0.0; self.m2 = 0.0 + self.s1x = 0.0; self.s1y = 0.0; self.s1z = 0.0 + self.s2x = 0.0; self.s2y = 0.0; self.s2z = 0.0 + self.lambda1 = 0.0; self.lambda2 = 0.0 + self.eccentricity = 0.0; self.meanPerAno = 0.0 + self.eos_table_index = 0.0 + self.dist = 0.0 + self.phi = 0.0 + self.theta = 0.0 + self.fref = 20.0 + + +class _FakeLal(object): + """Two SI scale constants matched to lal.MSUN_SI / lal.PC_SI.""" + MSUN_SI = 1.98892e30 + PC_SI = 3.085677581e16 + + +def test_grid_write_read_roundtrip_with_units(): + """Write a P_list as a hyperpipeline grid, read it back, and verify + that mass values round-trip via the on-disk solar-mass <-> kg + conversion declared in PARAM_DISK_TO_SI.""" + fake_lal = _FakeLal() + P_list = [] + for m1_solar, m2_solar, s1z, ecc in [ + (30.0, 25.0, 0.4, 0.1), + (15.0, 12.0, -0.3, 0.05), + (1.4, 1.35, 0.0, 0.0), + ]: + P = _FakeP() + P.m1 = m1_solar * fake_lal.MSUN_SI # kg in P + P.m2 = m2_solar * fake_lal.MSUN_SI + P.s1z = s1z + P.eccentricity = ecc + P_list.append(P) + + cols = hpio.build_column_list(use_eccentricity=True) + fname = tempfile.NamedTemporaryFile("w", suffix=".dat", delete=False).name + try: + hpio.write_grid_from_P_list(fname, P_list, cols, lal_module=fake_lal, + lalsimutils_module=None) + # On-disk file should have m1 in solar masses. + arr, hdr = hpio.read_table(fname) + assert hdr == cols + np.testing.assert_allclose(arr["m1"], [30.0, 15.0, 1.4]) + np.testing.assert_allclose(arr["m2"], [25.0, 12.0, 1.35]) + np.testing.assert_allclose(arr["eccentricity"], [0.1, 0.05, 0.0]) + # lnL/sigma_lnL filled with 0 (no values supplied). + np.testing.assert_allclose(arr["lnL"], [0.0, 0.0, 0.0]) + + # Now read back into fresh P objects and verify mass is in kg. + P_back, hdr_back = hpio.read_grid_to_P_list( + fname, P_factory=_FakeP, lal_module=fake_lal, + valid_params={"m1","m2","s1x","s1y","s1z","s2x","s2y","s2z", + "lambda1","lambda2","eccentricity","meanPerAno", + "eos_table_index","dist"}) + assert len(P_back) == 3 + for orig, back in zip(P_list, P_back): + assert abs(orig.m1 - back.m1) < 1e10, (orig.m1, back.m1) # kg + assert abs(orig.m2 - back.m2) < 1e10, (orig.m2, back.m2) + assert abs(orig.s1z - back.s1z) < 1e-9 + assert abs(orig.eccentricity - back.eccentricity) < 1e-9 + finally: + os.unlink(fname) + print("test_grid_write_read_roundtrip_with_units: OK") + + +def test_grid_distance_unit_conversion(): + """Distance round-trips: stored in Mpc on disk, in metres in P.""" + fake_lal = _FakeLal() + P = _FakeP() + P.m1 = 30.0 * fake_lal.MSUN_SI + P.m2 = 25.0 * fake_lal.MSUN_SI + P.dist = 500.0 * 1e6 * fake_lal.PC_SI # 500 Mpc in metres + cols = ("lnL", "sigma_lnL", "m1", "m2", + "a1x", "a1y", "a1z", "a2x", "a2y", "a2z", "dist") + fname = tempfile.NamedTemporaryFile("w", suffix=".dat", delete=False).name + try: + hpio.write_grid_from_P_list(fname, [P], cols, lal_module=fake_lal) + arr, _ = hpio.read_table(fname) + np.testing.assert_allclose(arr["dist"], [500.0]) # Mpc on disk + # Read back into a P + P_back_list, _ = hpio.read_grid_to_P_list( + fname, P_factory=_FakeP, lal_module=fake_lal, + valid_params={"m1","m2","dist","a1x","a1y","a1z","a2x","a2y","a2z"}) + # dist should be back in metres; tolerance is generous because + # 1 Mpc = ~3e22 m so absolute differences are huge. + rel = abs(P.dist - P_back_list[0].dist) / P.dist + assert rel < 1e-12, rel + finally: + os.unlink(fname) + print("test_grid_distance_unit_conversion: OK") + + +def test_grid_auto_suffix_append(): + """write_grid_from_P_list mirrors ChooseWaveformParams_array_to_xml's + auto-suffix-append behaviour: callers can pass a basename without + extension and the writer adds '.dat' for them. This is what makes + the BasicIteration call sites (which pass `overlap-grid-N` to both + the legacy XML writer and the hyperpipeline writer) work uniformly.""" + P = _FakeP() + P.m1 = 30.0 * _FakeLal.MSUN_SI + P.m2 = 25.0 * _FakeLal.MSUN_SI + cols = hpio.build_column_list() + base = tempfile.NamedTemporaryFile("w", suffix="", delete=False).name + try: + hpio.write_grid_from_P_list(base, [P], cols, lal_module=_FakeLal()) + # The actual file should exist with a .dat suffix appended. + assert os.path.exists(base + ".dat"), "auto-suffix not appended" + assert not os.path.exists(base) or os.path.getsize(base) == 0, \ + "writer wrote to the un-suffixed path" + arr, _ = hpio.read_table(base + ".dat") + np.testing.assert_allclose(arr["m1"], [30.0]) + finally: + for f in (base, base + ".dat"): + if os.path.exists(f): + os.unlink(f) + # And: when caller already includes .dat, it's not double-appended. + base2 = tempfile.NamedTemporaryFile("w", suffix=".dat", delete=False).name + try: + hpio.write_grid_from_P_list(base2, [P], cols, lal_module=_FakeLal()) + assert not os.path.exists(base2 + ".dat"), "double-appended .dat" + arr, _ = hpio.read_table(base2) + np.testing.assert_allclose(arr["m1"], [30.0]) + finally: + for f in (base2, base2 + ".dat"): + if os.path.exists(f): + os.unlink(f) + print("test_grid_auto_suffix_append: OK") + + +def test_column_alias_bridge(): + """The on-disk column name 'a1x' bridges to ChooseWaveformParams.s1x.""" + P = _FakeP() + P.s1x = 0.42; P.s1y = -0.17; P.s2z = 0.6 + cols = hpio.build_column_list() + fname = tempfile.NamedTemporaryFile("w", suffix=".dat", delete=False).name + try: + hpio.write_grid_from_P_list(fname, [P], cols, lal_module=_FakeLal()) + arr, _ = hpio.read_table(fname) + # On disk: a1x carries the s1x value + np.testing.assert_allclose(arr["a1x"], [0.42]) + np.testing.assert_allclose(arr["a1y"], [-0.17]) + np.testing.assert_allclose(arr["a2z"], [0.6]) + # Round-trip back into a P + P_back, _ = hpio.read_grid_to_P_list( + fname, P_factory=_FakeP, lal_module=_FakeLal(), + valid_params={"s1x","s1y","s1z","s2x","s2y","s2z","m1","m2"}) + assert abs(P_back[0].s1x - 0.42) < 1e-12 + assert abs(P_back[0].s2z - 0.6) < 1e-12 + finally: + os.unlink(fname) + print("test_column_alias_bridge: OK") + + +def test_sky_column_alias_bridge(): + """Sky columns use explicit ecliptic names on disk and theta/phi on P.""" + P = _FakeP() + P.phi = 1.3 + P.theta = -0.2 + cols = hpio.build_column_list(use_sky=True) + fname = tempfile.NamedTemporaryFile("w", suffix=".dat", delete=False).name + try: + hpio.write_grid_from_P_list(fname, [P], cols, lal_module=_FakeLal()) + arr, hdr = hpio.read_table(fname) + assert hdr == cols + np.testing.assert_allclose(arr["ecliptic_longitude"], [1.3]) + np.testing.assert_allclose(arr["ecliptic_latitude"], [-0.2]) + P_back, _ = hpio.read_grid_to_P_list( + fname, P_factory=_FakeP, lal_module=_FakeLal(), + valid_params={"m1", "m2", "s1x", "s1y", "s1z", "s2x", "s2y", + "s2z", "phi", "theta"}) + assert abs(P_back[0].phi - 1.3) < 1e-12 + assert abs(P_back[0].theta + 0.2) < 1e-12 + finally: + os.unlink(fname) + print("test_sky_column_alias_bridge: OK") + + +def test_grid_no_lal_module_passthrough(): + """When lal_module is None, no unit conversion happens (raw passthrough).""" + P = _FakeP() + P.m1 = 30.0 # NOT scaled + P.m2 = 25.0 + cols = hpio.build_column_list() + fname = tempfile.NamedTemporaryFile("w", suffix=".dat", delete=False).name + try: + hpio.write_grid_from_P_list(fname, [P], cols, lal_module=None) + arr, _ = hpio.read_table(fname) + np.testing.assert_allclose(arr["m1"], [30.0]) # passthrough + finally: + os.unlink(fname) + print("test_grid_no_lal_module_passthrough: OK") + + +def test_consolidate_drops_high_sigma(): + cols = hpio.build_column_list() + f1 = tempfile.NamedTemporaryFile("w", suffix=".dat", delete=False).name + f2 = tempfile.NamedTemporaryFile("w", suffix=".dat", delete=False).name + try: + hpio.write_row(f1, cols, [-1.0, 0.01, 30, 25, 0,0,0, 0,0,0]) + hpio.write_row(f2, cols, [-2.0, 1.5, 30, 25, 0,0,0, 0,0,0]) # sigma>0.9 + arr, hdr = hpio.read_many([f1, f2]) + out, _ = hpio.consolidate(arr, hdr, sigma_cut=0.9) + assert len(out) == 1 + np.testing.assert_allclose(out["lnL"][0], -1.0) + finally: + for f in (f1, f2): + os.unlink(f) + print("test_consolidate_drops_high_sigma: OK") + + +if __name__ == "__main__": + test_default_roundtrip() + test_eccentricity_columns() + test_tides_with_eos_index() + test_sky_columns() + test_to_legacy_dat_default() + test_legacy_column_indices_consistency() + test_sniff_distinguishes_legacy() + test_sniff_recognizes_new_format() + test_env_flag() + test_concatenated_shards() + test_read_many_skips_empties_and_mismatches() + test_consolidate_weighted_average() + test_consolidate_drops_high_sigma() + test_grid_write_read_roundtrip_with_units() + test_grid_distance_unit_conversion() + test_grid_auto_suffix_append() + test_column_alias_bridge() + test_sky_column_alias_bridge() + test_grid_no_lal_module_passthrough() + print("\nAll hyperpipeline_io tests passed.") diff --git a/MonteCarloMarginalizeCode/Code/test/test_lisa_helper_contract.py b/MonteCarloMarginalizeCode/Code/test/test_lisa_helper_contract.py new file mode 100644 index 000000000..9e99ea097 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lisa_helper_contract.py @@ -0,0 +1,132 @@ +#!/usr/bin/env python +"""Contract tests for the standalone LISA CEPP helper.""" + +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__), "..", "..", "..")) +HELPER = os.path.join( + REPO_ROOT, "MonteCarloMarginalizeCode", "Code", "bin", "helper_LISA_Events.py" +) +CEPP = os.path.join( + REPO_ROOT, + "MonteCarloMarginalizeCode", + "Code", + "bin", + "create_event_parameter_pipeline_BasicIteration", +) + + +def _run_helper(tmp_path, *extra_args): + cmd = [ + sys.executable, + HELPER, + "--working-directory", + os.fspath(tmp_path), + "--zero-likelihood", + ] + cmd.extend(extra_args) + subprocess.run(cmd, check=True) + + +def test_lisa_helper_writes_cepp_contract_files(tmp_path): + _run_helper(tmp_path) + + expected = { + "proposed-grid.dat", + "args_ile.txt", + "args_cip_list.txt", + "args_test.txt", + "helper_transfer_files.txt", + "command-cepp-lisa.sh", + } + 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 "ecliptic_longitude" in columns + assert "ecliptic_latitude" in columns + + ile_args = (tmp_path / "args_ile.txt").read_text() + assert ile_args.startswith("X ") + assert "--LISA" in ile_args + assert "--h5-frame-FD" in ile_args + assert "--zero-likelihood" 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 + assert "--sim-grid" not in ile_args + assert "--output-file" not in ile_args + + 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 "--fname" not in cip_args + + test_args = (tmp_path / "args_test.txt").read_text() + assert test_args.startswith("X ") + assert "--always-succeed" in test_args + + cepp_command = (tmp_path / "command-cepp-lisa.sh").read_text() + assert "RIFT_HYPERPIPELINE_FORMAT=1" in cepp_command + assert "integrate_likelihood_extrinsic_batchmode_lisa" in cepp_command + assert os.fspath(tmp_path / "proposed-grid.dat") in cepp_command + + +def test_lisa_helper_bundle_renders_basic_cepp_dag(tmp_path): + _run_helper( + tmp_path, + "--grid-size", + "1", + "--n-iterations", + "1", + "--n-samples-per-job", + "1", + ) + + env = os.environ.copy() + env["RIFT_HYPERPIPELINE_FORMAT"] = "1" + env["PATH"] = os.path.dirname(HELPER) + os.pathsep + env.get("PATH", "") + cmd = [ + sys.executable, + CEPP, + "--ile-n-events-to-analyze", + "1", + "--input-grid", + os.fspath(tmp_path / "proposed-grid.dat"), + "--ile-exe", + HELPER.replace("helper_LISA_Events.py", "integrate_likelihood_extrinsic_batchmode_lisa"), + "--ile-args", + os.fspath(tmp_path / "args_ile.txt"), + "--cip-args-list", + os.fspath(tmp_path / "args_cip_list.txt"), + "--test-args", + os.fspath(tmp_path / "args_test.txt"), + "--working-directory", + os.fspath(tmp_path), + "--n-iterations", + "1", + "--n-samples-per-job", + "1", + "--n-copies", + "1", + "--request-memory-ILE", + "1024", + "--request-memory-CIP", + "1024", + ] + try: + subprocess.run(cmd, check=True, env=env, cwd=os.fspath(tmp_path)) + except (subprocess.CalledProcessError, FileNotFoundError) as exc: + pytest.skip("CEPP DAG render unavailable in this environment: {}".format(exc)) + + assert (tmp_path / "ILE.sub").exists() + assert (tmp_path / "CIP.sub").exists() + assert any(path.suffix == ".dag" for path in tmp_path.iterdir()) diff --git a/MonteCarloMarginalizeCode/Code/test/test_lisa_lalsimutils_compat.py b/MonteCarloMarginalizeCode/Code/test/test_lisa_lalsimutils_compat.py new file mode 100644 index 000000000..179caed2c --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lisa_lalsimutils_compat.py @@ -0,0 +1,64 @@ +#!/usr/bin/env python +"""Smoke tests for LISA lalsimutils compatibility helpers.""" + +import os +import tempfile + +import h5py +import numpy as np + +import RIFT.LISA.lalsimutils_compat as compat +import RIFT.lalsimutils as lalsimutils + + +def _write_cache(directory, h5_path): + cache_path = os.path.join(directory, "frames.cache") + with open(cache_path, "w") as cache_file: + cache_file.write(f"A A 0 1 file://localhost{h5_path}\n") + return cache_path + + +def test_frame_h5_to_hoff_reads_frequency_series(): + with tempfile.TemporaryDirectory() as tmpdir: + h5_path = os.path.join(tmpdir, "fd.h5") + with h5py.File(h5_path, "w") as h5_file: + h5_file.attrs["epoch"] = 0.0 + h5_file.attrs["f0"] = 0.0 + h5_file.attrs["deltaF"] = 0.25 + h5_file.attrs["length"] = 4 + h5_file.create_dataset("data", data=np.array([1, 2, 3, 4], dtype=np.complex128)) + + hoff = compat.frame_h5_to_hoff(_write_cache(tmpdir, h5_path), "A", verbose=False) + + assert hoff.data.length == 4 + assert hoff.deltaF == 0.25 + np.testing.assert_array_equal(hoff.data.data, np.array([1, 2, 3, 4], dtype=np.complex128)) + + +def test_h5_frame_to_non_herm_hoff_reads_time_series(): + with tempfile.TemporaryDirectory() as tmpdir: + h5_path = os.path.join(tmpdir, "td.h5") + with h5py.File(h5_path, "w") as h5_file: + h5_file.attrs["epoch"] = 0.0 + h5_file.attrs["f0"] = 0.0 + h5_file.attrs["deltaT"] = 0.5 + h5_file.attrs["length"] = 4 + h5_file.create_dataset("data", data=np.array([1, 0, 0, 0], dtype=np.float64)) + + hoff = compat.frame_data_to_non_herm_hoff( + _write_cache(tmpdir, h5_path), + "A:TD", + TDlen=-1, + h5_frame=True, + verbose=False, + ) + + assert hoff.data.length == 4 + assert hoff.deltaF == 0.5 + np.testing.assert_allclose(hoff.data.data, 0.5 * np.ones(4, dtype=np.complex128)) + + +def test_lisa_print_params_installer_is_local_to_compat_call(): + compat.install_choose_waveform_print_params_lisa() + + assert lalsimutils.ChooseWaveformParams.print_params_lisa is compat.print_params_lisa diff --git a/MonteCarloMarginalizeCode/Code/test/test_lisa_operational_synthetic.py b/MonteCarloMarginalizeCode/Code/test/test_lisa_operational_synthetic.py new file mode 100644 index 000000000..c8e7da8f8 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lisa_operational_synthetic.py @@ -0,0 +1,138 @@ +#!/usr/bin/env python +"""Operational smoke test for the LISA likelihood path.""" + +import os + +import lal +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.likelihood import factored_likelihood_LISA + + +def _synthetic_lisa_params(): + P = lalsimutils.ChooseWaveformParams() + P.m1 = 1.0e5 * lal.MSUN_SI + P.m2 = 8.0e4 * lal.MSUN_SI + P.s1z = 0.1 + P.s2z = -0.05 + P.dist = 1.0e9 * lal.PC_SI + P.fmin = 1.0e-3 + P.fref = 5.0e-3 + P.fmax = 0.125 + P.deltaT = 4.0 + P.deltaF = 1.0 / 4096.0 + P.approx = lalsim.IMRPhenomD + P.theta = 0.3 + P.phi = 1.0 + P.psi = 0.2 + P.incl = 0.4 + P.phiref = 0.1 + return P + + +def _write_cache(tmp_path): + cache_path = tmp_path / "lisa.cache" + rows = [ + ("A", "A-fake_strain-1000000-10000.h5"), + ("E", "E-fake_strain-1000000-10000.h5"), + ("T", "T-fake_strain-1000000-10000.h5"), + ] + cache_path.write_text( + "".join( + f"{channel} {channel} 0 1 file://localhost{os.fspath(tmp_path / filename)}\n" + for channel, filename in rows + ) + ) + return cache_path + + +def _flat_psd_like(data_dict): + psd_dict = {} + for channel, data in data_dict.items(): + psd = lal.CreateREAL8FrequencySeries( + "PSD", + lal.LIGOTimeGPS(0), + 0, + data.deltaF, + lalsimutils.lsu_HertzUnit, + data.data.length, + ) + psd.data.data[:] = 1.0e-40 + psd_dict[channel] = psd + return psd_dict + + +def _evaluate_lisa_lnL(rholms, cross_terms, modes, P, psi=None, inclination=None): + return factored_likelihood_LISA.FactoredLogLikelihoodAlignedSpinLISA( + rholms, + cross_terms, + P.theta, + P.phi, + np.array([P.psi if psi is None else psi]), + np.array([P.incl if inclination is None else inclination]), + np.array([P.phiref]), + np.array([P.dist]), + modes, + P.dist, + )[0] + + +def test_synthetic_lisa_tdi_precompute_and_likelihood(tmp_path): + P = _synthetic_lisa_params() + modes = [(2, 2)] + + hlms = lisa_lalsimutils_compat.hlmoff_for_LISA(P, Lmax=2, modes=modes) + generated_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(generated_data, os.fspath(tmp_path)) + + cache_path = _write_cache(tmp_path) + data_dict = { + channel: lisa_lalsimutils_compat.frame_h5_to_hoff( + os.fspath(cache_path), channel, verbose=False + ) + for channel in ("A", "E", "T") + } + psd_dict = _flat_psd_like(data_dict) + + _, cross_terms, _, rholms, _, _ = factored_likelihood_LISA.PrecomputeAlignedSpinLISA( + 0.0, + P.fref, + 8.0, + hlms, + None, + data_dict, + psd_dict, + P.fmin, + 0.5 / P.deltaT, + P.fmax, + P.deltaT, + P.theta, + P.phi, + analyticPSD_Q=False, + inv_spec_trunc_Q=False, + T_spec=0.0, + ) + + modes_array = np.array(list(hlms.keys())) + lnL_at_injection = _evaluate_lisa_lnL(rholms, cross_terms, modes_array, P) + lnL_offset = _evaluate_lisa_lnL( + rholms, cross_terms, modes_array, P, psi=P.psi + 0.8, inclination=P.incl + 0.5 + ) + + assert np.isfinite(lnL_at_injection) + assert np.isfinite(lnL_offset) + assert lnL_at_injection > lnL_offset diff --git a/MonteCarloMarginalizeCode/Code/test/test_lisa_response_import.py b/MonteCarloMarginalizeCode/Code/test/test_lisa_response_import.py new file mode 100644 index 000000000..8887b8cc2 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_lisa_response_import.py @@ -0,0 +1,23 @@ +"""Import smoke tests for dormant LISA support modules. + +These tests deliberately avoid waveform generation and data-file access. The +goal is to keep the parachuted LISA response code import-safe for CI while the +larger likelihood and pipeline integration lands in later stages. +""" + +import numpy as np + + +def test_lisa_response_import_and_pure_helper(): + import RIFT.LISA.response.LISA_response as lisa_response + + idx = lisa_response.get_closest_index(np.array([0.0, 2.0, 5.0]), 1.7) + assert idx == 1 + + +def test_lisa_factored_likelihood_import(): + import RIFT.likelihood.factored_likelihood_LISA as lisa_likelihood + + assert hasattr(lisa_likelihood, "PrecomputeAlignedSpinLISA") + assert hasattr(lisa_likelihood, "FactoredLogLikelihoodAlignedSpinLISA") + diff --git a/requirements.txt b/requirements.txt index 31d799627..b0da2ba37 100644 --- a/requirements.txt +++ b/requirements.txt @@ -22,4 +22,4 @@ ligo.skymap gwpy natsort scipy -htcondor[monitor] +#htcondor[monitor]