Stabilization needs C

Programmer la stabilisation d'un tas de sable sur une grille m x n, c'est très facile en C ; on commence par créer une matrice avec deux lignes et deux colonnes de plus, et on initialise la matrice à zéro par précaution:

int** c_grid_new (int m, int n) {
    int m2 = m + 2, n2 = n + 2;
    int **g = malloc(m2 * sizeof(int *));
    int i, j;

    for (i = 0; i < m2; i++) {
        g[i] = malloc(n2 * sizeof(int *));
        for (j = 0; j < n2; j++)
            g[i][j] = 0;
    }
    return g;
}

Voici alors le code pour parcourir la matrice, ébouler tous les sommets instables, et retourner le nombre d'éboulements:

long c_grid_fire_unstable (int **c, int m, int n) {
    /*
     * c est une matrice (m+2)x(n+2)
     */
    int i, j;
    long u = 0;

    for (i = 1; i <= m; i++)
        for (j = 1; j <= n; j++)
            if (c[i][j] > 3) {
                u += 1;
                c[i][j] -= 4;
                c[i-1][j] += 1;
                c[i+1][j] += 1;
                c[i][j-1] += 1;
                c[i][j+1] += 1;
            }
    return u;
}

Noter qu'un sommet peut devenir instable pendant l'exécution de cette fonction, suite à l'éboulement de l'un de ses prédécesseurs, et être alors éboulé; autrement dit le nombre d'éboulements déclenchés par cette fonction (et retourné par elle) peut être supérieur au nombre de sommets initialement instables. Heureusement ceci n'a aucune influence sur le résultat de la stabilisation, car l'ordre des éboulements est sans importance:

long c_grid_stabilize (int **c, int m, int n) {
    int u;
    long f = 0;
    while (1) {
        u = c_grid_fire_unstable (c, m, n);
        if (u == 0)
            return f;
        f += u;
    }
}

Il ne reste enfin qu'à libérer la mémoire allouée pour créer la matrice:

void c_grid_free (int **g, int m, int n) {
    int m2 = m + 2, n2 = n + 2, i;

    for (i = 0; i < m2; i++)
        free (g[i]);
    free(g);
}

On peut ajouter une fonction main écrite avec les conventions habituelles, pour faire quelques tests et vérifier ainsi qu'on n'a pas laissé traîner de bug stupide. La seule question qui peut sembler délicate est l'utilisation de ces fonctions C depuis Sage, mais en fait c'est très simple. On commence par insérer quelques déclarations magiques Cython:

cdef extern from "libstab.c":
    int** c_grid_new (int m, int n)
    void c_grid_free (int **g, int m, int n)
    long c_grid_stabilize (int **c, int m, int n)

avant la définition de la classe qui va utiliser ces fonctions (libstab.c est évidemment le nom que j'ai choisi pour le fichier qui contient les fonctions C). Après quoi on peut écrire la méthode de stabilisation d'une grille comme suit:

class ConfigGrid (Configuration):
    [...]        
    def stabilize(self):
        cdef int i, j
        cdef int m = self.sandpile.nbrows, n = self.sandpile.nbcols
        cdef int **c = c_grid_new (m, n)

        for i in range (1, m+1):
            for j in range (1, n+1):
                c[i][j] = self[i,j]
        f = c_grid_stabilize (c, m, n)
        for i in range (1, m+1):
            for j in range (1, n+1):
                self[i,j] = c[i][j]
        c_grid_free (c, m, n) 
        return f

Cython permet de mélanger des objets C (définis par cdef) et des objets Python; en particulier l'initialisation de la matrice C avec la configuration Python (qui est essentiellement un dictionnaire) est effectuée miraculeusement par l'instruction:

c[i][j] = self[i,j]

A gauche on a un élément d'un tableau de tableaux C, i et j sont des indices, tandis qu'à droite on un élément d'un dictionnaire Python, i,j est une clef ! On opère de même pour récupérer le résultat du calcul.

La suite est en principe compliquée, puisqu'il faut convertir le fichier jbspc.pyx (dans lequel figurent toutes les définitions de classes et de méthodes, y compris les quelques lignes écrites en Cython) en un fichier C, puis compiler le résultat sans oublier d'indiquer à gcc d'utiliser le fichier libstab.o et les bibliothèques Sage standard (même si pour le moment nous ne les utilisons pas). Grâce à Dieu, incarné par William Stein et ses apôtres, l'instruction load de Sage effectue tout le travail automatiquement:

sage: load /home/betrema/Eclipse/Sage/jbspc.pyx
Compiling /home/betrema/Eclipse/Sage/jbspc.pyx...
sage: c=ConfigGrid(200,200,4)
sage: time c.stabilize()
95629488
Time: CPU 1.49 s, Wall: 1.50 s

Ouf, l'ordinateur est enfin utilisé efficacement ! En résumé nous avons utilisé trois outils pour effectuer les calculs de stabilisation d'un tas de sable sur une grille :

1. module standard sandpile distribué avec Sage ;
2. module Python pur ;
3. module Cython avec code C.

En passant du module 1 au module 2 on gagne un facteur 100 en efficacité ; essentiellement parce que le module standard définit un tas de sable comme sous-classe d'un graphe, ce qui est logique mais se révèle très lent. Il faut aussi éviter le piège des entiers Sage, ici on a besoin des entiers machine (int) et non du type Integer.

En passant du module 2 au module 3 on gagne un facteur 300 : je n'ai pas donné les chiffres exacts, mais stabiliser une grille 200x200 avec le module 2 demande environ 450 secondes sur ma machine.

Au total on a donc diminué le temps de calcul d'un facteur 30000 (trente mille), ce qui me rend un peu perplexe sur la politique de "review" des modules Sage (voir le ticket 10618), même si le code C du module 3 ne traite que les grilles. Nous allons pouvoir désormais examiner les résultats des calculs pour des configurations de tailles raisonnables, il ne reste qu'à utiliser les modules de dessin (plot) de Sage, j'en parlerai bientôt.

Note. Pour voir le programme C qui résulte de la compilation du fichier Cython, entrer la commande :

sage -cython jbspc.pyx

et le résultat se trouve dans jbspc.c .