stable version of dist marg, and time-and-phase marg. Need to fix dist

and time marg before moving onto dist and tiem nad phase marg

examples/injection_examples/basic_tutorial.py

@@ -19,13 +19,13 @@ 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(170)
+np.random.seed(2157034243)
 
 # 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),
 # spins of both black holes (a, tilt, phi), etc.
 injection_parameters = dict(mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0, phi_12=1.7, phi_jl=0.3,
-                            luminosity_distance=4000., iota=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413,
+                            luminosity_distance=2000., iota=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413,
                             waveform_approximant='IMRPhenomPv2', reference_frequency=50., ra=1.375, dec=-1.2108)
 
 # Create the waveform_generator using a LAL BinaryBlackHole source function

tupak/likelihood.py

@@ -85,7 +85,7 @@ class GravitationalWaveTransient(Likelihood):
         if self.phase_marginalization:
             self.bessel_function_interped = None
             self.setup_phase_marginalization()
-            prior['psi'] = 0
+            prior['phase'] = 0
 
     @property
     def prior(self):
@@ -141,25 +141,24 @@ class GravitationalWaveTransient(Likelihood):
         if self.time_marginalization:
 
             delta_tc = self.waveform_generator.sampling_frequency
-            #print (1./delta_tc, (interferometer.duration) )
-            if self.distance_marginalization:
 
-                rho_opt_1 = optimal_snr_squared.real * \
-                             self.waveform_generator.parameters['luminosity_distance'] ** 2
-                rho_mf_1_tc_array = matched_filter_snr_squared_tc_array.real * \
-                            self.waveform_generator.parameters['luminosity_distance']
+            if self.distance_marginalization:
 
-                ## TRY time marg before dist marg:
-                #time_marged_log_l_dist_array = logsumexp(rho_mf_1_tc_array, b=delta_tc)
-                #log_l = self.interp_dist_margd_likelihood(rho_mf_1_tc_array, rho_opt_1)[0]
+                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 * \
+                            self.waveform_generator.parameters['luminosity_distance'] \
+                            / self.ref_dist
 
-                dist_marged_log_l_tc_array = self.interp_dist_margd_likelihood(rho_mf_1_tc_array, rho_opt_1)[0]
+                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)
+
             elif self.phase_marginalization:
 
-                phase_marged_log_l_tc_array = self.bessel_function_interped(abs(matched_filter_snr_squared_tc_array))
-                log_l = logsumexp(phase_marged_log_l_tc_array, axis=0, b=delta_tc)
+                log_l = logsumexp( self.bessel_function_interped(abs(matched_filter_snr_squared_tc_array)), b=delta_tc)\
+                                                    - optimal_snr_squared/2.
 
             else:
 
@@ -168,12 +167,14 @@ class GravitationalWaveTransient(Likelihood):
         elif self.distance_marginalization:
             #TODO: order should be:  (phase-> dist-> tc)
 
-            rho_opt_1 = optimal_snr_squared.real * \
-                         self.waveform_generator.parameters['luminosity_distance'] ** 2
-            rho_mf_1 = matched_filter_snr_squared.real * \
-                        self.waveform_generator.parameters['luminosity_distance']
+            rho_opt_ref = optimal_snr_squared.real * \
+                         self.waveform_generator.parameters['luminosity_distance'] ** 2 \
+                        / self.ref_dist**2.
+            rho_mf_ref = matched_filter_snr_squared.real * \
+                        self.waveform_generator.parameters['luminosity_distance'] \
+                        / self.ref_dist
 
-            log_l=self.interp_dist_margd_likelihood(rho_mf_1, rho_opt_1)[0]
+            log_l=self.interp_dist_margd_loglikelihood(rho_mf_ref, rho_opt_ref)[0]
 
         elif self.phase_marginalization:
             matched_filter_snr_squared = self.bessel_function_interped(abs(matched_filter_snr_squared))
@@ -198,48 +199,53 @@ class GravitationalWaveTransient(Likelihood):
                                               for distance in self.distance_array])
 
         ### Make the lookup table ###
+        self.ref_dist = 1000 # 1000 Mpc
+        self.dist_margd_loglikelihood_array = np.zeros((400,800))
 
-        self.dist_margd_likelihood_array = np.zeros((1000,1000))
+        self.rho_opt_ref_array = np.logspace(-3,4,self.dist_margd_loglikelihood_array.shape[0]) # optimal filter snr at fiducial distance of 1 Mpc
+        self.rho_mf_ref_array = np.hstack( (-np.logspace(2,-3,self.dist_margd_loglikelihood_array.shape[1]/2), np.logspace(-3,4,self.dist_margd_loglikelihood_array.shape[1]/2)) ) # matched filter snr at fiducial distance of 1 Mpc
 
-        self.rho_opt_1_array = np.linspace(5e5,7e8,self.dist_margd_likelihood_array.shape[0]) # optimal filter snr at fiducial distance of 1 Mpc
-        self.rho_mf_1_array = np.linspace(-45000,2e6,self.dist_margd_likelihood_array.shape[1])  # matched filter snr at fiducial distance of 1 Mpc
 
+        for ii, rho_opt_ref in enumerate(self.rho_opt_ref_array):
 
-        for ii, rho_opt_1 in enumerate(self.rho_opt_1_array):
+            for jj, rho_mf_ref in enumerate(self.rho_mf_ref_array):
 
-            for jj, rho_mf_1 in enumerate(self.rho_mf_1_array):
-
-                optimal_snr_squared_array = rho_opt_1 \
+                optimal_snr_squared_array = rho_opt_ref * self.ref_dist**2. \
                                             / self.distance_array ** 2
 
-                matched_filter_snr_squared_array = rho_mf_1 \
+                matched_filter_snr_squared_array = rho_mf_ref * self.ref_dist \
                                                 / self.distance_array
 
-                self.dist_margd_likelihood_array[ii][jj] = \
+                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)
 
         log_norm = logsumexp(0 / self.distance_array  - 0/self.distance_array**2. / 2, b=self.distance_prior_array * self.delta_distance)
-        self.dist_margd_likelihood_array -= log_norm
+        self.dist_margd_loglikelihood_array -= log_norm
 
 
-        self.interp_dist_margd_likelihood = interp2d(self.rho_mf_1_array, self.rho_opt_1_array, self.dist_margd_likelihood_array)
+        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_likelihood_array",self.dist_margd_likelihood_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 'psi' not in self.prior.keys() or not isinstance(self.prior['psi'], tupak.prior.Prior):
-            logging.info('No prior provided for polarization, using default prior.')
-            self.prior['psi'] = tupak.prior.create_default_prior('psi')
+        if 'phase' not in self.prior.keys() or not isinstance(self.prior['phase'], tupak.prior.Prior):
+            logging.info('No prior provided for phase at coalescence, using default prior.')
+            self.prior['phase'] = tupak.prior.create_default_prior('phase')
         self.bessel_function_interped = interp1d(np.linspace(0, 1e6, int(1e5)),
                                                  np.log([i0e(snr) for snr in np.linspace(0, 1e6, int(1e5))])
                                                  + np.linspace(0, 1e6, int(1e5)),
                                                  bounds_error=False, fill_value=-np.inf)
 
 
+
 class BasicGravitationalWaveTransient(Likelihood):
     """ A basic gravitational wave transient likelihood