[WIP] Add quick plotting of adaptive grid

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"]

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

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