Merge branch 'detector_response' into 'master'

move around the detector response to the signal into a new function

See merge request Monash/peyote!15

peyote/detector.py

@@ -74,36 +74,45 @@ class Interferometer:
         detector_tensor = 0.5 * (np.einsum('i,j->ij', self.x, self.x) - np.einsum('i,j->ij', self.y, self.y))
         return detector_tensor
 
-    def inject_signal(self, waveform_generator):
+    def get_detector_response(self, waveform_polarizations, parameters):
         """
-        Inject a signal into noise.
-
-        Adds the requested signal to self.data
-
-        :param waveform_generator: waveform_generator type
+        Get the detector response for a particular waveform
 
+        :param waveform_polarizations: dict, polarizations of the waveform
+        :param parameters: dict, parameters describing position and time of arrival of the signal
+        :return: detector_response: signal observed in the interferometer
         """
-        signal = waveform_generator.frequency_domain_strain()
-
-        for mode in signal.keys():
+        signal = {}
+        for mode in waveform_polarizations.keys():
             det_response = self.antenna_response(
-                waveform_generator.parameters['ra'],
-                waveform_generator.parameters['dec'],
-                waveform_generator.parameters['geocent_time'],
-                waveform_generator.parameters['psi'], mode)
+                parameters['ra'],
+                parameters['dec'],
+                parameters['geocent_time'],
+                parameters['psi'], mode)
 
-            signal[mode] *= det_response
+            signal[mode] = waveform_polarizations[mode] * det_response
         signal_ifo = sum(signal.values())
 
         time_shift = self.time_delay_from_geocenter(
-            waveform_generator.parameters['ra'],
-            waveform_generator.parameters['dec'],
-            waveform_generator.parameters['geocent_time'])
-        signal_ifo = signal_ifo * np.exp(-1j * 2 * np.pi * (time_shift
-                                                            + waveform_generator.parameters['geocent_time'])
-                                         * waveform_generator.frequency_array)
-
-        self.data += signal_ifo
+            parameters['ra'],
+            parameters['dec'],
+            parameters['geocent_time'])
+
+        signal_ifo = signal_ifo * np.exp(-1j * 2 * np.pi * (time_shift + parameters['geocent_time'])
+                                         * self.frequency_array)
+
+        return signal_ifo
+
+    def inject_signal(self, waveform_polarizations, parameters):
+        """
+        Inject a signal into noise.
+
+        Adds the requested signal to self.data
+
+        :param waveform_polarizations: dict, polarizations of the waveform
+        :param parameters: dict, parameters describing position and time of arrival of the signal
+        """
+        self.data += self.get_detector_response(waveform_polarizations, parameters)
 
     def unit_vector_along_arm(self, arm):
         """

peyote/likelihood.py

@@ -22,27 +22,6 @@ class Likelihood(object):
         else:
             self.__noise_log_likelihood = noise_log_likelihood
 
-    def get_interferometer_signal(self, waveform_polarizations, interferometer):
-        h = []
-        for mode in waveform_polarizations:
-            det_response = interferometer.antenna_response(
-                self.waveform_generator.parameters['ra'],
-                self.waveform_generator.parameters['dec'],
-                self.waveform_generator.parameters['geocent_time'],
-                self.waveform_generator.parameters['psi'], mode)
-            h.append(waveform_polarizations[mode] * det_response)
-        signal = np.sum(h, axis=0)
-
-        time_shift = interferometer.time_delay_from_geocenter(
-            self.waveform_generator.parameters['ra'],
-            self.waveform_generator.parameters['dec'],
-            self.waveform_generator.parameters['geocent_time'])
-        signal = signal * np.exp(-1j * 2 * np.pi * (time_shift
-                                                    + self.waveform_generator.parameters['geocent_time'])
-                                 * self.waveform_generator.frequency_array)
-
-        return signal
-
     def log_likelihood(self):
         log_l = 0
         waveform_polarizations = self.waveform_generator.frequency_domain_strain()
@@ -51,7 +30,7 @@ class Likelihood(object):
         return log_l.real
 
     def log_likelihood_interferometer(self, waveform_polarizations, interferometer):
-        signal_ifo = self.get_interferometer_signal(waveform_polarizations, interferometer)
+        signal_ifo = interferometer.get_detector_response(waveform_polarizations, self.waveform_generator.parameters)
 
         log_l = - 2. / self.waveform_generator.time_duration * np.vdot(interferometer.data - signal_ifo,
                                                                        (interferometer.data - signal_ifo)

peyote/utils.py

+from __future__ import division
 import logging
 import numpy as np
 from astropy.time import Time

tutorials/BasicTutorial.py

@@ -10,7 +10,7 @@ peyote.utils.setup_logger()
 time_duration = 1.
 sampling_frequency = 4096.
 
-simulation_parameters = dict(
+injection_parameters = dict(
     mass_1=36.,
     mass_2=29.,
     spin11=0,
@@ -29,19 +29,19 @@ simulation_parameters = dict(
     geocent_time=1126259642.413,
     psi=2.659
 )
+
 # Create the waveform_generator using a LAL BinaryBlackHole source function
 waveform_generator = peyote.waveform_generator.WaveformGenerator(
     sampling_frequency=sampling_frequency,
     time_duration=time_duration,
     source_model=peyote.source.lal_binary_black_hole,
-    parameters=simulation_parameters)
-
+    parameters=injection_parameters)
 hf_signal = waveform_generator.frequency_domain_strain()
 
 sampling_parameters = peyote.parameter.Parameter.\
-    parse_floats_to_parameters(simulation_parameters)
+   parse_floats_to_parameters(injection_parameters)
 
-#sampling_parameters = peyote.parameter.Parameter.\
+# sampling_parameters = peyote.parameter.Parameter.\
 #    parse_keys_to_parameters(simulation_parameters.keys())
 
 
@@ -49,15 +49,21 @@ sampling_parameters = peyote.parameter.Parameter.\
 H1 = peyote.detector.H1
 H1.set_data(sampling_frequency=sampling_frequency, duration=time_duration,
             from_power_spectral_density=True)
-H1.inject_signal(waveform_generator)
+H1.inject_signal(waveform_polarizations=hf_signal, parameters=injection_parameters)
 
 # Simulate the data in L1
 L1 = peyote.detector.L1
 L1.set_data(sampling_frequency=sampling_frequency, duration=time_duration,
             from_power_spectral_density=True)
-L1.inject_signal(waveform_generator)
+L1.inject_signal(waveform_polarizations=hf_signal, parameters=injection_parameters)
+
+# Simulate the data in V1
+V1 = peyote.detector.V1
+V1.set_data(sampling_frequency=sampling_frequency, duration=time_duration,
+            from_power_spectral_density=True)
+V1.inject_signal(waveform_polarizations=hf_signal, parameters=injection_parameters)
 
-IFOs = [H1, L1]
+IFOs = [H1, L1, V1]
 
 # Plot the noise and signal
 fig, ax = plt.subplots()
@@ -77,12 +83,12 @@ likelihood = peyote.likelihood.Likelihood(IFOs, waveform_generator)
 # New way way of doing it, still not perfect
 sampling_parameters['mass_1'].prior = peyote.prior.Uniform(lower=35, upper=37)
 sampling_parameters['luminosity_distance'].prior = peyote.prior.Uniform(lower=30, upper=200)
-# sampling_parameters['geocent_time'].prior = peyote.prior.Uniform(lower=1126259642.413 - 0.1,
-#                                                                    upper=1126259642.413 + 0.1)
+#sampling_parameters["geocent_time"].prior = peyote.prior.Uniform(lower=injection_parameters["geocent_time"] - 0.1,
+#                                                                  upper=injection_parameters["geocent_time"]+0.1)
 
 result = peyote.sampler.run_sampler(likelihood, priors=sampling_parameters, sampler='nestle', verbose=True)
 print(result)
-truths = [simulation_parameters[x] for x in result.search_parameter_keys]
+truths = [injection_parameters[x] for x in result.search_parameter_keys]
 
 
 fig = corner.corner(result.samples, truths=truths, labels=result.search_parameter_keys)

tutorials/make_standard_data.py

@@ -41,7 +41,7 @@ IFO = peyote.detector.H1
 IFO.set_data(
     from_power_spectral_density=True, sampling_frequency=sampling_frequency,
     duration=time_duration)
-IFO.inject_signal(waveform_generator)
+IFO.inject_signal(waveform_polarizations=waveform_generator.frequency_domain_strain(), parameters=simulation_parameters)
 hf_signal_and_noise = IFO.data
 frequencies = peyote.utils.create_fequency_series(
     sampling_frequency=sampling_frequency, duration=time_duration)