From f23eaada321d983d05ec1c593fb8571bfc8990bc Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Fri, 3 Jul 2026 08:09:50 -0400 Subject: [PATCH 1/2] Add cubic time interpolation for NoLoop likelihood --- .../Code/RIFT/likelihood/Q_inner_product.py | 56 ++++++++++++ .../RIFT/likelihood/cuda_Q_inner_product.cu | 64 +++++++++++++ .../RIFT/likelihood/factored_likelihood.py | 91 ++++++++++++++++--- .../Code/test/test_noloop_time_interp.py | 45 +++++++++ 4 files changed, 244 insertions(+), 12 deletions(-) create mode 100644 MonteCarloMarginalizeCode/Code/test/test_noloop_time_interp.py diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/Q_inner_product.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/Q_inner_product.py index 0cb5f5a4b..42b8cd02e 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/Q_inner_product.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/Q_inner_product.py @@ -52,3 +52,59 @@ def Q_inner_product_cupy(Q, A, start_indices, window_size): ) return out + + +def Q_inner_product_cubic_cupy(Q, A, start_indices, fractional_offsets, window_size): + """Cubic-interpolated Q inner product for fractional detector-time offsets. + + ``start_indices`` are the integer floor indices of the first requested time + sample and ``fractional_offsets`` are the corresponding fractional parts in + [0, 1). The CUDA kernel uses a four-point cubic Lagrange stencil along the + time axis and zero extension outside the precomputed Q buffer. + """ + num_time_points, num_lms = Q.shape + num_extrinsic_samples, _ = A.shape + + assert not cupy.isfortran(Q) + assert not cupy.isfortran(A) + + out = cupy.empty( + (num_extrinsic_samples, window_size), + dtype=cupy.complex128, + order="C", + ) + + global _cuda_code + if _cuda_code is None: + path = os.path.join(os.path.dirname(__file__), 'cuda_Q_inner_product.cu') + if not (os.path.isfile(path)): + path = os.path.join(os.path.split(os.path.dirname(__file__))[0], 'cuda_Q_inner_product.cu') + with open(path, 'r') as f: + _cuda_code = f.read() + Q_prod_fn = cupy.RawKernel(_cuda_code, "Q_inner_cubic") + else: + Q_prod_fn = cupy.RawKernel(_cuda_code, "Q_inner_cubic") + + float_prec = 16 + # The cubic stencil uses more registers than the nearest-neighbor kernel, so + # use a conservative default block shape to keep launches portable across + # older GPUs. Tune upward on newer cards with RIFT_Q_CUBIC_THREADS_X/Y. + num_threads_x = int(os.environ.get("RIFT_Q_CUBIC_THREADS_X", "4")) + num_threads_y = int(os.environ.get("RIFT_Q_CUBIC_THREADS_Y", "128")) + block_size = num_threads_x, num_threads_y, 0 + grid_size = ( + (num_extrinsic_samples+num_threads_x-1)//num_threads_x, + 0, + 0, + ) + args = ( + Q, A, start_indices, fractional_offsets, window_size, + num_time_points, num_extrinsic_samples, num_lms, + out, + ) + Q_prod_fn( + grid_size, block_size, args, + shared_mem=cupy.int32(0), + ) + + return out diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/cuda_Q_inner_product.cu b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/cuda_Q_inner_product.cu index d472d996a..8cbe32db3 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/cuda_Q_inner_product.cu +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/cuda_Q_inner_product.cu @@ -47,4 +47,68 @@ extern "C" { } } // if } // Q_inner + + __global__ void Q_inner_cubic( + const double2 * Q, const double2 * A, + const int * index_start, + const double * fractional_offset, + int window_size, + int num_time_points, + int num_extrinsic_samples, + int num_lms, + double2 * out + ){ + size_t sample_idx = threadIdx.x + blockDim.x*blockIdx.x; + size_t t_idx = threadIdx.y + blockDim.y * blockIdx.y; + + if (sample_idx < num_extrinsic_samples) { + int i_first_time = index_start[sample_idx]; + double u = fractional_offset[sample_idx]; + double w_m1 = -u * (u - 1.0) * (u - 2.0) / 6.0; + double w_0 = (u + 1.0) * (u - 1.0) * (u - 2.0) / 2.0; + double w_p1 = -(u + 1.0) * u * (u - 2.0) / 2.0; + double w_p2 = (u + 1.0) * u * (u - 1.0) / 6.0; + + for (size_t i_time = t_idx; i_time < window_size; i_time+=blockDim.y) { + size_t i_output = sample_idx*window_size + i_time; + int q_time = i_first_time + (int)i_time; + double out_re = 0.0; + double out_im = 0.0; + + for (size_t i_lm = 0; i_lm < num_lms; ++i_lm) { + double q_re = 0.0; + double q_im = 0.0; + int q_m1 = q_time - 1; + int q_0 = q_time; + int q_p1 = q_time + 1; + int q_p2 = q_time + 2; + if (q_m1 >= 0 && q_m1 < num_time_points) { + double2 q = Q[((size_t)q_m1)*num_lms + i_lm]; + q_re += w_m1 * q.x; + q_im += w_m1 * q.y; + } + if (q_0 >= 0 && q_0 < num_time_points) { + double2 q = Q[((size_t)q_0)*num_lms + i_lm]; + q_re += w_0 * q.x; + q_im += w_0 * q.y; + } + if (q_p1 >= 0 && q_p1 < num_time_points) { + double2 q = Q[((size_t)q_p1)*num_lms + i_lm]; + q_re += w_p1 * q.x; + q_im += w_p1 * q.y; + } + if (q_p2 >= 0 && q_p2 < num_time_points) { + double2 q = Q[((size_t)q_p2)*num_lms + i_lm]; + q_re += w_p2 * q.x; + q_im += w_p2 * q.y; + } + double2 a = A[sample_idx*num_lms + i_lm]; + out_re += a.x*q_re - a.y*q_im; + out_im += a.x*q_im + a.y*q_re; + } + + out[i_output] = make_double2(out_re, out_im); + } + } + } // Q_inner_cubic } // extern diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood.py index 1e9b6df03..a5f87f089 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood.py @@ -1805,7 +1805,53 @@ def _factored_lnL_helper(kappa_sq, rho_sq): return kappa_sq - 0.5 * rho_sq -def DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop(tvals, P_vec, lookupNKDict, rholmsArrayDict, ctUArrayDict,ctVArrayDict,epochDict,Lmax=2,array_output=False,xpy=np, loglikelihood=_factored_lnL_helper,return_lnLt=False,phase_marginalization=False): +def _cubic_Q_window_numpy(Q_block, start_indices, fractional_offsets, npts): + """Return cubic-interpolated Q windows with zero extension. + + Q_block has shape (n_time, n_lm). The returned array has shape + (n_extrinsic, npts, n_lm), matching the CPU fallback layout used by NoLoop. + The interpolation is a local four-point cubic Lagrange stencil; at integer + offsets it reproduces the original samples exactly. + """ + npts_extrinsic = len(start_indices) + n_lms_det = Q_block.shape[1] + Qlms = np.zeros((npts_extrinsic, npts, n_lms_det), dtype=np.complex128) + tgrid = np.arange(npts) + n_time = Q_block.shape[0] + for i in range(npts_extrinsic): + idxs = int(start_indices[i]) + tgrid + u = float(fractional_offsets[i]) + u2 = u*u + u3 = u2*u + weights = ( + -u * (u - 1.0) * (u - 2.0) / 6.0, + (u + 1.0) * (u - 1.0) * (u - 2.0) / 2.0, + -(u + 1.0) * u * (u - 2.0) / 2.0, + (u + 1.0) * u * (u - 1.0) / 6.0, + ) + for offset, weight in zip((-1, 0, 1, 2), weights): + idxs_here = idxs + offset + valid = (idxs_here >= 0) & (idxs_here < n_time) + if np.any(valid): + Qlms[i, valid] += weight * Q_block[idxs_here[valid]] + return Qlms + + +def _nearest_Q_window_numpy(Q_block, start_indices, npts, xpy=np): + """Return nearest-grid Q windows with zero extension.""" + npts_extrinsic = len(start_indices) + n_lms_det = Q_block.shape[1] + Qlms = xpy.zeros((npts_extrinsic, npts, n_lms_det), dtype=np.complex128) + tgrid = np.arange(npts) + n_time = Q_block.shape[0] + for i in range(npts_extrinsic): + idxs = int(start_indices[i]) + tgrid + valid = (idxs >= 0) & (idxs < n_time) + Qlms[i, valid] = Q_block[idxs[valid]] + return Qlms + + +def DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop(tvals, P_vec, lookupNKDict, rholmsArrayDict, ctUArrayDict,ctVArrayDict,epochDict,Lmax=2,array_output=False,xpy=np, loglikelihood=_factored_lnL_helper,return_lnLt=False,phase_marginalization=False,time_interp='nearest'): """ DiscreteFactoredLogLikelihoodViaArray uses the array-ized data structures to compute the log likelihood, either as an array vs time *or* marginalized in time. @@ -1813,9 +1859,19 @@ def DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop(tvals, P_vec, lookupNKDic The timeseries quantities are computed via discrete shifts of an existing grid Note 'P' must have the *sampling rate* set to correctly interpret the event time. Note arguments passed are NOW ARRAYS, in contrast to similar function which does not have 'Vector' postfix + + time_interp : {'nearest', 'cubic'} + Detector-time sampling convention for the data term. 'nearest' + preserves the historical NoLoop integer-bin gather. 'cubic' evaluates + the precomputed Q_lm time series at the fractional detector arrival time + using a four-sample cubic Lagrange stencil, with zero extension outside + the precomputed buffer. """ global distMpcRef + if time_interp not in ('nearest', 'cubic'): + raise ValueError("time_interp must be 'nearest' or 'cubic'") + detectors = rholmsArrayDict.keys() npts = len(tvals) npts_extrinsic = len(P_vec.phi) @@ -1910,7 +1966,13 @@ def DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop(tvals, P_vec, lookupNKDic ) tfirst = t_det + tvals[0] - ifirst = (xpy.rint((tfirst) / deltaT) + 0.5).astype(np.int32) # C uses 32 bit integers : be careful + sample_first = tfirst / deltaT + if time_interp == 'nearest': + ifirst = (xpy.rint(sample_first) + 0.5).astype(np.int32) # C uses 32 bit integers : be careful + frac_first = None + else: + ifirst = xpy.floor(sample_first).astype(np.int32) + frac_first = (sample_first - xpy.floor(sample_first)).astype(np.float64) # ilast = ifirst + npts @@ -1983,15 +2045,22 @@ def DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop(tvals, P_vec, lookupNKDic # Shape Q = (npts_time_full, nlms) # Shape A=FY_conj = (npts_extrinsic, nlms) # shape result = (npts_extrinsic, npts_time_*window* = npts) - Q_prod_result = Q_inner_product.Q_inner_product_cupy( - Q, FY_conj, - ifirst, npts, - ) + if time_interp == 'nearest': + Q_prod_result = Q_inner_product.Q_inner_product_cupy( + Q, FY_conj, + ifirst, npts, + ) + else: + Q_prod_result = Q_inner_product.Q_inner_product_cubic_cupy( + Q, FY_conj, + ifirst, frac_first, npts, + ) else: - # Use old code completely unchanged ... very wasteful on memory management! - Qlms = xpy.empty((npts_extrinsic, npts, n_lms), dtype=np.complex128) - for i in range(npts_extrinsic): - Qlms[i] = rholmsArrayDict[det][...,ifirst[i]:(ifirst[i]+npts)].T + Q_block = rholmsArrayDict[det].T + if time_interp == 'nearest': + Qlms = _nearest_Q_window_numpy(Q_block, ifirst, npts, xpy=xpy) + else: + Qlms = _cubic_Q_window_numpy(Q_block, ifirst, frac_first, npts) if phase_marginalization: Qlms[:, :, 1] = xpy.conj(Qlms[:, :, 1]) @@ -2157,5 +2226,3 @@ def ComputeYlmsArrayVector(lookupNK, theta, phi): Ylms[indx] = lalylm(theta, phi, s, l, m) return Ylms - - diff --git a/MonteCarloMarginalizeCode/Code/test/test_noloop_time_interp.py b/MonteCarloMarginalizeCode/Code/test/test_noloop_time_interp.py new file mode 100644 index 000000000..334155d99 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/test_noloop_time_interp.py @@ -0,0 +1,45 @@ +import os + +os.environ.setdefault("RIFT_LOWLATENCY", "1") + +import numpy as np + +from RIFT.likelihood.factored_likelihood import ( + _cubic_Q_window_numpy, + _nearest_Q_window_numpy, +) + + +def test_cubic_q_window_matches_nearest_at_integer_samples(): + q = (np.arange(40) + 1j * np.arange(40, 80)).reshape(20, 2) + start = np.asarray([2, 7, 12], dtype=np.int32) + + nearest = _nearest_Q_window_numpy(q, start, 5) + cubic = _cubic_Q_window_numpy(q, start, np.zeros(len(start)), 5) + + assert np.allclose(cubic, nearest) + + +def test_cubic_q_window_reproduces_linear_midpoints(): + q = (np.arange(40) + 1j * np.arange(40, 80)).reshape(20, 2) + start = np.asarray([3], dtype=np.int32) + + cubic = _cubic_Q_window_numpy(q, start, np.asarray([0.5]), 4) + expected = 0.5 * (q[3:7] + q[4:8]) + + assert np.allclose(cubic[0], expected) + + +def test_cubic_q_window_reproduces_cubic_polynomial(): + times = np.arange(30, dtype=float) + values = times**3 - 2.0 * times**2 + 0.5 * times - 7.0 + q = np.column_stack([values, values + 1j * values]) + start = np.asarray([10], dtype=np.int32) + frac = np.asarray([0.37]) + + cubic = _cubic_Q_window_numpy(q, start, frac, 3) + target_times = start[0] + frac[0] + np.arange(3) + expected_values = target_times**3 - 2.0 * target_times**2 + 0.5 * target_times - 7.0 + expected = np.column_stack([expected_values, expected_values + 1j * expected_values]) + + assert np.allclose(cubic[0], expected) From 512c0012267572b78d33c681a5775170fd3948e7 Mon Sep 17 00:00:00 2001 From: Richard O'Shaughnessy Date: Fri, 3 Jul 2026 08:45:25 -0400 Subject: [PATCH 2/2] Wire cubic NoLoop time interpolation into ILE --- .../Code/bin/integrate_likelihood_extrinsic | 29 ++++++++++++--- .../integrate_likelihood_extrinsic_batchmode | 35 +++++++++++++++---- .../Code/bin/util_RIFT_pseudo_pipe.py | 3 ++ 3 files changed, 57 insertions(+), 10 deletions(-) diff --git a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic index 347e0bce3..e59523312 100755 --- a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic +++ b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic @@ -205,7 +205,7 @@ integration_params.add_option("--adapt-weight-exponent", type=float, default=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("--interpolate-time", default=False,help="If using the maintained NoLoop likelihood, evaluate Q_lm at fractional detector times using cubic interpolation instead of nearest sample bins. Accepts truthy values such as True/1/yes. (Default=false)") integration_params.add_option("--d-prior",default='Euclidean' ,type=str,help="Distance prior for dL. Options are dL^2 (Euclidean) 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.") @@ -252,7 +252,26 @@ for pin_param in LIKELIHOOD_PINNABLE_PARAMS: pinnable.add_option(option, type=float, help="Pin the value of %s." % pin_param) optp.add_option_group(pinnable) -opts, args = optp.parse_args() +def _normalize_interpolate_time_argv(argv): + out = [] + i = 0 + while i < len(argv): + out.append(argv[i]) + if argv[i] == "--interpolate-time" and (i + 1 == len(argv) or argv[i + 1].startswith("--")): + out.append("True") + i += 1 + return out + +opts, args = optp.parse_args(_normalize_interpolate_time_argv(sys.argv[1:])) + +def _truthy_option(value): + if isinstance(value, bool): + return value + if value is None: + return False + return str(value).strip().lower() in ("1", "true", "t", "yes", "y", "on") + +opts._noloop_time_interp = "cubic" if _truthy_option(opts.interpolate_time) else "nearest" 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]") @@ -1180,7 +1199,8 @@ else: # Sum over time for every point in other extrinsic params P.phi = xpy_default.mod(P.phi, 2*np.pi) lnL = factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop(tvals, - P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max,xpy=xpy_default) + P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max,xpy=xpy_default, + time_interp=opts._noloop_time_interp) nEvals +=len(right_ascension) return identity_convert(xpy_default.exp(lnL-manual_avoid_overflow_logarithm)) else: @@ -1261,7 +1281,8 @@ else: # Sum over time for every point in other extrinsic params P.phi = xpy_default.mod(P.phi, 2*np.pi) lnL = factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop(tvals, - P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max,xpy=xpy_default, loglikelihood=distmarg_loglikelihood) + P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max,xpy=xpy_default, loglikelihood=distmarg_loglikelihood, + time_interp=opts._noloop_time_interp) nEvals +=len(right_ascension) return identity_convert(xpy_default.exp(lnL-manual_avoid_overflow_logarithm)) diff --git a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode index 1ad5a917b..95cba9c99 100755 --- a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode +++ b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode @@ -271,7 +271,7 @@ integration_params.add_option("--adapt-weight-exponent", type=float, default=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("--interpolate-time", default=False,help="If using the maintained NoLoop likelihood, evaluate Q_lm at fractional detector times using cubic interpolation instead of nearest sample bins. Accepts truthy values such as True/1/yes. (Default=false)") integration_params.add_option("--d-prior",default='Euclidean' ,type=str,help="Distance prior for dL. Options are dL^2 (Euclidean), 'pseudo_cosmo', and 'cosmo' and 'cosmo_sourceframe' .") integration_params.add_option("--d-prior-redshift", action='store_true', help="If true, distance prior is computed in redshift. This option MAY be enforced for 'cosmo' sampling") 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.") @@ -329,7 +329,26 @@ for pin_param in LIKELIHOOD_PINNABLE_PARAMS: pinnable.add_option(option, type=float, help="Pin the value of %s." % pin_param) optp.add_option_group(pinnable) -opts, args = optp.parse_args() +def _normalize_interpolate_time_argv(argv): + out = [] + i = 0 + while i < len(argv): + out.append(argv[i]) + if argv[i] == "--interpolate-time" and (i + 1 == len(argv) or argv[i + 1].startswith("--")): + out.append("True") + i += 1 + return out + +opts, args = optp.parse_args(_normalize_interpolate_time_argv(sys.argv[1:])) + +def _truthy_option(value): + if isinstance(value, bool): + return value + if value is None: + return False + return str(value).strip().lower() in ("1", "true", "t", "yes", "y", "on") + +opts._noloop_time_interp = "cubic" if _truthy_option(opts.interpolate_time) else "nearest" # cosmo d prior tools for interpolation: not used normally, but set if needed final_scipy_interpolate=None @@ -1334,7 +1353,8 @@ def resample_samples(my_samples, 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) + P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max,xpy=xpy_default,return_lnLt=True, + time_interp=opts._noloop_time_interp) lnLt = identity_convert(lnLt) # back to CPU. Note we have removed offsets if opts.zero_likelihood: @@ -1620,7 +1640,8 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t lnL = factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop(tvals, - P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max,xpy=xpy_default) + P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max,xpy=xpy_default, + time_interp=opts._noloop_time_interp) # 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 @@ -1723,7 +1744,8 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t 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) + P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max,xpy=xpy_default, loglikelihood=distmarg_loglikelihood, phase_marginalization=True, + time_interp=opts._noloop_time_interp) # 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 @@ -1766,7 +1788,8 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t 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) + P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max,xpy=xpy_default, loglikelihood=distmarg_loglikelihood, + time_interp=opts._noloop_time_interp) # 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 diff --git a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py index 4c07b204b..36e6e7acd 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py @@ -211,6 +211,7 @@ def unsafe_parse_arg_string_dict(my_argstr): parser.add_argument("--internal-n-evaluations-per-iteration",default=None,type=int,help="Number of ILE evaluation points per iteration, if not set then pipeline selects experience-based default. Each ILE worker will do a fraction of this total workload.") parser.add_argument("--add-extrinsic",action='store_true') parser.add_argument("--add-extrinsic-time-resampling",action='store_true',help="adds the time resampling option. Only deployed for vectorized calculations (which should be all that end-users can access)") +parser.add_argument("--internal-ile-interpolate-time",action='store_true',help="Pass --interpolate-time True to ILE, enabling cubic interpolation of Q_lm at fractional detector arrival times in the maintained NoLoop likelihood.") parser.add_argument("--batch-extrinsic",action='store_true') parser.add_argument("--fmin",default=20,type=int,help="Mininum frequency for integration. template minimum frequency (we hope) so all modes resolved at this frequency") # should be 23 for the BNS parser.add_argument("--fmin-template",default=None,type=float,help="Mininum frequency for template. If provided, then overrides automated settings for fmin-template = fmin/Lmax") # should be 23 for the BNS @@ -934,6 +935,8 @@ def unsafe_parse_arg_string_dict(my_argstr): # - requested or # - AC + not freezeadapt line += " --force-reset-all " +if opts.internal_ile_interpolate_time: + line += " --interpolate-time True " if not(opts.manual_extra_ile_args is None): line += " {} ".format(opts.manual_extra_ile_args) # embed with space on each side, avoid collisions if '--declination ' in opts.manual_extra_ile_args: # if we are pinning dec, we aren't using a cosine coordinate. Don't mess up.