Commit 8bde891a authored by John Douglas Veitch's avatar John Douglas Veitch
Browse files

Merge branch 'o3_refactor' into 'master'

O3 refactor

See merge request !329
parents 367b097f 1a9375df
Pipeline #23005 failed with stages
in 231 minutes and 6 seconds
......@@ -146,33 +146,44 @@ def factored_log_likelihood(extr_params, rholms_intp, crossTerms, Lmax):
psi = extr_params.psi
dist = extr_params.dist
# N.B.: The Ylms are a function of - phiref b/c we are passively rotating
# the source frame, rather than actively rotating the binary.
# Said another way, the m^th harmonic of the waveform should transform as
# e^{- i m phiref}, but the Ylms go as e^{+ i m phiref}, so we must give
# - phiref as an argument so Y_lm h_lm has the proper phiref dependence
# FIXME: Strictly speaking, this should be inside the detector loop because
# there *could* be different l,m pairs for different detectors. This never
# happens in practice, so it's pulled out here, and we use the first
# detector as a reference.
Ylms = compute_spherical_harmonics(Lmax, incl, -phiref, rholms_intp[rholms_intp.keys()[0]])
lnL = 0.
for det in detectors:
CT = crossTerms[det]
# This is the GPS time at the detector
t_det = compute_arrival_time_at_detector(det, RA, DEC, tref)
F = complex_antenna_factor(det, RA, DEC, psi, t_det)
det_rholms = {} # rholms evaluated at time at detector
for key in rholms_intp[det]:
func = rholms_intp[det][key]
det_rholms[key] = func(float(t_det))
lnL += single_detector_log_likelihood(det_rholms, CT, Ylms, F, dist)
return lnL
# use EXTREMELY many bits
i = 0
_lnL = np.zeros(np.shape(extr_params.phi), dtype=np.float128)
for RA, DEC, phiref, incl, psi, tref, dist in zip(extr_params.phi, \
extr_params.theta, extr_params.phiref, \
extr_params.incl, extr_params.psi, \
extr_params.tref, extr_params.dist):
# N.B.: The Ylms are a function of - phiref b/c we are passively rotating
# the source frame, rather than actively rotating the binary.
# Said another way, the m^th harmonic of the waveform should transform as
# e^{- i m phiref}, but the Ylms go as e^{+ i m phiref}, so we must give
# - phiref as an argument so Y_lm h_lm has the proper phiref dependence
# FIXME: Strictly speaking, this should be inside the detector loop because
# there *could* be different l,m pairs for different detectors. This never
# happens in practice, so it's pulled out here, and we use the first
# detector as a reference.
Ylms = compute_spherical_harmonics(Lmax, incl, -phiref, rholms_intp[rholms_intp.keys()[0]])
lnL = 0.
for det in detectors:
CT = crossTerms[det]
# This is the GPS time at the detector
t_det = compute_arrival_time_at_detector(det, RA, DEC, tref)
F = complex_antenna_factor(det, RA, DEC, psi, t_det)
det_rholms = {} # rholms evaluated at time at detector
for key in rholms_intp[det]:
func = rholms_intp[det][key]
det_rholms[key] = func(float(t_det))
lnL += single_detector_log_likelihood(det_rholms, CT, Ylms, F, dist)
_lnL[i] = lnL
i += 1
return _lnL
def factored_log_likelihood_time_marginalized(tvals, extr_params, rholms_intp, rholms, crossTerms, det_epochs, Lmax, interpolate=False):
"""
......@@ -201,44 +212,56 @@ def factored_log_likelihood_time_marginalized(tvals, extr_params, rholms_intp, r
psi = extr_params.psi
dist = extr_params.dist
# N.B.: The Ylms are a function of - phiref b/c we are passively rotating
# the source frame, rather than actively rotating the binary.
# Said another way, the m^th harmonic of the waveform should transform as
# e^{- i m phiref}, but the Ylms go as e^{+ i m phiref}, so we must give
# - phiref as an argument so Y_lm h_lm has the proper phiref dependence
# FIXME: Strictly speaking, this should be inside the detector loop because
# there *could* be different l,m pairs for different detectors. This never
# happens in practice, so it's pulled out here, and we use the first
# detector as a reference.
Ylms = compute_spherical_harmonics(Lmax, incl, -phiref, rholms[rholms.keys()[0]])
lnL = 0.
delta_t = tvals[1] - tvals[0]
for det in detectors:
CT = crossTerms[det]
F = complex_antenna_factor(det, RA, DEC, psi, tref)
rho_epoch = float(det_epochs[det])
# This is the GPS time at the detector
t_det = compute_arrival_time_at_detector(det, RA, DEC, tref)
det_rholms = {} # rholms evaluated at time at detector
if ( interpolate ):
# use the interpolating functions.
for key, func in rholms_intp[det].iteritems():
det_rholms[key] = func(float(t_det)+tvals)
else:
# do not interpolate, just use nearest neighbors.
for key, rhoTS in rholms[det].iteritems():
# PRB: these can be moved outside this loop to after t_det
tfirst = float(t_det)+tvals[0]
ifirst = int((tfirst - rho_epoch) / delta_t + 0.5)
ilast = ifirst + len(tvals)
det_rholms[key] = rhoTS[ifirst:ilast]
lnL += single_detector_log_likelihood(det_rholms, CT, Ylms, F, dist)
maxlnL = np.max(lnL)
return maxlnL + np.log(np.sum(np.exp(lnL - maxlnL)) * (tvals[1]-tvals[0]))
# use EXTREMELY many bits
i = 0
_lnL = np.zeros(np.shape(extr_params.phi), dtype=np.float128)
for RA, DEC, phiref, incl, psi, dist in zip(extr_params.phi, \
extr_params.theta, extr_params.phiref, \
extr_params.incl, extr_params.psi, \
extr_params.dist):
# N.B.: The Ylms are a function of - phiref b/c we are passively rotating
# the source frame, rather than actively rotating the binary.
# Said another way, the m^th harmonic of the waveform should transform as
# e^{- i m phiref}, but the Ylms go as e^{+ i m phiref}, so we must give
# - phiref as an argument so Y_lm h_lm has the proper phiref dependence
# FIXME: Strictly speaking, this should be inside the detector loop because
# there *could* be different l,m pairs for different detectors. This never
# happens in practice, so it's pulled out here, and we use the first
# detector as a reference.
Ylms = compute_spherical_harmonics(Lmax, incl, -phiref, rholms[rholms.keys()[0]])
lnL = 0.
delta_t = tvals[1] - tvals[0]
for det in detectors:
CT = crossTerms[det]
F = complex_antenna_factor(det, RA, DEC, psi, tref)
rho_epoch = float(det_epochs[det])
# This is the GPS time at the detector
t_det = compute_arrival_time_at_detector(det, RA, DEC, tref)
det_rholms = {} # rholms evaluated at time at detector
if ( interpolate ):
# use the interpolating functions.
for key, func in rholms_intp[det].iteritems():
det_rholms[key] = func(float(t_det)+tvals)
else:
# do not interpolate, just use nearest neighbors.
for key, rhoTS in rholms[det].iteritems():
# PRB: these can be moved outside this loop to after t_det
tfirst = float(t_det)+tvals[0]
ifirst = int((tfirst - rho_epoch) / delta_t + 0.5)
ilast = ifirst + len(tvals)
det_rholms[key] = rhoTS[ifirst:ilast]
lnL += single_detector_log_likelihood(det_rholms, CT, Ylms, F, dist)
maxlnL = np.max(lnL)
_lnL[i] = maxlnL + np.log(np.sum(np.exp(lnL - maxlnL)) * (tvals[1]-tvals[0]))
i += 1
return _lnL
def single_detector_log_likelihood(rholm_vals, crossTerms, Ylms, F, dist):
"""
......
......@@ -1047,6 +1047,8 @@ def hlmoft(P, Lmax=2, Fp=None, Fc=None):
hlms = lalsim.SimInspiralChooseTDModes( \
P.phiref, P.deltaT, \
P.m1, P.m2, \
P.spin1x, P.spin1y, P.spin1z, \
P.spin2x, P.spin2y, P.spin2z, \
P.fmin, P.fref, P.dist, \
extra_params, Lmax, P.approx)
except RuntimeError:
......@@ -1054,7 +1056,7 @@ def hlmoft(P, Lmax=2, Fp=None, Fc=None):
P.m1, P.m2, \
P.spin1x, P.spin1y, P.spin1z, \
P.spin2x, P.spin2y, P.spin2z, \
P.dist, P.incl, P.phiref, \
P.dist, P.phiref, \
P.psi, P.eccentricity, P.meanPerAno, \
P.deltaT, P.fmin, P.fref, \
extra_params, P.approx)
......
......@@ -465,7 +465,7 @@ class MCSampler(object):
raise NanOrInf("maxlnL = inf")
if show_evaluation_log:
print int(time.time()), ": ", self.ntotal, eff_samp, math.log(maxval), numpy.log(int_val1/self.ntotal), numpy.log(int_val1/self.ntotal)-maxlnL, numpy.sqrt(var*self.ntotal)/int_val1
print "{0:.3f} : {1:d} {2:.5f} {3:.2f} {4:.2f} {5:.2f} {6:.3f}".format(time.time(), self.ntotal, eff_samp, math.log(maxval), numpy.log(int_val1 / self.ntotal), numpy.log(int_val1 / self.ntotal) - maxlnL, numpy.sqrt(var * self.ntotal) / int_val1)
if (not convergence_tests) and self.ntotal >= nmin and self.ntotal >= nmax and neff != float("inf"):
print >>sys.stderr, "WARNING: User requested maximum number of samples reached... bailing."
......
......@@ -487,21 +487,16 @@ if not opts.time_marginalization:
# to it.
#
def likelihood_function(right_ascension, declination, t_ref, phi_orb, inclination, psi, distance):
# 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, declination,
t_ref, phi_orb, inclination, 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 * lal.PC_SI # luminosity distance
P.phi = right_ascension # right ascension
P.theta = declination # declination
P.tref = t_ref + fiducial_epoch # ref. time (rel to epoch for data taking)
P.phiref = phi_orb # ref. orbital phase
P.incl = inclination # inclination
P.psi = psi # polarization angle
P.dist = distance * 1.e6 * lal.PC_SI # luminosity distance
lnL[i] = factored_likelihood.factored_log_likelihood(P, rholms_intp, cross_terms, opts.l_max)
i+=1
lnL = factored_likelihood.factored_log_likelihood(P, rholms_intp, cross_terms, opts.l_max)
return numpy.exp(lnL)
......@@ -530,23 +525,16 @@ else: # Sum over time for every point in other extrinsic params
def likelihood_function(right_ascension, declination, phi_orb, inclination,
psi, distance):
# use EXTREMELY many bits
lnL = numpy.zeros(right_ascension.shape,dtype=numpy.float128)
# PRB: can we move this loop inside the factored_likelihood? It might help.
i = 0
# choose an array at the target sampling rate. P is inherited globally
for ph, th, phr, ic, ps, di in zip(right_ascension, declination,
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 * lal.PC_SI # luminosity distance
lnL[i] = factored_likelihood.factored_log_likelihood_time_marginalized(tvals, P, rholms_intp, rholms, cross_terms, det_epochs, opts.l_max,interpolate=opts.interpolate_time)
i+=1
P.phi = right_ascension # right ascension
P.theta = declination # declination
P.tref = fiducial_epoch # see 'tvals', above
P.phiref = phi_orb # ref. orbital phase
P.incl = inclination # inclination
P.psi = psi # polarization angle
P.dist = distance * 1.e6 * lal.PC_SI # luminosity distance
lnL = factored_likelihood.factored_log_likelihood_time_marginalized(tvals, P, rholms_intp, rholms, cross_terms, det_epochs, opts.l_max,interpolate=opts.interpolate_time)
return numpy.exp(lnL)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment