Define the problem as a mathematical programming model

Project.toml

@@ -7,7 +7,11 @@ version = "0.1.0"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+Juniper = "2ddba703-00a4-53a7-87a5-e8b9971dde84"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"

src/GardenOptim.jl

@@ -17,6 +17,7 @@ export update!, fillgardenrandomly!, gardenevolution!, outputgarden
 include("classification.jl")
 include("loaddata.jl")
 include("mcmc.jl")
+include("optim.jl")
 
 "Save the garden to a CSV file."
 function outputgarden(garden::Matrix{Int}, plants::Vector{Symbol})

src/optim.jl

@@ -0,0 +1,68 @@
+using JuMP
+using Ipopt
+using Juniper
+
+function neighbourindices(mask::Matrix, idx::Int)::Vector{Int}
+    if mask[idx] == 0
+        return []
+    end
+
+    m, n = size(mask)
+    j, i = divrem(idx - 1, m)
+    i += 1
+    j += 1
+    indices = [(i, j-1), (i, j+1), (i-1, j), (i+1, j)]
+    # cells filled with 0 are not part of the garden
+    cartesianindices = [
+        CartesianIndex(k, l) for (k, l) in indices
+        if 0 < k <= m && 0 < l <= n && mask[k, l] != 0
+    ]
+    # convert to linear indices
+    LinearIndices(mask)[cartesianindices]
+end
+
+function neighbourmatrix(mask::Matrix)::Matrix{Bool}
+    N = length(mask)
+    d = zeros(N, N)
+    for i in 1:N
+        for j in neighbourindices(mask, i)
+            d[i, j] = 1
+        end
+    end
+    d
+end
+
+function definemodel(plantcounts::Vector, mask::Matrix, costs::Matrix)
+    N = length(mask)
+    Q = size(costs, 1)
+
+    optimizer = Juniper.Optimizer
+    nl_solver = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
+
+    model = Model(optimizer_with_attributes(optimizer, "nl_solver" => nl_solver))
+
+    @variable(model, x[1:N, 1:Q], Bin)
+
+    @NLobjective(
+        model,
+        Min,
+        sum(
+            costs[q, r] * x[i, q] * x[j, r]
+            for i = 1:N, q = 1:Q, r = 1:Q for j = neighbourindices(mask, i)
+        )
+    )
+
+    for i = 1:N
+        @constraint(model, sum(x[i, q] for q = 1:Q) == 1)
+        if mask[i] == 0
+            for q = 1:Q
+                @constraint(model, x[i, q] == 0)
+            end
+        end
+    end
+
+    for q = 1:Q
+        @constraint(model, sum(x[i, q] for i = 1:N) >= plantcounts[q])
+    end
+    model
+end

test/optim.jl

@@ -0,0 +1,15 @@
+using GardenOptim
+using LinearAlgebra
+using Test
+
+@testset "optim" begin
+    mask = zeros(Bool, 5, 7)
+    mask[[1, 2, 6, 7]] .= 1
+    @test GardenOptim.neighbourindices(mask, 1) == [6, 2]
+    @test GardenOptim.neighbourindices(mask, 2) == [7, 1]
+    @test GardenOptim.neighbourindices(mask, 3) == []
+    neighbourcount = sum([length(GardenOptim.neighbourindices(mask, i)) for i = 1:length(mask)])
+    @test sum(GardenOptim.neighbourmatrix(mask)) == neighbourcount
+    @test issymmetric(GardenOptim.neighbourmatrix(mask))
+end
+

test/runtests.jl

@@ -3,4 +3,5 @@ using Test
 @testset "GardenOptim.jl" begin
     include("classification.jl")
     include("mcmc.jl")
+    include("optim.jl")
 end