Aller au contenu

Sac à dos et algorithmes probabilistes

Ce TP en C accompagne le chapitre Algorithmes d'approximation, algorithmes probabilistes et le TD sur le sac à dos.

Implémentation du sac à dos

On va supposer que l'instance du sac à dos est décrite par :

  • int n le nombre d'éléments (\(n \ge 1\)) ;
  • int *w le tableau des poids (strictement positifs) ;
  • int *v le tableau des valeurs (strictement positives) ;
  • int W le poids max.

Cas fractionnaire

Question 1. Écrire une fonction de prototype

float *resout_sacados_frac(int n, int *w, int *v, int W);

qui alloue et renvoie un tableau contenant la solution gloutonne du sac à dos fractionnaire.

On pourra utiliser le tri à bulle suivant qui place dans la variable ordre un tableau des indices des objets dans l'ordre décroissant de l'efficacité.

int *ordre = malloc(sizeof(int) * n);
for(int i = 0; i < n; i++)
    ordre[i] = i;
for(int i = 0; i < n; i++)
{
    for(int j = 0; j < n-1; j++)
    {
        float ej = (float)v[ordre[j]]
            / (float)w[ordre[j]];
        float ejn = (float)v[ordre[j+1]]
            / (float)w[ordre[j+1]];
        if (ej < ejn) {
            int temp = ordre[j];
            ordre[j] = ordre[j+1];
            ordre[j+1] = temp;
        }
    }
}
Corrigé 
float *resout_sacados_frac(int n, int *w, int *v, int W)
{
    int *ordre = malloc(sizeof(int) * n);
    for(int i = 0; i < n; i++)
        ordre[i] = i;
    for(int i = 0; i < n; i++)
    {
        for(int j = 0; j < n-1; j++)
        {
            float ej = (float)v[ordre[j]]
                / (float)w[ordre[j]];
            float ejn = (float)v[ordre[j+1]]
                / (float)w[ordre[j+1]];
            if (ej < ejn) {
                int temp = ordre[j];
                ordre[j] = ordre[j+1];
                ordre[j+1] = temp;
            }
        }
    }

    float *f = malloc(sizeof(float) * n);
    for(int i = 0; i < n; i++)
        f[i] = 0.0;

    int o = 0;
    while(o < n && W >= w[ordre[o]]) {
        W -= w[ordre[o]];
        f[ordre[o]] = 1.0;
        o++;
    }

    if (o < n)
        f[ordre[o]] = (float)W/(float)w[ordre[o]];

    free(ordre);

    return f;
}

Programmation dynamique

Question 2. Écrire une fonction de prototype

bool *resout_sacados(int n, int *w, int *v, int W);

qui alloue et renvoie un tableau contenant la solution exacte par programmation dynamique sous la forme d'un tableau de \(n\) booléens indiquant si on a pris le \(i\)-ème élément. On utilisera le tableau \(A(i,p)\) du TD (plus petit poids d'un sous-ensemble des \(i\) premiers objets de valeur exactement \(p\)).

Corrigé

Le tableau \(A\) est stocké à plat : la case \((i, p)\) est à l'indice p + i*(n*V+1). La valeur W + 1 joue le rôle de \(\infty\) (seules les cases \(\le W\) nous intéressent). Le tableau used mémorise si l'objet \(i\) est utilisé dans le minimum, pour reconstruire la solution.

bool *resout_sacados(int n, int *w, int *v, int W)
{
    bool *f = malloc(sizeof(bool) * n);
    for(int i = 0; i < n; i++)
        f[i] = false;

    int V = v[0];
    for (int i = 1; i < n; i++)
        if(v[i] > V)
            V = v[i];

    int *tab = malloc(sizeof(int) * (n+1)*(1+n*V));
    bool *used = malloc(sizeof(bool) * (n+1)*(1+n*V));
    for(int i = 0; i < (n+1)*(n*V+1); i++)
    {
        used[i] = false;
        tab[i] = W + 1;
    }
    tab[0] = 0;
    for(int i = 1; i <= n; i++) {
        tab[i * (n*V+1)] = 0;
        for(int p = 1; p <= n * V; p++) {
            int val = tab[p + (i-1)*(n*V+1)];
            if (p >= v[i-1]) {
                int valw = w[i-1] + tab[(i-1)*(n*V+1)+p-v[i-1]];
                if (valw < val) {
                    used[p + i*(n*V+1)] = true;
                    val = valw;
                }
            }
            tab[p + i*(n*V+1)] = val;
        }
    }

    // plus grande valeur atteignable avec un poids <= W
    int p = n * V;
    while(p > 0 && tab[n*(n*V+1) + p] > W) p--;

    // reconstruction de la solution
    int i = n;
    while(p != 0) {
        if (used[p+i*(n*V+1)]) {
            f[i-1] = true;
            p = p - v[i-1];
        }
        i--;
    }

    free(tab);
    free(used);

    return f;
}

Schéma d'approximation (FPTAS)

Question 3. Écrire une fonction de prototype

bool *approx_sacados(int n, int *w, int *v, int W, float eps);

qui renvoie la solution de l'algorithme d'approximation du TD pour la précision \(\epsilon\).

Corrigé

On suppose, comme dans le TD, que chaque objet rentre seul dans le sac (sinon on l'élimine au préalable).

bool *approx_sacados(int n, int *w, int *v, int W, float eps)
{
    int V = 0;
    for (int i = 0; i < n; i++)
    {
        if (v[i] > V)
            V = v[i];
    }
    float factor = (V * eps) / n;
    int *va = malloc(sizeof(int) * n);
    for (int i = 0; i < n; i++)
    {
        va[i] = (int)(v[i] / factor);
    }
    bool *f = resout_sacados(n, w, va, W);
    free(va);
    return f;
}

Des tests

int n = 10;
int W = 842;
int w[] = {415,241,555,842,841,59,626,535,715,541};
int v[] = {515,341,655,942,941,159,726,635,815,641};

Qui doit renvoyer

// pour fractionnaire
{ 1.0, 1.0, 0.0, 0.0, 0.0,
  1.0, 0.0, 0.237, 0.0, 0.0 }
// pour la programmation dynamique
// objets 1, 5, 9 pour une valeur 1141
0 1 0 0 0 1 0 0 0 1

Pour \(\epsilon = 0{,}5\), l'approximation renvoie par exemple les objets 1, 5 et 7, de valeur \(1135 \ge 0{,}5 \times 1141\).

On pourra aussi tester resout_sacados sur de petites instances aléatoires en comparant avec une recherche exhaustive sur les \(2^n\) sous-ensembles, et vérifier que approx_sacados renvoie une solution admissible de valeur \(\ge (1 - \epsilon) OPT\).

On peut trouver une calculatrice en ligne ici pour vérifier les résultats : https://augustineaykara.github.io/Knapsack-Calculator/

Fonctionnement d'un générateur aléatoire

On va créer ici notre propre générateur aléatoire. Il s'agit d'une famille de suites, comme peuvent l'être les suites de Syracuse, qui dépend de la valeur de \(u_0\), appelée le germe du générateur.

Voici le générateur RANDU utilisé par IBM dans les années 1960, de très mauvaise qualité, qui a induit en erreur des centaines de scientifiques :

\[u_{n+1} \equiv (65539 u_n) \mod 2^{31}\]

avec \(u_0\) un nombre impair.

Note : on peut écrire x & 0x7fffffff pour calculer le modulo \(2^{31}\) d'un entier positif.

Ici, ce générateur ne renvoie que des nombres impairs. On en déduit un générateur de float dans \(]0;1[\) par division : \(f_n = \frac{u_n}{2^{31}}\).

Question 4. On utilisera une variable globale de type int pour contenir le dernier terme calculé de la suite.

Écrire des fonctions de prototype

void srand_ibm(int seed);
int rand_ibm();

qui initialise le germe (on rajoute 1 s'il est pair) pour srand_ibm et renvoie le prochain terme de la suite (et change donc la valeur de la variable globale).

Corrigé

Le produit 65539 * u dépasse la capacité d'un int : en C, un dépassement sur un entier signé est un comportement indéfini. On calcule donc en unsigned int, où les calculs se font modulo \(2^{32}\), avant de réduire modulo \(2^{31}\).

int u = 1;
void srand_ibm(int seed)
{
    u = seed + 1 - seed % 2;
}

int rand_ibm()
{
    u = (int)((65539u * (unsigned int)u) & 0x7fffffffu);
    return u;
}

On en déduit, grâce à la fonction ldexpf(f, n) de math.h qui calcule \(f \times 2^{n}\), une fonction randu renvoyant un flottant :

float randu()
{
    return ldexpf(rand_ibm(), -31);
}

On peut voir que ce générateur est mauvais en traçant les points \((f_k, f_{k+1}, f_{k+2})\) formés de trois valeurs successives : ils se répartissent sur 15 plans parallèles.

Triplets successifs (f_k, f_k+1, f_k+2) de randu : ils se répartissent sur 15 plans parallèles. Triplets successifs (f_k, f_k+1, f_k+2) de randu : ils se répartissent sur 15 plans parallèles.
Triplets successifs $(f_k, f_{k+1}, f_{k+2})$ de randu : ils se répartissent sur 15 plans parallèles.

Question 5. Quelle est la proportion de nombres \(< \frac{1}{2}\) renvoyée par ce générateur pour un grand nombre de tirages ?

Corrigé

Expérimentalement, on obtient une proportion très proche de \(\frac{1}{2}\) (par exemple \(0{,}5005\) sur \(10^6\) tirages depuis le germe 1) : ce test unidimensionnel ne détecte pas le défaut du générateur. Le défaut apparaît sur les triplets : comme \(65539 = 2^{16} + 3\), on a \(65539^2 = 2^{32} + 6 \cdot 2^{16} + 9 \equiv 6 \times 65539 - 9 \pmod{2^{31}}\), donc

\[u_{k+2} \equiv 6 u_{k+1} - 9 u_k \pmod{2^{31}}\]

Ainsi \(9 f_k - 6 f_{k+1} + f_{k+2}\) est toujours un entier, compris entre \(-5\) et \(9\) : d'où les 15 plans de la figure.

À partir de maintenant, on va utiliser srand et rand donnés par la bibliothèque standard du C.

À noter qu'il y a de meilleures fonctions disponibles avec srand48 et lrand48.

La fonction rand renvoie une valeur entre 0 et RAND_MAX, qui vaut sûrement \(2^{31}-1\) dans votre implémentation.

Question 6. Écrire une fonction int randint(int a, int b) qui renvoie un entier uniformément entre \(a\) et \(b-1\), ici, on supposera que \(0 \le a < b \le\) RAND_MAX.

Corrigé

La solution naïve a + rand() % (b - a) n'est pas exactement uniforme quand \(b - a\) ne divise pas RAND_MAX + 1 : les petits restes sont un peu plus fréquents. On rejette donc les tirages au-delà du plus grand multiple de \(b - a\) (moins un) inférieur ou égal à RAND_MAX. En notant \(N\) = RAND_MAX + 1 et \(k = b - a\), il y a \(r = N \bmod k\) valeurs rejetées, avec \(r \le k - 1\) et \(r \le N - k\) (au moins un bloc de \(k\) valeurs est conservé), donc \(r < \frac{N}{2}\) : chaque tirage est accepté avec probabilité \(> \frac{1}{2}\) et le nombre moyen de tirages est inférieur à 2. C'est un algorithme de Las Vegas.

int randint(int a, int b)
{
    int n = b - a;
    // (RAND_MAX % n + 1) % n vaut (RAND_MAX + 1) % n sans débordement
    int limite = RAND_MAX - (RAND_MAX % n + 1) % n;
    int r = rand();
    while (r > limite)
        r = rand();
    return a + r % n;
}

Karger contracte, Karger remplace

On va implémenter l'algorithme de Karger vu en classe.

Pour cela, on considère une représentation des graphes différente de celle qu'on a vue en classe.

struct edge {
    int x; // arête x -- y
    int y;
    int id; // identifiant de l'arête
};
typedef struct edge edge;

struct multigraph {
    edge *edges;
    int nedges;
};
typedef struct multigraph multigraph;

Un multigraphe est alors un tableau d'arêtes et on identifie un sommet à un entier. L'identifiant permettra de suivre l'arête après les différentes contractions. On pourra utiliser l'indice de l'arête dans le graphe de départ.

Question 7. Écrire une fonction de prototype

multigraph *contract(multigraph *g, edge e);

qui renvoie le multigraphe obtenu en contractant l'arête \(e\). On va allouer et créer un nouveau multigraphe.

On pourra choisir \(\min(x,y)\) pour le nouveau sommet.

Corrigé

On renomme le sommet \(\max(x, y)\) en \(\min(x, y)\) et on supprime les boucles ainsi créées (dont \(e\) elle-même).

multigraph *contract(multigraph *g, edge e)
{
    int garde = e.x < e.y ? e.x : e.y;
    int supprime = e.x < e.y ? e.y : e.x;
    multigraph *h = malloc(sizeof(multigraph));
    h->edges = malloc(sizeof(edge) * g->nedges);
    h->nedges = 0;
    for (int i = 0; i < g->nedges; i++) {
        edge f = g->edges[i];
        if (f.x == supprime) f.x = garde;
        if (f.y == supprime) f.y = garde;
        if (f.x != f.y) {
            h->edges[h->nedges] = f;
            h->nedges++;
        }
    }
    return h;
}

Question 8. Écrire une fonction de prototype

multigraph *karger_coupure(multigraph *g);

qui renvoie une coupure obtenue par une passe de l'algorithme de Karger.

Corrigé

On suppose le graphe connexe. Un multigraphe connexe sans boucle a exactement deux sommets si et seulement si toutes ses arêtes relient la même paire de sommets (avec au moins trois sommets, il y a au moins deux paires reliées), ce qui évite de maintenir le nombre de sommets.

bool deux_sommets(multigraph *g)
{
    for (int i = 1; i < g->nedges; i++) {
        edge f = g->edges[i];
        edge e = g->edges[0];
        if (!((f.x == e.x && f.y == e.y) || (f.x == e.y && f.y == e.x)))
            return false;
    }
    return true;
}

multigraph *karger_coupure(multigraph *g)
{
    // copie de g pour ne pas le libérer
    multigraph *h = malloc(sizeof(multigraph));
    h->nedges = g->nedges;
    h->edges = malloc(sizeof(edge) * g->nedges);
    for (int i = 0; i < g->nedges; i++)
        h->edges[i] = g->edges[i];
    while (!deux_sommets(h)) {
        edge e = h->edges[randint(0, h->nedges)];
        multigraph *h2 = contract(h, e);
        free(h->edges);
        free(h);
        h = h2;
    }
    return h;
}

void free_multigraph(multigraph *g)
{
    free(g->edges);
    free(g);
}

Les champs id des arêtes restantes donnent les arêtes de la coupure dans le graphe de départ. Chaque contraction est en \(O(m)\) et il y en a \(n - 2\) : une passe est en \(O(nm)\) (plus le test deux_sommets, lui aussi en \(O(m)\) par tour).

Question 9. Testez votre fonction sur le graphe exemple donné en classe. Estimez le nombre d'essais avant d'obtenir une coupure minimale.

On peut réaliser d'autres tests avec des recollements de deux graphes complets.

Corrigé

Avec deux graphes complets \(K_5\) (sommets \(0\) à \(4\) et \(5\) à \(9\)) reliés par les deux arêtes \(\{0, 5\}\) et \(\{1, 6\}\), la coupure minimale est formée de ces deux arêtes. Sur \(10\,000\) passes, on l'obtient environ \(2\,600\) fois, soit une probabilité d'environ \(0{,}26\) et donc en moyenne un peu moins de 4 essais (loi géométrique), bien mieux que la borne \(\frac{2}{n(n-1)} = \frac{1}{45}\) du cours, qui donnerait 45 essais.