Aller au contenu

Algorithme de Ford-Fulkerson en C

Ce TP accompagne les parties Couplage maximal dans un graphe biparti du chapitre Algorithmique avancée des graphes. Il reprend en C le TP Couplage avant de passer aux flots.

Tous les programmes de ce TP utilisent stdio.h, stdlib.h, string.h, stdbool.h et, pour la dernière partie, limits.h.

Manipulations de base sur les graphes

On pourra sauter cette partie si on est à l'aise sur ces notions et partir du fichier tp_couplage_1.c.

On va considérer le type suivant pour des graphes avec listes chaînées d'adjacence :

struct edge_list {
    int dst;
    struct edge_list *next;
};
typedef struct edge_list edge_list;

struct graph {
    int nvertices;
    bool directed;
    edge_list **edge_lists;
};
typedef struct graph graph;

Question 1. Écrire une fonction de prototype

graph *graph_create(int nvertices, bool directed);

qui crée un graphe vide à nvertices sommets et dont le caractère orienté est donné par le booléen directed.

Corrigé 
graph *graph_create(int nvertices, bool directed) {
    graph *g = malloc(sizeof(graph));

    g->nvertices = nvertices;
    g->edge_lists = malloc(sizeof(edge_list *) * nvertices);
    for(int i = 0; i < nvertices; i++)
        g->edge_lists[i] = NULL;
    g->directed = directed;

    return g;
}

Question 2. Écrire des fonctions de prototypes

void edge_list_free(edge_list *l);
void graph_free(graph *g);

permettant de libérer la mémoire occupée par un graphe.

Corrigé 
void edge_list_free(edge_list *l)
{
    while(l != NULL)
    {
        edge_list *t = l->next;
        free(l);
        l = t;
    }
}

void graph_free(graph *g)
{
    for(int i = 0; i < g->nvertices; i++)
        edge_list_free(g->edge_lists[i]);
    free(g->edge_lists);
    free(g);
}

Question 3. Écrire une fonction de prototype

bool is_connected(graph *g, int x, int y);

qui teste s'il existe une arête \(x \rightarrow y\).

Corrigé 
bool is_connected(graph *g, int x, int y)
{
    edge_list *l = g->edge_lists[x];
    while(l != NULL)
    {
        if(l->dst == y)
            return true;
        l = l->next;
    }
    return false;
}

Question 4. Écrire une fonction de prototype

void add_edge(graph *g, int src, int dst);

qui ajoute une arête de src à dst dans le graphe. On ne fera aucune vérification quant au fait que l'arête soit déjà présente.

Corrigé 
void add_edge(graph *g, int src, int dst)
{
    edge_list *e = malloc(sizeof(edge_list));
    e->dst= dst;
    e->next = g->edge_lists[src];
    g->edge_lists[src] = e;
    if (!g->directed)
    {
        edge_list *e = malloc(sizeof(edge_list));
        e->dst= src;
        e->next = g->edge_lists[dst];
        g->edge_lists[dst] = e;
    }
}

Question 5. Écrire une fonction de prototype

void add_edges(graph *g, int *vtx_pairs, int npairs);

qui prend en entrée un tableau t de 2*npairs entiers et qui rajoute des arêtes de t[2i] à t[2i+1].

Corrigé 
void add_edges(graph *g, int *vtx_pairs, int npairs)
{
    for(int i = 0; i < npairs*2; i+=2)
    {
        add_edge(g, vtx_pairs[i], vtx_pairs[i+1]);
    }
}

Question 6. Écrire une fonction de prototype

void graph_print(graph *g);

qui affiche le graphe sous un format lisible. On s'en servira pour faire des vérifications.

Corrigé 
void graph_print(graph *g)
{
    for(int i = 0; i < g->nvertices; i++)
    {
        printf("%d : ", i);
        edge_list *l = g->edge_lists[i];
        while(l != NULL)
        {
            printf(" %d,", l->dst);
            l = l->next;
        }
        printf("\n");
    }
}

Question 7. Écrire une fonction de prototype

int *dfs_init(graph *g, int src);

qui renvoie un tableau pred valant -1 pour tous les sommets autres que src et valant src pour src.

Corrigé 
int *dfs_init(graph *g, int src)
{
    int *pred = malloc(sizeof(int) * g->nvertices);
    memset(pred, -1, sizeof(int) * g->nvertices);
    pred[src] = src;
    return pred;
}

memset remplit octet par octet : cela fonctionne pour -1 car tous les octets de l'entier -1 valent 0xFF, mais ce ne serait pas le cas pour une autre valeur comme 1.

Question 8. Écrire une fonction de prototype

void dfs(graph *g, int *pred, int x);

qui effectue une étape récursive d'un DFS en partant du sommet x et en remplissant le tableau de prédécesseurs pred : pred[y] = x si y est découvert par le sommet x.

Corrigé 
void dfs(graph *g, int *pred, int x)
{
    edge_list *l = g->edge_lists[x];
    while(l != NULL)
    {
        int y = l->dst;
        if (pred[y] < 0)
        {
            pred[y] = x;
            dfs(g, pred, y);
        }
        l = l->next;
    }
}

Question 9. Écrire une fonction de prototype

int *dfs_launch(graph *g, int src);

qui renvoie le tableau pred à l'issue d'un DFS depuis src.

Corrigé 
int *dfs_launch(graph *g, int src)
{
    int *pred = dfs_init(g, src);
    dfs(g, pred, src);
    return pred;
}

Retour sur les couplages de graphe biparti

On pourra sauter cette partie si on est à l'aise sur ces notions et partir du fichier tp_couplage_2.c.

On va reprendre les notations du TP précédent. Pour indiquer les deux parties des sommets d'un graphe biparti, on va appeler les premiers sommets des sources et les autres des buts. L'appartenance à un ensemble ou à l'autre est déterminée par un tableau de booléens sources.

graph *g1 = graph_create(5, false);
int vtx_pairs1[] = { 0, 3, 1, 3, 1, 4, 2, 3, 2, 4 };
bool sources1[] = { true, true, true, false, false };
add_edges(g1, vtx_pairs1, 5);

graph *g2 = graph_create(12, false);
int vtx_pairs2[] = { 0, 8, 0, 7, 2, 6, 2, 9, 3, 8, 4, 8, 4, 9, 5, 11 };
bool sources2[] = { true, false, true, true, true, true,
    false, false, false, false, false, false };
add_edges(g2, vtx_pairs2, 8);

Pour les couplages, on va considérer un tableau d'entiers matching où matching[x] = y si \(\{x,y\}\) est dans le couplage, et matching[x] = -1 si \(x\) est libre.

Question 10. Écrire une fonction de prototype

graph *residual_graph(graph *g, bool *sources, int *matching);

qui renvoie le graphe résiduel orienté associé au couplage du graphe biparti.

Corrigé 
graph *residual_graph(graph *g, bool *sources, int *matching)
{
    graph *gr = graph_create(g->nvertices, true);
    for(int i = 0; i < g->nvertices; i++)
    {
        edge_list *l = g->edge_lists[i];
        while(l != NULL)
        {
            if((sources[i] && matching[i] != l->dst)
               || (!sources[i] && matching[i] == l->dst))
                add_edge(gr, i, l->dst);
            l = l->next;
        }
    }
    return gr;
}

Question 11. Écrire une fonction de prototype

void swap_along_path(int *matching, int *pred, int x);

qui considère un sommet but libre x atteint par un DFS depuis un sommet source libre dans le graphe résiduel et bascule les arêtes du chemin correspondant en utilisant le tableau pred des prédécesseurs.

Corrigé 
void swap_along_path(int *matching, int *pred, int x)
{
    bool swap = true;
    while(pred[x] != x)
    {
        if(swap)
        {
            matching[x] = pred[x];
            matching[pred[x]] = x;
        }
        swap = !swap;
        x = pred[x];
    }
}

Question 12. Écrire une fonction de prototype

bool augment(graph *g, bool *sources, int *matching, int *pred);

qui considère un graphe biparti, un couplage et un tableau pred issu d'un DFS depuis une source libre et renvoie true après avoir basculé un chemin augmentant s'il en existe un, et false sinon.

Corrigé 
bool augment(graph *g, bool *sources, int *matching, int *pred)
{
    for(int i =0; i < g->nvertices; i++)
    {
        if(!sources[i] && matching[i] < 0 && pred[i] >= 0)
        {
            swap_along_path(matching, pred, i);
            return true;
        }
    }

    return false;
}

Question 13. Écrire une fonction de prototype

int *maximal_matching(graph *g, bool *sources);

qui renvoie un tableau matching décrivant un couplage maximal.

Corrigé 
int *maximal_matching(graph *g, bool *sources)
{
    int *matching = malloc(sizeof(int) * g->nvertices);
    memset(matching, -1, sizeof(int) * g->nvertices);
    for(int i = 0; i < g->nvertices; i++)
    {
        if(sources[i] && matching[i] < 0)
        {
            graph *gr = residual_graph(g, sources, matching);
            int *pred = dfs_launch(gr, i);
            augment(g, sources, matching, pred);
            free(pred);
            graph_free(gr);
        }
    }
    return matching;
}

Comme dans le TP en OCaml, il suffit de chercher une seule fois un chemin augmentant depuis chaque source (voir le cours).

Ford-Fulkerson

Note si vous vous sentez à l'aise, il peut être intéressant de directement essayer d'implémenter Ford-Fulkerson. On fera attention à préserver la complexité, ce qui implique de bien réfléchir à l'organisation de votre code. C'est un exercice délicat mais très instructif. Cette partie vous propose un choix d'implémentation parmi d'autres.

Pour pouvoir travailler sur l'algorithme de Ford-Fulkerson, il va être nécessaire de reporter des changements de flot le long d'un chemin dans le graphe résiduel vers le graphe initial. On va donc avoir besoin de considérer les arêtes dans un chemin et pas seulement les sommets donnés par la relation de prédécesseurs et on va avoir besoin d'accéder directement à une arête du graphe initial depuis une des deux arêtes créées dans le graphe résiduel.

C'est pour cela qu'on enrichit le type précédent :

struct edge_list {
    int src;
    int dst;
    int flow; // le flot sur l'arete
    int capacity; // la capacite de l'arete
    struct edge_list *ref; // NULL ou l'arete associee
    struct edge_list *next;
};
typedef struct edge_list edge_list;

struct graph {
    int s; /* graphe s-->*-->t */
    int t;
    int nvertices;
    edge_list **edge_lists;
};
typedef struct graph graph;

Question 14. Reprendre la fonction graph_create pour tenir compte de s et t.

graph *graph_create(int nvertices, int s, int t);

Ici, le graphe est forcément orienté.

Corrigé 
graph *graph_create(int nvertices, int s, int t) {
    graph *g = malloc(sizeof(graph));

    g->nvertices = nvertices;
    g->s = s;
    g->t = t;
    g->edge_lists = malloc(sizeof(edge_list *) * nvertices);
    for(int i = 0; i < nvertices; i++)
        g->edge_lists[i] = NULL;

    return g;
}

Les fonctions edge_list_free et graph_free sont inchangées.

Question 15. Reprendre la fonction add_edge pour tenir compte de cela, cette fonction affectera un flot nul et une référence NULL et renverra l'arête créée :

edge_list *add_edge(graph *g, int src, int dst, int capacity);
Corrigé 
edge_list *add_edge(graph *g, int src, int dst, int capacity)
{
    edge_list *e = malloc(sizeof(edge_list));
    e->src= src;
    e->dst= dst;
    e->next = g->edge_lists[src];
    g->edge_lists[src] = e;
    e->capacity = capacity;
    e->flow = 0;
    e->ref = NULL;

    return e;
}

Question 16. Reprendre la fonction graph_print pour afficher le flot.

void graph_print(graph *g);
Corrigé 
void graph_print(graph *g)
{
    for(int i = 0; i < g->nvertices; i++)
    {
        edge_list *l = g->edge_lists[i];
        while(l != NULL)
        {
            printf("%d ---%d/%d---> %d\n", i, l->flow, l->capacity, l->dst);
            l = l->next;
        }
    }
}

Question 17. Reprendre les fonctions de parcours pour qu'elles correspondent aux prototypes suivants :

edge_list **dfs_init(graph *g);
void dfs(graph *g, edge_list **pred, int x);
edge_list **dfs_launch(graph *g);

Tous les DFS démarreront de g->s.

Ici, le tableau des prédécesseurs indique l'arête empruntée pour remonter. Comme la source g->s du DFS aura une valeur NULL dans ce tableau, il est nécessaire de la tester explicitement pour ne pas la considérer à tort comme non visitée.

Corrigé 
edge_list **dfs_init(graph *g)
{
    edge_list **pred = malloc(sizeof(edge_list *) * g->nvertices);
    memset(pred, 0, sizeof(edge_list *) * g->nvertices);
    return pred;
}

void dfs(graph *g, edge_list **pred, int x)
{
    edge_list *l = g->edge_lists[x];
    while(l != NULL)
    {
        int y = l->dst;
        if (pred[y] == NULL && y != g->s)
        {
            pred[y] = l;
            dfs(g, pred, y);
        }
        l = l->next;
    }
}

edge_list **dfs_launch(graph *g)
{
    edge_list **pred = dfs_init(g);
    dfs(g, pred, g->s);
    return pred;
}

Question 18. Écrire une fonction de prototype

graph *residual_graph(graph *g);

qui renvoie le graphe résiduel associé au flot présent sur les arêtes de g.

On fera en sorte que les arêtes du graphe résiduel fassent référence aux arêtes de g dont elles proviennent. On va utiliser le champ capacity pour la capacité résiduelle et le champ flow aura la valeur 1 si on est dans l'orientation de g et -1 si on est dans l'orientation inverse (retour du flot).

Corrigé 
graph *residual_graph(graph *g)
{
    graph *gr = graph_create(g->nvertices, g->s, g->t);
    for(int i = 0; i < g->nvertices; i++)
    {
        edge_list *l = g->edge_lists[i];
        while(l != NULL)
        {
            if (l->flow > 0)
            {
                edge_list *e = add_edge(gr, l->dst, i, l->flow);
                e->ref = l;
                e->flow = -1;
            }
            if (l->flow < l->capacity)
            {
                edge_list *e = add_edge(gr, i, l->dst, l->capacity-l->flow);
                e->ref = l;
                e->flow = 1;
            }
            l = l->next;
        }
    }
    return gr;
}

Question 19. Écrire une fonction de prototype

void adjust_flow_along_path(edge_list **pred, int t);

qui, étant donné un tableau pred issu d'un DFS depuis la source du graphe de flot et le sommet t destination visité, va remonter le chemin augmentant pour ajuster le flot dans le graphe initial.

Il faudra deux parcours ici pour obtenir le minimum puis pour ajuster le flot.

Corrigé 
void adjust_flow_along_path(edge_list **pred, int t)
{
    int min_flow = INT_MAX;
    int y = t;
    while(pred[y] != NULL)
    {
        if(min_flow > pred[y]->capacity)
            min_flow = pred[y]->capacity;
        y = pred[y]->src;
    }
    y = t;
    while(pred[y] != NULL)
    {
        pred[y]->ref->flow += pred[y]->flow * min_flow;
        y = pred[y]->src;
    }
}

La remontée s'arrête en g->s, seul sommet visité dont le prédécesseur vaut NULL.

Question 20. Écrire une fonction de prototype

void maximal_flow(graph *g);

qui calcule un flot maximal depuis un flot nul.

Corrigé 
void maximal_flow(graph *g)
{
    bool has_augmented = true;
    while(has_augmented)
    {
        has_augmented = false;
        graph *gr = residual_graph(g);
        edge_list **pred = dfs_launch(gr);
        if (pred[g->t] != NULL)
        {
            has_augmented = true;
            adjust_flow_along_path(pred, g->t);
        }
        free(pred);
        graph_free(gr);
    }
}

Chaque itération coûte \(O(|S| + |A|)\) et augmente la valeur du flot d'au moins 1 : on retrouve la complexité \(O(|f^*|(|S| + |A|))\) du cours.

Question 21. Testez votre fonction sur le réseau suivant, où la source \(s\) est le sommet 0 et le puits \(t\) le sommet 5 (les étiquettes sont les capacités) :

int main(void)
{
    graph *g = graph_create(6, 0, 5);
    add_edge(g, 0, 1, 20);
    add_edge(g, 0, 2, 10);
    add_edge(g, 1, 2, 10);
    add_edge(g, 1, 3, 5);
    add_edge(g, 3, 2, 15);
    add_edge(g, 2, 4, 10);
    add_edge(g, 4, 3, 10);
    add_edge(g, 3, 5, 15);
    add_edge(g, 4, 5, 20);
    maximal_flow(g);
    graph_print(g);

    graph_free(g);
}
Corrigé

On obtient un flot de valeur 15, par exemple celui de la figure du cours :

0 ---0/10---> 2
0 ---15/20---> 1
1 ---5/5---> 3
1 ---10/10---> 2
2 ---10/10---> 4
3 ---10/15---> 5
3 ---0/15---> 2
4 ---5/20---> 5
4 ---5/10---> 3

Le flot obtenu dépend de l'ordre des listes d'adjacence, mais sa valeur est toujours 15, la capacité de la coupe \(\{0, 1, 2\}\).

Question 22. Assurez-vous de n'avoir aucune fuite mémoire.

Corrigé

On peut compiler avec gcc -Wall -Wextra -fsanitize=address -g : à la fin du programme, l'outil signale les blocs alloués non libérés. On peut aussi utiliser valgrind --leak-check=full ./a.out sous Linux.

Question 23. Appliquez cela au cas d'un graphe biparti pour retrouver la méthode vue précédemment.

Corrigé

On construit le réseau du cours : les sommets \(0, \dots, n-1\) du graphe biparti, une source \(s = n\) et un puits \(t = n + 1\), des arêtes de capacité 1 de \(s\) vers chaque source, de chaque but vers \(t\), et de la source vers le but pour chaque arête du graphe. La valeur d'un flot maximal est la taille d'un couplage maximal, et les arêtes entre sources et buts portant un flot 1 forment ce couplage.

int max_matching_size(int nvertices, bool *sources, int *vtx_pairs, int npairs)
{
    int s = nvertices, t = nvertices + 1;
    graph *g = graph_create(nvertices + 2, s, t);
    for(int i = 0; i < nvertices; i++)
    {
        if (sources[i])
            add_edge(g, s, i, 1);
        else
            add_edge(g, i, t, 1);
    }
    for(int i = 0; i < npairs; i++)
    {
        int x = vtx_pairs[2*i];
        int y = vtx_pairs[2*i+1];
        if (sources[x])
            add_edge(g, x, y, 1);
        else
            add_edge(g, y, x, 1);
    }
    maximal_flow(g);
    int value = 0;
    for(edge_list *l = g->edge_lists[s]; l != NULL; l = l->next)
        value += l->flow;
    graph_free(g);
    return value;
}

Sur le graphe g2 précédent, on obtient 5. Le graphe résiduel de ce réseau est exactement celui de la partie précédente, complété par \(s\) et \(t\) : chaque itération de Ford-Fulkerson bascule un chemin augmentant du couplage.

Remarque

Les corrigés de ce TP ont été testés, avec -fsanitize=address,undefined, en comparant sur des milliers de petits graphes aléatoires la valeur du flot obtenu avec la capacité minimale d'une coupe calculée en énumérant toutes les coupes, et la taille des couplages obtenus avec celle d'un couplage maximal calculé par énumération.