Dessiner un tas de sable (plot)

Pour afficher un tas de sable avec Sage il faut définir une méthode plot(), comme dans le module standard de David Perkinson. Ici je ne m'intéresse qu'aux tas de sable sur une grille rectangulaire, et je vais donc écrire une méthode pour la classe ConfigGrid :

from sage.plot.colors import Color
from sage.plot.plot import Graphics
from sage.plot.line import line
from sage.plot.text import text
from sage.plot.scatter_plot import scatter_plot
[...]
class ConfigGrid (Configuration):
    def __init__ (self, m, n, sable = 0):
        g = SandpileGrid (m, n)
        Configuration.__init__ (self, g, sable)
        
    def stabilize(self):
        [...]

    def plot (self, palette=0, cell=None, composition='', **kwds):

Le premier argument, palette, est probablement inattendu, mais j'aime bien pouvoir varier les palettes de couleurs, j'en reparlerai. Je souhaite afficher une configuration comme un ensemble de cellules, rondes a priori, avec pour chacune le nombre de grains de sable, indiqué à la fois par une couleur et par un chiffre (ou plusieurs). Le second argument, cell, désigne la taille de chaque cellule, et le troisième, composition, sert à spécifier des options, j'en reparlerai. Enfin la méthode admet des arguments complémentaires en nombre indéfini, passés sous la forme keyword=value, selon un protocole Python simple, agréable, et classique.

La méthode plot() retourne un objet graphique, appelé traditionnellement G, et constitué de primitives graphiques, qu'on ajoute une à une à G. Ici la primitive essentielle est fournie par la fonction scatter_plot() :

    def plot (self, palette=0, cell=None, composition='', **kwds):
        m = self.sandpile.nbrows
        n = self.sandpile.nbcols
        
        # ms = marker size
        dmax = max (m, n)
        if dmax < 10: ms = 500
        else:  ms = 4000 // dmax
        # l'argument l'emporte, s'il est fourni
        if cell != None: ms = cell
            
        G = Graphics()
        p = [[x, m-y-1] for y in range(m) for x in range(n)]
        c = [couleur(self[i+1, j+1]) for i in range(m) for j in range(n)]
        ec = 'black'    # edgecolor
        # scatter_plot : zorder = 5 par défaut
        # on peut modifier les options 'marker' et 'alpha'
        # exemple: plot (palette=3, marker='s')
        s = scatter_plot (p, markersize = ms,
            facecolor = c, edgecolor = ec, **kwds)
        G += s

La fonction scatter_plot() prend en argument une liste de points, qui seront représentés par des marques (disque, carré, croix, etc); il faut prendre garde au piège habituel : la ligne i correspond à l'ordonnée y renversée (la ligne 1 est en haut, etc), tandis que la colonne j correspond à l'abscisse x. La fonction admet ensuite des arguments optionnels : taille des marques, couleur, et couleur de leur bord. La spécification des couleurs est très souple, et peut être une liste : ici chaque couleur dépend du nombre de grains de sable et de la palette, je ne donnerai pas le code de la fonction couleur (il est très simple, il faut juste faire attention à utiliser les palettes de manière circulaire).

La fonction scatter_plot() admet deux autres arguments optionnels, indiqués dans le code ci-dessus, et récupérés automatiquement s'ils sont fournis lors de l'appel à la méthode plot(), grâce à **kwds. De son côté elle transmet d'éventuels arguments inconnus d'elle à la méthode show(), qui affiche un objet graphique.

Tracer les lignes de la grille est très simple, on utilise la primitive graphique line() :

        for i in range (m):
            G += line ([(-0.5, i),(n-0.5, i)])
        for j in range (n):
            G += line ([(j, -0.5),(j, m-0.5)])

Ces lignes passent sous les cellules dessinées par scatter_plot(), comme par magie : en fait l'ordre de superposition des composantes d'un objet graphique est déterminé par le paramètre zorder, qui vaut 0 par défaut pour line(), et 5 pour scatter_plot().

Il ne reste plus qu'à ajouter la représentation décimale de chaque altitude, si on le souhaite, en utilisant la primitive graphique text() :

        fs = 8 + ms // 50
        G += sum (text (str(self[i+1, j+1]), (j, m-i-1),
                        fontsize = fs, zorder = 10)
                  for i in range (m) for j in range (n))

Comme d'habitude il faut écrire soigneusement la correspondance entre le couple (ligne, colonne) de la configuration, et le couple (x, y) sur le dessin. Le paramètre zorder spécifie que le texte est visible, car placé au-dessus des marques dessinées par scatter_plot().

La fin du code de plot() se comprend d'elle-même :

        G.xmin(-1)
        G.xmax(n)
        G.ymin(-1)
        G.ymax(m)
        G.axes(False)
        G.set_aspect_ratio(1)
        return G

Le code de la méthode ConfigGrid.plot() est maintenant presque complet. Je ne donne pas les palettes, ce sont de simples listes de couleurs, que je suis allé chercher sur le site Color Brewer — le graphisme est un métier.

Quand on utilise cette méthode, on souhaite rapidement disposer de pas mal d'options : afficher ou non les lignes, les chiffres, etc. Le premier réflexe est d'ajouter des paramètres un à un, mais on est rapidement débordé. J'ai choisi d'ajouter un seul paramètre composition, qui est une chaîne de caractères, chaque option étant repérée par une lettre, d'où par exemple le fragment de code suivant :

        digits = dmax < 80
        if 'd' in composition:
            digits = not digits

et le code donné un peu plus haut pour afficher les altitudes avec la primitive text() est gardé par :

        if digits:

Je peux ainsi basculer le choix par défaut (chiffres affichés pour les grilles de côtés inférieurs à 80). Voici quelques exemples de ce que je peux obtenir :

c=ConfigGrid (6, 8, 4); c.stabilize(); c.plot()

c=ConfigGrid (100, 100, 12)
c.stabilize()
c.plot(palette=2, figsize=12)

L'option 'q' permet d'afficher seulement le quart nord-ouest de la grille :

c.plot(palette=2, composition='q', figsize=12)

L'option 'b' supprime le bord de chaque cellule, l'option 'd' supprime l'affichage des chiffres, l'option 'l' supprime l'affichage des lignes de la grille, et 's' désigne une marque carrée (square) :

c.plot(palette=4, cell=20, composition='bdlq', marker='s')

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 .

Stabiliser le sable

Puisqu'il est impossible de mener des calculs sur des graphes de grande taille avec les modules standard de Sage, essayons un développement ascendant (bottom-up), en construisant un module ad hoc, simplifié à l'extrême. On est loin par contre de repartir de zéro, car on peut pêcher sans vergogne des fragments du module sandpile (dit standard dans la suite de cet article) de David Perkinson, merci à lui.

Commençons par une classe que j'appellerai SandpileGraph (pour éviter les collisions avec la classe Sandpile de Perkinson; par contre je garderai autant que possible les noms des attributs):

class SandpileGraph (dict):
    def __init__ (self, g, sink):
        dict.__init__ (self, g)
        self.sink = sink

Une instance de cette classe sera donc un dictionnaire, avec un attribut supplémentaire, le puits. L'argument g doit être un dictionnaire qui représente un graphe, et l'argument sink doit être l'un des sommets de ce graphe; inutile de vérifier systématiquement la cohérence de ces arguments, ce serait beaucoup de travail pour pas grand-chose.

Si S est une instance de SandpileGraph, parcourir les sommets et leurs voisins s'écrit très simplement:

    for u in S:          # Une clef du dictionnaire est un sommet
        for v in S[u]:   # La valeur associée à une clef est la
            [...]        # liste des voisins

et on ajoute tout de suite deux méthodes à la classe:

    def out_degree (self, v):
        return len (self[v])
    
    def boundary (self):
        return self[self.sink]

Le bord est défini comme la liste des voisins du puits; attention, ce n'est pas la convention habituelle, selon laquelle aucune arête ne doit sortir du puits.

Dans le modèle du tas de sable, une configuration est une fonction qui associe à chaque sommet du graphe un entier naturel, censé représenter le nombre de grains de sable déposés sur ce sommet. On peut voir aussi une configuration comme une surface discrète, chaque entier représentant une altitude (ou une hauteur). Une configuration sera évidemment codée à son tour par un dictionnaire:

class Configuration (dict):
    def __init__ (self, S, sand):
        dict.__init__ (self, sand)
        self.sandpile = S

L'argument S doit être une instance de la classe SandpileGraph, et sand doit être un dictionnaire avec les mêmes clefs que S, et des valeurs entières. Comme d'habitude, je ne vérifie pas ces assertions sur les arguments, ce serait fastidieux et l'utilisateur est quand même censé savoir un minimum ce qu'il fait.

Il est commode de prévoir le cas où l'argument sand est un simple entier n, pour construire une configuration de hauteur uniforme, je modifie donc légèrement la méthode __init__ , et j'en profite pour ajouter l'attribut fired : c'est un dictionnaire supplémentaire, qui codera un script de mise à feu:

    def __init__ (self, S, sand = 0):
        if isinstance (sand, (int, Integer)):
            sand = S._dict_uniform (sand)
        dict.__init__ (self,sand)
        self.sandpile = S
        self.reset_fired()
        
    def reset_fired (self):
        # vivement la syntaxe 2.7
        # self.fired = {v:0 for v in self}
        f = {}
        for v in self:
            f[v] = 0
        self.fired = f
        

Je ne donne pas le code (évident ou presque, on y reviendra) de la méthode _dict_uniform ; ne pas oublier que Sage possède deux types d'entiers, int et Integer, on y reviendra aussi.

Python 2.7 est disponible depuis juillet 2010, et ajoute aux traditionnelles "list comprehensions" les "set and dictionary comprehensions". Aujourd'hui (octobre 2010), Sage est distribué avec la version 2.6 de Python, voir le ticket #9958, donc je ne peux pas employer la syntaxe naturelle pour initialiser l'attribut fired.

Attention à ne pas définir fired comme une configuration (ie une instance de la classe Configuration), ce serait pourtant logique du point de vue mathématique, mais cette configuration aurait à nouveau un attribut fired, et ainsi de suite, d'où une chaîne de références infinie. Nous sommes prêts à entrer dans le vif du sujet:

    def fire_vertex (self, v):
        self.fired[v] +=1
        S = self.sandpile
        self[v] -= S.out_degree(v)
        for u in S[v]:
            # pour éviter un test on accepte de
            # briser self[S.sink] = 0
            # if u != S.sink:
            self[u] += 1

    def unstable (self):
        S = self.sandpile
        # garantit que le puits ne sera pas dans la liste
        self [S.sink] = 0
        return [v for v in self if self[v] >= S.out_degree(v)]

    def fire_unstable (self):
        u = self.unstable()
        if not u:       # liste vide
            return False
        for v in u:
            self.fire_vertex(v)
        return True

    def stabilize (self, reset = True):
        if reset:
            self.reset_fired()
        while self.fire_unstable():
            pass
        return sum (self.fired[v] for v in self)

Je choisis d'effectuer les transformations (ie mises à feu, stabilisation) sur place, c'est le plus simple et le plus efficace; en prime chacune peut alors retourner un booléen ou un compteur.

Pour passer aux tests sur des grilles, il reste seulement à définir les sous-classes :

class SandpileGrid (SandpileGraph):
    def __init__ (self, m, n):
        self.nbrows = m
        self.nbcols = n
        g = sandpile_grid (m, n)
        SandpileGraph.__init__ (self, g, 'sink')
        
class ConfigGrid (Configuration):
    def __init__ (self, m, n, sable = 0):
        g = SandpileGrid (m, n)
        Configuration.__init__ (self, g, sable)

La fonction sandpile_grid construit le dictionnaire qui définit une grille avec m lignes et n colonnes (voir l'essentiel du code dans le message du 13 octobre). Premier test:

sage: delta=ConfigGrid(100,100,4)
sage: time delta.stabilize()
6109740
Time: CPU 109.24 s, Wall: 109.69 s

Dans le cas général l'altitude de chaque sommet de la configuration δ est égale au degré du sommet; ici on voit que chacun des 10000 sommets a été mis à feu plus de 600 fois, en moyenne, durant la stabilisation. Le temps d'exécution n'est pas terrible (presque deux minutes), mais il est meilleur qu'avec le module standard, et de loin:

sage: time g=grid_sandpile(100,100)
Time: CPU 16.08 s, Wall: 16.35 s
sage: delta=g.all_k_config(4)
sage: time sdelta=delta.stabilize()
Time: CPU 3250.75 s, Wall: 3263.01 s

Déjà 16 secondes pour construire la grille (grid_sandpile est une fonction du module standard, à ne pas confondre avec ma version, que j'ai appelée … sandpile_grid), et presque une heure pour stabiliser δ, argh ! Noter que le module standard n'effectue pas les transformations sur place: ci-dessus sdelta est une nouvelle configuration; c'est plus intuitif, mais moins efficace, en temps comme en mémoire.

Sage réserve d'autres surprises, le module que j'ai écrit est pour l'instant du pur Python (je ne fais appel à aucun module Sage), mais Sage effectue un prétraitement (preprocessing) avant de passer une instruction à l'interprète Python, et convertit en particulier les entiers littéraux, comme 4, en Integer(4). Integer est le type des entiers nobles, aussi grands qu'on veut, et munis à cet effet de méthodes sophistiquées; un objet de ce type est en fait une enveloppe (wrapper) pour un entier de la bibliothèque MPIR; int est le type des entiers du pauvre, ou entiers machine, codés sur 32 ou 64 bits. Il se trouve qu'un tas de sable n'utilise jamais d'entier plus grand qu'un entier machine, et que la différence de temps d'exécution est spectaculaire:

sage: delta_int=ConfigGrid(100,100,int(4))
sage: time delta_int.stabilize()
6109740
Time: CPU 50.65 s, Wall: 50.97 s

Le temps d'exécution a été divisé par plus de deux ! Avec le module standard le temps d'exécution passe de 3250 secondes à 3217, ce qui est moins spectaculaire :-) Pour éviter ce genre de plaisanterie, la méthode _dict_uniform de la classe SandpileGraph doit donc être écrite comme suit:

    def _dict_uniform (self, n):
        d = {}
        # forcer la conversion Integer -> int
        n = int(n)
        for v in self:
            d[v] = n
        d[self.sink] = 0
        return d

Essayons ensuite d'accélérer encore la stabilisation avec une méthode ad hoc, qui fait tout elle-même et ne calcule pas l'attribut fired (script de mise à feu):

    def quick_stabilize (self):
        S = self.sandpile
        ff = 0
        while True:
            f = 0
            self[S.sink] = 0
            for v in self:
                d = S.out_degree (v)
                if self[v] >= d:
                    f += 1
                    self[v] -= d
                    for u in S[v]:
                        self[u] += 1
            if f == 0:
                return ff
            ff += f                    

sage: delta=ConfigGrid(100,100,4)
sage: time delta.quick_stabilize()
6109740
Time: CPU 29.95 s, Wall: 30.10 s

Il n'y a plus grand-chose à gratter dans cette voie, sauf une amélioration magique de 25% du temps de calcul en faisant une copie du module (il s'appelle jbsp.py) avec le suffixe pyx, qui désigne un module écrit en Cython. On n'attend aucune amélioration sans exploiter les particularités de Cython, mais la magie fonctionne:

sage: load /home/betrema/Eclipse/Sage/jbsp.pyx
Compiling /home/betrema/Eclipse/Sage/jbsp.pyx...
sage: delta=ConfigGrid(100,100,4)
sage: time delta.quick_stabilize()
6109740
Time: CPU 23.39 s, Wall: 23.50 s

Le nombre de sommets de la grille est une fonction quadratique du côté, et l'expérience semble indiquer que la complexité de l'algorithme de stabilisation de la configuration δ sur une grille (exprimée par le nombre de sommets mis à feu) serait quadratique en fonction du nombre de sommets. Donc en multipliant le côté de la grille par 2, on multiplie approximativement le temps de calcul par 16, soit environ 6 minutes pour une grille 200x200, et 4 heures pour une grille 500x500, avec la solution qu'on vient d'exposer.

Ce sera tout pour aujourd'hui, c'est déjà bien long, on verra une autre fois comment écrire l'algorithme de stabilisation directement en C (dans le cas d'une grille ce n'est pas difficile) et comment utiliser cette fonction C avec Sage.