8-fold integral by using quadpy

[herr-mensch]

Hi Nico,

I accidentally found your project under the name quadpy, as I have been trying to
evaluate an eight-fold integral in sagemath using scipy. However, even for a simple
test integrand, i.e. 1, the computation is very-very slow, while Mathematica does it
over a tiny fraction of a second.

I decided to post this issue, because I don’t understand how to use quadpy
for the 8-fold integral with a scalar integrand, such that each integration variable
runs over a fixed interval, I mean I have smth like this in the Mathematica-language:

NIntegrate[ f[x1,x2,x3,x4,y1,y2,y3,y4] , {x1,-1,1}, {x2,-1,1}, {x3,-1,1}, {x4,-1,1}, {y1,-1,1}, {y2,-1,1}, {y3,-1,1}, {y4,-1,1}]
where f is a complex function.

Could you please explain how this integration can be realized by using quadpy?

[nschloe]

Alright, so there are basically two approaches:

• cn (n-cube).
• nquad (not implemented yet, see Quadpy should have dblquad. #254)

You can check out any scheme of quadpy.cn and try your luck. You'll get an approximation, mind you. Scipy has nquad which does an (arguably crude) adaptive integration. quadpy still needs to plug that hole.

[herr-mensch]

Thank you for your rapid response!

I should have said that I am very far from being an Informatiker. I mean I realised that I need to use smth like quadpy.cn in order to evaluate that integral, but the problem I experience is that I do not really understand how my mathematica example can be translated into the quadpy language... I would very much appreciate if you could provide that translation as I have little knowledge about this issue.

[nschloe]

Well, then maybe you need more explanations. What does the function look like in mathematica? Do you need the exact result of the integral or a floating-point approximation?

[herr-mensch]

Let me then be more concrete.

I already played with nquad from scipy by the evaluation of a 6-fold integral with the integrand 1 by running this
nquad(lambda x1,x2,x3,y1,y2,y3: 1, [[-1,+1], [-1,+1], [-1,+1], [-1,1], [-1,+1], [-1,+1]]). It takes of the order of 1 minute
which is very long for this ja trivial integral. What I would like to know how I can e.g. compute this same integral
but by using quadpy... but I do not understand how to implement that within quadpy.

It is quite okay to me if I have a floating-point approximation.

[nschloe]

What does your function look like? Most importantly: How smooth is it?

[herr-mensch]

The problem is not in that function, rather I don't understand how to implement
nquad(lambda x1,x2,x3,y1,y2,y3: 1, [[-1,+1], [-1,+1], [-1,+1], [-1,1], [-1,+1], [-1,+1]])
using your package quadpy.

[nschloe]

I'm just trying to find out if quadpy will do any good for you. We can work it out and you'll get some value, but it might not be accurate. It depends on the smoothness of your function.

[herr-mensch]

I first want to test with this simple example, because what I get up to now with either nested quad or nquad is very unsatisfactory even with that trivial integrand. You know I have a job running in a computer cluster already for ca. 5 days which must return the evaluation of 1 over a 8D cube with the length 1. I can of course compute this without any program. You see that smth is, therefore, deeply wrong with nested quad and nquad in this case.

If quadpy can compute this 8D integral with the integrand 1 during a small fraction of a second, then I will insert my integrand, i.e. f(x1,x2,x3,x4,y1,y2,y3,y4), which I actually need to evaluate. This complex function is quite good, but bulky.
In particular, it has two localisations, respectively, in x-space and y-space, and exponentially vanishes far away from where it has those two smooth enough peaks.

[nschloe]

So what about the main example in the readme? https://github.com/nschloe/quadpy/#n-cube-cn

[herr-mensch]

That's the thing I don't understand from this example how to implement a 8D integral.

For instance, I have x and y variables, but there is only "x[0]" in that example. How can this x[0] be related to my variables?
You see I just would like to have more details. These all might be trivial, but remember that I am not a computer scientist..

[nschloe]

If you're in 8D, you have x[0] to x[7].

[herr-mensch]

Ah, I see now! I was thinking that x[0] is not one-dimensional in that example..
It means I have f = f(x[0],..,x[7]), where x[i] \in (-1,1) for i = {0,..,7}.

You have there also this "quadpy.cn.ncube_points([0.0, 1.0], [0.1, 0.9], [-1.0, 1.0], [-1.0, -0.5])". How should it be in my case?

[nschloe]

I think

quadpy.cn.ncube_points([[-1,+1], [-1,+1], [-1,+1], [-1,1], [-1,+1], [-1,+1]])

[herr-mensch]

Thanks a lot! I will try this when I have quadpy installed. [The quadpy package should have been installed over this night on the computer cluster, but don't have it yet..]

[nschloe]

You can always

pip install quadpy

on any machine and play around with it.

[herr-mensch]

I have played by now several hours with quadpy which seems to be quite nice. However, I now have this issue

sage: import numpy
sage: import quadpy
sage: def f(x): return numpy.sqrt(x[0])
sage: scheme = quadpy.cn.dobrodeev_1978(2)
sage: scheme.integrate(lambda x: f(x), quadpy.cn.ncube_points([-10.0, 10.0], [-10.0, 10.0]))
/export/pc/sage/sage-9.0/src/bin/sage-ipython:1: RuntimeWarning: invalid value encountered in sqrt
#!/usr/bin/env sage-python
nan

Can you please tell why I have this problem?

Apparently, one needs to help numpy here by defining a new sqrt function [for real x]:
def new_sqrt(x): return numpy.exp((numpy.log(numpy.abs(x))+piI(numpy.abs(x)-(x))/(2*numpy.abs(x)))/2).
But is there anything better than this, which can work with quadpy?

[nschloe]

You're trying the take the sqrt of a negative number, hence the error.

[herr-mensch]

That was ja clear. The problem is apparently overcome if one instead uses numpy.sqrt(x[0]+I*0).

I have now another issue:
-- If mpmath is used, then the evaluation of my f(x[0],...x[7]) at a given point is very fast (\ll 1 sec).
-- if numpy is used, then the evaluation of my f(x) at the same point is very long (\sim 5 sec).
Clearly, numpy is not preferred by me in this case.

I get an error message if I use mpmath in the definition of my f(x) which I then try to integrate over a 8D cube using quadpy.
Can this problem be overcome anyhow within quadpy?

[nschloe]

Seems to be a problem in your function, not quadpy.

[nschloe]

Also, it's impossible for me to say anything if all the information I have is "it doens't work". You need to provide a minimal example that displays the issue (MWE).

[herr-mensch]

Here is the example:

sage: import mpmath
sage: import quadpy
sage: def f(x): return mpmath.exp(x[0])
sage: scheme = quadpy.cn.dobrodeev_1978(2)
sage: scheme.integrate(lambda x: f(x), quadpy.cn.ncube_points([-1.0, +1.0], [-1.0, +1.0]))

NotImplementedError                       Traceback (most recent call last)
<ipython-input-5-7f9c2dd5d151> in <module>()
----> 1 scheme.integrate(lambda x: f(x), quadpy.cn.ncube_points([-RealNumber('1.0'), +RealNumber('1.0')], [-RealNumber('1.0'), +RealNumber('1.0')]))

/export/pc/sage/sage-9.0/local/lib/python3.7/site-packages/quadpy/cn/_helpers.py in integrate(self, f, ncube, dot)
     17         detJ = get_detJ(self.points.T, ncube)
     18         ref_vol = 2 ** numpy.prod(len(ncube.shape) - 1)
---> 19         return ref_vol * dot(f(x) * abs(detJ), self.weights)
     20 
     21     def points_inside(self):

<ipython-input-5-7f9c2dd5d151> in <lambda>(x)
----> 1 scheme.integrate(lambda x: f(x), quadpy.cn.ncube_points([-RealNumber('1.0'), +RealNumber('1.0')], [-RealNumber('1.0'), +RealNumber('1.0')]))

<ipython-input-3-0c3a1834f964> in f(x)
----> 1 def f(x): return mpmath.exp(x[Integer(0)])

/export/pc/sage/sage-9.0/local/lib/python3.7/site-packages/sage/libs/mpmath/ext_main.pyx in sage.libs.mpmath.ext_main.Context._sage_exp (build/cythonized/sage/libs/mpmath/ext_main.c:15923)()
   1260             return rc
   1261         else:
-> 1262             raise NotImplementedError("unknown argument")
   1263 
   1264     def _sage_cos(ctx, x, **kwargs):

NotImplementedError: unknown argument
sage:

How can this be resolved?

[herr-mensch]

I don't know why some lines are larger and bold..

[herr-mensch]

But if I use numpy, then it works:

sage: import numpy
sage: import quadpy
sage: def f(x): return numpy.exp(x[0])
sage: scheme = quadpy.cn.dobrodeev_1978(2)
sage: scheme.integrate(lambda x: f(x), quadpy.cn.ncube_points([-1.0, +1.0], [-1.0, +1.0]))
4.700780178854748
sage:

[nschloe]

Please use triple backticks when posting code on GitHub (https://github.github.com/gfm/#fenced-code-blocks). I've edited your posts accordingly.

[nschloe]

To debug, just print some values to see where the error occurs. (That's what I would do.)

[herr-mensch]

I am not sure that I understand what you mean under "...print some values to see...".
Do you mean the values of my function? If affirmative, then

sage: N(f([t10,x10,y10,z10,t20,x20,y20,z20]))
-0.00429899797392792 + 0.00130472796050442*I
sage: N(f([t10,x10,y10+0.1,z10,t20,x20,y20,z20]))
-0.000713019262896895 - 0.00437659603528060*I
sage: N(f([t10,x10,y10,z10,t20,x20,y20+0.1,z20]))
-0.00419528008788525 + 0.00127344625071768*I
sage: N(f([t10,x10,y10,z10,t20,x20,y20+2.0,z20]))
-1.27135393547545e-6 + 3.87040475911282e-7*I
sage: N(f([t10,x10,y10+2.0,z10,t20,x20,y20,z20]))
-5.29283575737129e-10 - 3.44362062552914e-10*I
sage:

However, the problem is not in this function, rather one has to use numpy in order to take integrals
with the help of quadpy. I wonder if it is possible to use mpmath instead of numpy within quadpy..

[nschloe]

If your plan is to use mpmath to get "more precision", then I must disappoint you. Unless your function is a polynomial of low degree, the error in the result does not come from too few digits but from the general approximation error. In principle, quadpy may work with mpmath, but I can't see why one should do that.

[herr-mensch]

it is because if I use mpmath the above values of f(x[0],..,x[7]) are computed really fast, but if I use numpy, then the evaluation lasts ca. 5 sec. This is even worth than Mathematica in which I was able to reduce the evaluation time to ca. 1 sec.

But how could I use mpmath within quadpy?

[nschloe]

I don't know.

[herr-mensch]

okay, thanks.

Do you know how to get rid of this error message:

sage: scheme.integrate(lambda x: norm_1i(x), quadpy.cn.ncube_points([-10.0, +10.0], [-10.0, +10.0], [-10.0, +10.0]))
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-3-1dbdbc8e5ae0> in <module>()
----> 1 scheme.integrate(lambda x: norm_1i(x), quadpy.cn.ncube_points([-RealNumber('10.0'), +RealNumber('10.0')], [-RealNumber('10.0'), +RealNumber('10.0')], [-RealNumber('10.0'), +RealNumber('10.0')]))

/export/pc/sage/sage-9.0/local/lib/python3.7/site-packages/quadpy/cn/_helpers.py in integrate(self, f, ncube, dot)
     15         ncube = numpy.asarray(ncube)
     16         x = transform(self.points.T, ncube).T
---> 17         detJ = get_detJ(self.points.T, ncube)
     18         ref_vol = 2 ** numpy.prod(len(ncube.shape) - 1)
     19         return ref_vol * dot(f(x) * abs(detJ), self.weights)

/export/pc/sage/sage-9.0/local/lib/python3.7/site-packages/quadpy/cn/_helpers.py in get_detJ(xi, cube)
     93     J = numpy.array(J)
     94     J = numpy.moveaxis(J, (0, 1), (-2, -1))
---> 95     out = numpy.linalg.det(J)
     96     return out
     97 

/export/pc/sage/sage-9.0/local/lib/python3.7/site-packages/numpy/linalg/linalg.py in det(a)
   2091     t, result_t = _commonType(a)
   2092     signature = 'D->D' if isComplexType(t) else 'd->d'
-> 2093     r = _umath_linalg.det(a, signature=signature)
   2094     r = r.astype(result_t, copy=False)
   2095     return r

TypeError: No loop matching the specified signature and casting
was found for ufunc det
sage: 

[nschloe]

This looks like a wrong definition of the function. I'd be happy to help you further, please contact me privately for rates ([email protected]).

[nschloe]

Closing for lack of feedback. Feel free to reopen at any time!