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.
Aucun commentaire:
Enregistrer un commentaire