dimanche 23 octobre 2011

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.

lundi 17 octobre 2011

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

jeudi 13 octobre 2011

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!