Aller au contenu

Vertex Cover et TSP

On va reprendre ici certains algorithmes étudiés dans le cours pour réaliser des approximations.

MinVertexCover

On considère un graphe non orienté donné en OCaml uniquement par ses listes d'adjacence en identifiant un sommet à son indice.

let g = [|[1];[0;2];[1;3;4];[2;4;5;6];[2;3;5];[3;4;6];[3;5]|]

Question 1. Écrire une fonction aretes qui prend un graphe et renvoie la liste de ses arêtes sous la forme d'une liste de couples \((i,j)\) où \(i < j\).

Corrigé 
let aretes ladj =
    ladj
    |> Array.mapi (fun i l -> (* pour chaque sommet i et liste adj l *)
         l |> List.filter ((<)i) (* on ne garde que les voisins > i *)
           |> List.map (fun j -> (i,j))) (* on transforme les voisins en aretes *)
    |> Array.fold_left (@) [] (* on concatene le tableau de listes *)

Sur g, on obtient [(0,1); (1,2); (2,3); (2,4); (3,4); (3,5); (3,6); (4,5); (5,6)].

Question 2. Écrire une fonction retire_coincidentes qui prend une arête \(a\) et une liste d'arêtes \(l\) et renvoie la liste \(l' \subset l\) des arêtes qui ne partagent pas d'extrémités avec \(a\).

Corrigé 
let retire_coincidentes (x,y) =
    List.filter (fun (x',y') -> x <> x' && y <> y' && x <> y' && y <> x')

Question 3. En déduire une fonction min_vertex_cover_2a qui renvoie une couverture de sommets par la 2-approximation du cours.

Corrigé 
let min_vertex_cover_2a ladj =
     let rec aux a =
        match a with
        | [] -> []
        | (x,y)::q -> x::y::aux (retire_coincidentes (x,y) q)
     in aux (aretes ladj)

Question 4. Testez vos fonctions.

Corrigé

Sur g, on obtient [0; 1; 2; 3; 4; 5], soit 6 sommets, alors qu'une couverture optimale en a 4 (par exemple \(\{1, 3, 4, 5\}\) ; il en faut au moins un pour l'arête \(\{0,1\}\) et au moins trois parmi \(2, \dots, 6\) pour couvrir les triangles \(2,3,4\), \(3,4,5\) et \(3,5,6\) et l'arête \(\{5, 6\}\)). On a bien \(6 \le 2 \times 4\).

Pour tester plus systématiquement, on peut comparer sur des graphes aléatoires à une recherche exhaustive de la couverture optimale (on teste les \(2^n\) sous-ensembles de sommets pour \(n \le 10\)).

Union-find

Question 5. Reprendre la définition de la structure Union-Find en OCaml (voir le cours).

On pourra par exemple considérer le type :

type 'a union_find = {
    values_to_index : ('a, int) Hashtbl.t;
    parent : int array;
    rank : int array
}

où values_to_index permet d'accéder aux éléments par leur valeur.

On pourra alors définir des fonctions :

(* crée un union_find *)
let make_sets (l:'a list) : 'a union_find = ...
(* renvoie l'indice de la racine de l'ensemble contenant x *)
let find (uf:'a union_find) (x:'a) : int = ...
(* fusionne deux ensembles donnés par leur racine *)
let union (uf:'a union_find) (i:int) (j:int) = ...
Corrigé 
(* crée un union_find *)
let make_sets (l:'a list) : 'a union_find =
    let n = List.length l in
    let h = Hashtbl.create n in
    List.iteri (fun i a -> Hashtbl.add h a i) l;
    { values_to_index = h; parent = Array.init n Fun.id; rank = Array.make n 0 }

(* renvoie l'indice de la racine de l'ensemble contenant x *)
let find (uf:'a union_find) (x:'a) : int =
    let rec remonte i =
        if uf.parent.(i) = i
        then i
        else remonte uf.parent.(i)
    in remonte (Hashtbl.find uf.values_to_index x)

(* fusionne deux ensembles donnés par leur racine *)
let union (uf:'a union_find) (i:int) (j:int) =
    if i <> j
    then begin
        let ri = uf.rank.(i) in
        let rj = uf.rank.(j) in
        if ri = rj
        then begin
            uf.parent.(j) <- i;
            uf.rank.(i) <- ri + 1
        end else if ri < rj then uf.parent.(i) <- j else uf.parent.(j) <- i
    end

let size uf = Array.length uf.parent
let mem uf x = Hashtbl.mem uf.values_to_index x

Application à 3SatAllOrNothing

On rappelle le problème 3SatAllOrNothing (voir les exercices de calculabilité) : il s'agit étant donné une formule de 3-SAT de déterminer s'il existe une valuation telle que chaque clause contient soit trois littéraux vrais, soit trois littéraux faux.

Pour résoudre ce problème, on va définir une relation d'équivalence sur les littéraux : la plus petite relation d'équivalence telle que \(l \equiv l'\) dès que \(l\) et \(l'\) sont dans une même clause de la formule. Deux littéraux équivalents doivent avoir la même valeur. Il suffira alors de vérifier dans un second temps qu'on peut donner une valeur à chaque classe de façon cohérente : les classes de \(x\) et de \(\neg x\) doivent recevoir des valeurs opposées, ce qui revient à 2-colorier le graphe dont les sommets sont les classes et dont les arêtes relient la classe de \(x\) à celle de \(\neg x\).

On considère le type suivant :

type lit = Pos of char | Neg of char
type sat3 = (lit * lit * lit) list

Question 6. Résoudre le problème en utilisant une structure union_find.

Corrigé 
let litteraux f =
    List.map (fun (i,j,k) -> [i;j;k]) f
    |> List.concat
    |> List.sort_uniq Stdlib.compare

let variables f =
   litteraux f
   |> List.map (function Pos x | Neg x -> x)
   |> List.sort_uniq Stdlib.compare

let construit_classes f =
    let uf = make_sets (litteraux f) in
    List.iter (fun (l1, l2, l3) ->
           union uf (find uf l1) (find uf l2);
           union uf (find uf l1) (find uf l3)) f;
    uf

let deux_coloriable ladj =
   let n = Array.length ladj in
   let echec = ref false in
   let couleurs = Array.make n (-1) in
   let rec parcours couleurs col x =
       if couleurs.(x) = -1
       then begin
          couleurs.(x) <- col;
          List.iter (fun j -> parcours couleurs ((col+1) mod 2) j) ladj.(x)
       end else if couleurs.(x) <> col then echec := true
   in
   for i = 0 to n - 1 do
       if couleurs.(i) = -1
       then parcours couleurs 0 i
   done;
   not !echec, couleurs

let resout f =
   let var = variables f in
   let classes = construit_classes f in
   let ladj = Array.make (size classes) [] in
   List.iter (fun v ->
       if mem classes (Pos v) && mem classes (Neg v)
       then begin
           let pi = find classes (Pos v) and ni = find classes (Neg v) in
           if not (List.mem ni ladj.(pi)) then ladj.(pi) <- ni :: ladj.(pi);
           if not (List.mem pi ladj.(ni)) then ladj.(ni) <- pi :: ladj.(ni)
       end) var;
    deux_coloriable ladj

resout [ (Pos 'a', Pos 'b', Pos 'c');
    (Pos 'a', Pos 'b', Neg 'd');
    (Pos 'd', Neg 'b', Neg 'c')]

La première composante du résultat vaut true ici. Si \(x\) et \(\neg x\) sont dans la même classe, la classe est reliée à elle-même et le 2-coloriage échoue, comme attendu. Le tableau couleurs (indexé par les racines) donne une valuation : par exemple, les littéraux de couleur \(0\) sont vrais.

TSP métrique

Le problème du TSP métrique consiste, étant donné un graphe complet \(G = (S,A)\) non orienté et une fonction de distance \(d : S^2 \rightarrow \mathbb{R}_+\) vérifiant l'inégalité triangulaire

\[ \forall x, y, z \in S, d(x,y) \le d(x,z) + d(z,y) \]

à déterminer la longueur minimale d'un cycle hamiltonien de \(G\).

Ce problème est NP-difficile. On va donner une 2-approximation en considérant l'algorithme suivant :

  • on construit un arbre couvrant de poids minimal (par exemple, avec Kruskal) ;
  • on considère le cycle induit par le préordre d'un parcours de l'arbre couvrant.

Question 7. On vous demande ici de mettre en place cette approximation pour résoudre l'instance suivante :

let points = [|
 (525,187); (295,946); (364,787); (533,5); (203,862); (625,11); (824,675);
 (369,58); (308,785); (726,766); (155,22); (197,879); (881,599); (862,953);
 (328,391); (983,744); (426,695); (478,583); (821,911); (79,60); (665,272);
 (191,773); (114,281); (800,567); (222,693); (846,961); (215,740); (353,769);
 (750,726); (908,909); (309,658); (867,649); (592,106); (16,637); (54,796);
 (449,649); (980,654); (594,819); (764,101); (921,307); (928,647); (738,196);
 (550,578); (478,489); (449,755); (189,415); (433,545); (967,478); (149,416);
 (144,117); (298,972); (854,144)
 |]
Corrigé 
type 'a arbre = Noeud of 'a * 'a arbre list

let arbre valeurs ladj racine =
    let visite = Array.make (Array.length ladj) false in
    let rec construit x =
        visite.(x) <- true;
        Noeud(valeurs.(x),
            List.map construit (List.filter (fun y -> not visite.(y)) ladj.(x)))
    in construit !racine

let rec preordre n =
    let Noeud(x,l) = n in
    x :: List.concat (List.map preordre l)

let kruskal dist points =
    let n = Array.length points in
    let aretes = ref [] in
    for i = 0 to n-1 do
       for j = i+1 to n-1 do
           aretes := (i, j, dist points.(i) points.(j))::!aretes
       done
    done;
    let aretes = List.sort (fun (_,_,a) (_,_,b) -> Stdlib.compare a b) !aretes in
    let ladj = Array.make n [] in
    let uf = make_sets (List.init n Fun.id) in
    List.iter (fun (i,j,_) ->
        let ri = find uf i and rj = find uf j in
        if ri <> rj
        then begin
            union uf ri rj;
            ladj.(i) <- j :: ladj.(i);
            ladj.(j) <- i :: ladj.(j)
        end) aretes;
    let racine = ref 0 in
    while uf.parent.(!racine) <> !racine do
        incr racine
    done;
    arbre points ladj racine

let longueur dist l =
    let p0 = List.hd l in
    let rec aux l =
       match l with
       | [] -> 0.
       | [p] -> dist p0 p
       | p1::p2::q -> dist p1 p2 +. aux (p2::q)
    in aux l

let dist_euclidienne (x,y) (x',y') =
   sqrt (float_of_int ((x-x')*(x-x') + (y-y')*(y-y')))

let approx_tsp dist =
    longueur dist (preordre (kruskal dist points))

approx_tsp dist_euclidienne

On obtient un cycle de longueur environ \(7521{,}27\), pour un arbre couvrant minimal de poids environ \(4809{,}02\) : la longueur optimale est donc comprise entre \(4809\) et \(7521\).

Question 8. Même question, mais cette fois-ci, on considère que les points sont sur un tore. Concrètement, cela signifie qu'ils sont sur un carré de dimension \(\range{0}{999}^2\) dont les bords gauche/droite et haut/bas sont identifiés.

Note : il s'agit d'un challenge présent sur le site http://caesum.com/game pour avoir le droit d'être anonyme sur ses forums. Le record sur le site est de 5107,36 le 9 novembre 2001.

Corrigé 
let dist_tore (x,y) (x',y') =
    let dx = min (min (abs (x-x')) (abs (x-x'-1000))) (abs (1000+x-x')) in
    let dy = min (min (abs (y-y')) (abs (y-y'-1000))) (abs (1000+y-y')) in
    sqrt (float_of_int (dx * dx + dy * dy))

approx_tsp dist_tore

On obtient environ \(6556{,}18\) (arbre couvrant minimal de poids environ \(4508{,}06\)), assez loin du record.