working (dist,time,phase) marg

examples/injection_examples/basic_tutorial.py

@@ -19,7 +19,7 @@ label = 'basic_tutorial'
 tupak.utils.setup_logger(outdir=outdir, label=label)
 
 # Set up a random seed for result reproducibility.  This is optional!
-np.random.seed(2157034243)
+np.random.seed(88170235)
 
 # We are going to inject a binary black hole waveform.  We first establish a dictionary of parameters that
 # includes all of the different waveform parameters, including masses of the two black holes (mass_1, mass_2),
@@ -55,7 +55,7 @@ for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl','psi', 'ra', 'd
 priors['luminosity_distance'] = tupak.prior.create_default_prior(name='luminosity_distance')
 
 # Initialise the likelihood by passing in the interferometer data (IFOs) and the waveoform generator
-likelihood = tupak.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator,time_marginalization=False, distance_marginalization=True, prior=priors)
+likelihood = tupak.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator,time_marginalization=False, phase_marginalization=False, distance_marginalization=False, prior=priors)
 
 # Run sampler.  In this case we're going to use the `dynesty` sampler
 result = tupak.sampler.run_sampler(likelihood=likelihood, priors=priors, sampler='dynesty', npoints=1000,

tupak/detector.py

@@ -69,6 +69,7 @@ class Interferometer(object):
         self.frequency_array = []
         self.sampling_frequency = None
         self.duration = None
+        self.time_marginalization = False
 
     @property
     def minimum_frequency(self):
@@ -226,7 +227,7 @@ class Interferometer(object):
             det_response = self.antenna_response(
                 parameters['ra'],
                 parameters['dec'],
-                parameters['geocent_time'],
+                self.epoch,#parameters['geocent_time'],
                 parameters['psi'], mode)
 
             signal[mode] = waveform_polarizations[mode] * det_response
@@ -237,9 +238,14 @@ class Interferometer(object):
         time_shift = self.time_delay_from_geocenter(
             parameters['ra'],
             parameters['dec'],
-            parameters['geocent_time'])
+            self.epoch)#parameters['geocent_time'])
+
+        if self.time_marginalization:
+            dt = time_shift # when marginalizing over time we only care about relative time shifts between detectors and marginalized over
+                            # all candidate coalescence times
+        else:
+            dt = self.epoch - (parameters['geocent_time'] - time_shift)
 
-        dt = self.epoch - (parameters['geocent_time'] - time_shift)
         signal_ifo = signal_ifo * np.exp(
                 -1j * 2 * np.pi * dt * self.frequency_array)
 
@@ -688,7 +694,7 @@ def get_interferometer_with_fake_noise_and_injection(
     interferometer = get_empty_interferometer(name)
     interferometer.set_data(
         sampling_frequency=sampling_frequency, duration=time_duration,
-        from_power_spectral_density=True)
+        from_power_spectral_density=True, epoch=(injection_parameters['geocent_time']+2)-time_duration)
     interferometer.inject_signal(
         waveform_polarizations=injection_polarizations,
         parameters=injection_parameters)

tupak/likelihood.py

@@ -72,9 +72,6 @@ class GravitationalWaveTransient(Likelihood):
         self.phase_marginalization = phase_marginalization
         self.prior = prior
 
-        #if self.time_marginalization:
-        #    prior['geocent_time'] = 0 #FIXME: should be fixed to injection value; probably easiest to pass in a delta function prior
-
         if self.distance_marginalization:
             self.distance_array = np.array([])
             self.delta_distance = 0
@@ -127,45 +124,57 @@ class GravitationalWaveTransient(Likelihood):
 
             if self.time_marginalization:
 
+                interferometer.time_marginalization = self.time_marginalization
+
                 if matched_filter_snr_squared_tc_array is None:
 
                     matched_filter_snr_squared_tc_array = 4. * (1./interferometer.duration) \
                                                             * np.fft.ifft( signal_ifo.conjugate()[0:-1] * interferometer.data[0:-1] \
-                                                            / interferometer.power_spectral_density_array[0:-1]).real * len(interferometer.data[0:-1])
+                                                            / interferometer.power_spectral_density_array[0:-1]) * len(interferometer.data[0:-1])
                 else:
 
                     matched_filter_snr_squared_tc_array += 4. * (1./interferometer.duration) \
                                                         * np.fft.ifft( signal_ifo.conjugate()[0:-1] * interferometer.data[0:-1] \
-                                                        / interferometer.power_spectral_density_array[0:-1] ).real * len(interferometer.data[0:-1])
+                                                        / interferometer.power_spectral_density_array[0:-1] ) * len(interferometer.data[0:-1])
 
         if self.time_marginalization:
 
-            delta_tc = self.waveform_generator.sampling_frequency
+            delta_tc = 1./self.waveform_generator.sampling_frequency
+            tc_log_norm = np.log(interferometer.duration*delta_tc)
 
             if self.distance_marginalization:
 
                 rho_opt_ref = optimal_snr_squared.real * \
                              self.waveform_generator.parameters['luminosity_distance'] ** 2 \
                             / self.ref_dist**2.
-                rho_mf_ref_tc_array = matched_filter_snr_squared_tc_array.real * \
+                rho_mf_ref_tc_array = matched_filter_snr_squared_tc_array * \
                             self.waveform_generator.parameters['luminosity_distance'] \
                             / self.ref_dist
 
-                dist_marged_log_l_tc_array = self.interp_dist_margd_loglikelihood(rho_mf_ref_tc_array, rho_opt_ref)[0]
 
-                log_l = logsumexp(dist_marged_log_l_tc_array, axis=0, b=delta_tc)
+                if self.phase_marginalization:
+
+                    phase_margd_rho_mf_tc_array = self.bessel_function_interped(abs(rho_mf_ref_tc_array))
+
+                    dist_marged_log_l_tc_array = self.interp_dist_margd_loglikelihood(phase_margd_rho_mf_tc_array, rho_opt_ref)
+                    log_l = logsumexp(dist_marged_log_l_tc_array, axis=0, b=delta_tc) - tc_log_norm
+
+                else:
+
+                    dist_marged_log_l_tc_array = self.interp_dist_margd_loglikelihood(rho_mf_ref_tc_array.real, rho_opt_ref)
+
+                    log_l = logsumexp(dist_marged_log_l_tc_array, axis=0, b=delta_tc) - dist_margd_loglikelihood_array
 
             elif self.phase_marginalization:
 
                 log_l = logsumexp( self.bessel_function_interped(abs(matched_filter_snr_squared_tc_array)), b=delta_tc)\
-                                                    - optimal_snr_squared/2.
+                                                    - optimal_snr_squared/2. - tc_log_norm
 
             else:
 
-                log_l = logsumexp(matched_filter_snr_squared_tc_array, axis=0, b=delta_tc) - optimal_snr_squared / 2
+                log_l = logsumexp(matched_filter_snr_squared_tc_array.real, axis=0, b=delta_tc) - optimal_snr_squared / 2. - tc_log_norm
 
         elif self.distance_marginalization:
-            #TODO: order should be:  (phase-> dist-> tc)
 
             rho_opt_ref = optimal_snr_squared.real * \
                          self.waveform_generator.parameters['luminosity_distance'] ** 2 \
@@ -219,21 +228,15 @@ class GravitationalWaveTransient(Likelihood):
                 self.dist_margd_loglikelihood_array[ii][jj] = \
                                                     logsumexp(matched_filter_snr_squared_array - \
                                                      optimal_snr_squared_array / 2, \
-                                                     b=self.distance_prior_array * self.delta_distance)
+                                                     b= self.distance_prior_array * self.delta_distance)
 
-        log_norm = logsumexp(0 / self.distance_array  - 0/self.distance_array**2. / 2, b=self.distance_prior_array * self.delta_distance)
+        log_norm = logsumexp(0. / self.distance_array  - 0./self.distance_array**2.,\
+                            b=self.distance_prior_array* self.delta_distance)
         self.dist_margd_loglikelihood_array -= log_norm
 
 
         self.interp_dist_margd_loglikelihood = interp2d(self.rho_mf_ref_array, self.rho_opt_ref_array, self.dist_margd_loglikelihood_array)
 
-        #np.save("dist_margd_loglikelihood_array",self.dist_margd_loglikelihood_array)
-        #sys.exit()
-
-    '''def setupt_time_marginalization(self):
-        if 'geocent_time' not in self.prior.keys() or not isinstance(self.prior['geocent_time'], tupak.prior.Prior):
-            logging.info('No prior provided for geocent_time, using default prior.')
-            self.prior['geocent_time'] ='''
 
     def setup_phase_marginalization(self):
         if 'phase' not in self.prior.keys() or not isinstance(self.prior['phase'], tupak.prior.Prior):

tupak/prior.py

@@ -418,9 +418,9 @@ def create_default_prior(name):
         Default prior distribution for that parameter, if unknown None is returned.
     """
     default_priors = {
-        'mass_1': Uniform(name=name, minimum=5, maximum=100),
-        'mass_2': Uniform(name=name, minimum=5, maximum=100),
-        'chirp_mass': Uniform(name=name, minimum=5, maximum=100),
+        'mass_1': Uniform(name=name, minimum=20, maximum=100),
+        'mass_2': Uniform(name=name, minimum=20, maximum=100),
+        'chirp_mass': Uniform(name=name, minimum=25, maximum=100),
         'total_mass': Uniform(name=name, minimum=10, maximum=200),
         'mass_ratio': Uniform(name=name, minimum=0.125, maximum=1),
         'symmetric_mass_ratio': Uniform(name=name, minimum=8 / 81, maximum=0.25),