[Open Source][C++] C++11 et analyse numérique
Bonjour,
Description courte
-------------------
J'ai développé une petite bibliothèque (un fichier .hpp) en C++11 qui permet de manipuler de manières plus agréable des tableaux de données multidimensionnels, données qui sont sensées représenter des valeurs discrétisées d'un problème continu (par exemple en mathématiques, ou très couramment en physique).
Le fichier est auto-documenté, mais on peut générer une documentation Doxygen si besoin. Si le concept vous intéresse, je vous invite à lire la description en détail ci dessous, pour voir ce quoi cette bibliothèque est capable, et si elle peut vous être utile (j'ai du mal à résumer convenablement ce que j'ai voulu faire sans trop faire appel à du vocabulaire de physicien).
Le tout est distribué sous licence zlib, et téléchargeable ici : [dernière version].
Note : pour compiler, ce petit bout de programme nécessite certaines fonctionnalités du C++11 :
- les templates variadiques
- range-based for loop
- le mot clé default pour les fonctions membres
- les listes d'initialisations
- le mot clé auto
- le mot clé decltype
- les typedefs de template (on peut s'en passer en définissant la variable pré-processeur NO_TEMPLATE_TYPEDEF)
Toutes ces fonctionnalités sont implémentées dans gcc-4.7, et je crois que seule la dernière manque à gcc-4.6.
Description longue
-------------------
1] Introduction
Il existe de nombreux domaines (par exemple en recherche, et je pense tout particulièrement à la physique puisqu'il s'agit de mon propre domaine) où l'on a des équations mathématiques à résoudre, mais où tout traitement analytique est impossible. Soit parce que le volume de donnée est trop important, soit parce que les équations en jeu sont trop complexes. Dans ces cas là, on fait appel aux ordinateurs, et on utilise des méthodes numériques.
Dans beaucoup de cas, cela revient à considérer une version discrète d'un problème continu. Par exemple, si l'on étudie la déformation d'une corde de longueur L donnée, on choisit N points sur la corde dont on va étudier la trajectoire, et considérer que la corde est uniquement constituée de ces N points.
Le résultat réel (que l'on obtiendrait si l'on pouvait faire le traitement analytique) est obtenu dans la limite où N est infini (c'est la définition du continu), mais cette limite est numériquement inatteignable à cause de problèmes de précision. On se contente donc en général de prendre N suffisamment grand. Dans ce qui suit, N est appelé le nombre d'échantillons (number of samples pour Shakespear).
La manière la plus simple de traiter ce genre de problème en C++ est de créer un simple tableau de taille N :
Code:
double position_corde[N];
où chaque élément du tableau représente un des N points considérés (ici il s'agit d'un type double, donc ce serait par exemple la déformation transversale de la corde). Puis, on fait évoluer la corde dans le temps. Par exemple :
Code:
1
2
3
4
5
| size_t wrap(int i, size_t N) { return i < 0 ? i%N + N : i%N; }
for (size_t t = 0; t < 10; ++t)
for (size_t i = 0; i < N; ++i)
position_corde[i] = position_corde[wrap(i-1, N)]; |
À chaque pas de temps t, un point i de la corde prend la position du point précédent.
Cette approche, techniquement très simple, présente cependant quelques inconvénients. En premier lieu, il faut "se souvenir" qu'au point i correspond la position réelle i*L/N (et c'est le cas le plus simple, où on attribue la position 0 au premier point). Ensuite, il existe certaines opérations mathématiques (comme la dérivation) qui agissent sur un point et son entourage (typiquement les premiers, voir seconds voisins). Si leur implémentation est simple pour une fonction à un paramètre (représentée donc par un tableau unidimensionnel), elle l'est nettement moins pour une fonction à plusieurs paramètres (comme, par exemple, si notre corde étais un tapis...). Enfin, c'est une approche largement héritée du C, qui ne fait aucun usage des forces du C++.
J'ai donc voulu voir à quel point on pouvait rendre plus agréable l'utilisation de ces concepts grâce au C++, et en particulier grâce à sa toute dernière version, le C++11. J'ai donc mis au point une toute petite bibliothèque, constituée d'un seul fichier (eh oui, tout est template...).
2] Fonctions à une variable
Le résultat tient en une seule classe, nommée utils::numerical_function<T>. C'est une classe template : T correspond au type de donnée contenue (par exemple double dans l'exemple précédent). En interne, elle contient simplement un std::vector, qui contient lui-même les données.
La première chose que j'ai voulu implémenter, c'est une syntaxe proche des fonctions usuelles. Par exemple, position_corde(x) au lieu de position_corde[i], où x est le paramètre réel, équivalent à i*L/N. Pour que ça puisse fonctionner, il faut bien sûr que la fonction stocke avec elle l'information nécessaire, c'est à dire "de combien d'unités sont espacés deux points adjacents, et quel est la position du premier point". Il faut aussi gérer les positions qui sortent du cadre imposé : dans le cas de la corde, que faire si x devient négatif, ou plus grand que L ? Doit-on engueuler l'utilisateur qui ne sait pas ce qu'il fait, laisser planter le programme, ou bien doit-on considérer que le point N+1 est en réalité le point 0 (c'est ce qu'on a fait dans le premier exemple, en prenant i%N) ? Cela revient à choisir les conditions aux limites, comme on dit en physique.
Ensuite, quel point choisir si aucun entier i ne correspond exactement à x*N/L ? Doit on prendre le plus proche ? Ou faire une interpolation quelconque entre les deux points les plus proches ? C'est le choix de l'interpolation.
La première étape quand on construit une telle fonction va donc être de spécifier ces informations :
Code:
1
2
3
4
5
6
7
8
9
| // On créé d'abord l'objet, en spécifiant combien de points on veut traiter
utils::numerical_function<double> position_corde(N);
// ... on choisis le mode d'interpolation (INTERPOLATION_LINEAR
// correspond à une interpolation linéaire, et INTERPOLATION_NEAREST
// choisit le point le plus proche)
position_corde.set_interpolation_method(utils::INTERPOLATION_LINEAR);
// ... puis on précise que l'intervalle est de [0, L], et les conditions
// aux limites (boundaries en anglais) sont périodiques
position_corde.set_range(utils::range(0, L), utils::BOUNDARY_PERIODIC); |
Par défaut :
- l'interpolation est INTERPOLATION_NEAREST (c'est la méthode la plus rapide),
- les conditions aux limites sont "laisser planter le programme si on dépasse". (c'est ce qui se fait de plus performant, si on sait ce qu'on fait, c'est à dire qu'on est sûr de ne jamais dépasser),
- et l'intervalle, si l'utilisateur oublie de préciser, est arbitrairement choisi à [0,1].
Il est possible de changer le nombre d'échantillons N par la suite grâce à la méthode resize(N);.
Une fois que tout ça est fait, on peut utiliser pleinement la fonction. Sauf que...
La syntaxe position_corde(x) (équivalente à position_corde.get_value(x)) est belle et pratique si on veut lire les données contenues dans la fonction. Si en revanche on souhaite écrire de nouvelles choses dedans, on est obligé de travailler avec les indices. J'ai donc mis à disposition une autre méthode pour accéder aux données, via la surcharge de l'opérateur []. De fait, position_corde[i] a strictement la même signification que dans l'exemple précédent (c'est équivalent à position_corde.get_data(i)).
Enfin, pour se simplifier un peu la vie, j'ai aussi ajouté une méthode index_to_param(i), qui permet d'obtenir le x correspondant à l'indice i fourni (et param_to_index(x) fait l'inverse).
Avec tout ça, on peut donc faire exactement la même chose qu'avec la syntaxe "C", mais on a aussi la possibilité de travailler avec les grandeurs réelles (le x) plutôt qu'avec les indices si on le souhaite. Qui plus est, la fonction porte avec elle les informations nécessaires pour convertir un indice en une valeur réelle : on a donc pas besoin de "se souvenir" de quoi que ce soit. Il suffit de tout spécifier au début.
3] Fonctions à plusieurs variables
On a donc traité le premier problème de l'approche "C".
Maintenant, qu'en est-il des fonctions à plusieurs variables ? Je ne voulais pas m'amuser à faire une classe pour la fonctions à une variable, puis une autre pour les fonctions à deux variables et ainsi de suite. J'ai donc opté pour l'un des plus beaux ajouts du C++11 : les templates variadiques.
En réalité, j'ai un peu menti : la fonction est de type utils::numerical_function<T, P>, où P est le nombre de variables de la fonction (1 correspond par exemple à la corde, 2 au tapis, etc. Par défaut, P vaut 1).
Dans le cas P=2, avec N1 échantillons pour la première varible et N2 pour la seconde, et pour un tapis rectangulaire de taille L1 par L2, l'initialisation de la fonction se fait de la manière suivante (grâce aux templates variadiques) :
Code:
1
2
3
4
5
6
7
8
9
10
11
12
13
| // On créé d'abord l'objet, en spécifiant combien de points on veut traiter
utils::numerical_function<double, 2> position_tapis(N1, N2);
// ... on choisis le mode d'interpolation (INTERPOLATION_LINEAR
// correspond à une interpolation linéaire, et INTERPOLATION_NEAREST
// choisit le point le plus proche)
position_tapis.set_interpolation_method(utils::INTERPOLATION_LINEAR);
// ... puis on précise que l'intervalle est de [0, L1] pour la première variable
// et [0, L2] pour la seconde, et les conditions aux limites (boundaries en
// anglais) sont périodiques pour les deux variables
position_tapis.set_range(
utils::range(0, L1), utils::BOUNDARY_PERIODIC,
utils::range(0, L2), utils::BOUNDARY_PERIODIC
); |
Pas beaucoup plus compliqué !
On peut alors accéder aux valeurs en utilisant la notation "fonction" avec les paramètres réels position_tapis(x, y) ou la notation "tableau" avec les indices position_tapis[i][j].
Dans ce contexte, il faut préciser pour les fonctions index_to_param et param_to_index à quel paramètre on fait référence. En effet, l'écart entre deux points i et i+1 n'est pas nécessairement le même qu'entre j et j+1. Dans le cas où l'on a plusieurs variables, on doit alors utiliser la syntaxe index_to_param(p, i), où p est un entier qui donne le rang du paramètre concerné (c'est à dire 0 pour le premier paramètre, 1 pour le second, etc. L'ordre est celui dans lequel les paramètres sont donnés lors d'un appel à la fonction : f(x1, x2, x3, ...)).
En résumé, là encore, on a les même fonctionnalités qu'avec la méthode "C", auxquelles s'ajoute le confort supplémentaire apporté par le C++, quand on le souhaite.
4] Les fonctionnalités supplémentaires
Mais on peut faire plus... D'abord, sur l'accès aux données. J'ai ajouté deux classes, qui sont je pense familières à la majorité d'entre vous : iterator et const_iterator. Elles permettent de parcourir l'ensemble des données contenues par une fonction, de la même manière qu'on parcourt un conteneur standard (en particulier, les fonctions begin() et end() sont disponibles). Cela permet de tirer profit d'un autre ajout du C++11 : les range-based for loops (je n'ai pas pu imaginer de traduction potable...) :
Code:
1
2
| for (double& value : position_tapis)
value = // ... |
... qui serait l'équivalent de :
Code:
1
2
3
| for (size_t i = 0; i < N1; ++i)
for (size_t j = 0; j < N2; ++j)
position_tapis[i][j] = // ... |
Ce n'est pas toujours très pratique, puisqu'on n'a pas d'information sur la "position" de value dans la boucle, mais ça a son intérêt dans certaines situations.
Qui plus est, il y a certaines opérations que l'on est souvent amené à faire, et qui peuvent être pénibles sur ces grands ensembles de données.
Par exemple, mettre toutes les valeurs à zéro (ou n'importe quelle autre nombre). On peut le faire ici de deux manières : soit avec la fonction fill(v), où v est un nombre, soit avec l'opérateur d'affectation : position_tapis = v (en interne, l'opérateur d'affectation redirige vers la fonction fill).
On peut aussi vouloir ajouter un même nombre à toutes les valeurs : on peut le faire avec shift(v), ou avec les opérateurs += et -=.
On peut encore vouloir multiplier toutes les valeurs par un même nombre, ce qui se fait grâce à scale(v), ou avec les opérateurs *= et /=.
Et cetera.
Il est également possible d'utiliser les opérateurs numériques usuels f+v, f*v, etc, mais attention aux performances (en principe, pour N opérations, N copies sont effectuées, si le compilateur est assez malin pour faire de la RVO).
Pour des traitement plus complexes, contenant un grand nombre d'opérations, il peut s'avérer plus efficace de ne faire qu'une seule traversée de toutes les données, et donc d'implémenter le traitement "à la main", à l'intérieur d'une boucle.
On peut aussi utiliser la fonction fill autrement, pour affecter des valeurs différentes à différents points. Typiquement, c'est utile si l'on veut initialiser les données à une fonction connue, par exemple cos(x). Il existe plusieurs manières de le faire, car la fonction fill attend un paramètre qui dispose de l'opérateur () : ce peut être une fonction C++, un foncteur, une std::function, une lambda, ...
Cet opérateur () doit prendre P arguments, qui sont sensés correspondre aux P paramètres de notre fonction. Par exemple pour une fonction à un paramètre :
Code:
1
2
3
4
5
6
| utils::numerical_function<double, 1> corde(N);
corde.set_range(utils::range(0, 3.14159), utils::BOUNDARY_PERIODIC);
// Totalement équivalent :
corde.fill(&cos);
corde = cos;
corde = [](double x) { return cos(x); }; |
et pour une fonction à deux paramètres :
Code:
1
2
3
4
5
6
7
8
9
| utils::numerical_function<double, 2> tapis(N1, N2);
tapis.set_range(
utils::range(-1, 1), utils::BOUNDARY_PERIODIC,
utils::range(-1, 1), utils::BOUNDARY_PERIODIC
);
// Totalement équivalent :
tapis.fill(&atan2);
tapis = atan2;
tapis = [](double x, double y) { return atan2(x, y); }; |
Dernière fonctionnalité de cette classe : les opérations de dérivations et intégrations. Numériquement, on peut calculer la dérivée d'une fonction en un point avec une bonne approximation en mélangeant savamment 5 valeurs proches du point considéré (typiquement : le point lui même et ses 4 plus proches voisins). J'ai donc mis à disposition la méthode derivate(p, x, y, ...), où p est le rang du paramètre par rapport auquel on souhaite dériver (comme dans index_to_param), et x, y, ... sont les coordonnées du point où l'on souhaite calculer la dérivée. Une alternative est derivate_i(p, i, j, ...), où le point considéré est donné directement par ses indices : c'est plus performant, mais parfois moins pratique. Aussi, j'ai implémenté la méthode get_derivative(p=0), qui retourne une nouvelle fonction qui contient la dérivée de la précédente par rapport au paramètre p en chaque point :
Code:
1
2
3
4
| utils::numerical_function<double> f(N), f_prime;
f.set_range(utils::range(0, 3.14159));
f = cos;
f_prime = f.get_derivative(); |
Les données contenues dans f_prime doivent correspondre à la dérivée du cosinus, c'est à dire à -sin(x).
Il est également possible d'intégrer la fonction. Numériquement, cela correspond tout simplement à faire une somme... Si on intègre k fois une fonction à P arguments (avec P>=k), on obtient une nouvelle fonction avec P-k arguments (si P=k, on obtient simplement une seule valeur). Pour ce faire, on utilise la méthode integrate(p, begin, end, ...) où là encore p est le rang du paramètre par rapport auquel on veut intégrer, et begin et end forment l'intervalle sur lequel on souhaite intégrer. Les autre paramètres sont d'autres triplets du même genre, si l'on veut intégrer plusieurs variables à la fois (c'est plus performant que de le faire une par une). À noter cependant : si plusieurs variables sont intégrées en même temps, le code suppose que l'utilisateur a fourni les arguments de integrate par ordre de p croissant.
Comme pour derivate, il existe une alternative integrate_i pour laquelle l'intervalle d'intégration est spécifié directement en indices (c'est également plus performant).
Pour clore cette longue présentation, je vais rapidement présenter la classe printer, qui permet de créer une sortie formatée pour les données d'une fonction. L'exemple le plus simple est la sauvegarde de ces données dans un fichier texte, où les différents arguments sont placés sur différentes colonnes (la valeur correspondante de la fonction étant sur la dernière d'entre elles) :
Code:
1
2
3
4
5
6
7
| utils::numerical_function<double> f(40);
f.set_range(utils::range(0, 3.14159));
f = cos;
utils::numerical_function<double>::printer p("data.txt");
p.add(f);
p.print(); |
À titre d'exemple, j'ai aussi implémenté une sortie vers le logiciel de visualisation gnuplot, qui permet de sauvegarder (sans interaction avec l'utilisateur) le graphe d'une fonction dans un fichier .png.
5] Le mot de la fin
Voilà pour cette présentation. J'ai trouvé que c'était un exercice agréable pour apprendre à manipuler ces templates variadiques. J'ai pu commencer doucement avec des fonction simples comme resize() ou set_range(), pour finir sur des comportements bien plus complexes, par exemple pour integrate().
Au niveau des performances, je n'ai pas encore mis au point de benchmark qui permette de juger de la pénalité induite par tout ce "sucre syntaxique". Pour une fonction à un seul argument, je pense que l'on doit s'approcher très près des performances "C" (si l'on active les optimisations du compilateur). En revanche, il doit être possible d'améliorer les performances pour un nombre d'argument supérieur à 1. En particulier, les données sont actuellement stockées dans un vector de vector de vector de ... etc. Il me semble que les performances pourraient être meilleures en implémentant ça comme un gros tableau unidimensionnel, mais il faudrait faire des tests pour être sûr. Le problème, c'est que ce petit changement bouleverse complètement l'implémentation que j'ai faite d'une partie des fonctionnalités présentées ci-dessus. Donc à voir...
Pour le moment, ça fonctionne :)
