#Installation, Required first time write:
#]
#add JuMP#v0.18.5
#add Clp
#add MathProgBase
using JuMP      #load the package JuMP
using Clp       #load the package Clp (an open linear-programming solver)
#using Gurobi   #The commercial optimizer Gurobi requires installation

#The ? can be put infront of commands, variables, and functions to get more information.
#Note that is applied on the main object, hence to ask about an element in an array do:
#element = array[5]
#?element

#Build the model and get variables and constraints back (see intro_mod.jl)
include("intro_mod.jl")
m, x, nutrition_demands = build_diet_model("intro_dat.jl")
print(m) # prints the model instance

setsolver(m, ClpSolver(LogLevel = 1))
#setsolver(m, GurobiSolver())
solve(m)

println("z =  ", getobjectivevalue(m))   		# display the optimal solution
println("x =  ", getvalue(x))
println("Foods in solution: ", foods[[getvalue(x[i]) > 0 for i in I]])
println("reduced cost =  ", getdual(x))
protein_demand = nutrition_demands[findfirst(nutrients .== "protein")]
println("protein dual =  ", protein_demand)

# You can always define aid functions to simply your life, as below
# Moreover, it's good practice to place this functions in a seperate file
# and use include("lp_util_functions.jl"), to keep the code structured.
"""
 Gets the current slack of the constraint, for feasible solution it's always positive.
 For double sided inequalities it's the least slack that is given.
"""
function getslack(constraint::ConstraintRef)::Float64  # If you dont want you dont have to specify types
  lin_constr = LinearConstraint(constraint)
  row_val = getvalue(lin_constr.terms)
  return min(lin_constr.ub - row_val, row_val - lin_constr.lb)
end
fat_demand = nutrition_demands[findfirst(nutrients .== "fat")]
println("fat slack =  ", getslack(fat_demand))


# Note the level of sodium, modify the model to restrict it e.g.
amount_of_sodium = @expression(m, sum(N[i,4]*x[i] for i in I))
sodium_constarint = @constraint(m, amount_of_sodium <= 2000 )
solve(m)
println("amount of sodium", getvalue(amount_of_sodium))

# And modify the constraint
JuMP.setRHS(sodium_constarint, 1500)
status = solve(m)
println("Solve status = ", status)
# Multiple times
JuMP.setRHS(sodium_constarint, 2000)
status = solve(m)

#To modify the objective, change c then change the objective by calling:
#@objective(m, Min, sum(c[i]*x[i] for i in I))

# Available through the backend MathProgBase, its possible to get the basis status (on LP problems)
using MathProgBase
"""
Get the basis status of a variable, i.e., if it's in the optmial basis or at is lower or upper bound
"""
function getbasisstatus(var::Variable)
  col_basis, row_basis = MathProgBase.getbasis(internalmodel(var.m))
  return col_basis[var.col]
end
println("Basis status of ", foods[1], " = ", getbasisstatus(x[1]))
println("Basis status of ", foods[2], " = ", getbasisstatus(x[2]))
"""
Get the basis status of a constraint, i.e., if the assosiated slack in the basis or the constraint is at its lower or upper bound
"""
function getbasisstatus(constraint::ConstraintRef)
  col_basis, row_basis = MathProgBase.getbasis(internalmodel(constraint.m))
  return row_basis[constraint.idx]
end
println("Basis status of fat demand = ", getbasisstatus(fat_demand))
println("Basis status of protein demand = ", getbasisstatus(protein_demand))