Merge branch 'custom-interpolation' into 'master'

Add custom interpolation

See merge request !917

bilby/core/utils.py

@@ -925,7 +925,10 @@ class Counter(object):
 class UnsortedInterp2d(interp2d):
 
     def __call__(self, x, y, dx=0, dy=0, assume_sorted=False):
-        """  Wrapper to scipy.interpolate.interp2d which preserves the input ordering.
+        """  Modified version of the interp2d call method.
+
+        This avoids the outer product that is done when two numpy
+        arrays are passed.
 
         Parameters
         ----------
@@ -942,8 +945,34 @@ class UnsortedInterp2d(interp2d):
         array_like: See superclass
 
         """
-        unsorted_idxs = np.argsort(np.argsort(x))
-        return super(UnsortedInterp2d, self).__call__(x, y, dx=dx, dy=dy, assume_sorted=False)[unsorted_idxs]
+        from scipy.interpolate.dfitpack import bispeu
+        out_of_bounds_x = (x < self.x_min) | (x > self.x_max)
+        out_of_bounds_y = (y < self.y_min) | (y > self.y_max)
+        bad = out_of_bounds_x | out_of_bounds_y
+        if isinstance(x, float) and isinstance(y, float):
+            if bad:
+                output = self.fill_value
+                ier = 0
+            else:
+                output, ier = bispeu(*self.tck, x, y)
+        else:
+            if isinstance(x, np.ndarray):
+                output = np.zeros_like(x)
+                x_ = x[~bad]
+            else:
+                x_ = x * np.ones_like(y)
+            if isinstance(y, np.ndarray):
+                output = np.zeros_like(y)
+                y_ = y[~bad]
+            else:
+                y_ = y * np.ones_like(x)
+            output[bad] = self.fill_value
+            output[~bad], ier = bispeu(*self.tck, x_, y_)
+        if ier == 10:
+            raise ValueError("Invalid input data")
+        elif ier:
+            raise TypeError("An error occurred")
+        return output
 
 
 #  Instantiate the default argument parser at runtime