Comments (9)
This issue is closed. The Heaviside function cannot be integrated by this method, as we agree.
from quadpy.
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)
from quadpy.
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 =========================
from quadpy.
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").
from quadpy.
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.
from quadpy.
What quadpy version are you using?
from quadpy.
from quadpy.
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.
from quadpy.
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.
from quadpy.
Related Issues (20)
- The source code is not available HOT 1
- Status of package and licensing HOT 4
- Quadpy is now asking for a sigma license HOT 18
- Can't import HOT 8
- Orthopy version no longer supported HOT 9
- version of orthopy no longer supported HOT 5
- Return points and values for integrate_adaptive HOT 1
- Importing quadpy causes python to quit without any error message HOT 4
- Issue upgrading quadpy HOT 7
- Unable to find valid license for Sigma. HOT 3
- Recurrent issue with sigma license HOT 5
- NameError: name 'c2' is not defined HOT 3
- Quadpy RuntimeError HOT 7
- Adaptive multivariate integration over hypercubes HOT 9
- Integrate over multiple limits at once using Quadpy HOT 8
- newest pip version not importable HOT 1
- quadpy.c1.gauss_patterson(5) problem HOT 2
- License Key Requirement for Quadpy Library HOT 1
- Quadrature results inconsistent with scipy.integrate.quad and mpmath.quad HOT 3
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from quadpy.