Repository: sigma-py/quadpy (language: python)

Issue:
Inconsistency in c2 integration?
Thank you very much for the great package. I am very new to it, it might be that I do not have the right understanding. I am testing it for applications in a numerical PDE environment. I just look at two 2D test function supplied in **quadpy/examples_test.py** with minimal changes. The function **test_t2()** returns the expected results: the integral value on a single triangle is the same as for identically copied triangles. The function **test_c2()** has a different behavior: the integral value on a single square is just the half of the value returned for identically copied squares. The value for the single square is right. Is that okay?

Here are the two functions (quadpy 0.16.5, python 3.7):
```python
import numpy
import quadpy

def test_t2():
    val = quadpy.t2.schemes["cubtri"]().integrate(
        lambda x: numpy.exp(x[0]), numpy.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]])
    )
    print('single triangle:', val)
    
    val = quadpy.t2.schemes["cubtri"]().integrate(
        lambda x: [numpy.exp(x[0]), numpy.exp(x[1])],
        numpy.stack(
            [
                [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]],
                [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]],
                [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]],
                [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]],
                [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]],
            ],
            axis=-2,
        ),
    )
    print('multiple triangles:\n', val)
    assert val.shape == (2, 5)


def test_c2():
    val= quadpy.c2.schemes["stroud_c2_5_4"]().integrate(
        lambda x: numpy.exp(x[0]),
        quadpy.c2.rectangle_points([0.0, 1.0], [0.0, 1.0]),
    )
    print('single square:', val)

    val = quadpy.c2.schemes["stroud_c2_5_4"]().integrate(
        lambda x: [numpy.exp(x[0]), numpy.exp(x[1])],
        numpy.stack(
            [
                quadpy.c2.rectangle_points([0.0, 1.0], [0.0, 1.0]),
                quadpy.c2.rectangle_points([0.0, 1.0], [0.0, 1.0]),
                quadpy.c2.rectangle_points([0.0, 1.0], [0.0, 1.0]),
                quadpy.c2.rectangle_points([0.0, 1.0], [0.0, 1.0]),
                quadpy.c2.rectangle_points([0.0, 1.0], [0.0, 1.0]),
            ],
            axis=-2,
        ),
    )
    print('multiple squares:\n', val)
    assert val.shape == (2, 5)


test_t2()
test_c2()
```
And the result is:
**single triangle:** 0.7182818284626384
**multiple triangles:**
 [[0.71828183 0.71828183 0.71828183 0.71828183 0.71828183]
 [0.71828183 0.71828183 0.71828183 0.71828183 0.71828183]]

**single square:** 1.718281004372522
**multiple squares:**
 [[3.43656201 3.43656201 3.43656201 3.43656201 3.43656201]
 [3.43656201 3.43656201 3.43656201 3.43656201 3.43656201]]

I guess the last values should also be 1.718....
Thank you very much for any help.



Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index fc3a2f28d..343f68582 100644
--- a/README.md
+++ b/README.md
@@ -51,12 +51,12 @@ pip install quadpy
 ```
 and do
 ```python
-import numpy
+import numpy as np
 import quadpy
 
 
 def f(x):
-    return numpy.sin(x) - x
+    return np.sin(x) - x
 
 
 val, err = quadpy.quad(f, 0.0, 6.0)
@@ -68,15 +68,15 @@ and "intervals" in spaces of arbitrary dimension.
 
 To integrate over a _triangle_, do
 ```python
-import numpy
+import numpy as np
 import quadpy
 
 
 def f(x):
-    return numpy.sin(x[0]) * numpy.sin(x[1])
+    return np.sin(x[0]) * np.sin(x[1])
 
 
-triangle = numpy.array([[0.0, 0.0], [1.0, 0.0], [0.7, 0.5]])
+triangle = np.array([[0.0, 0.0], [1.0, 0.0], [0.7, 0.5]])
 
 # get a "good" scheme of degree 10
 scheme = quadpy.t2.get_good_scheme(10)
@@ -111,7 +111,7 @@ domains at once, you can provide them all in one `integrate()` call, e.g.,
 <!--exdown-skip-->
 ```python
 # shape (3, 5, 2), i.e., (corners, num_triangles, xy_coords)
-triangles = numpy.stack(
+triangles = np.stack(
     [
         [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]],
         [[1.2, 0.6], [1.3, 0.7], [1.4, 0.8]],
@@ -125,7 +125,7 @@ triangles = numpy.stack(
 The same goes for functions with vectorized output, e.g.,
 ```python
 def f(x):
-    return [numpy.sin(x[0]), numpy.sin(x[1])]
+    return [np.sin(x[0]), np.sin(x[1])]
 ```
 
 More examples under [test/examples_test.py](test/examples_test.py).
@@ -162,12 +162,12 @@ for how to generate Gauss formulas for your own weight functions.
 
 Example:
 ```python
-import numpy
+import numpy as np
 import quadpy
 
 scheme = quadpy.c1.gauss_patterson(5)
 scheme.show()
-val = scheme.integrate(lambda x: numpy.exp(x), [0.0, 1.0])
+val = scheme.integrate(lambda x: np.exp(x), [0.0, 1.0])
 ```
 
 ### 1D half-space with weight function exp(-r) (_E<sub>1</sub><sup>r</sup>_)
@@ -206,12 +206,12 @@ val = scheme.integrate(lambda x: x ** 2)
 
 Example:
 ```python
-import numpy
+import numpy as np
 import quadpy
 
 scheme = quadpy.u2.get_good_scheme(7)
 scheme.show()
-val = scheme.integrate(lambda x: numpy.exp(x[0]), [0.0, 0.0], 1.0)
+val = scheme.integrate(lambda x: np.exp(x[0]), [0.0, 0.0], 1.0)
 ```
 
 ### Triangle (_T<sub>2</sub>_)
@@ -254,12 +254,12 @@ Apart from the classical centroid, vertex, and seven-point schemes we have
 
 Example:
 ```python
-import numpy
+import numpy as np
 import quadpy
 
 scheme = quadpy.t2.get_good_scheme(12)
 scheme.show()
-val = scheme.integrate(lambda x: numpy.exp(x[0]), [[0.0, 0.0], [1.0, 0.0], [0.5, 0.7]])
+val = scheme.integrate(lambda x: np.exp(x[0]), [[0.0, 0.0], [1.0, 0.0], [0.5, 0.7]])
 ```
 
 ### Disk (_S<sub>2</sub>_)
@@ -285,12 +285,12 @@ val = scheme.integrate(lambda x: numpy.exp(x[0]), [[0.0, 0.0], [1.0, 0.0], [0.5,
 
 Example:
 ```python
-import numpy
+import numpy as np
 import quadpy
 
 scheme = quadpy.s2.get_good_scheme(6)
 scheme.show()
-val = scheme.integrate(lambda x: numpy.exp(x[0]), [0.0, 0.0], 1.0)
+val = scheme.integrate(lambda x: np.exp(x[0]), [0.0, 0.0], 1.0)
 ```
 
 ### Quadrilateral (_C<sub>2</sub>_)
@@ -324,12 +324,12 @@ val = scheme.integrate(lambda x: numpy.exp(x[0]), [0.0, 0.0], 1.0)
 
 Example:
 ```python
-import numpy
+import numpy as np
 import quadpy
 
 scheme = quadpy.c2.get_good_scheme(7)
 val = scheme.integrate(
-    lambda x: numpy.exp(x[0]),
+    lambda x: np.exp(x[0]),
     [[[0.0, 0.0], [1.0, 0.0]], [[0.0, 1.0], [1.0, 1.0]]],
 )
 ```
@@ -398,23 +398,23 @@ val = scheme.integrate(lambda x: x[0] ** 2)
 
 Example:
 ```python
-import numpy
+import numpy as np
 import quadpy
 
 scheme = quadpy.u3.get_good_scheme(19)
 # scheme.show()
-val = scheme.integrate(lambda x: numpy.exp(x[0]), [0.0, 0.0, 0.0], 1.0)
+val = scheme.integrate(lambda x: np.exp(x[0]), [0.0, 0.0, 0.0], 1.0)
 ```
 Integration on the sphere can also be done for functions defined in spherical
 coordinates:
 ```python
-import numpy
+import numpy as np
 import quadpy
 
 
 def f(theta_phi):
     theta, phi = theta_phi
-    return numpy.sin(phi) ** 2 * numpy.sin(theta)
+    return np.sin(phi) ** 2 * np.sin(theta)
 
 
 scheme = quadpy.u3.get_good_scheme(19)
@@ -432,12 +432,12 @@ val = scheme.integrate_spherical(f)
 
 Example:
 ```python
-import numpy
+import numpy as np
 import quadpy
 
 scheme = quadpy.s3.get_good_scheme(4)
 # scheme.show()
-val = scheme.integrate(lambda x: numpy.exp(x[0]), [0.0, 0.0, 0.0], 1.0)
+val = scheme.integrate(lambda x: np.exp(x[0]), [0.0, 0.0, 0.0], 1.0)
 ```
 
 
@@ -466,13 +466,13 @@ val = scheme.integrate(lambda x: numpy.exp(x[0]), [0.0, 0.0, 0.0], 1.0)
 
 Example:
 ```python
-import numpy
+import numpy as np
 import quadpy
 
 scheme = quadpy.t3.get_good_scheme(5)
 # scheme.show()
 val = scheme.integrate(
-    lambda x: numpy.exp(x[0]),
+    lambda x: np.exp(x[0]),
     [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.5, 0.7, 0.0], [0.3, 0.9, 1.0]],
 )
 ```
@@ -494,13 +494,13 @@ val = scheme.integrate(
 
 Example:
 ```python
-import numpy
+import numpy as np
 import quadpy
 
 scheme = quadpy.c3.product(quadpy.c1.newton_cotes_closed(3))
 # scheme.show()
 val = scheme.integrate(
-    lambda x: numpy.exp(x[0]),
+    lambda x: np.exp(x[0]),
     quadpy.c3.cube_points([0.0, 1.0], [-0.3, 0.4], [1.0, 2.1]),
 )
 ```
@@ -512,13 +512,13 @@ val = scheme.integrate(
 
 Example:
 ```python
-import numpy
+import numpy as np
 import quadpy
 
 scheme = quadpy.p3.felippa_5()
 
 val = scheme.integrate(
-    lambda x: numpy.exp(x[0]),
+    lambda x: np.exp(x[0]),
     [
         [0.0, 0.0, 0.0],
         [1.0, 0.0, 0.0],
@@ -538,12 +538,12 @@ val = scheme.integrate(
 
 Example:
 ```python
-import numpy
+import numpy as np
 import quadpy
 
 scheme = quadpy.w3.felippa_3()
 val = scheme.integrate(
-    lambda x: numpy.exp(x[0]),
+    lambda x: np.exp(x[0]),
     [
         [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.5, 0.7, 0.0]],
         [[0.0, 0.0, 1.0], [1.0, 0.0, 1.0], [0.5, 0.7, 1.0]],
@@ -598,14 +598,14 @@ val = scheme.integrate(lambda x: x[0] ** 2)
 
 Example:
 ```python
-import numpy
+import numpy as np
 import quadpy
 
 dim = 4
 scheme = quadpy.tn.grundmann_moeller(dim, 3)
 val = scheme.integrate(
-    lambda x: numpy.exp(x[0]),
-    numpy.array(
+    lambda x: np.exp(x[0]),
+    np.array(
         [
             [0.0, 0.0, 0.0, 0.0],
             [1.0, 2.0, 0.0, 0.0],
@@ -627,12 +627,12 @@ val = scheme.integrate(
 
 Example:
 ```python
-import numpy
+import numpy as np
 import quadpy
 
 dim = 4
 scheme = quadpy.un.dobrodeev_1978(dim)
-val = scheme.integrate(lambda x: numpy.exp(x[0]), numpy.zeros(dim), 1.0)
+val = scheme.integrate(lambda x: np.exp(x[0]), np.zeros(dim), 1.0)
 ```
 
 
@@ -650,12 +650,12 @@ val = scheme.integrate(lambda x: numpy.exp(x[0]), numpy.zeros(dim), 1.0)
 
 Example:
 ```python
-import numpy
+import numpy as np
 import quadpy
 
 dim = 4
 scheme = quadpy.sn.dobrodeev_1970(dim)
-val = scheme.integrate(lambda x: numpy.exp(x[0]), numpy.zeros(dim), 1.0)
+val = scheme.integrate(lambda x: np.exp(x[0]), np.zeros(dim), 1.0)
 ```
 
 ### n-Cube (_C<sub>n</sub>_)
@@ -675,13 +675,13 @@ val = scheme.integrate(lambda x: numpy.exp(x[0]), numpy.zeros(dim), 1.0)
 
 Example:
 ```python
-import numpy
+import numpy as np
 import quadpy
 
 dim = 4
 scheme = quadpy.cn.stroud_cn_3_3(dim)
 val = scheme.integrate(
-    lambda x: numpy.exp(x[0]),
+    lambda x: np.exp(x[0]),
     quadpy.cn.ncube_points([0.0, 1.0], [0.1, 0.9], [-1.0, 1.0], [-1.0, -0.5]),
 )
 ```
diff --git a/quadpy/cn/_helpers.py b/quadpy/cn/_helpers.py
index bf6e4817b..6c2ca7a65 100644
--- a/quadpy/cn/_helpers.py
+++ b/quadpy/cn/_helpers.py
@@ -26,6 +26,16 @@ def __init__(
 
     def integrate(self, f: Callable, ncube, dot=np.dot):
         ncube = np.asarray(ncube)
+        if (
+            ncube.shape[: self.dim] != tuple(self.dim * [2])
+            or ncube.shape[-1] < self.dim
+        ):
+            expected = ", ".join(str(i) for i in self.dim * [2])
+            string = ", ".join(str(val) for val in ncube.shape)
+            raise QuadpyError(
+                f"Wrong domain shape. Expected ({expected}, ..., n >= {self.dim}), "
+                f"got ({string})."
+            )
         x = transform(self.points, ncube).T
         detJ = get_detJ(self.points, ncube)
 
@@ -36,14 +46,14 @@ def integrate(self, f: Callable, ncube, dot=np.dot):
                 f"Wrong return value shape {fx.shape}. " f"Expected (..., {string})."
             )
 
-        ref_vol = 2 ** np.prod(len(ncube.shape) - 1)
+        ref_vol = 2 ** np.prod(self.dim)
         return ref_vol * dot(fx * abs(detJ), self.weights)
 
-    def points_inside(self) -> bool:
+    def points_inside(self):
         """Are all points strictly inside the domain?"""
         return np.all((-1 < self.points) & (self.points < 1))
 
-    def points_inside_or_boundary(self) -> bool:
+    def points_inside_or_boundary(self):
         """Are all points inside the domain or on its boundary?"""
         return np.all((-1 <= self.points) & (self.points <= 1))
 
diff --git a/quadpy/enr/_helpers.py b/quadpy/enr/_helpers.py
index a79b09cc2..d0dc12644 100644
--- a/quadpy/enr/_helpers.py
+++ b/quadpy/enr/_helpers.py
@@ -1,6 +1,9 @@
+from typing import Callable
+
 import ndim
 import numpy as np
 
+from .._exception import QuadpyError
 from ..helpers import QuadratureScheme
 
 
@@ -11,7 +14,17 @@ def __init__(self, name, dim, weights, points, degree, source, tol=1.0e-14):
         assert points.shape[0] == dim
         super().__init__(name, weights, points, degree, source, tol)
 
-    def integrate(self, f, dot=np.dot):
+    def integrate(self, f: Callable, dot=np.dot):
         flt = np.vectorize(float)
         ref_vol = ndim.enr.volume(self.dim)
-        return ref_vol * dot(f(flt(self.points)), flt(self.weights))
+
+        x = flt(self.points)
+        fx = np.asarray(f(x))
+
+        assert len(x.shape) == 2
+        if len(x.shape) == 2 and fx.shape[-1] != x.shape[1]:
+            raise QuadpyError(
+                f"Wrong return value shape {fx.shape}. Expected (:, {x.shape[1]})."
+            )
+
+        return ref_vol * dot(fx, flt(self.weights))
diff --git a/quadpy/enr2/_helpers.py b/quadpy/enr2/_helpers.py
index 0457fd2a3..65259495b 100644
--- a/quadpy/enr2/_helpers.py
+++ b/quadpy/enr2/_helpers.py
@@ -1,6 +1,7 @@
 import ndim
 import numpy as np
 
+from .._exception import QuadpyError
 from ..helpers import QuadratureScheme
 
 
@@ -14,4 +15,13 @@ def __init__(self, name, dim, weights, points, degree, source, tol=2.7e-14):
     def integrate(self, f, dot=np.dot):
         flt = np.vectorize(float)
         ref_vol = ndim.enr2.volume(self.dim, "physicists")
-        return ref_vol * dot(f(flt(self.points)), flt(self.weights))
+
+        x = flt(self.points)
+        fx = np.asarray(f(x))
+
+        assert len(x.shape) == 2
+        if len(x.shape) == 2 and fx.shape[-1] != x.shape[1]:
+            raise QuadpyError(
+                f"Wrong return value shape {fx.shape}. Expected (:, {x.shape[1]})."
+            )
+        return ref_vol * dot(fx, flt(self.weights))
diff --git a/quadpy/tn/_helpers.py b/quadpy/tn/_helpers.py
index 4b7b9a5d5..b3b6d5e51 100644
--- a/quadpy/tn/_helpers.py
+++ b/quadpy/tn/_helpers.py
@@ -1,4 +1,5 @@
 import math
+from typing import Callable
 
 import numpy as np
 import sympy
@@ -21,14 +22,16 @@ def points_inside(self):
     def points_inside_or_boundary(self):
         return np.all((0 <= self.points) & (self.points <= 1))
 
-    def integrate(self, f, simplex, dot=np.dot):
+    def integrate(self, f: Callable, simplex, dot=np.dot):
         flt = np.vectorize(float)
         simplex = np.asarray(simplex)
-        if simplex.shape[0] != self.dim + 1:
+        if simplex.shape[0] != self.dim + 1 or simplex.shape[-1] < self.dim:
             string = ", ".join(str(val) for val in simplex.shape)
             raise QuadpyError(
-                f"Wrong domain shape. Expected ({self.dim + 1}, ...), got ({string})."
+                f"Wrong domain shape. Expected ({self.dim + 1}, ..., n >= {self.dim}), "
+                f"got ({string})."
             )
+
         x = transform(flt(self.points), simplex.T)
         vol = get_vol(simplex)
 
diff --git a/setup.cfg b/setup.cfg
index 1153a0d1b..0893c592b 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,6 +1,6 @@
 [metadata]
 name = quadpy
-version = 0.16.5
+version = 0.16.6
 author = Nico Schlömer
 author_email = nico.schloemer@gmail.com
 description = Numerical integration, quadrature for various domains


Test patch:
diff --git a/test/test_cn.py b/test/test_cn.py
index 41ec61c91..46219dbd0 100644
--- a/test/test_cn.py
+++ b/test/test_cn.py
@@ -67,7 +67,24 @@ def test_scheme(scheme):
     )
 
 
+# https://github.com/nschloe/quadpy/issues/401
+def test_multiple_cubes_volume():
+    val = quadpy.c2.schemes["stroud_c2_5_4"]().integrate(
+        lambda x: x[0],
+        np.stack(
+            [
+                quadpy.c2.rectangle_points([0.0, 1.0], [0.0, 1.0]),
+                quadpy.c2.rectangle_points([0.0, 1.0], [0.0, 1.0]),
+                quadpy.c2.rectangle_points([0.0, 1.0], [0.0, 1.0]),
+                quadpy.c2.rectangle_points([0.0, 1.0], [0.0, 1.0]),
+                quadpy.c2.rectangle_points([0.0, 1.0], [0.0, 1.0]),
+            ],
+            axis=-2,
+        ),
+    )
+    assert val.shape == (5,)
+    assert np.all(np.abs(val - 0.5) < 1.0e-13)
+
+
 if __name__ == "__main__":
-    n_ = 4
-    scheme_ = quadpy.cn.Stroud(n_, "Cn 7-1")
-    test_scheme(scheme_, 1.0e-14)
+    test_multiple_cubes_volume()