Random Apollonian Packing II : Simulations en deux dimensions Pierre Auclair

Nous avons vu dans un précédent billet comment dessiner la fractale du cercle d'Apollonius, un empilement de sphères déterministe en deux dimensions. Aujourd'hui nous allons découvrir ensemble la version aléatoire de cette fractrale, le Random Apollonian Packing.

Mise en place

Commençons par mettre en place tous les outils dont nous allons avoir besoin pour travailler. Au départ, nous avons une boite rectangulaire de hauteur et largeur finies.

import numpy as np

WIDTH, HEIGHT = 1, 1
coin = np.array([WIDTH, HEIGHT])

Notre boite est un rectangle à deux dimensions qui comprend tous les points tels que $x\in [0, \text{WIDTH}]$ et $y\in[0, \text{HEIGHT}]$. Par commodité, on a défini un vecteur coin qui représente le coin supérieur droit. Cette boite va contenir un ensemble de cercles déterminés par leurs centres et leurs rayons. On initialise des listes pour les contenir.

centres, rayons = np.array([]), np.array([])

Pour visualiser ce qui se passe dans cette boite, on définit une fonction qui utilise matplotlib. Pour mieux voir les insertions, le dernier cercle inséré sera dessiné en rouge et tous les autres en bleu.

import matplotlib.pyplot as plt

# Permet d'utiliser les polices LateX dans les figures
plt.rc("text", usetex=True)
plt.rc("font", family="serif")

def plot_configuration(ax, centres, rayons, color_last=False):
    """"Dessine la boite et les cercles qu'elle contient

    Parameters
    ==========
    - ax: pyplot subplot
        Axes utilisés pour le dessin
    - centres: numpy array
        Centres des cercles de dimensions (n_cercles, 2)
    - rayons: float
        Rayons des cercles de dimensions (n_cercles, )

    Returns
    =======
    None
    """

    if color_last:
        # Le dernier cercle insere sera rouge
        color = "red"
    else:
        color= "blue"

    for centre, rayon in zip(centres, rayons):        
        circle = plt.Circle(
            (centre[0], centre[1]), # Centre
            rayon, # Rayon
            fill=False, # Vide
            edgecolor=color, # Couleur du contour
            linewidth=2e-1, # Epaisseur du contour
        )
        color = "blue"
        ax.add_patch(circle)

    # Restriction a notre boite
    ax.set_xlim(0, WIDTH)
    ax.set_ylim(0, HEIGHT)

    # Que les cercles restent des cercles
    ax.set_aspect('equal', 'box')

Méthode 1 : naive

Le processus pour ajouter un cercle est simple :

On sélectionne un point au hasard uniformément dans la boite
On vérifie que ce point n'est pas contenu dans un cercle déjà tracé
Ce point devient le centre d'un cercle qui grandit jusqu'à toucher le bord de la boite ou un autre cercle

Et on répète cette procédure à l'infini. L'espace disponible dans notre boite diminue tandis que les rayons des cercles insérés vont devenir de plus en plus petits. La méthode naive la plus simple consiste à calculer la distance avec tous les autres cercles déjà insérés. En Python "naif", cela donne

def add_circle(centres, rayons):
    """"Ajoute un cercle dans la boite.

    Parameters
    ==========
    - centres: numpy array
        Centres des cercles de dimensions (n_cercles, 2)
    - rayons: numpy array
        Rayons des cercles de dimensions (n_cercles, )

    Returns
    =======
    tuple of numpy arrays
    """

    rayon_max = 0

    # Tant que l'on est dans un autre cercle, on retire un nouveau point
    while rayon_max <= 0:

        # Tir aleatoire d'un nouveau point
        nouv_centre = np.random.uniform(high=(WIDTH, HEIGHT))

        # Rayon maximum autorise pour rentrer dans la boite
        rayon_max = np.min([nouv_centre, coin - nouv_centre])

        # On calcule la distance aux autres cercles
        # Si rayon_max <= 0 alors notre point est compris dans un autre cercle
        for i in range(len(centres)):
            rayon_max = min(rayon_max, np.linalg.norm(nouv_centre - centres[i]) - rayons[i])

    # On insère le nouveau cercle
    # Les details ne sont pas importants
    if len(centres) == 0:
        centres = np.array([nouv_centre])
        rayons = np.array([rayon_max])
    else:
        centres = np.insert(centres, 0, nouv_centre, axis=0)
        rayons = np.insert(rayons, 0, rayon_max, axis=0)

    return centres, rayons

Visualisons tout ça, étape par étape

# Nouvelle figure
fig = plt.figure(figsize=(12 * WIDTH, 6 * HEIGHT))

centres, rayons = np.array([]), np.array([])
for i in range(10):
    # On prend un sous-axe par insertion
    ax = fig.add_subplot(2, 5, 1 + i)

    # On ajoute un cercle
    centres, rayons = add_circle(centres, rayons)

    # On dessine
    plot_configuration(ax, centres, rayons, color_last=True)
    ax.set_title(f"i = {1 + i}")

fig.tight_layout()
plt.show()

Et voila, on a déjà obtenu nos premières figures ! On peut augmenter le nombre de cercles jusqu'à 1000 pour voir ce qu'on obtient.

centres, rayons = np.array([]), np.array([])
for i in range(1000):
    centres, rayons = add_circle(centres, rayons)

fig, ax = plt.subplots(figsize=(6 * WIDTH, 6 * HEIGHT))

plot_configuration(ax, centres, rayons)

plt.show()

Magnifique, n'est-ce pas ? Cependant, cette méthode est très lente et peu optimisée. Il me faut plusieurs secondes pour 1000 cercles, et on peut faire beaucoup mieux...

Méthode 2 : accélération avec `numpy` et la vectorisation

Ce qui est très lent avec cette première méthode naive, c'est qu'elle utilise Python ! En particulier, lorqu'on utilise une boucle for pour calculer la distance avec tous les autres cercles. Quand on commence à connaitre Python, on apprend à éviter d'utiliser ces boucles qui sont très inefficaces et on apprend à utiliser les librairies écrites en C comme numpy et la vectorisation qui offrent un énorme gain en termes de performances.

Par exemple, on peut remplacer cette boucle for par un seul appel à la fonction np.linalg.norm. En effet cette fonction est vectorisable, c'est à dire qu'elle peut s'appliquer directement sur des tableaux de données, au lieu d'être appelée pour chacune des cases. Le gain est énorme pour un changement minimal, voyez plutôt

def add_circle_2(centres, rayons):
    """"Ajoute un cercle dans la boite.

    Parameters
    ==========
    - centres: numpy array
        Centres des cercles de dimensions (n_cercles, 2)
    - rayons: numpy array
        Rayons des cercles de dimensions (n_cercles, )

    Returns
    =======
    tuple of numpy arrays
    """

    rayon_max = 0

    # Tant que l'on est dans un autre cercle, on retire un nouveau point
    while rayon_max <= 0:

        # Tir aleatoire d'un nouveau point
        nouv_centre = np.random.uniform(high=(WIDTH, HEIGHT))

        # Rayon maximum pour rentrer dans la boite
        rayon_max = np.min([nouv_centre, coin - nouv_centre])

        # On s'assure qu'il y a d'autres cercles
        if len(centres) > 0:

            # Un seul appel a cette fonction, au lieu d'une boucle for
            rayon_max = min(
                    np.min(
                        np.linalg.norm(nouv_centre - centres, axis=1)
                        - rayons
                    ),
                    rayon_max,
            )


    # On insère le nouveau cercle
    if len(centres) == 0:
        centres = np.array([nouv_centre])
        rayons = np.array([rayon_max])
    else:
        centres = np.insert(centres, 0, nouv_centre, axis=0)
        rayons = np.insert(rayons, 0, rayon_max, axis=0)

    return centres, rayons

On n'a changé que quelques lignes. Comparons les performances de ces deux méthodes pour insérer 1000 cercles

import time

t0 = time.time()

# On reinitialise
centres, rayons = np.array([]), np.array([])
for i in range(1000):
    centres, rayons = add_circle(centres, rayons)

print(f"Methode 1 : {time.time() - t0} secondes")

t0 = time.time()

# On reinitialise
centres, rayons = np.array([]), np.array([])
for i in range(1000):
    centres, rayons = add_circle_2(centres, rayons)

print(f"Methode 2 : {time.time() - t0} secondes")

Methode 1 : 10.739489316940308 secondes
Methode 2 : 0.26102399826049805 secondes

Derrière son apparente simplicité d'utilisation, Python est en réalité un langage très subtile si on veut en tirer du code performant. Il y a bien souvent plusieurs manières d'effectuer une tache, mais avec une grande diversité de temps d'execution.

Dans un premier, il faut toujours faire attention à la manière dont on écrit des boucles. Python pénalise fortement l'utilisation des boucles for et les appels aux fonctions. Ça vaut toujours le coup de les remplacer par des appels groupés à des fonctions vectorisées, d'autant plus si elles font appels à des librairies C ou Fortran.

Avec la vectorisation, nous avons grandement gagné en termes de performances mais la complexité de notre algorithme reste la même. La $(n + 1)$-ième insertion nécessite $n$ comparaisons pour chaque tirage aléatoire, avec un taux de rejet qui augmente au fur et à mesure de la simulation. En effet, l'espace vide diminue rapidement et la probabilité de tirer un point en dehors des cercles déjà tracés diminue de même.

Pour optimiser la complexité de l'algorithme d'insertion, j'ai proposé d'utiliser une variation sur le concept des arbres KD.

Méthode 3 : Arbres de partionnement

L'idée de base, c'est de considérer un sous-ensemble de cercles, disons $\mathbb{A}$, au lieu de prendre tous ceux contenus dans notre boite. On peut construire une mini-boite $\mathbb{B}$ qui contient tous les cercles de $\mathbb{A}$, donc $\mathbb{A} \subset \mathbb{B}$. Ainsi, pour un point $\vec{x}$ quelconque, la distance de $\vec{x}$ à mes cercles peut être bornée inférieurement très facilement $$d(\vec{x}, \mathbb{A}) \geq d(\vec{x}, \mathbb{B})$$ Dit autrement, en calculant une seule distance entre $\vec{x}$ et ma boite $\mathbb{B}$, je peux décider que tous les cercles de $\mathbb{A}$ sont trop loins pour m'intéresser, sans avoir besoin de calculer la distance avec chacun de ces cercles. De sacrées économies !

On va donc partitionner notre boite et notre ensemble de cercles en un arbre de sous-boites et de sous-ensembles de cercles. L'intérêt d'une telle approche est que dans certains cas, il ne sera pas nécessaire de lister les cercles contenus dans une sous-boite pour savoir que ceux-ci sont trop loin. En théorie, on peut espérer que l'insertion se fasse en $\mathcal{O}(\ln n)$ au lieu de $\mathcal{O}(n)$ comme dans l'algorithme naif.

Par soucis de pédagogie, nous allons construire notre arbre directement en Python. C'est a priori une mauvaise idée, car Python va pénaliser ce type implémentation, et il serait plus malin d'utiliser une librairie pour ça. Mais ce sera plus clair comme ceci et on verra que le gain de cette méthode est déjà impressionnant.

La cellule de base de notre arbre est un objet avec un certain nombre de méthodes :

__init__ pour initialiser une sous-boite de taille donnée
dist_from_box pour calculer la distance d'un point à notre boite. C'est là tout l'intérêt de la procédure, si cette distance est trop grande, pas besoin de regarder les cercles que notre boite renferme !
find_collision va regarder le rayon maximal autorisé en considérant les cercles dans la boite courante et ses enfants
add_circle une méthode pour ajouter un nouveau cercle autorisé à notre structure en arbre
split qui permet de transformer une feuille en un arbre composé de deux branches
plot pour faire de la visualisation

THRESHOLD = 100


class Noeud:

    def __init__(self, x_min, x_max, y_min, y_max):

        self.centres = np.array([])
        self.rayons = np.array([])

        self.leaf = True

        self.x_min = x_min
        self.x_max = x_max
        self.y_min = y_min
        self.y_max = y_max

    def dist_from_box(self, x, y):

        # S'il est dans la boite, alors 0
        if (
            (self.x_min < x)
            and (x < self.x_max)
            and (self.y_min < y)
            and (y < self.y_max)
        ):
            distance = 0
        # S'il est au dessus ou en dessous
        elif (self.x_min < x) and (x < self.x_max):
            distance = min(abs(y - self.y_max), abs(self.y_min - y))
        # S'il est a gauche ou a droite
        elif (self.y_min < y) and (y < self.y_max):
            distance = min(abs(x - self.x_max), abs(self.x_min - x))
        else:
            distance = min(
                np.sqrt((self.x_min - x) ** 2 + (self.y_min - y) ** 2),
                np.sqrt((self.x_min - x) ** 2 + (self.y_max - y) ** 2),
                np.sqrt((self.x_max - x) ** 2 + (self.y_min - y) ** 2),
                np.sqrt((self.x_max - x) ** 2 + (self.y_max - y) ** 2),
            )
        return distance


    def find_collision(self, nouv_centre, rayon_max):

        if (self.centres.size > 0) and (rayon_max > 0):

            rayon_max = min(
                np.min(
                    np.linalg.norm(nouv_centre - self.centres, axis=1) - self.rayons
                ),
                rayon_max,
            )


        if (not self.leaf) and (rayon_max > 0):

            distance_1 = self.child_1.dist_from_box(
                nouv_centre[0],
                nouv_centre[1]
            )
            distance_2 = self.child_2.dist_from_box(
                nouv_centre[0],
                nouv_centre[1]
            )

            # On commence avec l'enfant le plus proche
            if distance_1 < distance_2:

                if rayon_max > distance_1:
                    rayon_max = min(
                        rayon_max, self.child_1.find_collision(nouv_centre, rayon_max)
                    )
                if rayon_max > distance_2:
                    rayon_max = min(
                        rayon_max, self.child_2.find_collision(nouv_centre, rayon_max)
                    )
            else:
                if rayon_max > distance_2:
                    rayon_max = min(
                        rayon_max, self.child_2.find_collision(nouv_centre, rayon_max)
                    )
                if rayon_max > distance_1:
                    rayon_max = min(
                        rayon_max, self.child_1.find_collision(nouv_centre, rayon_max)
                    )

        return rayon_max

    def add_circle(self, centre, rayon):

        if self.leaf:
            # On insère le nouveau cercle
            if len(self.centres) == 0:
                self.centres = np.array([centre])
                self.rayons = np.array([rayon])
            else:
                self.centres = np.concatenate((self.centres, [centre]))
                self.rayons = np.concatenate((self.rayons, [rayon]))
            if len(self.centres) > THRESHOLD:
                self.split()
        else:
            position = centre[self.split_direction]

            # Ajouter a gauche ou a droite
            if position + rayon < self.split_value_1:
                self.child_1.add_circle(centre, rayon)
            # Add to the upper/right part
            elif position - rayon > self.split_value_2:
                self.child_2.add_circle(centre, rayon)
            else:
                # On insère le nouveau cercle
                if len(self.centres) == 0:
                    self.centres = np.array([centre])
                    self.rayons = np.array([rayon])
                else:
                    self.centres = np.concatenate((self.centres, [centre]))
                    self.rayons = np.concatenate((self.rayons, [rayon]))

    def split(self):

        self.split_direction = np.argmax(
            [self.x_max - self.x_min, self.y_max - self.y_min]
        )

        left_most = np.sort(
            [
                self.centres[indice][self.split_direction] + self.rayons[indice]
                for indice in range(len(self.centres))
            ]
        )
        right_most = np.sort(
            [
                self.centres[indice][self.split_direction] - self.rayons[indice]
                for indice in range(len(self.centres))
            ]
        )[::-1]

        # split_index = np.argmax(right_most < left_most)

        self.split_value_1 = left_most[len(self.rayons) // 2]
        self.split_value_2 = right_most[len(self.rayons) // 2]

        # Split along the x-axis
        if self.split_direction == 0:
            # Left child
            self.child_1 = Noeud(
                self.x_min,
                self.split_value_1,
                self.y_min,
                self.y_max,
            )

            # Right child
            self.child_2 = Noeud(
                self.split_value_2,
                self.x_max,
                self.y_min,
                self.y_max,
            )
        # Split along the y-axis
        else:
            # Lower child
            self.child_1 = Noeud(
                self.x_min,
                self.x_max,
                self.y_min,
                self.split_value_1,
            )

            # Upper child
            self.child_2 = Noeud(
                self.x_min,
                self.x_max,
                self.split_value_2,
                self.y_max,
            )

        indices_1 = np.array(
            [
                indice
                for indice in range(len(self.centres))
                if self.centres[indice][self.split_direction] + self.rayons[indice]
                <= self.split_value_1
            ],
            dtype=np.int32,
        )
        indices_2 = np.array(
            [
                indice
                for indice in range(len(self.centres))
                if (self.centres[indice][self.split_direction] - self.rayons[indice]
                >= self.split_value_2)
                and (indice not in indices_1)
            ],
            dtype=np.int32,
        )
        indices_0 = np.array(
            [
                indice
                for indice in range(len(self.centres))
                if (indice not in indices_1) and (indice not in indices_2)
            ],
            dtype=np.int32,
        )

        self.leaf = False
        self.child_1.centres = self.centres[indices_1]
        self.child_2.centres = self.centres[indices_2]
        self.child_1.rayons = self.rayons[indices_1]
        self.child_2.rayons = self.rayons[indices_2]
        self.centres = self.centres[indices_0]
        self.rayons = self.rayons[indices_0]

    def plot(self, ax, recursive=False):

        ax.plot(
            [self.x_min, self.x_max, self.x_max, self.x_min, self.x_min],
            [self.y_min, self.y_min, self.y_max, self.y_max, self.y_min]
        )

        for centre, rayon in zip(self.centres, self.rayons):
            circle = plt.Circle(
                (centre[0], centre[1]),
                rayon,
                fill=False,
                edgecolor="blue",
                linewidth=1e-1,
            )
            ax.add_patch(circle)


        if recursive and not self.leaf:
            self.child_1.plot(ax, recursive=recursive)
            self.child_2.plot(ax, recursive=recursive)

        ax.set_xlim(0, WIDTH)
        ax.set_ylim(0, HEIGHT)
        ax.set_aspect("equal", "box")

J'ai peu commenté le code, mais je peux détailler certaines de ses subtilités :

Comme on va le voir, ce code est bien plus rapide que les méthodes précédentes et on commence à être limité par le temps d'enregistrement des cercles. On se retrouve face à un choix à faire. L'utilisation de numpy.array permet de vectoriser le calcul de distances mais il est inefficace pour ajouter des éléments. La raison en est que c'est un tableau de données contigues dans la mémoire, comme en C. Rapide à lire mais pas fait pour ajouter des éléments. Il se trouve que, de manière contre-intuitive, il est plus rapide de faire un np.concatenate qu'un np.insert. Encore une subtilité Pythonique
Autre choix à effectuer, la taille des feuilles de notre arbre. Si les feuilles contiennent beaucoup de cercles, on peut profiter de la vectorisation mais on perd l'intérêt de l'algorithme. À l'inverse, Python va pénaliser les appels récursifs aux méthodes et les appels à des fonctions non vectorisées. Ce choix de taille est le paramètre THRESHOLD.

Notre troisième méthode d'insertion fait donc appel à cette structure en arbre. L'insertion se fait de manière très similaire aux précédentes

def add_circle_3(base):
    """"Ajoute un cercle dans la boite.

    Parameters
    ==========
    - base: Noeud
        Arbre de partitionnement

    Returns
    =======
    base
    """

    rayon_max = 0
    coin = np.array([WIDTH, HEIGHT])

    # Tant que l'on est dans un autre cercle, on retire un nouveau point
    while rayon_max <= 0:

        # Tir aleatoire d'un nouveau point
        nouv_centre = np.random.uniform(high=(WIDTH, HEIGHT))            

        # Rayon maximum pour rentrer dans la boite
        rayon_max = np.min([nouv_centre, coin - nouv_centre])

        # Comparaison avec tous les cercles de l'arbre
        rayon_max = base.find_collision(nouv_centre, rayon_max)


    # On insère le nouveau cercle
    base.add_circle(nouv_centre, rayon_max)

    return base

Notre structure en arbre est donc composée d'une multitude de boites. Pour avoir un aperçu de cette décomposition, on peut regarder ce qui se passe après le premier appel à split, c'est-à-dire lors la première division.

Le premier enfant (panneau du milieu) va contenir tous les cercles localisés à gauche de la boite, le second enfant (panneau de droite) tous les cercles situés à droite et la base (panneau de gauche) gardera en mémoire tous ceux qui ne peuvent pas rentrer dans l'une de ces deux sous-boites. Soit qu'ils sont trop gros soit qu'ils sont pile à la frontière.

base = Noeud(0, WIDTH, 0, HEIGHT)

for i in range(101):
    base = add_circle_3(base)

fig = plt.figure(figsize=(18 * WIDTH, 6 * HEIGHT))

ax = fig.add_subplot(1, 3, 1)
base.plot(ax)

ax = fig.add_subplot(1, 3, 2)
base.child_1.plot(ax, recursive=True)

ax = fig.add_subplot(1, 3, 3)
base.child_2.plot(ax, recursive=True)

plt.show()

Évaluation des performances

Comparons ces trois méthodes en regardant jusqu'où elles arrivent en une minute

# Methode 1
t0 = time.time()
centres, rayons = np.array([]), np.array([])
time_1 = []
while time.time() - t0 < 60:
    centres, rayons = add_circle(centres, rayons)
    time_1.append(time.time() - t0)

# Methode 2
t0 = time.time()
centres, rayons = np.array([]), np.array([])
time_2 = []
while time.time() - t0 < 60:
    centres, rayons = add_circle_2(centres, rayons)
    time_2.append(time.time() - t0)

# Methode 3
t0 = time.time()
base = Noeud(0, WIDTH, 0, HEIGHT)
time_3 = []
while time.time() - t0 < 60:
    base = add_circle_3(base)
    time_3.append(time.time() - t0)

# Figure
fig, ax = plt.subplots()

ax.plot(time_1, label="Methode 1")
ax.plot(time_2, label="Methode 2")
ax.plot(time_3, label="Methode 3")

ax.legend()
ax.set_xlabel("Nombre de cercles")
ax.set_ylabel("Temps en secondes")

plt.show()

Ce graphique parle de lui-même. S'il y a un léger coût à l'entrée pour mettre en place et gérer cette structure de données, le gain en performances est flagrant !

Méthode 4 : Conversion dans un langage compilé

Python est excellent langage de prototypage, il est très facile d'accès et d'utilisation, et permet de tester rapidement différents algorithmes. Néanmoins, c'est un langage très lent quoi qu'il en soit, et il est bien souvent inévitable de transcrire notre prototype dans un langage compilé performant comme C, C++ ou Fortran ou bien d'utiliser des librairies comme numba qui s'en chargent pour nous.

À titre de comparaison, le code de la méthode 3 écrit en Fortran et executé sur la même machine (je vous épargne le code cette fois-ci)

n_range, t_range = np.loadtxt("fortran.txt").transpose()
n_range = n_range[t_range < 60]
t_range = t_range[t_range < 60]

# Figure
fig, ax = plt.subplots()

ax.plot(time_3, label="Python")
ax.plot(n_range, t_range, label="Fortran")

ax.legend()
ax.set_xlabel("Nombre de cercles")
ax.set_ylabel("Temps en secondes")

plt.show()

Conclusion et la suite

Dans ce billet, nous avons vu la version aléatoire du cercle d'Apollonius. Ça a été une bonne excuse pour discuter des limitations du langage Python, de vectorisation, de structures de données. J'espère vous avoir convaincu que Python est un très bon langage de prototypage, mais qu'il ne faut pas pour autant négliger des langages plus bas niveaux comme C, C++ ou Fortran qui deviennent parfois nécessaires si la performance est un critère important pour vous.

Dans le prochain et dernier billet de cette série, j'aimerais discuter des propriétés particulières de cette fractale, et présenter certains résultats de ma publication dans Physical Review E.

Je vous laisse avec une visualisation d'une simulation contenant un million de cercles. Vous pouvez y zoomer à l'envie et parcourir tous ses recoins !

Random Apollonian Packing II : Simulations en deux dimensions