using LinearAlgebra
using Alpine
using JuMP
using Cbc, Ipopt, Pavito, Juniper

radius = [25, 25, 7];
N = length(radius);
L = 100;
W = 60;

pavito_solver=PavitoSolver(mip_solver=CbcSolver(logLevel=0), cont_solver=IpoptSolver(print_level=0), mip_solver_drives=false, log_level=0);
juniper_solver=JuniperSolver(IpoptSolver(print_level=0,sb="yes"), mip_solver = CbcSolver(logLevel=0))

#m = Model(solver = juniper_solver)
# For global optmium (works better if we have CPLEX instead of Pavito or Cbc)
m = Model(solver = AlpineSolver(presolve_bt=0, nlp_solver = IpoptSolver(print_level=0, sb="yes"), mip_solver=pavito_solver, minlp_solver = juniper_solver))


@variable(m, delta[1:N], Bin );# switch variable of the circle 'i', whether it exists or not

@variable(m, 0 <= x[i=1:N] <= (L - radius[i]) );# x coordinate of center
@variable(m, 0 <= y[i=1:N] <= (W - radius[i]) );# y coordinate of center

@objective(m, Min, L*W - π*dot(delta,radius.^2) );


for i=1:N-1
    for j = i+1:N
        # non-convex constraint
        @NLconstraint(m, (x[i] - x[j])^2 + (y[i] - y[j])^2 >= ((radius[i] + radius[j])^2)*(delta[i] + delta[j] - 1) );
        #@NLconstraint(m, (x[i] - x[j])^2 + (y[i] - y[j])^2 >= ((radius[i] + radius[j])^2)*(delta[i] * delta[j]) );
    end
end

print(m)
status = solve(m)

println("Objective value: ", getobjectivevalue(m))
println("x = ", getvalue(x))
println("y = ", getvalue(y))
println("delta = ", getvalue(delta))