Merge pull request #29 from numericalEFT/mc-tools

Mc tools

src/grid/interpolate.jl

@@ -864,4 +864,50 @@ function differentiate1D(::CompositeDifferentiate, data, xgrid, x)
     return differentiate1D(view(data, head:tail), xgrid.subgrids[i], x)
 end
 
+# locate and volume for monte carlo
+
+"""
+    function locate(grid, x)
+
+return the index of grid point closest to x.
+Useful for Monte Carlo algorithm when variable x is continuous
+while histogram is stored on grid.
+
+#Arguments:
+- grid: one-dimensional grid of x
+- x: point to locate
+"""
+function locate(grid::AbstractGrid, x)
+    @assert x >= grid.bound[1] && x <= grid.bound[2]
+    i = floor(grid, x)
+    return abs(x-grid[i])<abs(x-grid[i+1]) ? i : (i+1)
+end
+
+"""
+    function volume(grid, i)
+
+return the volume of grid point i.
+The volume is defined as the length/area/volume/... of histogram bar
+represented by grid point i.
+In 1D grids of this package, it is defined as the length of interval
+between (grid[i-1]+grid[i])/2 and (grid[i]+grid[i+1])/2, and for edge points
+one side is replaced by boundary points.
+When index i is omitted, the length of the whole grid is returned.
+It is guaranteed that volume(grid)==sum(volume(grid, i) for i in 1:length(grid)).
+
+#Arguments:
+- grid: one-dimensional grid
+- i: index of grid point
+"""
+function volume(grid::AbstractGrid, i)
+    if i != 1 && i != length(grid)
+        return (grid[i+1]-grid[i-1])/2
+    elseif i == 1
+        return (grid[i+1]+grid[i])/2 - grid.bound[1]
+    else
+        return grid.bound[2]-(grid[i]+grid[i-1])/2
+    end
+end
+volume(grid::AbstractGrid) = grid.bound[2] - grid.bound[1]
+
 end

test/mc.jl

@@ -0,0 +1,82 @@
+@testset "MC" begin
+    # testing locate and volume functions provided for monte carlo
+    rng = MersenneTwister(1234)
+
+    function test_locate_volume(grid; δ = 1e-6)
+        vol = 0.0
+        for (ti, t) in enumerate(grid.grid)
+            # test locate
+            @test ti == Interp.locate(grid, t)
+            @test ti == Interp.locate(grid, (ti == 1) ? t : t - δ)
+            @test ti == Interp.locate(grid, (ti == length(grid)) ? t : t + δ)
+            vol += Interp.volume(grid, ti)
+        end
+
+        @test vol ≈ Interp.volume(grid)
+    end
+
+    @testset "Locate and Volume" begin
+        Ng, Ns = 10, 4
+        β = π
+
+        tgrid = SimpleGrid.Uniform{Float64}([0.0, β], Ng)
+        test_locate_volume(tgrid)
+
+        tgrid = SimpleGrid.BaryCheb{Float64}([0.0, β], Ng)
+        test_locate_volume(tgrid)
+
+        tgrid = SimpleGrid.Log{Float64}([0.0, β], Ng, 0.001, true)
+        test_locate_volume(tgrid)
+
+        tgrid = CompositeGrid.LogDensedGrid(:cheb, [0.0, β],[β/2,], Ng, 0.001, Ns)
+        test_locate_volume(tgrid)
+
+    end
+
+    function test_mc_histogram(grid; Npp=1e6, f = (x -> x))
+        # println("testing:$(typeof(grid))")
+        Ng = length(grid)
+        a, b = grid.bound[1], grid.bound[2]
+        Nmc = Npp * Ng
+
+        hist = zeros(Ng)
+        for i in 1:Nmc
+            x = rand(rng) * (b - a) + a
+            ind = Interp.locate(grid, x)
+            vol = Interp.volume(grid, ind)
+            hist[ind] += f(x) / vol * Interp.volume(grid)
+        end
+
+        for (ti, t) in enumerate(grid.grid)
+            vi = Interp.volume(grid, ti)
+            vall = Interp.volume(grid)
+            if ti != 1 && ti != Ng
+                @test isapprox(hist[ti] / Nmc, f(t),
+                               rtol=5 / sqrt(Nmc*vi/vall), atol = 5*vi/vall)
+            else
+                # edge points has extra error because grid point is not at center of interval
+                @test isapprox(hist[ti] / Nmc, f(t),
+                               rtol=5 / sqrt(Nmc*vi/vall), atol = 5*vi/vall)
+            end
+        end
+    end
+
+    @testset "MC histogram" begin
+        Ng, Nl, Ns = 16, 8, 4
+        β = π
+
+        tgrid = SimpleGrid.Uniform{Float64}([0.0, β], Ng)
+        test_mc_histogram(tgrid)
+
+        tgrid = SimpleGrid.BaryCheb{Float64}([0.0, β], Nl)
+        test_mc_histogram(tgrid)
+
+        tgrid = SimpleGrid.Log{Float64}([0.0, β], Nl, 0.1, true)
+        test_mc_histogram(tgrid)
+
+        tgrid = CompositeGrid.LogDensedGrid(:cheb, [0.0, β],[β/2,], Nl, 0.1, Ns)
+        test_mc_histogram(tgrid)
+
+    end
+
+end

test/runtests.jl

@@ -6,6 +6,7 @@ if isempty(ARGS)
     include("grid.jl")
     include("interpolate.jl")
     include("io.jl")
+    include("mc.jl")
 else
     include(ARGS[1])
 end