Optimisation code Rust

Comment aller sur la Lune rapidement?

a marqué ce sujet comme résolu.

Bonjour,

J’étudie la trajectoire des missions Apollo (pour potentiellement en faire un article) et j’en profite pour apprendre le Rust.

Un peu de contexte : j’ai implémenté un algorithme itératif (je me perdais dans la version analytique) qui calcule la force d’attraction de la Terre et de la Lune pour en déduire l’accéleration de mon vaisseau. J’intègre deux fois pour en déduire la position.

En pratique ca donne ça (cliquez pour agrandir) : figure_1.png1.

J’ai donc un algorithme qui fonctionne, mais il est (relativement) lent, et j’aimerais l’optimiser. Le but est non seulement d’avoir quelque chose le plus rapide possible, mais tant qu’on y est, faire un petit benchmark (et un autre article, la liste s’allonge vite !) : j’ai commencé par l’implémenter en Python, mais vu la lenteur je l’ai re-implémenté en Rust et en C.

Sans plus attendre, voici les version C et Rust (Python restera loin derrière dans tous les cas), ainsi qu’un petit Makefile pour compiler et exécuter tout ça. Normalement le code est simple à comprendre, mais je peux détailler si besoin.

Actuellement, j’obtiens 1050–1100 ms pour C et 1400–1500 ms pour Rust (et ne parlons pas de Python, on a un facteur 50 ou 100 !) ( :o Rust est plus lent ?!) (je parle de la partie "sim" uniquement, mais si vous voulez optimiser aussi la partie export, pourquoi pas ;) ).

Ma requête est donc : quelles sont les pistes à explorer pour optimiser les version C et Rust (mais plutot Rust que je connais moins) ? Comme l’interet est d’apprendre Rust, si vous avez des remarques sur la forme, ou quoi que ce soit d’autre, je suis ouvert à tout :)

Merci d’avance


  1. Tout est à l’échelle ; distances en mètres ; repère orthonormé centré sur le barycentre Terre-Lune, en rotation avec le système, ce qui donne aux corps une position fixe. Pour info l’altitude initiale est de 185 km à comparer au 6378 km du rayon de la Terre et ~400 000 km de distance Terre-Lune.

Je n’ai pas de très bonnes notions de Rust, mais en jetant un coup d’œil rapide aux deux versions, j’ai vu quelques points qui peuvent améliorer les deux versions.

Pour les optimisations "simples" qui ne changent pas fondamentalement la manière dont le code fonctionne:

  • Il faudrait que tu simplifies un maximum tes calculs quand c’est possible. Par exemple, dans le calcul de l’accélération, tu calcules la norme d’un vecteur (en utilisant une racine carré) et tu la met au carré directement après. Le compilateur ne peux pas optimiser ça parce que tu travailles avec des flottants et que le résultat n’est pas forcément le même que si simplifie avant. Dans la même ligne, tu multiplies par la masse d’un objet pour diviser par cette même masse juste après.
  • Énormément de calcules sont dupliqués à chaque ticks et bénéficieraient à être sauvegardés. Par exemple, G*earth_mass*moon_mass est constant tout au long de l’exécution du programme. Il n’est pas nécessaire de recalculer ça en permanence. Même au sein d’un seul tick, je ne pense pas que le compilateur soit capable de mettre en cache tous les résultats intermédiaires pour les réutiliser plus tard (parce qu’il y a un nombre limité de registres).

Cependant, le plus gros des gains peuvent être tirés en changeant la manière dont ton code fonctionne. Ton système Terre-Lune n’est pas affecté par le vaisseau, ce qui est une simplification parfaitement valide. C’est un système extrêmement simple vu que les deux tournent à vitesse constante autour de leur barycentre. Il serait bien plus efficace de déterminer les paramètres d’orbite (vitesse, distance et barycentre) au début du programme et de calculer la position de la Terre et de la Lune comme ça à chaque tick.

Le plus gros point pour la fin: tu utilises une méthode de simulation simple (méthode d’Euler) en calculant l’accélération, puis la vitesse et enfin la position. Le problème, c’est que cette méthode accumule les erreur de calculs assez vite, ce qui te force à utiliser un pas de temps très faible. Il existe d’autres méthodes, notamment Runge-Kutta qui n’est pas beaucoup plus compliqué, mais qui réduit très nettement l’accumulation d’erreurs et te permettra d’augmenter ton pas de temps. Il existe aussi des méthodes qui préservent les quantités physiques qui sont constantes (quantité de mouvement et quantité d’énergie), ce que la méthode d’Eurer ou de Runge-Kutta ne font pas. Je ne peux malheureusement pas t’en conseiller une vu que je ne les ai jamais utilisé.

Salut,

C’est du détail, mais comme je constate que tu utilises du C99 (au moins), tu peux simplifier l’allocation de ton tableau bidimensionnel comme suit.

Vector (*poss)[3] = malloc(sizeof(Vector[ts][3]));

if (poss == NULL)
    exit(43);

Cela te permet ainsi de n’appeler malloc() qu’une seule fois.

+1 -0

Python ne peut pas rivaliser avec C ou Rust pour du code pur, mais comme dit par @unidan, on peut remarquer où se trouve les goulots d’étranglement et optimiser avec divers outils, comme numpy, pandas, pypy puis par la suite cython, numba, pycuda, …

En règle générale, si le travail de recherche est fait correctement, on peut arriver à du 1,5x plus lent que du code C avec une vitesse d’écriture du code accélérée.

+0 -0

Aussi, autre point d’optimisation :

setbuf(f, NULL);

En désactivant la temporisation, tu augmentes le nombre d’appels système pour l’écriture des résultats dans ton fichier out.txt et, en conséquence, le temps d’exécution.

+1 -0

Merci pour vos réponses :)

setbuf(f, NULL);

En désactivant la temporisation, tu augmentes le nombre d’appels système pour l’écriture des résultats dans ton fichier out.txt et, en conséquence, le temps d’exécution.

Taurre

Ah zut, un copier coller malheureux, je pensais désactiver la bufferisation ligne par ligne. J’aurais du vérifier le man…

Merci @Berdes, comme ça fait la deuxième fois que j’entend parler de Runge-Kutta, je vais aller jeter un oeil, même si ca ne donne pas envie vue de l’extérieur ! Mais bon, ca risque de me reservir plus tard, ca va être un bon investissement.

Merci @Berdes, comme ça fait la deuxième fois que j’entend parler de Runge-Kutta, je vais aller jeter un oeil, même si ca ne donne pas envie vue de l’extérieur ! Mais bon, ca risque de me reservir plus tard, ca va être un bon investissement.

Tu ferais mieux de faire de l’Euler semi-implicite. RK n’est pas symplectique, et tu gagnes pas grand chose en terme de stabilité par rapport à du Euler explicite (sans compter que c’est pénible à coder et difficile à déboguer).

Ah zut, un copier coller malheureux, je pensais désactiver la bufferisation ligne par ligne. J’aurais du vérifier le man…

Nodraak

Par défaut, les fichiers sont normalement temporisés par blocs, il est donc rare d’avoir à modifier la temporisation de ce côté.

Sinon, je retire ce que j’ai dit concernant malloc() (le « C’est du détail ») : vu que visiblement tu l’appelles SIM_DURATION / SIM_DT fois (soit 2.592.000 fois), tu devrais à mon avis avoir une amélioration substantielle en ne l’appelant qu’une seule fois. ^^"

+1 -0

J’ai pas regardé le code dans le détail, mais les 2.592.000 tours de boucles annoncés par @Taurre me paraissent beaucoup. Si tu as vraiment autant de pas de temps, c’est beaucoup trop vue la tête de la courbe que tu obtiens, quelque milliers voire dizaines de milliers devraient suffire. Avant de micro-optimiser, essaye de réduire le nombre de pas de temps, tu vas gagner beaucoup plus.

+0 -0

J’ai pas regardé le code dans le détail, mais les 2.592.000 tours de boucles annoncés par @Taurre me paraissent beaucoup. Si tu as vraiment autant de pas de temps, c’est beaucoup trop vue la tête de la courbe que tu obtiens, quelque milliers voire dizaines de milliers devraient suffire. Avant de micro-optimiser, essaye de réduire le nombre de pas de temps, tu vas gagner beaucoup plus.

adri1

J’affiche volontairement (relativement) peu de point, mais pour la simulation j’ai besoin d’un pas assez petit1. Cela vaut surtout quand ma direction change beaucoup: un pas de 1h en orbit terrestre serait une catastrophe (pour comparaison, une orbite = ~1h30), mais durant le voyage (qui dure ~3–4 jours) en ligne (presque) droite, cela aurait une influence très faible.


  1. Apres tout est relatif. A priori 1 sec devrait largement suffire, mais pour le challenge, je voulais descendre à un pas aussi petit que possible.

Un pas de temps, ça se prend pas au pif en se disant "ça devrait suffire" ou "pour le challenge", ça se calcule (regarde du côté des critères de stabilité, c’est facile à comprendre pour du Euler explicite). En l’occurence, tu as sûrement intérêt à avoir un pas de temps adaptatif puisque tu as des différences importantes d’accélération au cours du trajet. C’est probablement le point critique en terme de perfs de ton programme.

+2 -0

Re! Petite mise à jour concernant mon aventure.

Les code

Je ferai probablement un petit repo git, mais en attendant, voici les dernières versions : C, Rust et Python. Pour rappel, les anciennes versions : C et Rust. (Et le script d’exécution bricolé pour exécuter tout ça facilement). (Et le script de benchmark, attention le code est très moche).

Améliorations implémentées

  • Je fixe la fréquence du CPU pour éviter que le délai d’augmentation de la fréquence ou la baisse de fréquence automatique (pour éviter la surchauffe) impacte le benchmark (j’ai un laptop qui a tendance à avoir facilement chaud): sudo bash -c "echo 1200000 > /sys/devices/system/cpu/cpu0/cpufreq/scaling_max_freq".
  • En plus des flags --release (Rust) et -02 ( C ), j’ai ajouté -Ofast pour C et -OO pour Python (pour ce dernier, pas sur que ça change grand chose ^^ ).
  • Je n’exporte plus "que" toutes les 50 sec (EXPORT_DT), au lieu d’un export à chaque step (on obtient 13 k points, soit 500 km ou 5 deg d’écart en orbite autour de la Terre, une précision largement suffisante pour l’affichage).
  • Méthode d’intégration (Euler) optimisée. J’ai pu facilement implémenter RK4, grâce à https://github.com/TheHappyKoala/Harmony-of-the-Spheres/blob/0c450c032fb1fd3cfcb1d32b7b76e27205fdde9f/src/js/Physics/Integrators/RK4.js.

Résultats

  • 5 exécutions, moyenne des 3 valeurs medianes, les durées sont en ms
  • A noter que pour éviter de passer 30 mn à attendre Python, j’ai augmenté le pas de simulation (facteur 100 pour Python et 10 pour Pypy). La valeur indiquée tient compte de ce facteur, pour faciliter la comparaison.
Langage \ algo v1 (euler-naive) euler-naive euler-optim rk4 
Python - 373 900 334 200 1 978 500
Pypy - 13 190 11 780 972 100
C 3 015 3 552 448 2 012
Rust 4 039 3 944 352 3 328

Plusieurs observations ou remarques, en vrac :

  • RK4 n’est pas comparable directement, j’en parle plus bas.
  • Fun fact pour Python, numpy prend 1 bonne seconde à importer et est 75% plus lent que math. A éviter donc si ce n’est pas nécessaire.
  • Pour la partie export, en réduisant le nombre de valeurs exportées (de 4 / sec à 1 / 50 sec), on gagne les perfs attendues (un facteur ~ 200): de 5 sec en Rust et 35 en C à 50–200 ms.
  • L’optimisation d’Euler permet de gagner un facteur ~ 10 avec C et Rust, mais pas Python. J’avoue que je ne l’explique pas. Je vais probablement étudier l’ASM C/Rust et le bytecode Python.
  • A retenir : C et Rust sont similaires, Pypy/Rust = facteur 3–33, Python/Rust = facteur 100–1000, Python/PyPy = facteur 30.

RK4 vs Euler

Autant il est facile de comparer mes algorithmes euler-naive et euler-optim, autant comparer RK4 l’est moins. Le gain en performances n’est pas évident, il dépend de la précision (step). Parfois RK4 est plus avantageux, parfois moins avantageux.

Je vous met un graph de la précision en fonction du temps de simulation. Les deux courbes (Euler et RK4) se croisent à un moment (vers sim_duration = 700 ms et precision = 2500 km). Si le sujet vous intéresse, j’ai un article en cours d’écriture que je publierais sur mon blog (soyez pas pressés hein, je vise une publication dans 1–3 mois).

(Attention, les axes sont logarithmiques)

Précision en fonction du temps de simulation - axes log-log
Précision en fonction du temps de simulation - axes log-log

Prochaines étapes / Pistes à explorer

  • Calculer la positions de la Terre et de la Lune avec le barycentre et la vitesse angulaire : pas sûr de réellement gagner des perfs, car on remplace des multiplications/divisions par des cos/sin (matrice de rotation).
  • Utiliser un pas de temps dt variable : on peu probablement gagner un facteur entre 1 et 10, donc à explorer en priorité. Par contre je n’ai pas compris comment déterminer le critère de stabilité "facile pour Euler explicite" : ma seconde dérivée (acc) est complexe, et donc la position ne peut pas être exprimée facilement avec une exponentielle. Tu pourrais m’éclairer @adri1 ? Je pourrais bricoler un truc à ma sauce, mais autant faire les choses bien.
  • S’amuser avec l’ASM C/Rust/Pypy et le bytecode Python pour comprendre les optimisations et comparer.
  • Quand le code sera finalisé, faire un benchmark exhaustif puis écrire les articles de blog
  • Quelqu’un veut se chauffer pour me faire la version Java, Javascript ou n’importe quel autre langage ? Les contributions sont bienvenues :)
  • En plus des flags --release (Rust) et -02 ( C ), j’ai ajouté -Ofast pour C et -OO pour Python (pour ce dernier, pas sur que ça change grand chose ^^ ).
Nodraak

Les options -O2 et -Ofast se contredisent. -Ofast revient à préciser -O3, mais en autorisant en plus le compilateur à ne pas respecter le standard.

Sinon, j’aurai été curieux de voir la différence en supprimant l’allocation dynamique via une boucle, comme je l’avais suggéré plus haut.

+0 -0

Les options -O2 et -Ofast se contredisent. -Ofast revient à préciser -O3, mais en autorisant en plus le compilateur à ne pas respecter le standard.

Ah zut, effectivement, j’aurais du retirer -O2. Encore une correction à apporter pour la prochaine itération.

Sinon, j’aurai été curieux de voir la différence en supprimant l’allocation dynamique via une boucle, comme je l’avais suggéré plus haut.

Taurre

Pardon, j’avais ignoré ce point, car le malloc est en dehors de la boucle de simulation, qui est la seule partie que je trouver pertinente à chronometrer. Je viens de mesurer la durée du malloc, ~ 1.900 ms vs 0.015 ms.

Pardon, j’avais ignoré ce point, car le malloc est en dehors de la boucle de simulation, qui est la seule partie que je trouver pertinente à chronometrer. Je viens de mesurer la durée du malloc, ~ 1.900 ms vs 0.015 ms.

Nodraak

Pas de soucis, merci pour avoir pris la mesure (surtout si vite). ^^
Cela fait déjà une belle différence. :)

+0 -0

Utiliser un pas de temps dt variable : on peu probablement gagner un facteur entre 1 et 10, donc à explorer en priorité. Par contre je n’ai pas compris comment déterminer le critère de stabilité "facile pour Euler explicite" : ma seconde dérivée (acc) est complexe, et donc la position ne peut pas être exprimée facilement avec une exponentielle. Tu pourrais m’éclairer @adri1 ? Je pourrais bricoler un truc à ma sauce, mais autant faire les choses bien.

J’aurais dit plutôt que tu peux gagner un facteur 100 perso. Effectivement, la détermination du critère de stabilité est plus délicate que ce que je pensais à vue de nez. Pas sûr que ce soit possible, même, en fait, comme les équations sont fortement non linéaires. J’ai fouillé sans succès. Une possibilité assez classique lorsque le pas de temps est difficile à déterminer analytiquement est de laisser le pas de temps évoluer par lui même en suivant l’évolution de la solution. Le principe est là.

Sinon ton implémentation de Runge-Kutta n’est pas bonne, ton v devrait être calculé avec la vitesse actualisée, pas avec l’ancienne. Il faudrait d’ailleurs vérifier que tu as le droit de télescoper deux RK sur des équations de premier ordre pour obtenir un RK sur une équation de second ordre…

En regardant le code C, je me suis demandé si le compilateur était effectivement capable de déduire assez d’informations pour être capable d’inliner Body_compute_acc_optim et de correctement réduire la boucle en une ou deux itérations suivant le cas. Un petit coup d’œil à l’assembleur généré indique que ce n’est pas le cas: https://godbolt.org/z/Yv_-kl. J’ai donc décidé de faire une version C++ en (sur)utilisant les templates pour m’assurer que le compilateur possède toutes les informations nécessaires pour bien optimiser le code.

J’en ai aussi profiter pour corriger la manière dont la simulation est faite. Comme @adri1 l’a indiqué, il n’est pas correcte de télescoper deux Euler (ou deux RK) sur du premier ordre pour obtenir l’équivalent au deuxième ordre. La manière correcte est de subtilement transformer le problème du deuxième ordre en un problème du premier ordre. Cette transformation se fait très facilement en établissant que l’état du système est le couple (position, vitesse). La dérivée deviens alors (vitesse, accélération). Les méthodes d’Euler et RK s’appliquent alors de manière triviale. De manière intéressante, le fait que tu utilises la nouvelle vitesse pour mettre à jour la position (en télescopant deux intégrations) corrige en grosse partie l’accumulation d’erreur typique de la méthode d’Euler. C’est probablement une des raisons pour lesquels RK4 te donne peux de gains en vitesse/erreur. Avec mon code, RK avec dt=0.1s est plus précis qu’Euler avec dt=0.001s. Au passage, quels valeurs as-tu utilisé pour calculer l’erreur?

Au final, le code C++ que j’ai fait est à peu près deux fois plus rapide que ton code C (avec euler-optim) en utilisant clang -O2 (-O3 n’apporte aucun gain), ce qui confirme mon intuition que le compilateur manque d’informations pour faire toutes les optimisations possible sur ton code C. J’ai aussi séparé ce qui correspond à du code réutilisable (dans le namespace anonyme) et ce qui est du code spécifique pour le problème à résoudre (en dehors du namespace anonyme).

Pour référence, RK4 me donne les valeurs suivantes à la fin de la simulation (dt=0.0001s pour 10m de simulation): 863950 -4666739 -7228 379617140 587951 -57658129 -333453994. Diminuer le pas de temps ne change quasiment rien aux valeurs: passer de dt=0.001s à dt=0.0001s ne change la position du vaisseau que de 0.7km.

Edit: Je viens juste de réaliser que je n’ai pas expliqué l’utilisation de absl::StrFormat. C’est une fonction de la lib abseil qui donne des performances très nettement supérieur à la version originale C++ ou même fprintf. Chez moi, la version C++ native prends 60ms, fprintf 45ms et absl::StrFormat 10ms.

+0 -0
Connectez-vous pour pouvoir poster un message.
Connexion

Pas encore membre ?

Créez un compte en une minute pour profiter pleinement de toutes les fonctionnalités de Zeste de Savoir. Ici, tout est gratuit et sans publicité.
Créer un compte