#************************************************************************
# Beer Flow Routing, "Mathematical Programming Modelling" (42112)
# Intro definitions
using JuMP
using HiGHS
#************************************************************************

#************************************************************************
# Data
Bars=["A","B","C","D","E","F","G"]
B=length(Bars)
X_coor=[3 9 8 2 6 10 11]
Y_coor=[10 10 6 5 2 2 7]
Distance = zeros(Float64,B,B)
for b1 in 1:B
    for b2 in 1:B
        Distance[b1,b2]= sqrt( (X_coor[b1]-X_coor[b2])*(X_coor[b1]-X_coor[b2]) + (Y_coor[b1]-Y_coor[b2])*(Y_coor[b1]-Y_coor[b2]) )
    end
end

Network=zeros(Int64,B,B)

# A-B
Network[1,2]=1
Network[2,1]=1

# A-C
Network[1,3]=1
Network[3,1]=1

# A-D
Network[1,4]=1
Network[4,1]=1

# B-C
Network[2,3]=1
Network[3,2]=1

# C-E
Network[3,5]=1
Network[5,3]=1

# C-F
Network[3,6]=1
Network[6,3]=1

# F-G
Network[6,7]=1
Network[7,6]=1

Demand=zeros(Int64,B,B)

# A  
Demand[1,2]=70

# B 
Demand[2,1]=50

# C
Demand[1,3]=80
Demand[2,3]=90

# D
Demand[1,4]=20
Demand[2,4]=50

# E
Demand[1,5]=80
Demand[2,5]=10

# F
Demand[1,6]=20
Demand[2,6]=100

# G
Demand[1,7]=40
Demand[2,7]=50

PipeCap=500

#************************************************************************

#************************************************************************
# Model
BN = Model(HiGHS.Optimizer)

# flow from node k to node l on edge from node i to node j
@variable(BN, x[k=1:B,l=1:B,i=1:B,j=1:B] >= 0)

# design binary variable
@variable(BN, y[k=1:B,i=1:B,j=1:B], Bin)
for k=1:B
    if sum(Demand[k,l] for l=1:B)==0
        for i=1:B
            for j=1:B
                fix(y[k,i,j],0; force = true)
            end
        end
    end
end

for k=1:B
    for i=1:B
        fix(y[k,i,i],0; force = true)
    end
end

# Minimize path cost
@objective(BN,
           #Min, pc
           Min, sum( Distance[i,j]*y[k,i,j] for k=1:B,i=1:B,j=1:B)
           )

# balance constraint
@constraint(BN, [i=1:B,k=1:B,l=1:B],
            sum( x[k,l,j,i] for j=1:B) - sum( x[k,l,i,j] for j=1:B) ==
            (i==k ? -Demand[k,l] : 0) + (i==l ? Demand[k,l] : 0)
            )

# network design constraint
@constraint(BN, [k=1:B,i=1:B,j=1:B],
            sum( x[k,l,i,j] for l=1:B) <= PipeCap*y[k,i,j]
            )
#************************************************************************

#************************************************************************
# solve
optimize!(BN)
println("Termination status: $(termination_status(BN))")
#************************************************************************

#************************************************************************
# Report results
if termination_status(BN) == MOI.OPTIMAL
    println("RESULTS:")
    println("BN result: $(objective_value(BN))")

    for k=1:B
        for i=1:B
            for j=1:B
                if value(y[k,i,j]) > 0.001
                    println("beer: $(k) y:  $(Bars[i]) -> $(Bars[j]) val: $(value(y[k,i,j]))  Distance: $(Distance[i,j])")   
                end
            end
        end
    end
else
  println("  No solution")
end
#************************************************************************