Voici à quoi ressemble un fpl (fully packed loop) typique (tiré au sort avec équiprobabilité, via un algorithme coupling from the past). Les points rouges (resp. bleus) marquent les +1 (resp. -1) de la matrice à signes alternants correspondante.

En voici un autre, inscrit cette fois dans un carré de côté 128; les circuits ne sont pas dessinés, seuls apparaissent les chemins, qui relient deux points sur les bords du carré. Les cinq chemins les plus longs sont coloriés en rouge (chemin le plus long, de taille 1669), vert (deux chemins de tailles 1116 et 544) ou bleu.

Voici enfin une matrice à signes alternants de même taille (toujours tirée au sort avec équiprobabilité), dont chaque élément est colorié selon le modèle à six sommets. Les points rouges et bleus sont, comme précédemment, les +1 et les -1, tandis que les autres sont coloriés comme leurs frères des quatre coins.

Les calculs ont été effectués avec un module Julia.

Promotion and growth diagrams

Here is the picture associated with dual evacuation (with an svg version):

This time the skew tableau t is depicted by the top row, and next rows, downward, result from successive dual evacuations of elements of t, selected in decreasing order.

The local commutativity rules enable us to build the whole family of diagrams, from the top row and from the diagonal, where all partitions are equal to the outer shape of t. The first column depicts the tableau ε^*(t), and symmetry between rows and columns yields a 2d-proof that dual evacuation ε^* is an involution.

This dual version of the preceding post would be of little interest, if it was not a preamble to a much nicer formula involving promotion. This operator ∂ is defined like evacuation of the minimal element, except that we put in the cell freed by this evacuation, a new entry larger than any other. The following picture shows p successive applications of ∂ to a tableau t (bottom row), where p is the size of t, i.e., the number of entries:

(here is an svg version).

Because the local rules that allow us to build the picture from the bottom row and from one of the diagonals are always the same, we discover that the left half of the picture, bordered by the central column, is exactly the same as the picture associated to evacuation, published in the preceding post; and on the right half of the picture we recognize a dual evacuation, applied to the central column (beware that diagrams in this part of the picture are different from those located at the start of this post, because top row is different). Hence we get a very nice 2d-proof of the following formula, again due to Schützenberger:

∂ ^p = ε ε^*

where the operators are composed in the "right" order, i.e. first ε then ε^*. The tableau ε(t) appears as the central column on this picture.

"Jeu de taquin" and growth diagrams

Growth diagrams yield a 2d-view of transformations involved in "jeu de taquin", a creation of Marcel-Paul Schützenberger, see Promotion and Evacuation by Richard P. Stanley, 2009. Here is a 2d-proof that evacuation (as defined by Schützenberger) is an involution:

This picture is constructed using a personal little Sage module, about skew tableaux. The standard module sage.combinat.skew_tableau is, as usual for enumerative combinatorics, useless. Here is an svg version of the picture.

Each row depicts a skew tableau as a sequence of growing diagrams: each diagram depicts a partition, and differs from the preceeding (in the same row) by adding a single cell, depicted by a blue square. Likewise each column is made (downward) by increasing diagrams, where the red square denotes the added cell. Green squares denote cells which are "simultaneously blue and red" (I know, this way of mixing colors is definitely wrong).

Here the bottom row is a given skew tableau t, and successive rows, bottom-up, result from successive evacuations of elements of t, as defined by Schützenberger.

The key observation is the local commutativity rule, which is equivalent to the definition of a taquin slide: starting from a diagram A, we may go right, where is the diagram B, then go down, where is D; or we may first go down, where is C, and then right, where is D. If blue and red squares involved are independent (i.e. they don't form a domino), these "two steps" walks differ by the order in which cells are added. Otherwise these walks are identical, D has an unique green square, that is a "free cell" during a "jeu de taquin" slide.

These rules allow to reconstruct the whole family of diagrams, from the initial tableau t, depicted by the bottom row, and from the diagonal, where all partitions are equal to the inner shape of t. The last column depicts the tableau ε(t) where ε denotes the evacuation operator, and because local commutativity rules are invariant by exchanging rows and columns, the bottom row depicts as well the tableau that results of applying ε to the last column. This is the expected 2d-proof that ε is an involution.

For the record let's look at the two bottom rows of the picture, to analyze their relationship with "jeu de taquin" sliding. The skew tableau t (bottom row) is the following left one, and the evacuation of entry 1 yield the next four tableaux, where a question mark ? denotes the "free cell" during sliding of the entries 1, 2, 5, 9. The last of these tableaux, on the right, results from this sliding:

* * 1 3 10    * * ? 3 10    * * 2 3 10    * * 2 3 10    * * 2 3 10
* * 2 5 11    * * 2 5 11    * * ? 5 11    * * 5 ? 11    * * 5 9 11
* 4 7 9       * 4 7 9       * 4 7 9       * 4 7 9       * 4 7 ?
6 8 12        6 8 12        6 8 12        6 8 12        6 8 12

Compare with the two bottom rows of the main picture, reproduced here for convenience:

In the bottom row, diagrams 1, 2, 5, 9 (with numbering starting at 0) have a green cell, in the same position as the question mark identifying free cell in "jeu de taquin"; other diagrams have a red cell, in the same position as the preceeding green one. It should be clear that on the other hand, this description of green and red cells results precisely from the local commutativity rules, and yields the top row among these two.

Note: jeu de taquin, promotion, evacuation and growth diagrams make sense for any finite poset, see Schützenberger and the cited article by Richard Stanley. This post is about the seminal case, where the underlying poset is N x N, and finite ideals are partitions, which form the Young lattice.

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.

Premières images de l'identité

Avant de voir comment effectuer rapidement les calculs avec Sage, voici tout de suite quelques résultats de petite taille (cliquer sur une image pour afficher dans un nouvel onglet la version originale en haute définition):

Grille 50x50

Grille 100x100

Grille 200x200

Pour les grilles suivantes, de plus grande taille, les codes couleurs sont: bleu = 0, vert = 1, jaune = 2, marron = 3. On ne voit que le quart nord-ouest de la grille, sinon les images sont de trop grande taille (et les autres quarts sont identiques par symétrie).

Grille 400x400

Grille 512x512

Grille 1000x1000

Sage et les tas de sable, speed dating

En ce moment je cherche à comprendre la structure de l'élément neutre (identity) du groupe du tas de sable, d'abord dans le cas standard d'une grille rectangulaire. Je calcule des exemples avec Sage, qui dispose depuis la version 4.7.1 d'un module sandpile écrit par David Perkinson (des versions préliminaires de ce module étaient disponibles depuis qq années en téléchargement séparé).

Au début étaient les graphes. Qu'est-ce qu'un graphe pour Sage ? Bien sûr il faut commencer par regarder la documentation (manuel de référence de Sage, section Graph Theory), mais ensuite j'aime bien savoir qui est qui, et pour cela il suffit d'interroger Python, qui a des capacités exceptionnelles d'introspection:

sage: sage.graphs
<module 'sage.graphs' from '[...]/sage/graphs/__init__.pyc'>

(je remplace par [...] toute la partie superflue du chemin, qui dépend de l'installation). Merci Python, la réponse indique que sage.graphs est un dossier, considéré comme un module car il comporte un fichier __init__.py (et sa version 'compilée' .pyc). Continuons à interroger Python, qui possède une fonction intégrée (builtin) dir très intelligente (on peut l'appliquer à n'importe quel objet):

sage: dir(sage.graphs)
['bipartite_graph', 'chrompoly', 'digraph', 'generic_graph', 'graph', 'graph_coloring', 'graph_database', etc.]

(j'ai largement élagué la liste fournie en réponse par dir et le 'etc.' est évidemment de moi). Sage possède un mécanisme célèbre pour obtenir une variante interactive de cette liste: il suffit d'appuyer sur la touche de tabulation après avoir entré sage.graphs. (y compris le point final).

sage: sage.graphs.digraph
<module 'sage.graphs.digraph' from '[...]/sage/graphs/digraph.pyc'>

Cette fois ce n'est plus un dossier, mais un module Python ordinaire. La fonction dir (ou la touche de tabulation sous Sage) permet à nouveau d'explorer cet objet:

sage: dir(sage.graphs.digraph)
['DiGraph', 'GenericGraph', etc.]

sage: DiGraph
<class 'sage.graphs.digraph.DiGraph'>

Ici apparaît une particularité (contestable) de Sage, DiGraph est un nom connu, avec des milliers et des milliers d'autres, importés automatiquement dès qu'on lance une session.

sage: dir (DiGraph)
['add_cycle', 'add_edge', 'add_edges', 'add_path','add_vertex', 'add_vertices', 'adjacency_matrix', 'breadth_first_search', 'distance', 'is_planar', 'is_regular', 'is_strongly_connected', 'is_subgraph', 'is_tree', 'kirchhoff_matrix', 'laplacian_matrix', 'latex_options', 'shortest_path', etc.]

Cette fois la liste originale comporte 329 attributs, dont la plupart sont des méthodes applicables à un graphe dirigé (c'est-à-dire à une instance de la classe DiGraph); j'en ai sélectionné une vingtaine à titre d'exemples.

Tout ceci est magnifique, sauf pour celui qui se pose la question stupide: ai-je besoin de cette usine à gaz et, au fond, comment peut-on représenter un graphe ?

Python possède une structure de données parfaitement adaptée, le dictionnaire, qui permet de définir très simplement une fonction (au sens mathématique, mapping) sans aucune restriction, ou presque, sur les objets du domaine. Revenons à notre objectif initial, étudier les tas de sable définis sur une grille rectangulaire:

1. Quels sont les sommets naturels de cette grille ? Les couples d'entiers.
2. Il faut leur ajouter, c'est la particularité du modèle, un sommet spécial, le puits, relié aux sommets du bord de la grille. Comment le coder ? (-1,-1) est possible, mais les sommets du graphe, c'est-à-dire les éléments du domaine du dictionnaire (en jargon, les clefs) ne sont pas forcément de même type. On peut donc coder le puits par un simple entier, ou par une chaîne de caractères, si, si ! Au travail:

def grille(m,n):
    g = {}
    # corners first
    g[1,1] = [(1,2), (2,1), 'sink', 'sink']
    g[m,1] = [(m-1,1), (m,2), 'sink', 'sink']
    g[1,n] = [(1,n-1), (2,n), 'sink', 'sink']
    g[m,n] = [(m-1,n), (m,n-1), 'sink', 'sink']
    g['sink'] = [(1,1), (1,1), (m,1), (m,1),\
                 (1,n), (1,n), (m,n), (m,n)]

La première instruction indique que g est un dictionnaire. La notation g[i,j] est une abréviation pour g[(i,j)], la syntaxe Python est intelligente, dès qu'il n'y a pas d'ambiguïté. Ne pas croire que g serait un tableau à deux dimensions (array), au sens informatique: g est le codage par un dictionnaire d'un tel tableau, et ce codage permet, entre autres, d'avoir un sommet alien (une chaîne de caractères).

On l'aura compris, g associe à chaque sommet la liste de ses voisins, et les coins de la grille sont reliés deux fois au puits. Noter que ce graphe non orienté est représenté par un graphe orienté: si v apparaît dans la liste des voisins de u, alors u apparaît dans la liste des voisins de v; la taille de l'espace alloué au graphe en mémoire double, c'est le prix à payer pour la simplicité de la représentation. Noter aussi que d'habitude aucune arête n'est issue du puits, ce qui est dommage pour le burning algorithm. Continuons avec un bord:

    # top edge
    for j in range(2,n):
        g[1,j] = [(1,j-1), (1,j+1), (2,j), 'sink']
        g['sink'].append((1,j))

Sautons le code similaire pour les trois autres bords, et traitons les sommets intérieurs de la grille:

    # inner vertices
    for i in range(2,m):
        for j in range(2,n):
            g[i,j] =[(i-1,j), (i+1,j), (i,j-1), (i,j+1)]
    return g

C'est un peu long à écrire, mais très simple et efficace. Bien que ce soit du pur Python, avec à la fin une double boucle qui souffre de la lenteur de son exécution par un interpréteur, une grille 200x200 est construite en moins d'une seconde, et une grille 1000x1000 en 35 secondes, sur mon ordinateur équipé d'un Intel Core i3 530. Par contre:

sage: time sage.sandpiles.sandpile.grid_sandpile(200,200)
Digraph on 40001 vertices
Time: CPU 241.39 s, Wall: 242.85 s

Difficile à croire, non ? Et ce n'est pas la faute de David Perkinson (ou très peu :-), son module utilise logiquement la classe DiGraph, et c'est elle qui est d'une lenteur catastrophique, bien qu'un dictionnaire du type de celui construit ci-dessus, fasse évidemment partie des formats d'entrée acceptables pour construire une instance de DiGraph:

sage: time DiGraph(grille(200,200))
Multi-digraph on 40001 vertices
Time: CPU 224.88 s, Wall: 225.76 s

Moralité, il n'existe pas d'usine à gaz gratuite !

Mea culpa, novembre 2014: ces 225 secondes nécessaires pour transformer un dictionnaire en graphe ne révèlent pas que les modules Sage sur les graphes entreraient dans la catégorie usine à gaz, comme je l'ai cru en 2011. La cause de ce temps de calcul prohibitif est beaucoup plus élémentaire: deux lignes de code défectueuses dans le module Sage, découvertes et corrigées en qq minutes par Nathann Cohen, un dimanche matin vers 6h!