How to call multi-dimensional vector value Vectorized Cubature numerical integration?

General Usage
,
1

As the topic, I want to use parallelized numerical integration. However, since the special receives of hcubature with vector-valued integrands, and I don’t know how to set each components of v respectively, there is a problem for calling higher dimensional case. Here are my attempts. First I use f(x,y)=x^2+y^2 as an 1D example:

using Cubature,BenchmarkTools
@btime hcubature(1,(x,v) -> v[1]= x[1]^2+x[2]^2, (0,0),(1,1))
  2.222 μs (73 allocations: 2.42 KiB)
([0.6666666666666665], [2.220446049250313e-16])
@btime hcubature(1,(x,v) -> v[1]= x[1]^2+x[2]^2, (0,0),(1,1))
  2.222 μs (73 allocations: 2.42 KiB)
([0.6666666666666665], [2.220446049250313e-16])
@btime hcubature_v((x,v) -> v[:]= x[1,:].^2+x[2,:].^2, (0,0),(1,1))
  1.020 μs (14 allocations: 1.38 KiB)
(0.6666666666666665, 2.220446049250313e-16)
@btime hcubature_v(1,(x,v) -> v[:]= x[1,:].^2+x[2,:].^2, (0,0),(1,1))
1.020 μs (14 allocations: 1.39 KiB)
([0.6666666666666665], [2.220446049250313e-16])

Both normal version and vectorized version works well.
Now let’s turn to a 2D example, f(x,y)=(x^2+y^2,x^3+y^3):

@btime hcubature(2,(x,v) -> v[:]=[x[1]^2+x[2]^2,x[1]^3+x[2]^3], (0,0),(1,1))
  2.756 μs (90 allocations: 3.78 KiB)
([0.6666666666666665, 0.4999999999999999], [2.220446049250313e-16, 1.1102230246251565e-16])
@btime hcubature_v(2,(x,v) -> v[1:2,:]=[x[1,:].^2+x[2,:].^2;x[1,:].^3+x[2,:].^3], (0,0),(1,1))
@btime hcubature_v(2,(x,v) -> v[:,:]=[x[1,:].^2+x[2,:].^2;x[1,:].^3+x[2,:].^3], (0,0),(1,1))

Unparallel integration is collect but vectorized versions return errors. I also tried let and to define a function as an integrand outside then use. But all of them failed. I think I should know how to set integrand (x,v) -> v[1,:],v[2,:]... respectively but not define as a list. Then I can find a way to use.

2

I find one method, by using begin and end.

using Cubature,BenchmarkTools
@btime hcubature(2,((x,v) -> begin v[1]=x[1]^2+x[2]^2;v[2]=x[1]^3+x[2]^3;v end), (0,0),(1,1))
2.244 μs (73 allocations: 2.45 KiB)
([0.6666666666666665, 0.4999999999999999], [2.220446049250313e-16, 1.1102230246251565e-16])
@btime hcubature_v(2,(x,v) -> begin v[1,:]=x[1,:].^2+x[2,:].^2; v[2,:]=x[1,:].^3+x[2,:].^3;v end, (0,0),(1,1))
  1.230 μs (19 allocations: 2.36 KiB)
([0.6666666666666665, 0.4999999999999999], [2.220446049250313e-16, 1.1102230246251565e-16])
1 Like
3

This won’t work — it just appends everything into a long 1d vector, because a slice like x[1,:] is a 1d array (“column vector”), not a 1-row matrix as in Matlab.

(Note also that for efficiency I would tend to try to work in-place in your integrands rather than generating lots of temporary arrays.)

Might be easier to write things out as explicit loops. See, for example, this test integrand which is used with hcubature_v.