start to split the likelihood into bitesize chunks

The current likelihood function does a lot. I think we could split it into neater chunks that will make it easier to use and adapt if we want to start analytically marginalising some parameters.

The current version shouldn't break anything but I feel like this should probably get some discussion before being merged.

@gregory.ashton @paul-lasky @moritz.huebner

peyote/likelihood.py

@@ -6,30 +6,36 @@ class Likelihood:
         self.interferometers = interferometers
         self.source = source
 
+    def get_interferometer_signal(self, waveform_polarizations, interferometer):
+        h = []
+        for mode in waveform_polarizations:
+            det_response = interferometer.antenna_response(
+                self.source.ra, self.source.dec,
+                self.source.geocent_time, self.source.psi, mode)
+
+            h.append(waveform_polarizations[mode] * det_response)
+        signal = np.sum(h, axis=0)
+
+        time_shift = interferometer.time_delay_from_geocenter(
+            self.source.ra, self.source.dec,
+            self.source.geocent_time)
+        signal *= np.exp(-1j * 2 * np.pi * time_shift * self.source.frequency_array)
+
+        return signal
+
     def log_likelihood(self):
         log_l = 0
         waveform_polarizations = self.source.frequency_domain_strain()
         for interferometer in self.interferometers:
-            h = []
-            for mode in waveform_polarizations:
-                det_response = interferometer.antenna_response(
-                    self.source.ra, self.source.dec,
-                    self.source.geocent_time, self.source.psi, mode)
-
-                h.append(waveform_polarizations[mode] * det_response)
-
-            signal_ifo = np.sum(h, axis=0)
-
-            time_shift = interferometer.time_delay_from_geocenter(
-                self.source.ra, self.source.dec,
-                self.source.geocent_time)
-            signal_ifo *= np.exp(-1j*2*np.pi*time_shift*self.source.frequency_array)
+            log_l += self.log_likelihood_interferometer(waveform_polarizations, interferometer)
+        return log_l.real
 
-            log_l -= 4. / self.source.time_duration * np.vdot(
-                interferometer.data - signal_ifo,
-                (interferometer.data - signal_ifo) / (
-                    interferometer.power_spectral_density_array))
+    def log_likelihood_interferometer(self, waveform_polarizations, interferometer):
+        signal_ifo = self.get_interferometer_signal(waveform_polarizations, interferometer)
 
+        log_l = - 4. / self.source.time_duration * np.vdot(interferometer.data - signal_ifo,
+                                                           (interferometer.data - signal_ifo)
+                                                           / interferometer.power_spectral_density_array)
         return log_l.real