### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# ╔═╡ 7845dc2e-5a8b-4a67-adb2-92df4b762df6
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add([
		Pkg.PackageSpec(name="InteractiveUtils"),
		Pkg.PackageSpec(name="JuMP"),
		Pkg.PackageSpec(name="GLPK"),
		Pkg.PackageSpec(name="XLSX"),
		Pkg.PackageSpec(name="DataFrames"),
		Pkg.PackageSpec(name="Plots")
	])
	using InteractiveUtils, JuMP, GLPK, XLSX, DataFrames, Plots
end

# ╔═╡ 61e397ee-9fa5-4925-8505-b0bc9d7df0c3
md"## Coverage-based location problems"

# ╔═╡ 3b4e2c9c-21f3-4e84-983d-97ee73675840
md" To see the location of your notebook type `pwd()`"

# ╔═╡ b00d14b1-7f36-47fe-9206-25f0738ae044
df_ambulance = DataFrame(XLSX.readtable("Data/dataset.xlsx", "ambulance")...)

# ╔═╡ da8b7342-15cb-4515-9e74-487814689540
df_demand = DataFrame(XLSX.readtable("Data/dataset.xlsx", "demand")...)

# ╔═╡ ab3ea4bc-c3bb-462c-9529-1c4895548d88
plt = scatter(df_ambulance.:x, df_ambulance."y", xlabel = "x coordinate", ylabel = "y coordinate", grid = false, label = "Ambulance")

# ╔═╡ 58cd468c-9b7a-4d14-84fd-79b72c6b1c25
scatter!(plt, df_demand.:x, df_demand."y", label = "Demand")

# ╔═╡ f35e03f4-2ba2-4471-b093-e66677057cf1
md"### Definition of a coverage"

# ╔═╡ c802dac0-d3bb-4e3a-88f0-2c81a1638b1b
begin
	for amb in eachrow(df_ambulance)
		r = 30
		th = Array(0:2*pi/100:2*pi+2*pi/100) # theta from 0 to 2pi 
		X = r*cos.(th) + ones(101,1)*amb.:x
		Y = r*sin.(th) + ones(101,1)*amb.:y
		plot!(plt, X,Y, c=:blue, legend=false)
	end
end

# ╔═╡ f60c5f40-395c-4e6f-ad41-c3b1c3fceaa2
plt

# ╔═╡ dc18dacb-9c5f-422d-9600-2494a7a693ba
md" ## How many ambulance should be located to cover all demands?"

# ╔═╡ cb3bc13c-68e5-40c7-932c-25f64933976e
md"""
### Step 1) Set the decision variable
"""

# ╔═╡ 07eea2bf-f092-4ec1-a82b-a43a63e9fb09
md"""
# 
$X_{i}: ~1~\text{if an ambulance is located at site}~i,~ 0~\text{otherwise}$
"""

# ╔═╡ aef7cb22-359c-4fe9-accd-7517f341cd4b
md"""
#### Step 2) Set the objective function
"""

# ╔═╡ 54088c2b-b4bc-4786-8a5a-e057c0ddc2f4
md"""
# 
$\min \sum_{i \in V}{X_i}$
"""

# ╔═╡ be9f41c6-64aa-4b97-a885-2c66a03ea8ce
md"""
#### Step 3) Set the constraint
"""

# ╔═╡ 766a9164-5959-4af6-ab6d-6da21f52cc9a
md"""
##### parameter
$a_{ij}:~ 1\text{if ambulance}~i~\text{can cover demand}~j,~ 0~\text{otherwise}$
threshold for coverage = 30 distance unit 

Use the Euclidean distance for the distance calculation
"""

# ╔═╡ 7ae345b1-49e9-4a6a-94d5-7db7ce5a82ef
md"""
$\sum_{i}{a_{ij}\cdot X_{i}} \geq 1~~\forall j \in D$
$X_{i} \in \{0,1\} ~~\forall i \in V$
"""

# ╔═╡ ddeb045f-028f-4abc-b0bb-433818b6c31e
md"#### Exercise 1: Implement the formulation and find the solution"

# ╔═╡ 3fe75bdc-121f-4e72-a5fe-8e0f7b8384fa
begin
	function getDistance(from, to)
		dist = sqrt((from.x-to.x)^2 + (from.y-to.y)^2)
	end
end

# ╔═╡ 17f93ee7-18c6-424f-ae99-e0490f17b46d
begin 
	nAmb = size(df_ambulance,1)
	nDemand = size(df_demand, 1)
	
	threshold = 30
	Aij = zeros(nAmb, nDemand)
	for i in 1:nAmb
		for j in 1:nDemand
			if getDistance(df_ambulance[i, :], df_demand[j, :]) <=30
				Aij[i,j] = 1
			end
		end
	end
	LocationModel = JuMP.Model()
	function setModel(m)
		@variable(m, X[1:nAmb], Bin)
		@objective(m, Min, sum(X[i] for i in 1:nAmb))
		for j in 1:nDemand
		 	@constraint(m, sum(Aij[i,j]*X[i] for i in 1:nAmb)>=1)
	    end
		#or put the for loop inside the macro
		#@constraint(m, [j in 1:nDemand], sum(Aij[i,j]*X[i] for i in 1:nAmb) >= 1)		
	end
	setModel(LocationModel)
	set_optimizer(LocationModel, GLPK.Optimizer)
	optimize!(LocationModel)
end

# ╔═╡ 746424da-9a97-483f-839d-9da9097a8dbb
termination_status(LocationModel)

# ╔═╡ 6f6ea59b-f020-4f94-8a11-2930158988f1
objective_value(LocationModel)

# ╔═╡ 93cef9e8-447d-46b8-b6cd-cb5d5f54262b
value.(LocationModel[:X])

# ╔═╡ e90579c8-4103-4972-9f6f-b788aacb7a9e
begin
	plt_sol = scatter(df_ambulance.:x, df_ambulance.:y, xlabel = "x coordinate", 							ylabel = "y coordinate", grid = false, label = "Ambulance")
	plt_sol = scatter!(plt_sol, df_demand.:x, df_demand."y", label = "Demand")
	for i in 1:nAmb
		if  value(LocationModel[:X][i]) == 1
			r = threshold
			th = Array(0:2*pi/100:2*pi+2*pi/100) # theta from 0 to 2pi 
			X = r*cos.(th) + ones(101,1)*df_ambulance.:x[i]
			Y = r*sin.(th) + ones(101,1)*df_ambulance.:y[i]
			plt_sol = plot!(plt_sol, X,Y, c=:blue, legend=false)
		end
	end
	plt_sol
end


# ╔═╡ 0b558788-94e9-43fb-9b1c-1c275a5bd2fe
md"-----
#### Exercise 2 : What if I can locate 3 ambulances only?
##### To where should I locate the ambulances to maximize the coverage?"

# ╔═╡ 3d1aae4f-c796-4fd9-a64f-59f128984d85
md"""
# 
$Y_{i}: ~1~\text{if demand}~j~\text{is covered by at least one ambulance},~ 0~\text{otherwise}$
"""

# ╔═╡ 37b04c56-adeb-402a-908b-604bee1c79bc
md"""
# 
$\max \sum_{j \in D}{Y_j}$
$\sum_{i}{a_{ij}\cdot X_{i}} \geq Y_j~~\forall j \in D$
$\sum_{i}{X_i} \leq p$
$X_{i} \in \{0,1\} ~~\forall i \in V$
$Y_{j} \in \{0,1\} ~~\forall j \in D$
"""

# ╔═╡ 4545b103-eb48-4507-910e-6b4b36b04249
begin
	MaximumCoverage = JuMP.Model()
	maxAmb = 3
	function setModel_forMaximization(m, maxAmb)
		@variable(m, X[1:nAmb], Bin)
		@variable(m, Y[1:nDemand], Bin)
		@objective(m, Max, sum(Y[j] for j in 1:nDemand))
		
		for j in 1:nDemand
			@constraint(m, sum(Aij[i,j]*X[i] for i in 1:nAmb)>=Y[j])
		end
		
		@constraint(m,  sum(X[i] for i in 1:nAmb) <= maxAmb)		
	end
	setModel_forMaximization(MaximumCoverage, maxAmb)
	set_optimizer(MaximumCoverage, GLPK.Optimizer)
	optimize!(MaximumCoverage)
end

# ╔═╡ ca438635-2b60-420d-a1e1-a483cf1a0d58
objective_value(MaximumCoverage)

# ╔═╡ 5bd776a0-44b6-4469-b05a-2a8dbf23a702
begin
	plt_sol_2 = scatter(df_ambulance.:x, df_ambulance.:y, xlabel = "x coordinate", 							ylabel = "y coordinate", grid = false, label = "Ambulance")
	scatter!(plt_sol_2, df_demand.:x, df_demand."y", label = "Demand")
	 for i in 1:nAmb
		if  value(MaximumCoverage[:X][i]) == 1
			r = threshold
			th = Array(0:2*pi/100:2*pi+2*pi/100) # theta from 0 to 2pi 
			X = r*cos.(th) + ones(101,1)*df_ambulance.:x[i]
			Y = r*sin.(th) + ones(101,1)*df_ambulance.:y[i]
			plot!(plt_sol_2, X,Y, c=:blue, legend=false)			
		end
	end
	plt_sol_2
end

# ╔═╡ Cell order:
# ╟─61e397ee-9fa5-4925-8505-b0bc9d7df0c3
# ╠═7845dc2e-5a8b-4a67-adb2-92df4b762df6
# ╟─3b4e2c9c-21f3-4e84-983d-97ee73675840
# ╠═b00d14b1-7f36-47fe-9206-25f0738ae044
# ╠═da8b7342-15cb-4515-9e74-487814689540
# ╠═ab3ea4bc-c3bb-462c-9529-1c4895548d88
# ╠═58cd468c-9b7a-4d14-84fd-79b72c6b1c25
# ╟─f35e03f4-2ba2-4471-b093-e66677057cf1
# ╠═c802dac0-d3bb-4e3a-88f0-2c81a1638b1b
# ╠═f60c5f40-395c-4e6f-ad41-c3b1c3fceaa2
# ╟─dc18dacb-9c5f-422d-9600-2494a7a693ba
# ╟─cb3bc13c-68e5-40c7-932c-25f64933976e
# ╟─07eea2bf-f092-4ec1-a82b-a43a63e9fb09
# ╟─aef7cb22-359c-4fe9-accd-7517f341cd4b
# ╟─54088c2b-b4bc-4786-8a5a-e057c0ddc2f4
# ╟─be9f41c6-64aa-4b97-a885-2c66a03ea8ce
# ╟─766a9164-5959-4af6-ab6d-6da21f52cc9a
# ╟─7ae345b1-49e9-4a6a-94d5-7db7ce5a82ef
# ╟─ddeb045f-028f-4abc-b0bb-433818b6c31e
# ╠═3fe75bdc-121f-4e72-a5fe-8e0f7b8384fa
# ╠═17f93ee7-18c6-424f-ae99-e0490f17b46d
# ╠═746424da-9a97-483f-839d-9da9097a8dbb
# ╠═6f6ea59b-f020-4f94-8a11-2930158988f1
# ╠═93cef9e8-447d-46b8-b6cd-cb5d5f54262b
# ╠═e90579c8-4103-4972-9f6f-b788aacb7a9e
# ╟─0b558788-94e9-43fb-9b1c-1c275a5bd2fe
# ╟─3d1aae4f-c796-4fd9-a64f-59f128984d85
# ╟─37b04c56-adeb-402a-908b-604bee1c79bc
# ╠═4545b103-eb48-4507-910e-6b4b36b04249
# ╠═ca438635-2b60-420d-a1e1-a483cf1a0d58
# ╠═5bd776a0-44b6-4469-b05a-2a8dbf23a702