sigma-py / quadpy Goto Github PK
View Code? Open in Web Editor NEW:triangular_ruler: Numerical integration (quadrature, cubature) in Python
:triangular_ruler: Numerical integration (quadrature, cubature) in Python
It is well known that the map from the monomials integrals to integration formulas is ill-conditioned (see, e.g., Gautschi who pointed it out in great detail for Gaussian rules). Hence, sanity tests in quadpy can be improved by not checking against monomial integrals, but integrals of orthogonal polynomials.
print(quadpy.line_segment.integrate(lambda x: scipy.exp(1j*x), [-1, 100],quadpy.line_segment.GaussKronrod(3)))
The code is for calculating integral of exp(j*x) from -1 to 100, however, it returned a wrong value:
(-11.219936196799733+10.78170013544869j)
According to wolframalpha, the result should be:
(0.335105 - 0.322017j)
But the result of the code provided by Nico Schlömer in stackoverflow is correct:
val = quadpy.line_segment.integrate(
lambda x: scipy.exp(1j*x),
[0, 1],
quadpy.line_segment.GaussKronrod(3)
)
print(val)
result:
(0.841470984808+0.459697694132j)
Am I missing something?
To this end, transform the coordinates of the quadrature schemes to theta, phi.
The nonlinear optimization problem from Witherden-Vincent/Papanicodopulos can be used.
Hi,
I'm using quadpy 0.12.5
Maybe I'm missing something, but I think there's an issue with the way the Gauss-Hermite integration is coded.
When looking at the wikipedia page, I see there should be a Pi^(-0.5) term added and a rescaling of the gauss-hermit points.
import quadpy
import numpy as np
from scipy import stats
from scipy import integrate
def f(x):
return x**2
#Manual Gauss-Hermite
points, weights = np.polynomial.hermite.hermgauss(5)
mu, sigma = 0, 1
rescaled_points = np.sqrt(2)*sigma*points + mu
manual_gauss_hermite = np.pi**(-0.5)*np.dot(f(rescaled_points), weights)
#Manual trapezoid
xvals = np.linspace(-10, 10, num=100)
yvals = f(xvals) * stats.norm.pdf(xvals)
trapezoid_int = integrate.trapz(yvals, xvals)
#Gauss-Hermite integration
quadpy_integration = quadpy.e1r2.integrate(f, quadpy.e1r2.GaussHermite(5))
print("Manual Gauss-Hermite: ", manual_gauss_hermite)
print("Trapezoid: ", trapezoid_int)
print("Quadpy Integration: ", quadpy_integration)
gives
Manual Gauss-Hermite: 0.9999999999999999
Trapezoid: 1.0
Quadpy Integration: 0.8862269254527577
So the trapezoid rule agrees with the rewrite of gauss-hermite, but not with the gauss-hermite procedure contained in quadpy
.
A. Grundmann and H.M. Moeller, Invariant integration formulas for the n-simplex by combinatorial methods, SIAM J. Numer. Anal. 15 (1978), 282–290.
See http://www.pmf.ni.ac.rs/pmf/publikacije/filomat/2016/30-4/30-4-28.pdf for details
I have an implementation using mpmath
, but it feels like I'm being very inefficient about it.
I recently installed quadpy on windows 8.1 with Anaconda3 distribution of Python and running the spyder script.
When i try to run the first example of the readme:
import numpy
import quadpy
def f(x):
return numpy.sin(x[0]) * numpy.sin(x[1])
triangle = numpy.array([[0.0, 0.0], [1.0, 0.0], [0.7, 0.5]])
val = quadpy.triangle.integrate(f, triangle, quadpy.triangle.Strang(9))`
I found the following text:
_Traceback (most recent call last):
File "<ipython-input-4-c9a12b6d6ed4>", line 2, in <module>
import quadpy
File "C:\Anaconda3\lib\site-packages\quadpy\__init__.py", line 24, in <module>
from . import line_segment
File "C:\Anaconda3\lib\site-packages\quadpy\line_segment\__init__.py", line 7, in <module>
from .gauss_kronrod import GaussKronrod
File "C:\Anaconda3\lib\site-packages\quadpy\line_segment\gauss_kronrod.py", line 6, in <module>
from orthopy import jacobi_recursion_coefficients, gauss_from_coefficients
ImportError: cannot import name 'jacobi_recursion_coefficients'_
I would appeciate any help.
Check out How and how not to check Gaussian quadrature formulae by Gautschi:
I only checked this for quadrilateral as it's the only one I use, but probably needs adjustment in other examples as well:
val = quadpy.quadrilateral.integrate(
lambda x: numpy.exp(x[0]),
quadpy.quadrilateral.rectangle_points([0.0, 1.0], [-0.3, 0.6]),
quadpy.quadrilateral.Stroud(6)
)
in the readme does no longer work as the indices refer to the indices in Stoud's publication (such as 'C2 7-3'
). Furthermore an explanation of index ordering would be nice, as this changed from the old versions of this package.
Like in: "Rectangle points are now given in [x0,x1][y0,y1]
order" (If I assume correctly)
In old versions it was [[x0,y0][x1,y1][x2,y2][x3,y3]]
; this might be confusing. At least for me it was.
See http://scicomp.stackexchange.com/a/26426/3980 for details.
CONSTRUCTING FULLY SYMMETRIC CUBATURE FORMULAE FOR THE SPHERE
http://pages.uoregon.edu/yuan/paper/LongVersion.pdf
S. Heo, Y. Xu, Constructing symmetric cubature formulae on a triangle, in: Z. Chen, Y. Li, C. A. Micchelli, Y. Xu (Eds.), Advances in computational mathematics: proceedings of the Guangzhou international symposium, Marcel Dekker, 1999, pp. 203–221.
https://books.google.de/books?id=wFRGGmxsxEEC&pg=PA203&lpg=PA203&dq=%22Constructing+symmetric+cubature+formulae+on+a+triangle%22&source=bl&ots=0izOzG3a2B&sig=qlpWb9imUf31o7iScseDwGFrbPY&hl=en&sa=X&redir_esc=y#v=onepage&q=%22Constructing%20symmetric%20cubature%20formulae%20on%20a%20triangle%22&f=false
It's already used incorrectly in multiple spots. See numpy/numpy#9473.
SIAM J. SCI. COMPUT. c 2014
Vol. 36, No. 1, pp. 267-288
SPECTRA OF MULTIPLICATION OPERATORS AS A NUMERICAL TOOL
B. VIOREANU AND V. ROKHLIN
Right now, the monomials are computed this way:
n
.x
.Recurrence relations of orthogonal polynomials provide a better idea, for example in 2D: The table
(0, 0)
(1, 0) (0, 1)
(2, 0) (1, 1) (0, 2)
... ... ...
can be constructed from top to bottom by simply multiplying elements from one row by x
(or y
) and appending one more element to the left or the right. Although recurrent (and hence not parallelizable), this should be way more efficient.
Hi Nico,
I wanted to use quadpy for my integrations. I tried the line_integral example given in the readme. It gives me a ValueError: einstein sum subscripts string contains too many subscripts for operand 0
.
I am not sure what is going wrong. Any ideas?
Thanks.
To be found in accupy.
Chuanmiao Chen, Michal Křížek, Liping Liu,
Numerical Integration over Pyramids,
Adv. Appl. Math. Mech., Vol. 5, No. 3, pp. 309-320, 2013,
https://doi.org/10.4208/aamm.12-m12110.
Right now, the basis for the approximation is the approximated representation of j
(with lambertw
). One might however also also be able to take the fully accurate cosh-representation of j
as a basis. A Newton-iteration is done anyways.
Use it everywhere.
Right now, numpy.vectorize(float)
must be called for each integration since the scheme may return the points and weights symbolically. This may take a while. To avoid this, add the keyword symbolic=True/False
to all schemes and compute the points and weights accordingly.
Thanks for this software!
I cannot create custom quadratures using the scheme proposed in the documentation.
See the section: Generating your own Gauss quadrature in three simple steps
I am using quadpy
0.11.4 and orthopy
0.5.3.
Lebedev integration can be tested against spherical harmonics only and is then already guaranteed to be exact for all Cartesian polynomials of the same degree. Test against spherical harmonics when they exact integral over the sphere is known (cf. https://math.stackexchange.com/questions/2377595/integral-of-spherical-harmonics-over-sphere).
http://www.sciencedirect.com/science/article/pii/S0898122115001224
Since the most recent release of orthopy (0.2.0, two days ago), calls to orthopy.recurrence_coefficients.jacobi() no longer accept the 'mode' keyword argument (and potentially more has stopped working, but this is the first thing I ran into):
File "/usr/local/lib/python2.7/dist-packages/quadpy/line_segment/gauss_kronrod.py", line 47, in __init__ orthopy.recurrence_coefficients.jacobi(length, a, b, mode='numpy') TypeError: jacobi() got an unexpected keyword argument 'mode'
Suggested fix: specify orthopy 0.1.7 as requirement in setup.py for the time being, then update quadpy to work with the newest release of orthopy.
Thanks!
Hi!
Thanks for a great module - it's very impressive!
I have a question. Currently the integration relies on the integrand being defined as in your example:
def f(x):
return np.sin(x[0]) + np.cos(x[1])
where the numpy vectorisation takes care of input vectors of shape x.shape = (2, N)
.
I am currently dealing with functions that assume the input is of shape x.shape = (N, 2)
, as in the following example:
def f(x):
return np.sin(x[:, 0]) + np.cos(x[:, 1])
yielding errors like
def integrate(f, simplex, scheme, dot=numpy.dot):
flt = numpy.vectorize(float)
x = transform(flt(scheme.points).T, simplex.T)
vol = get_vol(simplex)
> return vol * dot(f(x), flt(scheme.weights))
E ValueError: shapes (2,) and (7,) not aligned: 2 (dim 0) != 7 (dim 0)
Would it be possible to support functions that accept ''transposed'' input?
Hello, I have a complex-valued function which argument is a 3D point. I want to do an integration of it over a line segment (parametrized in t
). I get a ValueError: operands could not be broadcast together with shapes (3,) (1,15)
. Here is the sample code:
import numpy
from scipy import linalg
import quadpy
startPoint = numpy.array([0.5, 1.0, 1.2])
endPoint = numpy.array([-0.5, -1.0, 0.7])
length = linalg.norm(startPoint - endPoint)
gamma = 0.01 + 0.586j
def point(t):
return startPoint + (endPoint - startPoint)*t
def integrand(t):
pt = point(t)
r1 = linalg.norm(pt - startPoint)
r2 = linalg.norm(pt - endPoint)
nf = (r1 + r2 + length)/(r1 + r2 - length)
return numpy.log(nf)*numpy.exp(-gamma)
value, error_estimate = quadpy.line_segment.integrate_adaptive(integrand,[0,1],1e-10)
The error is raised because t
is passed as an array. If I try to vectorize the point
function using numpy.vectorize
, then the error changes to ValueError: setting an array element with a sequence.
. I can't think of a work around. Could you help me?
Hi,
Although not quite my expertise area, I think the code cannot generate some basic Lebedev quadratures.
Most simple of these being, just a single node at the origin.
Also:
{
"a3": [
[5.7735026918962575e-01]
],
"degree": 2
}
and:
{
"a2": [
[7.07106781186547524e-01]
],
"degree": 4
}
In both cases degree is just an arbitrary number that I have written in order to make it work with your code.
I came across this issue when I was trying to implement the algorithm described in https://aip.scitation.org/doi/full/10.1063/1.5008630.
I would like to use this library to find a sudo numerical Laplace transform of some discrete data. So I know that I would need a weight function of
Whight=lambda sigma, omega: scipy.exp(sigma+1j*omega)
but how would I construct the kernel around f
which is my input data 1D vector but besides that I have no idea on how to set this up with quadpy using the Quadrilateral or the RabinowitzRichter integration scheme or if its even doable
so for exsample by f
data may be given by
f=np.linspace(0,5,150)
as a basic exsample
Quad:
exp(-r^2):
triangle:
disk:
Actually test ncube integration and advertise it.
Hi, thanks for the great package!
Are there any resources for those of us who are looking to play around with different rules without diving deep into the literature? I don't understand the distinctions between many of the formulae (except of course the salient parts like order), and I'm not sure how to learn the distinction without combing through the original papers. Thanks!
Consider the code:
print(quadpy.tetrahedron.XiaoGimbutas(2).weights.size)
in version 0.11.4 the code returns "4"
in version 0.12.0 the code returns "96".
I stumbled upon your project in Stack Overflow when looking for a way to calculate integrals over arrays. Sadly it is poorly documented, so I have troubles finding out how to do the integration.
I have a function f(x: complex) -> np.ndarray
which returns a ndarray
of shape (N_l, N_e)
. How can I integrate this function? I tried
quadpy.line_segment.integrate_adaptive(
f,
[0, 1],
1e-8
)
as well as
quadpy.line_segment.integrate_adaptive(
lambda x: f(x).reshape(-1),
[0, 1],
1e-8
)
in both cases it complains over the mismatch of dimensions. Is there a way to calculate the integral?
All quadrature scheme classes do nothing else than setting the data in __init__
. Would we lose anything when converting all classes to functions that return the respective scheme?
Dear Nico,
first of all, many thanks for your job. Really impressive!
In case You may need it, at
http://www.math.unipd.it/~alvise/sets.html
I have several sets, mostly for interpolation, on your regions. But also a file for cubature over the square that might be helpful for Your job:
http://www.math.unipd.it/~alvise/POINTSETS/set_amr_square.m
All the best,
Alvise Sommariva
Kahan provides less digits and needs more time. However, by default quadpy use Kahan only https://github.com/nschloe/quadpy/blob/6d49cc07e96d644bdcf039fd66a541208ce96794/quadpy/helpers/misc.py#L50
Can probably be ported from C code below.
https://www.mathworks.com/matlabcentral/fileexchange/26800-xsum
The points and weights for many schemes are given symbolically; let the user benefit from that.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.