---
title: "Arbres quaternaires, compression et Barnes-Hut"
date: "2022-04-06"
patat:
  eval:
    tai:
      command: fish
      fragment: false
      replace: true
    ccc:
      command: fish
      fragment: false
      replace: true
  images:
    backend: auto
---

# Le cours précédent

## Questions


* Structure de données d'un arbre quaternaire?

. . .

```C
typedef struct _node {
    int info;
    struct _node *child[4];
} node;
```

. . .

* Dessin d'un noeud d'arbre quaternaire (avec correspondance `node`)?

. . .

```{.mermaid format=pdf width=400 loc=figs/}
graph TD;
    id0(("    "))-->|"child[0]"| id1(("info"));
    id0-->|"child[1]"| id2(("info"));
    id0-->|"child[2]"| id3(("info"));
    id0-->|"child[3]"| id4(("info"));
```

# Le cours précédent

## Questions

* Comment faire la symétrie d'axe horizontal?

. . .

```C
arbre symétrie(arbre)
    si !est_feuille(arbre)
        échanger(arbre.enfant[0], arbre.enfant[2])
        échanger(arbre.enfant[1], arbre.enfant[4])
        pour i de 0 à 3
            symétrie(arbre.enfant[i])
    retourne arbre
```

# Compression sans perte (1/N)

## Idée générale

* Regrouper les pixels par valeur

```
   SG=0  |  SD=1        SG=0   |  SD=1   
 21 | 12 | 4 |  4      21 | 12 |   4
  9 |  7 | 4 |  4       9 |  7 |  
-----------------  => ----------------- 
  1 |  1 | 0 | 31         1    |  0 | 31
  1 |  1 | 3 | 27              |  3 | 27
   IG=2  |  ID=3        IG=2   |  ID=3
```

* Comment faire?

# Compression sans perte (2/N)

## Que devient l'arbre suivant?

![](figs/quad_img_simple.svg)

. . . 

## Arbre compressé

![](figs/quad_img_simple_comp.svg)

# Compression sans perte (3/N)

* Si un noeud a tous ses enfants égaux:
    * Donner la valeur au noeud,
    * Supprimer les enfants.
* Remonter jusqu'à la racine.

## Écrire le pseudo-code (5min, matrix)

. . .

```C
rien compression_sans_pertes(arbre)
    si !est_feuille(arbre)
        pour i de 0 à 3
            compression_sans_pertes(arbre.enfant[i])
        si derniere_branche(arbre)
            valeur, toutes_égales = valeur_enfants(arbre)
            si toutes_egales
                arbre.info = valeur
                detruire_enfants(arbre)
```

# Compression sans perte (4/N)

\footnotesize

## Écrire le code C (5min, matrix)

. . .

```C
void lossless_compression(node *qt) {
    if (!is_leaf(qt)) {
        for (int i = 0; i < CHILDREN; i++) {
            lossless_compression(qt->child[i]);
        }
        if (is_last_branch(qt)) {
            int val = -1;
            if (last_value(qt, &val)) {
                qt->info = val;
                for (int i = 0; i < 4; ++i) {  
                    free(qt->child[i]);
                    qt->child[i] = NULL;
                }
            }
        }
    }
}
```

# Compression sans perte (5/N)

\footnotesize

```C
bool is_last_branch(node *qt) {
    for (int i = 0; i < 4; ++i) {
        if (!is_leaf(qt)) {
            return false;
        }
    }
    return true;
}
bool last_value(node *qt, int *val) {
    int info = qt->child[0];
    for (int i = 1; i < 4; ++i) {
        if (info != qt->child[i]) {
            return false;
        }
    }
    *val = info;
    return true;
}
```


# Compression avec perte (1/N)

## Idée générale

* Regrouper les pixels par valeur sous certaines conditions

```
   SG=0  |  SD=1        SG=0   |  SD=1   
 21 | 12 | 4 |  3      21 | 12 |    4
  9 |  7 | 4 |  4       9 |  7 |  
-----------------  => ------------------
  1 |  1 | 0 | 31         1    |  0 | 31
  2 |  1 | 3 | 27              |  3 | 27
   IG=2  |  ID=3        IG=2   |  ID=3
```

* On enlève si l'écart à la moyenne est "petit"?

# Compression avec perte (2/N)

## Que devient l'arbre suivant si l'écart est petit?

![](figs/quad_img_simple_variation.svg)

. . . 

## Arbre compressé

![](figs/quad_img_simple_comp_loss.svg)

# Compression avec perte (3/N)

## Comment mesurer l'écart à la moyenne?

. . .

* Avec l'écart-type

\begin{equation*}
\mu = \frac{1}{4}\sum_{i=0}^{3} p[i],\quad \sigma = \sqrt{\frac{1}{4}\sum_{i=0}^3 (\mu-p[i])
^2} = \sqrt{\frac{1}{4}\left(\sum_{i=0}^3p[i]^2\right)-\mu^2}
\end{equation*}

## Que devient l'algorithme?

. . .

* Si $\sigma<\theta$, $\theta$ est la **tolérance**:
    * Remplacer la valeur du pixel par la moyenne des enfants.
    * Remonter les valeurs dans l'arbre.

## Quelle influence de la valeur de $\theta$ sur la compression?

. . .

* Plus $\theta$ est grand, plus l'image sera compressée.

# Compression avec perte (4/N)

## Que devient l'arbre avec $\theta=0.5$?

![L'arbre original.](figs/quad_img_simple_variation.svg)

. . .

![Arbre compressé.](figs/quad_img_simple_comp_avg.svg)

# Compression avec perte (6/N)

## Modifications sur la structure de données?

. . .

* On stocke la moyenne, et la moyenne des carrés.

```C
struct noeud
    flottant moyenne, moyenne_carre
    node enfants[4]
```

* Comment on calcule `moyenne` et `moyenne_carre` sur chaque noeud (pseudo-code)?

# Calcul de la moyenne

## Pseudo-code (5min, matrix)

. . .

```C
rien moyenne(arbre) {
    si !est_feuille(arbre)
        pour enfant dans arbre.enfants
            moyenne(enfant)
        pour enfant dans arbre.enfants
            arbre.moyenne += enfant.moyenne
            arbre.moyenne_carre += enfant.moyenne_carre
        arbre.moyenne /= 4
        arbre.moyenne_carre /= 4
```

# La compression avec pertes

\footnotesize

## Pseudo-code (5min, matrix)

. . .

```C
rien compression_avec_pertes(arbre, theta)
    si !est_feuille(arbre)
        pour i de 0 à 3
            compression_avec_pertes(arbre.enfant[i])
        si derniere_branche(arbre)
            si racine(arbre.moyenne_carre - arbre.moyenne^2) < theta
                detruire_enfants(arbre)
```

## Le code en entier

```C
arbre = matrice_à_arbre(matrice)
moyenne(arbre)
compression_sans_pertes(arbre)
```

# La dynamique des corps célestes

## Slides très fortement inspirés du cours de J. Latt, Unige

## Simulation du problème à $N$-corps

* Prédiction du mouvement d'un grand nombre de corps célestes.
* Modélisation:
    * On se limite aux étoiles;
    * Chaque étoile est caractérisée par un point (coordonnées) et une masse;
    * On simule en deux dimensions.
    * Interactions uniquement par les lois de la gravitation Newtonienne (oui-oui c'est de la **physique**!).


# Les équations du mouvement

## Mouvement de la $i$-ème étoile

* Algorithme de Verlet ($t_{n+1}=t_n+\delta t$)

    $$
   \vec x_i(t_{n+1})= 2\vec x_i(t_n)-\vec x_i(t_{n-1})+\vec a_i(t_n)\delta t^2.
   $$

## Force de gravitation

* $\vec a_i(t_n)=\vec F_i/m_i$.
* Sur l'étoile $i$, la force résultante est donnée par

    $$
    \vec F_i=\sum_{j=1,j\neq i}^N \vec F_{ij}.
    $$
    avec
    $$
    \vec F_{ij}=\frac{G m_i m_j(\vec x_j-\vec x_i)}{||\vec x_j-\vec x_i||^3}.
    $$

# Algorithme du problème à $n$-corps

## Pseudo-code: structure de données

```C
struct étoile
    flottant m
    vec x, x_precedent, f
```

## Pseudo-code: itération temporelle

```C
rien iteration_temporelle(étoiles, dt)
    pour étoile_une dans étoiles
        étoile_une.f = 0
        pour étoile_deux dans étoiles
            si (étoile_un != étoile_deux)
                étoile_une.f += 
                    force(étoile_une, étoile_deux)
    pour étoile dans étoiles
        étoile.x, étoile.x_precedent = 
            verlet(étoile.x, étoile.x_precedent, 
                   étoile.f / étoile.m, dt)
```

# Algorithme du problème à $n$-corps

## Complexité

* Complexité de chacune des parties?

. . .

* $\mathcal{O}(N^2)$, $\mathcal{O}(N)$.

## En temps CPU pour **une itération**

\footnotesize

* Si temps pour $N=1$ on calcule en $1\mu s$:

+--------+-------+-------+-----------+
|     N  |  N^2  | t [s] | t [réel]  |
+--------+-------+-------+-----------+
|    10  |  10^2 | 1e-4  |           |
+--------+-------+-------+-----------+
|  10^4  |  10^8 | 1e+2  | ~1min     |
+--------+-------+-------+-----------+
|  10^6  | 10^12 | 1e+6  | ~11j      |
+--------+-------+-------+-----------+
|  10^9  | 10^18 | 1e+12 | ~30k ans  |
+--------+-------+-------+-----------+
|  10^11 | 10^22 | 1e+16 | ~300M ans |
+--------+-------+-------+-----------+

* Typiquement il y a des milliers-millions d'itérations.
* Il y a $10^{11}$ étoiles dans la galaxie.
* Houston we have a problem.

# Question

## Comment faire mieux, des idées?

. . .

* Si un groupe d'étoiles est suffisamment loin, on le modélise comme un corps unique situé en son centre de masse.
* Exemple: Si on simule plusieurs galaxies, on considère chaque galaxie comme un corps unique!
* Un arbre quaternaire est une structure parfaite pour regrouper les étoiles.

# Le cas à 10 corps

::: columns

:::: {.column width=50%}

## Illustration: le cas à 10 corps

![](figs/nbody_bare.png){width=60%}

::::

:::: {.column width=50%}

## Problématique

* On veut calculer la force sur $1$.

::::

:::

. . .


::: columns

:::: {.column width=50%}

## Illustration: le cas à 10 corps

![](figs/nbody_n2.png){width=60%}


::::

:::: {.column width=50%}

## Résultat

* Calcul et somme des forces venant des $9$ autre corps.

::::

:::

# Le cas à 10 corps

::: columns

:::: {.column width=50%}

## Réduction d'un groupe à un seul corps

![](figs/nbody_group.png){width=100%}

::::

:::: {.column width=50%}

## Idée

* On accélère le calcul en traitant un groupe comme un seul corps.
* Fonctionne uniquement si le groupe est assez loin.
* Autrement l'approximation est trop grossière.

::::

:::

# Solution: l'arbre quaternaire

## Corps célestes - arbre

![](figs/nbody_qt_withtree.png)

* On omet les noeuds vides pour éviter la surcharge.
* La numérotation est:
    * 0: ID
    * 1: SD
    * 2: IG
    * 3: SG

# Exemple d'insertion

::: columns

:::: {.column width=50%}

## Insertion corps 1

![](figs/corps1.png){width=100%}

::::

:::: {.column width=50%}

## Arbre, niveau 1

![](figs/arbre1.png){width=100%}

* Quadrant ID.
* La feuille est vide, on insère.

::::

:::

# Exemple d'insertion

::: columns

:::: {.column width=50%}

## Insertion corps 2

![](figs/corps2.png){width=100%}

::::

:::: {.column width=50%}

## Arbre, niveau 1

![](figs/arbre2.png){width=100%}

* Quadrant SD.
* La feuille est vide, on insère.

::::

:::

# Exemple d'insertion

::: columns

:::: {.column width=50%}

## Insertion corps 3 (1/N)

![](figs/corps3_1.png){width=100%}

::::

:::: {.column width=50%}

## Arbre, niveau 1

![](figs/arbre3_1.png){width=100%}

* Quadrant SD.
* La feuille est prise par 2.

::::

:::

# Exemple d'insertion

::: columns

:::: {.column width=50%}

## Insertion corps 3 (2/N)

![](figs/corps3_2.png){width=100%}

::::

:::: {.column width=50%}

## Arbre, niveau 2

![](figs/arbre3_2.png){width=100%}

* On crée un nouveau noeud.
* Deux corps dans le noeud ID.
* On crée un nouveau noeud.

::::

:::

# Exemple d'insertion

::: columns

:::: {.column width=50%}

## Insertion corps 3 (3/N)

![](figs/corps3_3.png){width=100%}

::::

:::: {.column width=50%}

## Arbre, niveau 3

![](figs/arbre3_3.png){width=100%}

* 2 va dans ID.
* 3 va dans SG.
* C'est des feuilles vides, tout va bien.

::::

:::

# Exemple d'insertion

::: columns

:::: {.column width=50%}

## Que fait-on avec les noeuds intérieurs?

* On les utilise pour:
    * stocker la masse totale;
    * stocker le centre de masse.

\begin{align}
m&=m_2+m_3,\\
\vec x &= \frac{m_2\vec x_2+m_3\vec x_3}{m}.
\end{align}

## Chaque feuille contient **une étoile**

::::

:::: {.column width=50%}

## Arbre

![](figs/arbre3_3.png){width=100%}

::::

:::

# Résumé

* Insertion du corps `c` dans le noeud `n` en partant de la racine.
* Si le noeud `n`
    * ne contient pas de corps, on y dépose `c`,
    * est interne, on met à jour masse et centre de masse. `c` est inséré récursivement dans le bon quadrant.
    * est externe, on subdivise `n`, on met à jour la masse et centre de masse, on insère récursivement les deux noeuds dans les quadrants appropriés.

## Remarque

* Il faut stocker les coordonnées des quadrants.
* Un noeud a un comportement différent s'il est interne ou externe.

# Algorithme d'insertion

## Structure de données

```C
struct node
    etoile e // externe: pour stocker 
    etoile sup_etoile // interne: pour stocker m, x
    quadrant q  // coordonnées du quadrant
    node enfants[4]
```

## Remarque: 

* On fait une simplification "moche": `sup_etoile` pourrait juste avoir une masse et une position.

# Algorithme d'insertion

\footnotesize

## Algorithme d'insertion, pseudo-code (15min, matrix)

. . .

```C
rien insertion_etoile(arbre, e)  
    si (!est_vide(arbre) && dans_le_quadrant(arbre.q, e.x)) {
        si (est_feuille(arbre))
            si (!contient_etoile(arbre))
                arbre.e = e
            sinon
                // on crée enfants et arbre.sup_etoile est initialisée
                subdivision_arbre(arbre, e) 
                pour enfant dans arbre.enfants
                    insertion_etoile(enfant, arbre.e)
                pour enfant dans arbre.enfants
                    insertion_etoile(enfant, e)
                destruction(arbre.e)
        sinon
            maj_masse_cdm(arbre.sup_etoile, e)
            pour enfant dans arbre.enfants
                insertion_etoile(enfant, e)
```

# Utilisation de l'arbre

* L'arbre est rempli: comment on calcule la force sur le corps 1?
* Parcours de l'arbre: 
    * si la distance entre 1 et le centre de masse est suffisante, on utilise la masse totale et centre de masse pour calculer la force.
    * sinon, on continue le parcours

# Calcul de la force

## Calcul de la force sur `1`

![](figs/force_1.png)

* Le cadrant ID ne contient que `1`, rien à faire.

# Calcul de la force

## Calcul de la force sur `1`

![](figs/force_2.png)

* Le cadrant SG ne contient `5` corps.

# Calcul de la force

## Calcul de la force sur `1`

![](figs/force_3.png)

* La distance entre `1` et le centre de masse de SG est `d`.

# Calcul de la force

## Calcul de la force sur `1`

![](figs/force_4.png)

* La distance entre `1` et le centre de masse de SG est `d`.
* Est-ce que `d` est assez grand?
* On va comparer avec la distance `d` avec la taille du quadrant `s`.

# Critère $\theta$

* On compare $d=||\vec x_1-\vec x_{cm}||$ avec $s$ la taille du quadrant.
* Le domain est assez éloigné si

    $$
    \frac{s}{d}<\theta,
    $$
* $\theta$ est la valeur de seuil.
* Une valeur typique est $\theta=0.5$, donc la condition devient

    $$
    d>2s.
    $$

# Calcul de la force

## Calcul de la force sur `1`

![](figs/force_4.png)

* Ici $d<2s$, domaine rejeté.
* ON descend dans l'arbre.

# Calcul de la force

## Calcul de la force sur `1`

![](figs/force_5.png)

* `s` est plus petit, mais....
* Cela ne suffit pas $d<2s$, domaine rejeté.

# Calcul de la force

## Calcul de la force sur `1`

![](figs/force_6.png)

* Les noeuds sont des feuilles, on calcule la force.
* On ajoute la force qu'ils exercent sur `1`.

# Algorithme pour le calcul de la force

Pour calculer la force sur un corps `c`, on parcourt l'arbre en commençant par la racine:

* Si le noeud `n` est une feuille et n'est pas `c`, on ajoute la force dûe à `n` sur `c`;
* Sinon si $s/d<\theta$, on traite `n` comme une feuille et on ajoute la force dûe à `n` sur `c`;
* Sinon on continue sur les enfants récursivement.


## Cotinuous notre exemple précédent!

# Calcul de la force

## Calcul de la force sur `1`

![](figs/force_7.png)

* Il y a deux corps dans le quadrant vert.
* Quel est le critère pour remplacer les étoiles par leur centre de masse?

. . .

* Et oui! $d>2s$ on peut remplacer les étoiles par leur centre de masse!

# Algorithme du calcul de force

## Écrire le psuedo-code du calcul de la force

\footnotesize

```C
rien maj_force_sur_etoile(arbre, e, theta)
    si est_vide(arbre)
        retourne

    si est_feuille(arbre) && contient_etoile(arbre) && dans_le_quadrant(arbre.q, e.x)
        maj_force(e, arbre.e)
    sinon si noeud_assez_loin(arbre, e, theta)
        maj_force(e, arbre.sup_etoile)
    sinon
        pour enfant dans enfants
            maj_force_sur_etoile(enfant, e, theta)
```
