the integrals with Heaviside function are very inaccurate

[sashamigdal]

The test below must return exactly 0; it returns 0.01 instead.

QuadPy.py::test_GroupIntegral PASSED                                     [100%]tol= 1e-14

quadpy: SphericalRestrictedIntegral = (0.012345679012345685+0j)
quadpy O(3) Restricted Integral 25.70 ms

import quadpy
import numpy as np
from scipy.linalg import eigh
from Timer import MTimer as Timer
#from ScipyQuad import SphericalRestrictedIntegral as SciPySRI

def test_Four():
    dim = 4
    scheme = quadpy.sn.dobrodeev_1970(dim)
    val = scheme.integrate(lambda x: np.ones(73), np.zeros(dim), 1.0)
    print(val)

def SphericalRestrictedIntegral(R):
    tt, W = eigh(R.imag)
    rr = np.array(R.real, dtype=float)
    RR = W.T @ rr @ W
    dim = 4
    scheme = quadpy.sn.dobrodeev_1970(dim)
    print("tol=", scheme.test_tolerance)

    vol = scheme.integrate(lambda x: np.ones(x.shape[1]), np.zeros(dim), 1.0)

    def funcC(x):
        XX = x[np.newaxis, :, :] * x[:, np.newaxis, :]
        imtrace = tt.dot(x ** 2)
        return np.exp(1j * np.trace(RR.dot(XX)) - imtrace) / vol * np.heaviside(-(imtrace)**2, 0.5)

    return scheme.integrate(funcC, np.zeros(dim), 1.0)

def test_GroupIntegral():
    R = np.array([np.random.normal() + 1j * np.random.normal() for _ in range(16)]).reshape((4, 4))
    R += R.T
    with Timer("quadpy O(3) Restricted Integral"):
        res = SphericalRestrictedIntegral(R)
        print("\nquadpy: SphericalRestrictedIntegral =", res)
    # with Timer("scipy O(3) Restricted Integral"):
    #     [rres, ires] = SciPySRI(R)
    #     print("\nscipy SphericalRestrictedIntegral =", rres[0] + 1j * ires[0] , "+/-" , rres[1] + 1j* ires[1])
    pass

[sashamigdal]

Here is the simplest version of the code displaying this bug. The integral computes correctly here if you remove assert. However, this assertion checks that all 4-vectors x[I,:] are on the unit sphere, and it fails.

def test_Four():
    dim = 4
    scheme = quadpy.sn.dobrodeev_1970(dim)
    def f(x):
        assert (np.min(np.sum(x ** 2, axis=0)) > 0.99)
        return np.ones(x.shape[1])
    val = scheme.integrate(f, np.zeros(dim), 1.0)*2/np.pi**2
    print(val)

[sashamigdal]

The "bug" was a result of my confusion: I used the method for a ball instead to]=fone for a sphere. Afyter fixing, the method works. Here is the test. with the output

def test_Four():
    dim = 4
    # scheme = quadpy.un.dobrodeev_1978(dim)
    scheme = quadpy.un.mysovskikh_2(dim)
    def f(x):
        assert (np.min(np.sum(x ** 2, axis=0)) > 0.99)
        return np.heaviside(x[0], 0.5)
    val = scheme.integrate(f, np.zeros(dim), 1.0)/(2*np.pi**2)
    print(val)

QuadPy.py::test_Four PASSED                                              [100%]0.5

========================= 1 passed, 1 warning in 1.60s =========================

[sashamigdal]

Here is a nontrivial test of the integration with the Heaviside function. This is a volume of randomly chosen half of the sphere S_3 (U_4). Naive symmetric points on a sphere, chosen in advance, would not interpolate across the random plane we used to cut the sphere.

So, why does it work??

def test_Four():
    dim = 4
    # scheme = quadpy.un.dobrodeev_1978(dim)
    scheme = quadpy.un.mysovskikh_2(dim)
    # print("\n",scheme.degree, "\n",scheme.points)
    # print(scheme.source)
    s = np.random.normal(size=dim)
    def f(x):
        assert x.shape[0] ==4
        assert (np.min(np.sum(x ** 2, axis=0)) > 0.99)
        return np.heaviside(s.dot(x), 0.5)
    val = scheme.integrate(f, np.zeros(dim), 1.0)/(2*np.pi**2)
    print(val)

However, in a less trivial example

def SphericalRestrictedIntegral(R):
    dim = 4
    scheme = quadpy.un.mysovskikh_2(dim)
    vol = scheme.integrate(lambda x: np.ones(x.shape[1]), np.zeros(dim), 1.0)
    def funcC(x):
        trc = np.sum(x*(R.dot(x)),axis=0)
        return np.exp(1j * trc) / vol * np.heaviside(trc.imag, 0.5)

    return scheme.integrate(funcC, np.zeros(dim), 1.0)

it produces results different from those produced by Mathematica and confirmed independently with "scipy. tplquad").

[sashamigdal]

Here are the results for the same integral over the unit sphere in 4 dimensions, with part of the sphere excluded by a Heaviside function

R = 
{{(-1.8865156212567178) + (-3.9332306782313617) I,(1.9313163965055087) + (1.945569085060555) I,(0.05165837220957231) + (0.16211997626152372) I,(0.6518737608730687) + (-0.13886204776087818) I},{(1.9313163965055087) + (1.945569085060555) I,(-2.4150290602528197) + (-0.4104850431708015) I,(0.5749968497969923) + (1.4490564635457646) I,(1.4261791614786743) + (1.6606551109080527) I},{(0.05165837220957231) + (0.16211997626152372) I,(0.5749968497969923) + (1.4490564635457646) I,(-2.696285384006767) + (-1.3987035224168607) I,(-1.727761829409315) + (-3.2702277051903845) I},{(0.6518737608730687) + (-0.13886204776087818) I,(1.4261791614786743) + (1.6606551109080527) I,(-1.727761829409315) + (-3.2702277051903845) I,(-1.4424998499797432) + (-1.1514209260936545) I}}
#########
results:

quadpy: SphericalRestrictedIntegral = (0.08611712973096886-0.06457595633916889j)
quadpy O(3) Restricted Integral 28.01 ms
scipy SphericalRestrictedIntegral = (0.05315358505360589-0.07107615617190838j) +/- (0.000926423291640542+0.0009594241392830782j)
scipy O(3) Restricted Integral 0.75 s

Mathematica Spherical Restricted Integral = Complex[0.053205332320443215, -0.07095290700221205]
Mathematica O(3) Restricted Integral 16.89 s

We see that Mathematica and Scipy agree with the requested preco=icion (3 digits for Scipy, 10 for Mathematica).

The "quadpy" results are significantly different., probably, the interpolation does not work in this case of the curved domain excluded from the sphere.

[nschloe]

What quadpy version are you using?

[sashamigdal]

Hi:

I installed quadpy 0.17.7, which is the latest version available to me.

Speaking about the inaccuracy with the Heaviside function, this is what is to be expected with the interpolation from fixed symmetric nodes as described in the original paper. The Heaviside function has a discontinuity in the integration region, whereas the integration method assumes differentiable functions, expandable in polynomials.

I will use your library for other cases because it offers the best speed.

On Tue, Jul 11, 2023 at 9:35 AM Nico Schlömer @.***> wrote:

What quadpy version are you using?

— Reply to this email directly, view it on GitHub https://github.com/sigma-py/quadpy/issues/478#issuecomment-1630847151, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAXUX3ZYKULCMYT4QVMAYBTXPVJAXANCNFSM6AAAAAA2DVWL6M . You are receiving this because you modified the open/close state.Message ID: @.***>

-- Sasha Migdal (609) 933 3896 cell

[nschloe]

The Heaviside function has a discontinuity in the integration region, whereas the integration method assumes differentiable functions, expandable in polynomials.

Indeed! To expand, quadpy is about Gaussian integration. Every scheme is accurate up to a certain degree, and the assumption is that it works well for integrands which aren't too wild. Discontinuities though are as wild as it gets -- Gaussian integration is not the right tool for the task. Instead, it is often possible to identify a subdomain in which the integrand is nonzero, and integrate over that. In the case of a characteristic function, such as Heaviside, evaluating the integral is really just asking for the area of a subdomain. There are much better tools for answering this question.

[nschloe]

Is this still a relevant issue? If yes, please post the code you'd like to see improved (with R) and I'll take a look.

[sashamigdal]

This issue is closed. The Heaviside function cannot be integrated by this method, as we agree.