JuliaMath / HCubature.jl

pure-Julia multidimensional h-adaptive integration
Other
153 stars 25 forks source link

Stubborn function #13

Closed yakir12 closed 5 years ago

yakir12 commented 5 years ago

I have a function that accepts a 2-dimentional vector and returns a scalar. If I sample only 100 points across both of its two domains I get this: screenshot_2018-11-30_09 32 24 Admittedly, there are only two regions where the function is not zero and they are kind of small. But no matter how small I set the rtol to, or how large I set the maxevals to, this still results in zero:

julia> hcubature(getsignal, [0, 0], [1, 0.5], rtol=10^-100, maxevals=10^100)
(0.0, 0.0)

and very quickly. Any ideas on what I can do?

giordano commented 5 years ago

Any ideas on what I can do?

Do you need to use this package or other packages are fine?

yakir12 commented 5 years ago

anything that works will be fine!

giordano commented 5 years ago

The divonne function in Cuba.jl has a peakfinder option to increase sampling around the given points. This is really useful when you have a function that is non-zero in very small regions and you know the position of those regions.

Honestly I've never used this feature in the julia interface but did use in the past directly with the C library, in these cases it's incredibly useful. Probably I can help you to do this in Cuba.jl as well.

Edit: probably the xgiven option should be sufficient, which is definitely easier to use

yakir12 commented 5 years ago

Great! I'll give it a solid try after lunch :)

yakir12 commented 5 years ago

Hmmm, seg fault:

julia> result = divonne((x, f) -> f = getsignal(x), 2, 1, nextra = 100)

signal (11): Segmentation fault
in expression starting at no file:0
unknown function (ip: 0xffffffffffffffff)
Iterate at /home/yakir/.julia/packages/Cuba/6hxNB/deps/usr/lib/libcuba.so (unknown line)
Integrate at /home/yakir/.julia/packages/Cuba/6hxNB/deps/usr/lib/libcuba.so (unknown line)
llDivonne at /home/yakir/.julia/packages/Cuba/6hxNB/deps/usr/lib/libcuba.so (unknown line)
dointegrate! at /home/yakir/.julia/packages/Cuba/6hxNB/src/divonne.jl:52 [inlined]
dointegrate at /home/yakir/.julia/packages/Cuba/6hxNB/src/Cuba.jl:185 [inlined]
#divonne#2 at /home/yakir/.julia/packages/Cuba/6hxNB/src/divonne.jl:145
#divonne at ./none:0
unknown function (ip: 0x7f3051f69193)
jl_fptr_trampoline at /usr/bin/../lib/libjulia.so.1 (unknown line)
jl_apply_generic at /usr/bin/../lib/libjulia.so.1 (unknown line)
unknown function (ip: 0x7f3070e53ab7)
unknown function (ip: 0x7f3070e53860)
unknown function (ip: 0x7f3070e54828)
unknown function (ip: 0x7f3070e54f04)
unknown function (ip: 0xfffffffffffffffe)
unknown function (ip: 0x7f304a9c109f)
unknown function (ip: 0xb)
unknown function (ip: 0x7f3070e553cc)
unknown function (ip: 0x7f306f87206c)
jl_toplevel_eval_in at /usr/bin/../lib/libjulia.so.1 (unknown line)
unknown function (ip: 0x7f306691c161)
jl_apply_generic at /usr/bin/../lib/libjulia.so.1 (unknown line)
unknown function (ip: 0x7f3066a8d6ba)
unknown function (ip: 0x7f3066a8d92b)
jl_apply_generic at /usr/bin/../lib/libjulia.so.1 (unknown line)
unknown function (ip: 0x7f306f858d8d)
unknown function (ip: 0xffffffffffffffff)
Allocations: 45132449 (Pool: 45126241; Big: 6208); GC: 91
Segmentation fault (core dumped)
stevengj commented 5 years ago

I added an initdiv keyword option that should address this: it allows you to force the function to be sampled on a finer grid initially. e.g. try checking out the master branch and passing initdiv=10.

yakir12 commented 5 years ago

Doing it now. Exciting! Thank you. Alternatively, you could use your Sobol to improve from a simple grid to a more uniform sampling.

stevengj commented 5 years ago

A simple grid is perfectly uniform. The purpose of Sobol sequences is to sample uniformly while lessening the curse of dimensionality. However, in HCubature we inevitably have the curse of dimensionality anyway, and we can only handle box subdomains.

(For a smooth function in low dimensions, HCubature should be much more efficient than any Monte–Carlo method.)

yakir12 commented 5 years ago

OK, I can confirm that, at least for me, it works! Thanks a lot!!!