From 3a2dc40ce5b9a7b368be95f9aa63fce06d20c970 Mon Sep 17 00:00:00 2001 From: Dimitri Lozeve Date: Sat, 22 Feb 2020 20:08:03 +0100 Subject: [PATCH 1/6] 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 From c5c13916792f02e22881c728144ce1956ee798c5 Mon Sep 17 00:00:00 2001 From: Dimitri Lozeve Date: Sun, 23 Feb 2020 14:39:12 +0100 Subject: [PATCH 2/6] Use CartesianIndices instead of tuples --- src/optim.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/optim.jl b/src/optim.jl index 03cbbeb..0f65f89 100644 --- a/src/optim.jl +++ b/src/optim.jl @@ -8,9 +8,7 @@ function neighbourindices(mask::Matrix, idx::Int)::Vector{Int} end m, n = size(mask) - j, i = divrem(idx - 1, m) - i += 1 - j += 1 + i, j = Tuple(CartesianIndices(mask)[idx]) indices = [(i, j-1), (i, j+1), (i-1, j), (i+1, j)] # cells filled with 0 are not part of the garden cartesianindices = [ From fc1d2eb7920c3940991069399c5fd2c3db49736c Mon Sep 17 00:00:00 2001 From: Dimitri Lozeve Date: Sun, 23 Feb 2020 22:18:10 +0100 Subject: [PATCH 3/6] Fix constraint to allow empty cells --- src/optim.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/optim.jl b/src/optim.jl index 0f65f89..340f3ed 100644 --- a/src/optim.jl +++ b/src/optim.jl @@ -51,7 +51,7 @@ function definemodel(plantcounts::Vector, mask::Matrix, costs::Matrix) ) for i = 1:N - @constraint(model, sum(x[i, q] for q = 1:Q) == 1) + @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) From 43b752649788fdfa39982ca9a5c53e4611780142 Mon Sep 17 00:00:00 2001 From: Dimitri Lozeve Date: Sun, 23 Feb 2020 22:20:10 +0100 Subject: [PATCH 4/6] Add constraints for fixed plants in garden --- src/optim.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/optim.jl b/src/optim.jl index 340f3ed..a5a4bef 100644 --- a/src/optim.jl +++ b/src/optim.jl @@ -30,7 +30,7 @@ function neighbourmatrix(mask::Matrix)::Matrix{Bool} d end -function definemodel(plantcounts::Vector, mask::Matrix, costs::Matrix) +function definemodel(plantcounts::Vector, garden::Matrix, mask::Matrix, costs::Matrix) N = length(mask) Q = size(costs, 1) @@ -57,6 +57,9 @@ function definemodel(plantcounts::Vector, mask::Matrix, costs::Matrix) @constraint(model, x[i, q] == 0) end end + if garden[i] != 0 + @constraint(model, x[i, garden[i]] == 1) + end end for q = 1:Q From fd5cffac9e0ae3b604cab1be79f7d61fcf49c146 Mon Sep 17 00:00:00 2001 From: Dimitri Lozeve Date: Sun, 23 Feb 2020 22:20:58 +0100 Subject: [PATCH 5/6] Switch to MOSEK solver --- Project.toml | 3 +-- src/optim.jl | 11 ++++------- 2 files changed, 5 insertions(+), 9 deletions(-) diff --git a/Project.toml b/Project.toml index dc07392..94c9ab9 100644 --- a/Project.toml +++ b/Project.toml @@ -7,12 +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" +MosekTools = "1ec41992-ff65-5c91-ac43-2df89e9693a4" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" Unicode = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" diff --git a/src/optim.jl b/src/optim.jl index a5a4bef..5b003b4 100644 --- a/src/optim.jl +++ b/src/optim.jl @@ -1,6 +1,5 @@ using JuMP -using Ipopt -using Juniper +using MosekTools function neighbourindices(mask::Matrix, idx::Int)::Vector{Int} if mask[idx] == 0 @@ -34,14 +33,12 @@ function definemodel(plantcounts::Vector, garden::Matrix, mask::Matrix, costs::M 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)) + optimizer = Mosek.Optimizer + model = Model(optimizer_with_attributes(optimizer, "QUIET" => false)) @variable(model, x[1:N, 1:Q], Bin) - @NLobjective( + @objective( model, Min, sum( From 3527a0ad0b21d17c60994b3eb851c995b41f4680 Mon Sep 17 00:00:00 2001 From: Dimitri Lozeve Date: Tue, 25 Feb 2020 21:24:00 +0100 Subject: [PATCH 6/6] Switch to CPLEX solver --- Project.toml | 2 +- src/optim.jl | 7 ++++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 94c9ab9..f529a88 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Dimitri Lozeve "] version = "0.1.0" [deps] +CPLEX = "a076750e-1247-5638-91d2-ce28b192dca0" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" @@ -11,7 +12,6 @@ JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" -MosekTools = "1ec41992-ff65-5c91-ac43-2df89e9693a4" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" Unicode = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" diff --git a/src/optim.jl b/src/optim.jl index 5b003b4..4474ba3 100644 --- a/src/optim.jl +++ b/src/optim.jl @@ -1,5 +1,5 @@ using JuMP -using MosekTools +using CPLEX function neighbourindices(mask::Matrix, idx::Int)::Vector{Int} if mask[idx] == 0 @@ -33,8 +33,9 @@ function definemodel(plantcounts::Vector, garden::Matrix, mask::Matrix, costs::M N = length(mask) Q = size(costs, 1) - optimizer = Mosek.Optimizer - model = Model(optimizer_with_attributes(optimizer, "QUIET" => false)) + optimizer = optimizer_with_attributes( + CPLEX.Optimizer, "CPX_PARAM_OPTIMALITYTARGET" => 2, "CPX_PARAM_EPINT" => 1e-8) + model = Model(optimizer) @variable(model, x[1:N, 1:Q], Bin)