diff --git a/Project.toml b/Project.toml index e0156ba..06a2623 100644 --- a/Project.toml +++ b/Project.toml @@ -5,15 +5,16 @@ version = "0.0.3" [deps] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [compat] -julia = "1" Distributions = "0.23, 0.24, 0.25" +julia = "1" [extras] -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] test = ["HCubature", "Test"] diff --git a/src/vegas.jl b/src/vegas.jl index 54c2671..ab58965 100644 --- a/src/vegas.jl +++ b/src/vegas.jl @@ -1,5 +1,6 @@ using Random using Distributions +using Plots struct VEGASResult{T1,T2,T3,T4,T5,T6} grid::T1 @@ -30,11 +31,11 @@ Defaults to ones(2) Kwargs: ------ - nbins: Number of bins in each dimension. -Defaults to 100. +Defaults to 1000. - ncalls: Number of function calls per iteration. -Defaults to 1000. +Defaults to 10000. - maxiter: Maximum number of iterations. -Defaults to 100. +Defaults to 10. - rtol: Relative tolerance required. Defaults to 1e-4. - atol: Absolute tolerance required. @@ -64,9 +65,9 @@ Computational Physics 27.2 (1978): 192-203. function vegas(func, a = [0.,0.], b = [1.,1.]; - maxiter = 100, - nbins = 100, - ncalls = 1000, + maxiter = 10, + nbins = 1000, + ncalls = 10000, rtol = 1e-4, atol = 1e-4, debug = false, @@ -350,3 +351,12 @@ function extract_from_bins!(m, optm, grid, res) dist end +function plot_grid(res::VEGASResult) + g = res.cumulative_grid + ndim = size(g,2) + ndim != 2 && throw(ErrorException("Can only print 2D grids.")) + p = plot(xlim = (g[1,1], g[end,1]), ylim = (g[1,2], g[end,2])) + vline!(g[:,1], c = 1) + hline!(g[:,2], c = 1) + p +end diff --git a/test/runtests.jl b/test/runtests.jl index f711661..d38eb5d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -56,3 +56,15 @@ v = vegas(f, zeros(3), fill(3.0, 3), nbins = 1000, ncalls = 10000).integral_esti end end +function ackley(x) + a, b, c = 20.0, -0.2, 2.0*π + len_recip = inv(length(x)) + sum_sqrs = zero(eltype(x)) + sum_cos = sum_sqrs + for i in x + sum_cos += cos(c*i) + sum_sqrs += i^2 + end + return (-a * exp(b * sqrt(len_recip*sum_sqrs)) - + exp(len_recip*sum_cos) + a + 2.71) +end