Optimisation code Rust

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) : ¹.

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 !) ( 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

Aussi, si tu veux optimiser un code sans changer l’algorithme, il faut profiler.

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.

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.

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.

Merci pour vos réponses

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.

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. ^^"

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.

J’affiche volontairement (relativement) peu de point, mais pour la simulation j’ai besoin d’un pas assez petit¹. 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.

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.

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.

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)

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

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.

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.

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.

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

À la vue du nombre de boucles dans ton code C, je pense que cyton + numpy peuvent apporter un gros plus et se rapprocher un peu plus de C et Rust sans les toucher, mais il sera nettement moins ridicule.

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.