@@ -57,6 +57,15 @@ GraphsOptim.min_cost_assignment!
 graph_matching
 ```
 
+## Coloring
+
+!!! danger "Work in progress"
+Come back later!
+
+```@docs
+minimum_coloring
+```
+
 ## Utils
 
 ```@docs

@@ -7,6 +7,7 @@ module GraphsOptim
 
 using Graphs: AbstractGraph, is_directed
 using Graphs: vertices, edges, nv, ne, src, dst, inneigrs, outneigrs
+using Graphs: Coloring, maximal_cliques
 using FillArrays: Zeros, Fill
 using HiGHS: HiGHS
 using JuMP: Model
@@ -21,10 +22,12 @@ using OptimalTransport: sinkhorn
 export min_cost_flow
 export min_cost_assignment
 export FAQ, GOAT, graph_matching
+export minimum_coloring
 
 include("utils.jl")
 include("flow.jl")
 include("assignment.jl")
 include("graph_matching.jl")
+include("coloring.jl")
 
 end

@@ -0,0 +1,41 @@
+"""
+minimum_coloring(
+g, max_nb_colors;
+optimizer
+)
+
+Finds a graph coloring using the smallest possible number of colors, assuming that it will not exceed `max_nb_colors`.
+
+Returns a vector of color indices.
+
+Beware: this is an NP-complete problem and so runtime can in the worst case increase exponentially in the size of the graph.
+"""
+function minimum_coloring(
+g::AbstractGraph, max_nb_colors::Integer; optimizer=HiGHS.Optimizer
+)
+model = Model(optimizer)
+
+@variable(model, x[1:nv(g), 1:max_nb_colors], Bin)
+@variable(model, y[1:max_nb_colors], Bin)
+
+@objective(model, Min, sum(y))
+
+for v in 1:nv(g)
+@constraint(model, sum(x[v, :]) == 1)
+end
+for v in 1:nv(g), c in 1:max_nb_colors
+@constraint(model, x[v, c] <= y[c])
+end
+for c in 1:max_nb_colors
+for e in edges(g)
+u, v = src(e), dst(e)
+@constraint(model, x[u, c] + x[v, c] <= 1)
+end
+end
+
+set_silent(model)
+optimize!(model)
+# return a vector of color indices
+c = [argmax(value.(x[v, :])) for v in 1:nv(g)]
+return c
+end

@@ -0,0 +1,38 @@
+using Graphs
+using GraphsOptim
+using IterTools
+using Test
+
+function queens_graph(n::Integer, m::Integer)
+g = SimpleGraph(n * m)
+for (ix, iy, jx, jy) in product(1:n, 1:m, 1:n, 1:m)
+dx = ix - jx
+dy = iy - jy
+if dx != 0 || dy != 0
+if dx == 0 || dy == 0 || dx == dy || dx == -dy
+add_edge!(g, (ix - 1) * m + iy, (jx - 1) * m + jy)
+end
+end
+end
+return g
+end
+
+queens_graph(n::Integer) = queens_graph(n, n)
+
+function kneser_graph(n::Integer, k::Integer)
+ss = collect(subsets(1:n, k))
+return SimpleGraph([isdisjoint(a, b) for a in ss, b in ss])
+end
+
+# https://oeis.org/A088202
+
+nb_colors(c) = length(unique(c))
+
+@test minimum_coloring(queens_graph(1), 10) |> nb_colors == 1
+@test minimum_coloring(queens_graph(2), 10) |> nb_colors == 4
+@test minimum_coloring(queens_graph(3), 10) |> nb_colors == 5
+@test minimum_coloring(queens_graph(4), 10) |> nb_colors == 5
+@test minimum_coloring(queens_graph(5), 10) |> nb_colors == 5
+@test minimum_coloring(queens_graph(6), 10) |> nb_colors == 7
+
+@test minimum_coloring(kneser_graph(11, 4), 10) |> nb_colors == 11 - 2 * 4 + 2

@@ -39,4 +39,8 @@ using Test
 @testset verbose = true "Graph matching" begin
 include("graph_matching.jl")
 end
+
+@testset verbose = true "Coloring" begin
+include("coloring.jl")
+end
 end;