bilby/core/prior.py

@@ -66,12 +66,17 @@ class PriorDict(OrderedDict):
         ----------
         filename: str
             Name of the file to be read in
+
+        Notes
+        -----
+        Lines beginning with '#' or empty lines will be ignored.
         """
 
-        prior = {}
+        comments = ['#', '\n']
+        prior = dict()
         with open(filename, 'r') as f:
             for line in f:
-                if line[0] == '#':
+                if line[0] in comments:
                     continue
                 elements = line.split('=')
                 key = elements[0].replace(' ', '')
@@ -240,6 +245,13 @@ class PriorDict(OrderedDict):
         """Empty redundancy test, should be overwritten in subclasses"""
         return False
 
+    def copy(self):
+        """
+        We have to overwrite the copy method as it fails due to the presence of
+        defaults.
+        """
+        return self.__class__(dictionary=OrderedDict(self))
+
 
 class PriorSet(PriorDict):
 
@@ -653,7 +665,7 @@ class PowerLaw(Prior):
             normalising = (1 + self.alpha) / (self.maximum ** (1 + self.alpha) -
                                               self.minimum ** (1 + self.alpha))
 
-        return (self.alpha * np.log(val) + np.log(normalising)) + np.log(1. * self.is_in_prior_range(val))
+        return (self.alpha * np.nan_to_num(np.log(val)) + np.log(normalising)) + np.log(1. * self.is_in_prior_range(val))
 
 
 class Uniform(Prior):

bilby/core/utils.py

@@ -9,6 +9,7 @@ import inspect
 import subprocess
 
 import numpy as np
+from scipy.interpolate import interp2d
 
 logger = logging.getLogger('bilby')
 
@@ -692,6 +693,19 @@ def run_commandline(cl, log_level=20, raise_error=True, return_output=True):
         process.communicate()
 
 
+class UnsortedInterp2d(interp2d):
+    """
+    Wrapper to scipy.interpolate.interp2d which preserves the input ordering.
+
+    See https://stackoverflow.com/questions/44941271/scipy-interp2d-returned-function-sorts-input-argument-automatically-and-undesira
+    for the implementation details.
+    """
+
+    def __call__(self, x, y, dx=0, dy=0):
+        unsorted_idxs = np.argsort(np.argsort(x))
+        return interp2d.__call__(self, x, y, dx=dx, dy=dy)[unsorted_idxs]
+
+
 #  Instantiate the default argument parser at runtime
 command_line_args, command_line_parser = set_up_command_line_arguments()
 #  Instantiate the default logging

bilby/gw/detector.py

@@ -1280,9 +1280,7 @@ class Interferometer(object):
         signal_ifo *= self.strain_data.frequency_mask
 
         time_shift = self.time_delay_from_geocenter(
-            parameters['ra'],
-            parameters['dec'],
-            self.strain_data.start_time)
+            parameters['ra'], parameters['dec'], parameters['geocent_time'])
         dt = parameters['geocent_time'] + time_shift - self.strain_data.start_time
 
         signal_ifo = signal_ifo * np.exp(

bilby/gw/likelihood.py

 from __future__ import division
 import numpy as np
-from scipy.interpolate import interp2d, interp1d
+from scipy.interpolate import interp1d
 
 try:
     from scipy.special import logsumexp
@@ -9,7 +9,7 @@ except ImportError:
 from scipy.special import i0e
 
 from ..core import likelihood
-from ..core.utils import logger
+from ..core.utils import logger, UnsortedInterp2d
 from ..core.prior import Prior, Uniform
 from .detector import InterferometerList
 from .prior import BBHPriorDict
@@ -182,13 +182,11 @@ class GravitationalWaveTransient(likelihood.Likelihood):
                 if self.phase_marginalization:
                     dist_marged_log_l_tc_array = self._interp_dist_margd_loglikelihood(
                         abs(rho_mf_ref_tc_array), rho_opt_ref)
-                    log_l = logsumexp(dist_marged_log_l_tc_array,
-                                      b=self.time_prior_array)
                 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,
-                                      b=self.time_prior_array)
+                log_l = logsumexp(dist_marged_log_l_tc_array,
+                                  b=self.time_prior_array)
             elif self.phase_marginalization:
                 log_l = logsumexp(self._bessel_function_interped(abs(
                     matched_filter_snr_squared_tc_array)),
@@ -249,8 +247,9 @@ class GravitationalWaveTransient(likelihood.Likelihood):
 
     def _setup_distance_marginalization(self):
         self._create_lookup_table()
-        self._interp_dist_margd_loglikelihood = interp2d(self._rho_mf_ref_array, self._rho_opt_ref_array,
-                                                         self._dist_margd_loglikelihood_array)
+        self._interp_dist_margd_loglikelihood = UnsortedInterp2d(
+            self._rho_mf_ref_array, self._rho_opt_ref_array,
+            self._dist_margd_loglikelihood_array)
 
     def _create_lookup_table(self):
         """ Make the lookup table """

test/gw_likelihood_test.py

@@ -43,7 +43,7 @@ class TestBasicGWTransient(unittest.TestCase):
         """Test log likelihood matches precomputed value"""
         self.likelihood.log_likelihood()
         self.assertAlmostEqual(self.likelihood.log_likelihood(),
-                               -4048.5113284032036, 3)
+                               -4054.2229111227016, 3)
 
     def test_log_likelihood_ratio(self):
         """Test log likelihood ratio returns the correct value"""
@@ -113,7 +113,7 @@ class TestGWTransient(unittest.TestCase):
         """Test log likelihood matches precomputed value"""
         self.likelihood.log_likelihood()
         self.assertAlmostEqual(self.likelihood.log_likelihood(),
-                               -4048.5113284032036, 3)
+                               -4054.2229111227016, 3)
 
     def test_log_likelihood_ratio(self):
         """Test log likelihood ratio returns the correct value"""