Certains d'entre-vous en ont entendu parler vendredi, je me suis décidé ce week-end à implémenter l'algorithme HashLife. C'est un algorithme subtil et intéressant pour calculer de nombreuses générations de l'automate cellulaire dit "jeu de la vie" (Conway's game of life). J'en avais déjà parlé sur progmod comme sujet d'atelier potentiel, mais vous aviez trouvé (sans doute à raison, je m'en rends compte maintenant) que c'est peut-être trop compliqué. Ensuite, et indépendamment, j'ai aussi essayé d'aider des gens implémentant cet algorithme et qui posaient des questions de performance, mais ce n'est pas agréable de triturer du code sans le comprendre car on ne connaît pas le principe d'ensemble.

J'ai donc écrit ce week-end une implémentation de HashLife en OCaml, en espérant la rendre claire et en faire, en même temps, une explication de l'algorithme. C'est une forme de literate programming en un sens, car le code est truffé de commentaire et est pensé pour se lire de haut en bas, sans avoir besoin de sauter d'un endroit à un autre comme souvent dans les programmes informatiques.

J'ai mis le code sur gitorious, et vous pouvez trouver en particulier l'unique fichier source ici:

https://gitorious.org/hashlife-in-literate-ocaml/hashlife-in-literate-ocaml/blobs/master/hashlife.ml

La coloration syntaxique est un peu foireuse (j'ai soumis un rapport de bug à ce sujet, on verra), et le code est assez long (un bon millier de lignes) car il y a pas mal de choses à expliquer.

J'ai bien aimé cet exercice de style, et j'espère avoir des retours. En particulier, je me demande s'il est vraiment possible de comprendre l'algorithme en lisant le code, ou si c'est une lubie qui paraît réaliste seulement à la personne qui écrit le code après avoir compris (ou pendant la compréhension).

J'ai joué le jeu : je connais le jeu de la vie, je l'ai implémenté plusieurs fois naïvement, mais je n'avais jamais entendu parler de HashLife. J'écris mes remarques au fur et à mesure que je lis le code, sans avoir peur d'avoir l'air bête. Du coup c'est pas forcément que c'est pas clair, juste que je me suis posé la question. Disons du literate code reading.

Je me suis posé pas mal de questions bêtes au début :

• « an innermost (n-2k)-by-(n-2k) square » la notation bizarre m'a obligé à vérifier que tu voulais bien dire un carré de côté n-2k.
• m1 ? What ? Ah, c'est vraiment un constante pour -1.
• ( copies the content of the subsquare of a map into a isolated map with just the right side length. ) <- J'ai du mal à voir ce que ça veut dire "right side length". Au début j'ai compris "on ne renvoie que le côté droit", mais visiblement le code fait autre chose. Je me suis retenu de voir comment cette fonction était utilisée. J'ai finalement compris que ça renvoit effectivement un sous-carré : lire le nom de la fonction et son code correctement, ça aide. Je me demande pourquoi on n'a pas min_pos_x ou min_pos_y, mais ça, j'imagine bien que je le comprendrai plus tard.
• Intéressant le () pour le test.
• 1 lsl rank fourbe.
• Sympathique le random pour le test.

Damn, les choses sérieuses commencent.

• Impression que c'est inconsistent : pour moi dimension 1 = 2^1 = 2, pas 4. Si c'est vraiment une erreur, c'est pas cool pour la compréhension. Je continue en considérant que tu voulais dire '2 individual cells'.
  • Each tree represent (sic) a square of side 2^n. We say say that its 'dimension' is n, or that it is a n-tree.
  • If the dimension is 1, it contains the actual boolean state of the 4 individual cells.
• J'ai essayé de comprendre les explications avec les carrés en ASCII en me prenant des exemples, mais j'ai du mal. Les carrés ont visiblement 64 cases, donc dimension 6 (2^6 = 64). Si je prends n = 4, ça fait un carré n+2. Maintenant, ce carré aurait quatre sous-arbres de dimension n+1=5, donc 32 cases. Mais comme 32 n'est pas un carré parfait, je ne vois pas comment on peut avoir un arbre de dimensions n+1. Soit, prenons un arbre de dimension n, ça fait 16, et ça marche (du coup le carré ligne 152 montre les sous-arbres n-2 dans mon modèle mental). Du coup le carré centre il est de dimension 4 non ? Vu qu'on a que 16 éléments. Mais bon y'a rien qui colle avec la suite, donc mon modèle mental doit être faux.
• En lisant le type de données plus bas je me rend compte que j'ai effectivement fait le mauvais choix, et qu'il fallait peut-être plutôt lire « Each tree represents a square of side 2^n+1. We say say that its 'dimension' is n, or that it is a n-tree. ».

Du coup c'est un peu chiant, je dois relire les explications en essayanat cette nouvelle possibilité. Bonne excuse pour aller au lit.

Mon premier message ici, pour commenter une implémentation sympa d'un algorithme sympa. D'abord, les explications qui accompagnent le texte et les choix de présentation du code : c'est une aide précieuse pour comprendre ce qui se passe... mais les non-initiés risquent de négliger des détails qui se révèlent importantes par la suite. Ça, tu n'y peux rien.

• fonction life au tout début : pas convaincu par l'écriture de la règle de naissance/survie
• fonction dimension : vérifier la cohérence de la structure, pourquoi pas, c'est rigoureux.
• fonction tree_of_map : pourquoi ne pas ajouter de cases vides quand on a un tableau dont les dimensions ne sont pas une puissance de 2 ? Bon, ça va casser la jolie vérification de map_of_tree (tree_of_map foo) = foo aussi...

Pourquoi avoir gardé des booléens pour les feuilles ? Pourquoi pas un seul entier et des opérations bit à bit ? Est-ce que ça aurait un impact sur les perfs (et/ou sur le GC) ? Faudrait que je bidouille ça dans un coin.

En tout cas merci pour le partage et ce code qui donne à réfléchir.

Cygal > Merci beaucoup pour tes remarques. Il reste quelques incohérences (et de nombreuses fautes d'anglais) dans mon explication et il faut bien des gens pour débugguer :-'

Je suis ravi d'apprendre tout ce que vous trouvez délicat, subtil, "fourbe" etc. C'est clairement des endroits où il faut rajouter des commentaires ou clarifier le code.

Pour ton deuxième bloc de questions, l'incompréhension vient d'une confusion entre le côté du carré, de taille 2^n, et sa surface 2^(2n): un carré de dimension 1 a un côté de taille 2^1 = 2, mais donc contient bien 2*2 = 4 cellules individuelles.

(Je trouve que le nom "Cell" comme constructeur d'un arbre de dimension 1 (qui contient... 4 cellules) est un peu pourri mais je n'avais pas trouvé mieux. En fait "Base" serait peut-être mieux, je vais sans doute changer quand je serai convaincu que c'est mieux.)

La bonne nouvelle pour toi, c'est que tu n'es qu'au début, et donc tu peux te permettre de relire l'ensemble :-'

Cacophrène > merci pour tes remarques aussi.

Pourquoi avoir gardé des booléens pour les feuilles ? Pourquoi pas un seul entier et des opérations bit à bit ? Est-ce que ça aurait un impact sur les perfs (et/ou sur le GC) ? Faudrait que je bidouille ça dans un coin.

Ça doit avoir un impact, mais j'ai commencé à écrire du code naïf et lisible (je veux dire dont l'algorithmique soit correcte mais dont le code soit clair au risque d'avoir un fort facteur multiplicatif).
Tu verras qu'à la fin, pour l'implémentation mémoizée, j'ai utilisé une astuce du genre. Elle est re-traduite en quadruplet de booléens pour faire le calcul, mais ça n'a pas d'impact sur les performances puisque les résultats sur les 1-tree sont très vite mémorisés et donc récupérés dans le cache. C'est la deuxième explication: certaines choses qui seraient malines en temps sont rendues inutiles par la mémoisation.

pas convaincu par l'écriture de la règle de naissance/survie

Pourrais-tu être un peu plus explicite sur ce qui te dérange ?

(Globalement le code est assez "piéton", avec un usage intensif de l'indentation comme suggestive de position spatiale, et moins de factorisation qu'on pourrait le souhaiter. Ici on pourrait utiliser une liste par exemple et faire un fold, mais du coup on perdrait l'intuition spatiale, donc ça me tente moins.)

fonction tree_of_map : pourquoi ne pas ajouter de cases vides quand on a un tableau dont les dimensions ne sont pas une puissance de 2 ? Bon, ça va casser la jolie vérification de map_of_tree (tree_of_map foo) = foo aussi...

Ça ne casserait pas la vérification tant que je teste seulement sur des "foo" de la bonne taille. Je n'ai pas fait cette complémentation à la puissance de 2 la plus proche parce que je n'en avais pas besoin à ce endroit (en fait j'arrivais de la direction inverse: j'ai du code subtil sur les trees que je veux tester, donc je les transforme en map où je fais tourner un algo plus clairement correct). C'est implémenté plus bas dans la partie "exemples" parce que ça devient utile sur des exemples concrets (chiant de se forcer à complémenter soi-même). Je vais rajouter un commentaire pour prévenir.

Tu verras qu'à la fin, pour l'implémentation mémoizée, j'ai utilisé une astuce du genre. Elle est re-traduite en quadruplet de booléens pour faire le calcul, mais ça n'a pas d'impact sur les performances puisque les résultats sur les 1-tree sont très vite mémorisés et donc récupérés dans le cache. C'est la deuxième explication: certaines choses qui seraient malines en temps sont rendues inutiles par la mémoisation.

D'accord pour le gain en temps qui est compensé par la mise en cache des résultats. Il restera malgré tout un gain en mémoire si l'on remplace un entier par un n-uplet.

Pourrais-tu être un peu plus explicite sur ce qui te dérange ?

Bon, c'est juste de la présentation. Dans cette fonction life, la notation if count = 2 then c else if count = 3 then true else false ne me paraît pas très intuitive parce qu'on a l'habitude de distinguer les cas cellule vivante / cellule morte. Du coup je m'attendais plutôt à un truc du genre (c && (count = 2 || count = 3)) || (not c && count = 2).

C'est implémenté plus bas dans la partie "exemples" parce que ça devient utile sur des exemples concrets (chiant de se forcer à complémenter soi-même). Je vais rajouter un commentaire pour prévenir.

Exact, je n'avais pas fait attention.

En fait je me suis basé sur la description du sujet de Jean-Cristophe Filliâtre:

Si deux de ses huit voisines sont allumées, elle ne change pas d’état. Si trois sont allumées, elle s’allume. Sinon, elle s’éteint

Pour les quadruplets de booléens, ils ne sont pas stockés en mémoire dans la version mémoisée, qui associe un résultat à l'identifiant de l'arbre correspondant.

Exact, je n'avais pas fait attention.

Non, ce n'est pas une question d'inattention: je voudrais qu'on puisse lire le code de haut en bas, donc il ne faut pas que le lecteur soit étonné par quelque chose qui s'explique ensuite. J'ai donc rajouté un commentaire à ce sujet dans le code au moment où tu t'es posé la question.

Je me suis prêté au jeu également, et bien que j'ai eu à revenir en arrière pour clarifier certaines choses par moments, l'algorithme est bien expliqué et le code assez clair dans l'ensemble (je précise que je n'avais jamais entendu parler de l'algorithme avant). La plupart des confusions qu'on peut avoir sont clarifiées en relisant attentivement.

Concernant le code, les fonctions c (dx, dy) = map.(x+dx).(y+dy), ainsi que set et m gagneraient à avoir des noms plus descriptifs, une contraction de relative_position ou quelque chose du genre. C'est plus verbeux, mais ça serait beaucoup plus clairs.

The idea is that our current technique splits a computation of 2^n steps (on a n+2-tree) into two computations of 2^(n+1) steps (on n+1-trees). For s ≤ 2^n, we could instead split into two subcomputations of 'a' and 'b' steps, with s = a + b. But we need to ensure that both a and b are smaller than 2^(n-1), because it is not possible to compute more than 2^(n-1) steps in the future starting from n+1-trees.

Ce n'est pas plutôt "into two computations of 2^(n-1) steps" ? Aussi il serait bien de préciser "steps in the future". À ma première lecture, par manque de concentration peut-être, en lisant computation et steps, j'ai cru que tu faisais allusion à la complexité algorithmique :-°

Concernant la séparation du nombre de steps ahead en a + b, les explications gagneraient à être un peu étendues et clarifiées. Lorsqu'on lit les explications sur le cas général, on a l'impression que le 2^(n-1) est une sorte de nombre spécial et qu'on doit forcément calculer la 2^(n-1)-ième génération des (n+1)-trees avant de passer à l'étape suivante de l'algorithme, qui calcule 2^(n-1) générations dans le futur aussi. Alors qu'en fait il s'agit simplement d'avoir une somme permettant d'atteindre la 2^n-ième génération, qui peut être choisie autrement (si j'ai compris correctement). Je ne sais pas ce que je dis est clair, mais je pense que comprendre ça plus tôt permettrait de mieux comprendre le cas s ≤ 2^n.

Sur une note plus positive, je trouve que l'organisation spatiale que tu donnes au code aide beaucoup à la compréhension, et que globalement c'est très clair.

PS : Je n'ai pas encore vu la partie sur la mémoization, j'ajouterais des remarques dessus si j'en ai par la suite.

Tu as tout à fait raison pour "two computations of 2^(n+1) steps", c'est bien 2^(n-1) qu'il faut lire ici. C'est une erreur d'inattention de ma part, j'ai corrigé.

Pour ton second paragraphe sur le rôle privilégié de la somme 2^(n-1)+2^(n-1), je pense qu'on peut effectivement dire qu'elle joue un rôle privilégié. En effet, vu la structure de données on est forcé, de toute façon, de renvoyer un arbre de taille 2^n, donc 2^(n-1) est le futur "maximal" qu'on peut calculer sur un arbre de cette taille là. Si on demande de calculer un nombre plus petit d'itérations, on ne renvoie pas pour autant un arbre plus grand (alors qu'on a assez d'information pour), donc d'une certaine manière on perd de l'information. Ça se voit si je demande le calcul sur 10 itérations: la technique naïve de calcul sur les matrices renverra une matrice avec seulement 10 cases en moins sur les bords, alors que les algorithmes par quadtree renverront un arbre de côté moitié dans tous les cas.

Pour éviter cet effet qui privilégie 2^(n-1) il faudrait une structure plus flexible qui peut contenir non pas 4 enfants, mais autant qu'il est nécessaire, et à la granularité nécessaire, pour couvrir la zone totale pouvant être calculée. Ça ne m'a pas l'air facile.

Plus généralement cette partie sur le calcul d'un nombre d'étapes quelconques je ne l'ai pas lue quelque part, le sujet de Filliâtre dont je suis parti demande essentiellement de deviner ça tout seul. Du coup je ne sais pas si c'est la meilleure solution: il y a peut-être des façons de faire ça qui n'ont pas cet inconvénient.

Dans tous les cas, je ne suis pas convaincu que parler de ce problème plus tôt aide à la compréhension. Je trouve que la première version de la fonction "result", avec les sous-arbres qui dansent et tout, est déjà suffisamment difficile à comprendre et qu'ajouter des considérations sur le nombre d'étape à calculer risquerait d'augmenter encore la difficulté. C'est pour ça que j'ai préféré faire les choses en deux temps.

Citation de gasche

Pour ton second paragraphe sur le rôle privilégié de la somme 2^(n-1)+2^(n-1), je pense qu'on peut effectivement dire qu'elle joue un rôle privilégié. En effet, vu la structure de données on est forcé, de toute façon, de renvoyer un arbre de taille 2^n, donc 2^(n-1) est le futur "maximal" qu'on peut calculer sur un arbre de cette taille là. Si on demande de calculer un nombre plus petit d'itérations, on ne renvoie pas pour autant un arbre plus grand (alors qu'on a assez d'information pour), donc d'une certaine manière on perd de l'information. Ça se voit si je demande le calcul sur 10 itérations: la technique naïve de calcul sur les matrices renverra une matrice avec seulement 10 cases en moins sur les bords, alors que les algorithmes par quadtree renverront un arbre de degré moitié dans tous les cas.

Effectivement je ne le voyais pas comme ça, c'est plus clair maintenant. :)