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!