From 3a2dc40ce5b9a7b368be95f9aa63fce06d20c970 Mon Sep 17 00:00:00 2001
From: Dimitri Lozeve <dimitri@lozeve.com>
Date: Sat, 22 Feb 2020 20:08:03 +0100
Subject: [PATCH] Define the problem as a mathematical programming model

---
 Project.toml       |  4 +++
 src/GardenOptim.jl |  1 +
 src/optim.jl       | 68 ++++++++++++++++++++++++++++++++++++++++++++++
 test/optim.jl      | 15 ++++++++++
 test/runtests.jl   |  1 +
 5 files changed, 89 insertions(+)
 create mode 100644 src/optim.jl
 create mode 100644 test/optim.jl

diff --git a/Project.toml b/Project.toml
index ae3cc07..dc07392 100644
--- a/Project.toml
+++ b/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"
diff --git a/src/GardenOptim.jl b/src/GardenOptim.jl
index 800f3d8..187f36c 100644
--- a/src/GardenOptim.jl
+++ b/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})
diff --git a/src/optim.jl b/src/optim.jl
new file mode 100644
index 0000000..03cbbeb
--- /dev/null
+++ b/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
diff --git a/test/optim.jl b/test/optim.jl
new file mode 100644
index 0000000..1ac5d3d
--- /dev/null
+++ b/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
+
diff --git a/test/runtests.jl b/test/runtests.jl
index fbbbf55..1e4d2f6 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,4 +3,5 @@ using Test
 @testset "GardenOptim.jl" begin
     include("classification.jl")
     include("mcmc.jl")
+    include("optim.jl")
 end