using JuMP, GLPK, Cbc

function graph_coloring(; fast = true, benchmark = false)
    m = Model(optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0))

    @variable(m, max_color, Int)
    @variable(m, c[1:49, 1:4], Bin)

    @objective(m, Min, max_color)

    for s in 1:49
        @constraint(m, sum(c[s, 1:4]) == 1)
    end

    if fast
        for s in 1:49
            @constraint(m, sum(c[s, i] * i for i in 1:4) <= max_color)
        end
    else
        for s in 1:49, i in 1:4
            @constraint(m, max_color >= i * c[s, i])
        end
    end

    for i in 1:4
        @constraint(m, c[1, i] + c[11, i] <= 1)
        @constraint(m, c[1, i] + c[8, i] <= 1)
        @constraint(m, c[11, i] + c[8, i] <= 1)
        @constraint(m, c[11, i] + c[22, i] <= 1)
        @constraint(m, c[11, i] + c[24, i] <= 1)
        @constraint(m, c[24, i] + c[22, i] <= 1)
        @constraint(m, c[24, i] + c[36, i] <= 1)
        @constraint(m, c[22, i] + c[8, i] <= 1)
        @constraint(m, c[22, i] + c[23, i] <= 1)
        @constraint(m, c[22, i] + c[36, i] <= 1)
        @constraint(m, c[8, i] + c[2, i] <= 1)
        @constraint(m, c[8, i] + c[6, i] <= 1)
        @constraint(m, c[8, i] + c[23, i] <= 1)
        @constraint(m, c[23, i] + c[6, i] <= 1)
        @constraint(m, c[23, i] + c[31, i] <= 1)
        @constraint(m, c[23, i] + c[41, i] <= 1)
        @constraint(m, c[23, i] + c[36, i] <= 1)
        @constraint(m, c[36, i] + c[31, i] <= 1)
        @constraint(m, c[36, i] + c[41, i] <= 1)
        @constraint(m, c[2, i] + c[4, i] <= 1)
        @constraint(m, c[2, i] + c[5, i] <= 1)
        @constraint(m, c[2, i] + c[6, i] <= 1)
        @constraint(m, c[6, i] + c[5, i] <= 1)
        @constraint(m, c[6, i] + c[15, i] <= 1)
        @constraint(m, c[6, i] + c[31, i] <= 1)
        @constraint(m, c[31, i] + c[15, i] <= 1)
        @constraint(m, c[31, i] + c[33, i] <= 1)
        @constraint(m, c[31, i] + c[37, i] <= 1)
        @constraint(m, c[31, i] + c[41, i] <= 1)
        @constraint(m, c[41, i] + c[37, i] <= 1)
        @constraint(m, c[41, i] + c[40, i] <= 1)
        @constraint(m, c[4, i] + c[10, i] <= 1)
        @constraint(m, c[4, i] + c[5, i] <= 1)
        @constraint(m, c[5, i] + c[10, i] <= 1)
        @constraint(m, c[5, i] + c[13, i] <= 1)
        @constraint(m, c[5, i] + c[15, i] <= 1)
        @constraint(m, c[15, i] + c[13, i] <= 1)
        @constraint(m, c[15, i] + c[35, i] <= 1)
        @constraint(m, c[15, i] + c[33, i] <= 1)
        @constraint(m, c[33, i] + c[35, i] <= 1)
        @constraint(m, c[33, i] + c[37, i] <= 1)
        @constraint(m, c[37, i] + c[46, i] <= 1)
        @constraint(m, c[37, i] + c[40, i] <= 1)
        @constraint(m, c[40, i] + c[46, i] <= 1)
        @constraint(m, c[40, i] + c[47, i] <= 1)
        @constraint(m, c[10, i] + c[7, i] <= 1)
        @constraint(m, c[10, i] + c[13, i] <= 1)
        @constraint(m, c[13, i] + c[7, i] <= 1)
        @constraint(m, c[13, i] + c[26, i] <= 1)
        @constraint(m, c[13, i] + c[35, i] <= 1)
        @constraint(m, c[35, i] + c[26, i] <= 1)
        @constraint(m, c[35, i] + c[32, i] <= 1)
        @constraint(m, c[35, i] + c[39, i] <= 1)
        @constraint(m, c[35, i] + c[46, i] <= 1)
        @constraint(m, c[46, i] + c[39, i] <= 1)
        @constraint(m, c[46, i] + c[43, i] <= 1)
        @constraint(m, c[46, i] + c[47, i] <= 1)
        @constraint(m, c[47, i] + c[43, i] <= 1)
        @constraint(m, c[7, i] + c[26, i] <= 1)
        @constraint(m, c[7, i] + c[49, i] <= 1)
        @constraint(m, c[26, i] + c[21, i] <= 1)
        @constraint(m, c[26, i] + c[32, i] <= 1)
        @constraint(m, c[32, i] + c[21, i] <= 1)
        @constraint(m, c[32, i] + c[25, i] <= 1)
        @constraint(m, c[32, i] + c[29, i] <= 1)
        @constraint(m, c[32, i] + c[34, i] <= 1)
        @constraint(m, c[32, i] + c[39, i] <= 1)
        @constraint(m, c[39, i] + c[34, i] <= 1)
        @constraint(m, c[39, i] + c[38, i] <= 1)
        @constraint(m, c[39, i] + c[44, i] <= 1)
        @constraint(m, c[39, i] + c[42, i] <= 1)
        @constraint(m, c[39, i] + c[43, i] <= 1)
        @constraint(m, c[43, i] + c[42, i] <= 1)
        @constraint(m, c[49, i] + c[21, i] <= 1)
        @constraint(m, c[49, i] + c[25, i] <= 1)
        @constraint(m, c[21, i] + c[25, i] <= 1)
        @constraint(m, c[42, i] + c[44, i] <= 1)
        @constraint(m, c[42, i] + c[48, i] <= 1)
        @constraint(m, c[25, i] + c[17, i] <= 1)
        @constraint(m, c[25, i] + c[29, i] <= 1)
        @constraint(m, c[3, i] + c[12, i] <= 1)
        @constraint(m, c[12, i] + c[9, i] <= 1)
        @constraint(m, c[12, i] + c[14, i] <= 1)
        @constraint(m, c[9, i] + c[14, i] <= 1)
        @constraint(m, c[9, i] + c[16, i] <= 1)
        @constraint(m, c[14, i] + c[16, i] <= 1)
        @constraint(m, c[14, i] + c[19, i] <= 1)
        @constraint(m, c[14, i] + c[18, i] <= 1)
        @constraint(m, c[18, i] + c[19, i] <= 1)
        @constraint(m, c[16, i] + c[18, i] <= 1)
        @constraint(m, c[16, i] + c[17, i] <= 1)
        @constraint(m, c[16, i] + c[20, i] <= 1)
        @constraint(m, c[17, i] + c[20, i] <= 1)
        @constraint(m, c[17, i] + c[28, i] <= 1)
        @constraint(m, c[17, i] + c[30, i] <= 1)
        @constraint(m, c[17, i] + c[29, i] <= 1)
        @constraint(m, c[20, i] + c[28, i] <= 1)
        @constraint(m, c[30, i] + c[27, i] <= 1)
        @constraint(m, c[30, i] + c[29, i] <= 1)
        @constraint(m, c[30, i] + c[34, i] <= 1)
        @constraint(m, c[27, i] + c[34, i] <= 1)
        @constraint(m, c[29, i] + c[34, i] <= 1)
        @constraint(m, c[34, i] + c[38, i] <= 1)
        @constraint(m, c[38, i] + c[45, i] <= 1)
        @constraint(m, c[45, i] + c[44, i] <= 1)
        @constraint(m, c[44, i] + c[48, i] <= 1)
    end


    optimize!(m)
    if !benchmark
        println(JuMP.objective_value(m))
        println(JuMP.termination_status(m))
        _, color_washington = findmax(JuMP.value.(c[1, :]))
        println("color_washington: ", color_washington)
    end
end
