add support for userdata and the corresponding tests #39

kunyuan · 2022-08-27T15:47:35Z

Hi,
Thank you for creating this fantastic package! I use it heavily for calculating Feynman diagrams in condensed matter physics.

For my applications, I need to provide the integral with a lot of userdata (parameters, the structure of diagrams, etc.). Using global constants makes the code much less flexible. So I need to add support of userdata to your package.

The implementation is straightforward: I pass a tuple of (func!, userdata) to the C function, then unpack it in the generic_integrand function which is on the julia side.

To keep the package back-compatible. The new userdata APi will be invoked only when the user provides the userdata keyword explicitly. I also didn't touch the existing core functions to avoid any possible regression.

The usage is shown as below:

# the previous API are still there
julia> @time result = vegas((x, f) -> f[1]=x[1])
  0.079806 seconds (1.10 M allocations: 37.326 MiB, 73.86% compilation time)
Component:
 1: 0.49999995061461283 ± 4.9696306524036384e-5 (prob.: 2.5596303730113235e-5)
Integrand evaluations: 232000
Number of subregions:  0
Note: The desired accuracy was reached

# Support userdata keyword.
julia> @time result = vegas((x, f, d) -> f[1]=x[1]+d, userdata=1)
  0.062823 seconds (229.48 k allocations: 10.900 MiB, 98.04% compilation time)
Component:
 1: 1.4999579838501624 ± 8.38583942796998e-5 (prob.: 0.012663461802755724)
Integrand evaluations: 13500
Number of subregions:  0
Note: The desired accuracy was reached

# Userdata can be any julia object.
julia> @time result = vegas((x, f, d) -> f[1]=x[1]+d[:a], userdata=Dict(:a=>1))
  0.062314 seconds (230.40 k allocations: 10.942 MiB, 97.97% compilation time)
Component:
 1: 1.4999579838501624 ± 8.38583942796998e-5 (prob.: 0.012663461802755724)
Integrand evaluations: 13500
Number of subregions:  0
Note: The desired accuracy was reached

An interesting observation of the above examples is that the new API with the userdata keyword significantly reduces the memory allocation. The packing and unpacking the integrand the userdata probably helps the julia compiler to learn more type information. If this is true, it probably makes sense to improve the old API to reduce memory allocation. But this could be for a separate PR.

A couple of new tests are added for the userdata keyword.

If everything looks fine to you, I could update the documentation, too.

giordano · 2022-08-27T15:53:35Z

Hi, thanks a lot for your contribution! However I'm a bit nervous about adding a field tagged with Any type. Did you see https://giordano.github.io/Cuba.jl/stable/#Passing-data-to-the-integrand-function-1? That should let you pass "user data" to the integrator, but without having to touch the code of the package

kunyuan · 2022-08-27T16:26:27Z

Hi, I have just made the Any type concrete.

Thank you for sending me this link. I have some concerns about that solution. Julia's official performance tips have mentioned that a "captured" value could cause type instability, which costs performance issue. See the following link. https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured.

Moreover, with the userdata keyword, the user will have more choice how to pass parameters to the integrand.

giordano · 2022-08-27T19:58:51Z

Ok, let's do it. Can you please update the documentation, including an example? Thanks!

kunyuan · 2022-08-28T02:05:02Z

Ok, let's do it. Can you please update the documentation, including an example? Thanks!

Sure, I will work on it.

kunyuan · 2022-08-28T21:41:17Z

Ok, let's do it. Can you please update the documentation, including an example? Thanks!

Hi, I have just updated index.md to include an example for the userdata keyword. Please take a look and feel free to modify and improve it. English is not my native language.

In addition, I need some help from you on the following things:

The julia doctests for the stochastic algorithms (vegas, suave and divonne) failed in my computer. The code generates slightly different answer but within the error bar. I guess the output of these algorithm depends on the random seed, and it is probably not proper to doctest them. Right now, I simply suppress the doctests right now. Is this okay?
I added a compatibility note. I put some reference to the next Cuba.jl release as version X.X.X. You may want to modify it.

Thank you very much!

giordano · 2022-08-28T21:46:03Z

I guess the output of these algorithm depends on the random seed, and it is probably not proper to doctest them.

The random seed is fixed to 0 by default:

Cuba.jl/src/Cuba.jl

Line 37 in d6e4736

const SEED = 0

I'm not convinced that's the reason.

I'd prefer to keep the doctests, at least they've been consistent for some time in CI.

kunyuan · 2022-08-28T23:49:27Z

I guess the output of these algorithm depends on the random seed, and it is probably not proper to doctest them.

The random seed is fixed to 0 by default:

Cuba.jl/src/Cuba.jl

Line 37 in d6e4736

const SEED = 0

I'm not convinced that's the reason.
I'd prefer to keep the doctests, at least they've been consistent for some time in CI.

I investigated this problem a little bit. I find that both vegas and divonne give the exactly the same results. However, suave has a small descrepency (but within the errorbar). The below is what I got:

julia> suave((x, f) -> f[1] = cos(x[1]))
Component:
 1: 0.8413748866950277 ± 7.772872640827611e-5 (prob.: 1.0)
Integrand evaluations: 23000
Number of subregions:  23
Note: The desired accuracy was reached

while according to the existing doc (the first example):

julia> suave((x, f) -> f[1] = cos(x[1]))
Component:
 1: 0.8413748866950329 ± 7.772872640815592e-5 (prob.: 1.0)
Integrand evaluations: 23000
Number of subregions:  23
Note: The desired accuracy was reached

The results here are calculated with your code in master branch, so tiny discrepency has nothing to do with the new userdata update.
Could it be caused by the hardware difference? I'm using Intel i9 CPU.

kunyuan · 2022-08-29T00:13:41Z

I guess the output of these algorithm depends on the random seed, and it is probably not proper to doctest them.

The random seed is fixed to 0 by default:

Cuba.jl/src/Cuba.jl

Line 37 in d6e4736

const SEED = 0

I'm not convinced that's the reason.
I'd prefer to keep the doctests, at least they've been consistent for some time in CI.

Anyway, I have added back the doctests. Although I am still getting tiny discrepency on my local computer, the CI doctests work just fine.

giordano

Before going on, in light of #40, would you agree to contribute to this package under the terms of the MIT license, instead of the current LGPL?

docs/src/index.md

kunyuan · 2022-09-15T14:16:45Z

Before going on, in light of #40, would you agree to contribute to this package under the terms of the MIT license, instead of the current LGPL?

Yes, I totally agree. Do you want to leave a comment below the #40 issue?

giordano · 2022-09-15T20:18:07Z

Yes, please! But before merging this PR I'd like to fix CI (see #41), and that's already proving to be complicated.

giordano · 2022-09-15T21:40:56Z

src/Cuba.jl

+function generic_integrand_userdata!(ndim::Cint, x_::Ptr{Cdouble}, ncomp::Cint,
+                            f_::Ptr{Cdouble}, func!_and_userdata, nvec::Cint)
+    func!, userdata = func!_and_userdata
+    # Get arrays from "x_" and "f_" pointers.
+    x = unsafe_wrap(Array, x_, (ndim, nvec))
+    f = unsafe_wrap(Array, f_, (ncomp, nvec))
+    func!(x, f, userdata)
+    return Cint(0)
+end


According to https://codecov.io/gh/giordano/Cuba.jl/compare/3ddf09bf2627745fab70fa90f0285350e5ac79d4...616a7a108c1c9bda05bbc0825aaf435c91a561af/diff this function is untested

In the new commit, I have added some tests for vectorized integration with userdata. It should be sufficient to cover the untested function.

giordano · 2022-09-22T16:02:09Z

Thanks a lot!

giordano · 2022-09-24T16:31:38Z

This is included in v2.3.0. Thanks again for your contribution!

add support for userdata and the corresponding tests

9c65279

make userdata type concrete

85663a5

kunyuan closed this Aug 27, 2022

kunyuan reopened this Aug 27, 2022

LittleBug added 2 commits August 27, 2022 23:38

working on index.html

cb564ea

add docs for userdata keyword in index.md

8b6b178

add back doctests

ffc567d

giordano reviewed Sep 15, 2022

View reviewed changes

docs/src/index.md Outdated Show resolved Hide resolved

LittleBug and others added 2 commits September 15, 2022 16:24

change version X.X.X to 2.3.0

1aa83b1

Merge branch 'master' into userdata

616a7a1

giordano reviewed Sep 15, 2022

View reviewed changes

LittleBug and others added 2 commits September 19, 2022 08:48

add tests for vectorization with userdata

fe398ac

Merge branch 'master' into userdata

caee8b5

giordano merged commit c58f22a into giordano:master Sep 22, 2022

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

add support for userdata and the corresponding tests #39

add support for userdata and the corresponding tests #39

kunyuan commented Aug 27, 2022

giordano commented Aug 27, 2022

kunyuan commented Aug 27, 2022 •

edited

Loading

giordano commented Aug 27, 2022

kunyuan commented Aug 28, 2022

kunyuan commented Aug 28, 2022

giordano commented Aug 28, 2022

kunyuan commented Aug 28, 2022 •

edited

Loading

kunyuan commented Aug 29, 2022

giordano left a comment

kunyuan commented Sep 15, 2022

giordano commented Sep 15, 2022

giordano Sep 15, 2022

kunyuan Sep 19, 2022

giordano commented Sep 22, 2022

giordano commented Sep 24, 2022

add support for userdata and the corresponding tests #39

add support for userdata and the corresponding tests #39

Conversation

kunyuan commented Aug 27, 2022

giordano commented Aug 27, 2022

kunyuan commented Aug 27, 2022 • edited Loading

giordano commented Aug 27, 2022

kunyuan commented Aug 28, 2022

kunyuan commented Aug 28, 2022

giordano commented Aug 28, 2022

kunyuan commented Aug 28, 2022 • edited Loading

kunyuan commented Aug 29, 2022

giordano left a comment

Choose a reason for hiding this comment

kunyuan commented Sep 15, 2022

giordano commented Sep 15, 2022

giordano Sep 15, 2022

Choose a reason for hiding this comment

kunyuan Sep 19, 2022

Choose a reason for hiding this comment

giordano commented Sep 22, 2022

giordano commented Sep 24, 2022

kunyuan commented Aug 27, 2022 •

edited

Loading

kunyuan commented Aug 28, 2022 •

edited

Loading