Add Metropolis-Hastings algorithm for Potts modelling

Project.toml

@@ -3,6 +3,9 @@ uuid = "a9e1c8f2-2681-4423-bf67-1649e7f3cb09"
 authors = ["Dimitri Lozeve <dimitri@lozeve.com>"]
 version = "0.1.0"
 
+[deps]
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
 [compat]
 julia = "1"
 

src/GardenOptim.jl

@@ -1,5 +1,50 @@
 module GardenOptim
 
-greet() = print("Hello World!")
+using Logging
+
+export update!
+
+function swap!(grid::Array{Int, 2}, i::Int, j::Int)
+    t = grid[i]
+    grid[i] = grid[j]
+    grid[j] = t
+    grid
+end
+
+function neighbours(grid::Array{Int, 2}, idx)
+    m = size(grid, 1)
+    j, i = divrem(idx - 1, m)
+    i += 1
+    j += 1
+    neighbourindices = [(i, j-1), (i, j+1), (i-1, j), (i+1, j)]
+    [grid[k, l] for (k, l) in neighbourindices if 0 < k <= m && 0 < l <= m]
+end
+
+function deltacost(grid::Array{Int, 2}, costs::Array{Float64, 2}, i::Int, j::Int)
+    cost = 0
+    for k in neighbours(grid, i)
+        cost += costs[k, grid[j]] - costs[k, grid[i]]
+    end
+    for k in neighbours(grid, j)
+        cost += costs[k, grid[i]] - costs[k, grid[j]]
+    end
+    cost
+end
+
+function update!(grid::Array{Int, 2}, costs::Array{Float64, 2}, beta::Float64 = 10.0)
+    N = length(grid)
+    i, j = 0, 0
+    while i == j
+        i, j = rand(1:N, 2)
+    end
+    d = deltacost(grid, costs, i, j)
+    @debug "cost difference $d"
+    if rand() < exp(- beta * d)
+        @debug "swapping indices $i and $j"
+        return swap!(grid, i, j)
+    end
+    grid
+end
+
 
 end # module

test/runtests.jl

@@ -2,5 +2,29 @@ using GardenOptim
 using Test
 
 @testset "GardenOptim.jl" begin
-    # Write your own tests here.
+    @testset "swap" begin
+        grid = ones(Int, 5, 5)
+        grid[2] = 5
+        grid[8] = -3
+        GardenOptim.swap!(grid, 2, 8)
+        @test grid[2] == -3
+        @test grid[8] == 5
+        @test grid[6] == 1
+    end
+    @testset "neighbours" begin
+        grid = ones(Int, 5, 5)
+        @test length(GardenOptim.neighbours(grid, 4)) == 3
+        @test length(GardenOptim.neighbours(grid, 5)) == 2
+        @test length(GardenOptim.neighbours(grid, 8)) == 4
+        @test length(GardenOptim.neighbours(grid, 25)) == 2
+        @test GardenOptim.neighbours(grid, 1) == [1, 1]
+    end
+    @testset "deltacost" begin
+        grid = ones(Int, 5, 5)
+        costs = [[1. 0. 2.]; [0. 1. -1.]; [2. -1. 1.]]
+        @test GardenOptim.deltacost(grid, costs, 3, 8) == 0
+        grid[3] = 2
+        grid[9] = 3
+        @test GardenOptim.deltacost(grid, costs, 3, 8) != 0
+    end
 end