BUG(core): fix extrapolation in trapz_product() (#105)

Fixes the incorrect extrapolation in `glass.core.array.trapz_product()`
by making sure the interval over which all functions are interpolated
telescopes down to the union of the input intervals. This caused a bug
in `glass.points.effective_bias()`.

Fixes: #104
Fixed: Incorrect extrapolation in `glass.core.array.trapz_product()`,
  causing a bug in `glass.points.effective_bias()`

glass/core/array.py

@@ -53,7 +53,8 @@ def trapz_product(f, *ff, axis=-1):
     '''trapezoidal rule for a product of functions'''
     x, _ = f
     for x_, _ in ff:
-        x = np.union1d(x, x_[(x_ > x[0]) & (x_ < x[-1])])
+        x = np.union1d(x[(x >= x_[0]) & (x <= x_[-1])],
+                       x_[(x_ >= x[0]) & (x_ <= x[-1])])
     y = np.interp(x, *f)
     for f_ in ff:
         y *= np.interp(x, *f_)

glass/core/test/test_array.py

@@ -90,3 +90,17 @@ def test_ndinterp():
                             [[[2.15], [2.25], [2.35], [2.45]],
                              [[2.45], [2.35], [2.25], [2.15]],
                              [[2.15], [2.45], [2.25], [2.35]]]], atol=1e-15)
+
+
+def test_trapz_product():
+    from glass.core.array import trapz_product
+
+    x1 = np.linspace(0, 2, 100)
+    f1 = np.full_like(x1, 2.0)
+
+    x2 = np.linspace(1, 2, 10)
+    f2 = np.full_like(x2, 0.5)
+
+    s = trapz_product((x1, f1), (x2, f2))
+
+    assert np.allclose(s, 1.0)