sigma-py / quadpy Goto Github PK

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

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?

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!

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.

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".

See http://scicomp.stackexchange.com/a/26426/3980 for details.

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?

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.

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

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.

Actually test ncube integration and advertise it.

https://arxiv.org/abs/1411.5631

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?

https://doi.org/10.1016/j.amc.2006.10.041

http://www.sciencedirect.com/science/article/pii/S0898122115001224

• triangle
• tet
• wedge
• pyramid
• quad
• hex

Quad:

• Waldron, Symmetries of linear functionals
• Cohen, Some integration formulae for symmetric functions of two variables
• H.J. Schmid, On cubature formulae with a minimal number of knots
• R. Franke, Obtaining cubatures for rectangles and other planar regions by using orthogonal polynomials
• Möller, Linear abhängige Punktfunktionale bei zweidimensionalen Interpolations- und Approximationsproblemen
• D.R. Hunkins, Cubatures of precision 2k and 2k+1 for hyperrectangles
• A. Haegemans and R. Piessens, Construction of cubature formulas of degree seven and nine symmetric planar regions, using orthogonal polynomials
• Franke, Orthogonal Polynomials and Approximate Multiple Integration
• R. Piessens and A. Haegemans, Cubature formulas of degree nine for symmetric planar regions
• P. Verlinden and R. Cools, The algebraic construction of a minimal cubature formula of degree 11 for the square

exp(-r^2):

• A. Haegemans, Tables of circularly symmetrical integration formulas of degree 2d-1 for two-dimensional circularly symmetrical regions

triangle:

• B. Griener and H.J. Schmid, An interactive tool to visualize common zeros of twodimensional polynomials

disk:

• Near-minimal cubature formulae on the disk, Brahim Benouahmane, Cuyt Annie, Irem Yaman

Right now, the monomials are computed this way:

1. Exponent pairs are generated for all n.
2. These vectors are concatenated.
3. The exponents are applied individually to an array 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.

http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf

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.

It's already used incorrectly in multiple spots. See numpy/numpy#9473.

Use it everywhere.

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

• Kahan: The local error is subtracted from the next element. 1 to 3 more valid digits, 2 to 9 times slower.
• Knuth: As if the sum is accumulated in a 128 bit float: about 15 more valid digits. 1.4 to 4 times slower. This is suitable for the most real world problems.

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.

• line segment
• triangle
• tetrahedron
• hexahedron
• quadrilateral
• pyramid
• wedge
• n-ball
• n-cube
• n-sphere
• sphere
• circle
• ball
• disk

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 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.

Check out How and how not to check Gaussian quadrature formulae by Gautschi:

• https://link.springer.com/article/10.1007/BF02218441,
• https://www.cs.purdue.edu/homes/wxg/selected_works/section_08/084.pdf.

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.

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.

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?

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).

To this end, transform the coordinates of the quadrature schemes to theta, phi.

The points and weights for many schemes are given symbolically; let the user benefit from that.

The nonlinear optimization problem from Witherden-Vincent/Papanicodopulos can be used.

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!

See sympy/sympy#13883 (comment).

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.

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?

• 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

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.

See https://stackoverflow.com/a/15858047/353337.

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.

• line segment
• triangle
• tetrahedron
• hexahedron
• quadrilateral
• pyramid
• wedge
• n-ball
• n-cube
• n-simplex
• n-sphere
• sphere
• circle
• ball
• disk
• e1r
• e2r
• e3r
• enr
• e1r2
• e2r2
• e3r2
• enr2

To be found in accupy.