Adding likelihood to multi_signal model for inference of LISA signals. (

gwastro#4582)

* Changes to likelihood for multiple modes

* Added multi mode likelihood for wfs that don't need detector response computed.

* Experimenting with multi signal likelihood and bbhx

* Futher likelihood testing for multiple modes

* Test

* Still testing likelihood

* Almost fixed the cross term likelihood

* Testing

* Final tests

* Multi signal likelihood for lisa bbhx cleanup

* Tidy before pr

* Added check for f_lower to set to -1 (pycbc default value) if f_lower not set in params.

* Update docs examples with delta_f value

* Undo hard coded delta_f

* Add error when f_lower not supplied

* Undo previous change from sequence. Added Exception when f_lower not provided in get_fd_det_waveform.

* Removed exception

* Correct indents

* Removed TODO from relbin. Removed delta_f as a required parameter when using the get_fd_det_waveform_sequence method

examples/inference/lisa_smbhb_inj/lisa_smbhb_relbin.ini

@@ -50,6 +50,7 @@ spin2z = 0.36905807298613247
 distance = 17758.367941273442
 inclination = 1.5970175301911231
 t_obs_start = 31536000
+f_lower = 1e-4
 ; Put LISA behind the Earth by ~20 degrees.
 t_offset = 7365189.431698299
 

examples/inference/lisa_smbhb_ldc/lisa_smbhb_relbin.ini

@@ -41,6 +41,7 @@ spin2z = 0.36905807298613247
 distance = 17758.367941273442
 inclination = 1.5970175301911231
 t_obs_start = 31536000
+f_lower = 1e-4
 ; Put LISA behind the Earth by ~20 degrees.
 t_offset = 7365189.431698299
 

pycbc/inference/models/brute_marg.py

@@ -200,7 +200,7 @@ def _apply_sky_point_rotation(self, pref, sky_num):
     def from_config(cls, cp, **kwargs):
         kwargs['config_object'] = cp
         return super(BruteLISASkyModesMarginalize, cls).from_config(
-            cp, 
+            cp,
             **kwargs
         )
 

pycbc/inference/models/relbin.py

@@ -39,7 +39,8 @@
 from .gaussian_noise import BaseGaussianNoise
 from .relbin_cpu import (likelihood_parts, likelihood_parts_v,
                          likelihood_parts_multi, likelihood_parts_multi_v,
-                         likelihood_parts_det, likelihood_parts_vector,
+                         likelihood_parts_det, likelihood_parts_det_multi,
+                         likelihood_parts_vector,
                          likelihood_parts_v_pol,
                          likelihood_parts_v_time,
                          likelihood_parts_v_pol_time,
@@ -221,7 +222,6 @@ def __init__(
             self.f[ifo] = numpy.array(d0.sample_frequencies)
             self.df[ifo] = d0.delta_f
             self.end_time[ifo] = float(d0.end_time)
-            # self.det[ifo] = Detector(ifo)
 
             # generate fiducial waveform
             f_lo = self.kmin[ifo] * self.df[ifo]
@@ -240,7 +240,7 @@ def __init__(
                                                     sample_points=fpoints,
                                                     **self.fid_params)
                 curr_wav = wave[ifo]
-                self.ta[ifo] = 0
+                self.ta[ifo] = 0.
             else:
                 fid_hp, fid_hc = get_fd_waveform_sequence(sample_points=fpoints,
                                                           **self.fid_params)
@@ -362,6 +362,7 @@ def setup_antenna(self, earth_rotation, mode, fedges):
             atimes = self.fid_params["tc"]
             if self.still_needs_det_response:
                 self.lik = likelihood_parts_det
+                self.mlik = likelihood_parts_det_multi
             else:
                 self.lik = likelihood_parts
                 self.mlik = likelihood_parts_multi
@@ -490,19 +491,33 @@ def multi_loglikelihood(self, models):
         if not hasattr(self, 'hihj'):
             self.calculate_hihjs(models)
 
-        # finally add in the lognl term from this model
-        for m1, m2 in itertools.combinations(models, 2):
-            for det in self.data:
-                a0, a1, fedge = self.hihj[(m1, m2)][det]
-
-                fp, fc, dtc, hp, hc, h00 = m1._current_wf_parts[det]
-                fp2, fc2, dtc2, hp2, hc2, h002 = m2._current_wf_parts[det]
-
-                h1h2 = self.mlik(fedge,
-                                 fp, fc, dtc, hp, hc, h00,
-                                 fp2, fc2, dtc2, hp2, hc2, h002,
-                                 a0, a1)
-                loglr += - h1h2.real # This is -0.5 * re(<h1|h2> + <h2|h1>)
+        if self.still_needs_det_response:
+            for m1, m2 in itertools.combinations(models, 2):
+                for det in self.data:
+                    a0, a1, fedge = self.hihj[(m1, m2)][det]
+
+                    dtc, channel, h00 = m1._current_wf_parts[det]
+                    dtc2, channel2, h002 = m2._current_wf_parts[det]
+
+                    c1c2 = self.mlik(fedge,
+                                     dtc, channel, h00,
+                                     dtc2, channel2, h002,
+                                     a0, a1)
+                    loglr += - c1c2.real # This is -0.5 * re(<h1|h2> + <h2|h1>)
+        else:
+            # finally add in the lognl term from this model
+            for m1, m2 in itertools.combinations(models, 2):
+                for det in self.data:
+                    a0, a1, fedge = self.hihj[(m1, m2)][det]
+
+                    fp, fc, dtc, hp, hc, h00 = m1._current_wf_parts[det]
+                    fp2, fc2, dtc2, hp2, hc2, h002 = m2._current_wf_parts[det]
+
+                    h1h2 = self.mlik(fedge,
+                                     fp, fc, dtc, hp, hc, h00,
+                                     fp2, fc2, dtc2, hp2, hc2, h002,
+                                     a0, a1)
+                    loglr += - h1h2.real # This is -0.5 * re(<h1|h2> + <h2|h1>)
         return loglr + self.lognl
 
     def _loglr(self):
@@ -541,10 +556,13 @@ def _loglr(self):
             # with detector response. Otherwise, skip detector response.
 
             if self.still_needs_det_response:
+                dtc = 0.
+
                 channel = wfs[ifo].numpy()
-                filter_i, norm_i = lik(freqs, 0.0, channel, h00,
+                filter_i, norm_i = lik(freqs, dtc, channel, h00,
                                        sdat['a0'], sdat['a1'],
                                        sdat['b0'], sdat['b1'])
+                self._current_wf_parts[ifo] = (dtc, channel, h00)
             else:
                 hp, hc = wfs[ifo]
                 det = self.det[ifo]

pycbc/inference/models/relbin_cpu.pyx

@@ -81,6 +81,37 @@ cpdef likelihood_parts_multi_v(double [::1] freqs,
         r0 = r0n
     return hd
 
+# Used for calculating the cross terms of
+# two signals when analyzing multiple signals
+# with no antenna response applied
+cpdef likelihood_parts_det_multi(double [::1] freqs,
+                     double dtc,
+                     double complex[::1] hp,
+                     double complex[::1] h00,
+                     double dtc2,
+                     double complex[::1] hp2,
+                     double complex[::1] h002,
+                     double complex[::1] a0,
+                     double complex[::1] a1,
+                     ) :
+    cdef size_t i
+    cdef double complex hd=0, r0, r0n, r1
+    cdef int N
+
+    N = freqs.shape[0]
+    for i in range(N):
+        r0n = conj((exp(-2.0j * 3.141592653 * dtc * freqs[i])
+               * (hp[i])) / h00[i])
+        r0n *= ((exp(-2.0j * 3.141592653 * dtc2 * freqs[i])
+               * (hp2[i])) / h002[i])
+        r1 = r0n - r0
+        if i > 0:
+            hd += a0[i-1] * r0 + a1[i-1] * r1
+
+        r0 = r0n
+    return hd
+
+
 # Standard likelihood
 cpdef likelihood_parts(double [::1] freqs,
                      double fp,
@@ -126,7 +157,7 @@ cpdef likelihood_parts_det(double [::1] freqs,
                      double [::1] b1,
                      ) :
     cdef size_t i
-    cdef double complex hd=0, r0, r0n, r1, x0, x1, x0n;
+    cdef double complex hd=0, r0, r0n, r1, x0=0, x1, x0n;
     cdef double hh=0
     cdef int N
 
@@ -147,6 +178,7 @@ cpdef likelihood_parts_det(double [::1] freqs,
         x0 = x0n
     return conj(hd), hh
 
+
 # Used where the antenna response may be frequency varying
 cpdef likelihood_parts_v(double [::1] freqs,
                      double[::1] fp,

pycbc/waveform/parameters.py

@@ -626,6 +626,9 @@ def docstr(self, prefix='', include_label=True):
 # behaviour
 td_required = ParameterList([f_lower, delta_t, approximant])
 fd_required = ParameterList([f_lower, delta_f, approximant])
+# The following is required for the FD sequence waveforms with detector
+# response already applied
+fd_det_sequence_required = ParameterList([f_lower, approximant])
 
 ####
 cbc_td_required = ParameterList([mass1, mass2, f_lower, delta_t, approximant])

pycbc/waveform/waveform.py

@@ -435,6 +435,7 @@ def props(obj, **kwargs):
     input_params = parse_mode_array(input_params)
     return input_params
 
+
 def check_args(args, required_args):
     """ check that required args are given """
     missing = []
@@ -547,16 +548,14 @@ def get_fd_det_waveform_sequence(template=None, **kwds):
         channels, values are FrequencySeries.
     """
     input_params = props(template, **kwds)
-    input_params['delta_f'] = -1
-    input_params['f_lower'] = -1
     if input_params['approximant'] not in fd_det_sequence:
         raise ValueError("Approximant %s not available" %
                             (input_params['approximant']))
     wav_gen = fd_det_sequence[input_params['approximant']]
     if hasattr(wav_gen, 'required'):
         required = wav_gen.required
     else:
-        required = parameters.fd_required
+        required = parameters.fd_det_sequence_required
     check_args(input_params, required)
     return wav_gen(**input_params)
 
@@ -734,8 +733,6 @@ def get_fd_det_waveform(template=None, **kwargs):
         domain. Keys are requested data channels, values are FrequencySeries.
     """
     input_params = props(template, **kwargs)
-    if 'f_lower' not in input_params:
-        input_params['f_lower'] = -1
     if input_params['approximant'] not in fd_det:
         raise ValueError("Approximant %s not available" %
                             (input_params['approximant']))