Utiliser GSL et Python dans un programme en C

Il m'arrive fréquemment de programmer en C pour les codes de calcul qui nécessitent des performances incompatibles avec des langages comme Python ou MatLab. Et pourtant, sur TangenteX.com, je ne propose plus que des codes en Python : pourquoi ?
D'abord parce que l'apprentissage du C est plus laborieux que celui de Python ou MatLab. Mais surtout parce que C est un langage généraliste et que la conception de codes de physique numérique nécessite l'usage de nombreuses librairies. Par exemple et surtout, C ne dispose pas de librairie graphique alors que nos codes doivent quasiment toujours afficher des courbes. C ne dispose pas non plus de librairie de calcul, nous sommes donc condamnés à ré-écrire les fonctions de calcul ou de résolution de systèmes ou autres EDO. C'est formateur, mais un peu pénible à force.
Certes, sur les plateformes professionnelles, nous disposons de toutes les librairies souhaitables, puissantes, bien documentées et très chères ! Mais comment faire sans ?
En réfléchissant au problème, j'ai trouvé une solution, sans doute pas la seule. J'aime bien la librairie graphique standard de Python, matplotlib. Et j'ai sous la main la librairie GSL, pour GNU Scientific Library, que j'ai déjà utilisé il y a longtemps avec Dev-ccp. Ces deux librairies remplissent la plupart de nos besoins en physique numérique, tant pour le graphique que pour le calcul. Voyons comment les utiliser dans un programme C.

Les pré-requis

J'ai travaillé simultanément sur deux systèmes différents : Linux Kubuntu version 16.04.2 et macOS Sierra 10.12.5. Mon objectif était de vérifier si l'exercice est plus facile dans un cas ou dans l'autre. La réponse est simple : c'est pareil ! Je n'ai pas tenté le coup sur Windows, mais cela doit être très semblable.

La librairie GSL

Que ce soit sur Mac ou sur Linux, il vous faudra tout d'abord installer la librairie GSL. La méthode dépend de votre environnement. Sous Kubuntu, le gestionnaire de package vous permet de faire ça très facilement. Sous macOS, vous trouverez dans Apps Store la méthode de chargement. GSL est aussi disponible sous Windows.
Pour information, j'utilise la version 2.3 de GSL, mais les versions antérieures devraient convenir également. La librairie GSL est gratuite, disponible sous licence GNU.

L'environnement de développement C

Sur Kubuntu ou macOS, j'utilise CodeBlocks comme environnement de développement en C et C++. Il est disponible gratuitement sous licence et facilement téléchargeable et installable dans les deux environnements. Il est aussi disponible sous Windows.
A noter que sous macOS, la version de CodeBlocks que j'utilise (v 13.12) est assez instable. Elle se plante de temps en temps sans prévenir, alors sauvegardez votre code souvent !

Le point le plus délicat est la configuration du linker afin qu'il retrouve les librairies que nous allons utiliser. Dans CodeBlocks, allez avec le menu à Settings/Compiler/Linker settings. Dans la fenêtre de gauche, ajoutez les chemins des librairies, en utilisant le bouton Add en bas de la fenêtre.
Sur macOS et avec mon installation, j'ai les deux chemins suivants:

Sur Kubuntu, j'ai les chemins suivants :

Bien sur, il vous faudra adapter ces chemins à votre propre installation.

L'environnement Python

L'environnement Python que vous utilisez habituellement conviendra parfaitement, toutes les librairies nécessaires sont installées par défaut. J'ai fait les essais avec Python 2.7. Je ne sais pas si c'est directement portable en 3.x.

L'application

L'objet du programme

Le principe du programme est simple et reprend les habitudes "ancestrales" : on fait le calcul, on stocke les résultats du calcul, on les exporte vers le logiciel de tracé graphique.
Pour illustrer mon propos, j'ai repris l'exemple "bateau" de la résolution de l'équation différentielle du pendule, dans sa version non linéaire. Nous allons donc retrouver les éléments de code dont nous avons l'habitude : la définition du système différentiel, la résolution du système et l'affichage de la courbe résultante. Mais les détails vont sans doute vous surprendre par rapport au code Python correspondant...

La partie GSL

Les #include

Pour pouvoir utiliser la librairie GSL, il faut d'abord définir les en-têtes de chaque fonction que nous allons appeler dans le code. Ces en-têtes sont listées dans les fichiers "header" attachés à chaque librarie. Dans notre cas, pour GSL, nous utiliserons les fichiers header suivants, que nous déclarerons avec l'instruction du pré-processeur C #include :

#include "gsl/gsl_errno.h"

#include "gsl/gsl_odeiv2.h"

#include "gsl/gsl_math.h"

La définition du système différentiel

A part la syntaxe particulière du C, vous ne trouverez rien d'anormal ici :
int SystemeDiff (double t, const double y[], double dydx[], void *params)
{
  double omega2  = * (double*) (params);
  dydx[0] = y[1];
  dydx[1] = - omega2*sin(y[0]);
  return GSL_SUCCESS;
}

La résolution de l'EDO

Là, vous allez regretter la concision et l'élégance du code Python ! La librairie GSL permet une grande maitrise des algorithmes de résolution des EDO, au prix d'une complexité de mise en oeuvre non négligeable. Je ne vais pas rentrer dans les détails, que vous retrouverez dans la doc en ligne de GSL, mais ici par exemple, je déclare vouloir utiliser l'algorithme RK4, sur un système différentiel de 2 équations, avec une erreur d'intégration de 1e-6 :

  const double eps_abs = 1e-6;
  const double eps_rel = 0.0;
  const gsl_odeiv2_step_type *T = gsl_odeiv2_step_rk4;
  gsl_odeiv2_step      *s = gsl_odeiv2_step_alloc(T,n_equations);
  gsl_odeiv2_control   *c = gsl_odeiv2_control_y_new(eps_abs,eps_rel);
  gsl_odeiv2_evolve    *e = gsl_odeiv2_evolve_alloc(n_equations);

puis j'indique à GSL où est défini le système différentiel à intégrer :
  gsl_odeiv2_system sys = {SystemeDiff, NULL, n_equations, &params };
et enfin la boucle de calcul, qui parcourt le temps entre 0 et tmax avec un pas égal à dt :
  fp = fopen(NomFile,"w+");
  while (t < tmax)
    {      
      status = gsl_odeiv2_evolve_apply_fixed_step(e,c,s,&sys,&t,dt,y);
      if (status != GSL_SUCCESS)
      {          
          printf("erreur calcul : %d\n", status);
          break;
      }      
      fprintf(fp,"%.5e;%.5e\n", t, y[0]);
    }
  fclose(fp);

Avant de lancer le calcul, j'ouvre un fichier texte dans lequel je vais stocker, pour chaque pas de calcul, la valeur de t et de \( \theta (t) \). La résolution pas à pas du système différentiel est assurée par la routine gsl_odeiv2_evolve_apply_fixed_step(e,c,s,&sys,&t,dt,y).

A l'issue de l'exécution de la boucle, je dispose d'un fichier qui contient toutes les informations nécessaires pour procéder à l'affichage de la courbe résultante.

Libération de la mémoire

GSL créé de nombreux objets en mémoire, qu'il vaut mieux supprimer sous peine d'avoir de sérieux problèmes (Python fait ça tout seul comme un grand). Les instructions suivantes sont destinées à supprimer ces objets :

gsl_odeiv2_evolve_free(e);
gsl_odeiv2_control_free(c);
gsl_odeiv2_step_free(s);

La partie Python

Les concepteurs de l'interpréteur Python ont eu la très bonne idée de produire une API (Application Programming Interface) qui permet d'utiliser les fonctions de l'interpréteur en les appelant depuis un programme écrit en un autre langage, en C par exemple.
C'est précisemment ce que nous allon faire en lançant depuis notre programme C l'interpréteur Python puis en lui demandant d'exécuter les instructions que nous lui indiquerons.

Les #include

Les en-têtes des fonctions C de l'API Python sont décrites dans le fichier Python.h, d'où l'instruction :

#include <Python.h>

L'appel des intructions Python

Je lance l'interpréteur Python avec les instructions :

  Py_SetProgramName(argv[0]);
  Py_Initialize();

La première instruction permet de récupérer l'environnement du programme et de ne pas se casser la tête avec les différents chemins utilisés par Python et autres variables d'environnement.
Maintenant, je passe à Python toutes les instructions que je veux qu'il exécute, et là vous êtes en terrain connu : 

PyRun_SimpleString("from numpy import loadtxt\n");
PyRun_SimpleString("from matplotlib.pyplot import plot, show,xlim\n");
PyRun_SimpleString("x,y = loadtxt('NumLab/RK4GSL/TempData.txt',unpack=True, delimiter=';')\n");
PyRun_SimpleString("xlim(min(x),max(x))\n");
PyRun_SimpleString("plot(x,y)\n");
PyRun_SimpleString("show()\n");

Enfin, je sors de l'interpréteur Python :

 Py_Finalize();

Le scriptPython est des plus standard, sauf peut-être la fonction loadtxt() qui permet de lire un fichier texte dont les données sont séparées par un caractère défini dans le paramètre delimiter, ici un ; . Très pratique cette instruction...
Bien sur, vous pouvez ajouter des instructions, pour ajouter un titre ou des libellés sur les axes. Vous pouvez aussi appeler directement le script Python en utilisant l'instruction PyRun_SimpleFile(), mais je trouve plus drôle de pouvoir modifier les instructions Python en ligne dans mon code C.

Le code source C

Le code source C du programme est disponible ici.

Les résultats

Voici la courbe obtenue pour un angle initial égal à 0.1 radian et une vitesse angulaire initiale égale à 0.0 radian/seconde. L'intégration est réalisée sur une période,  entre 0 et \( 2 \pi \), avec un pas temporel de 1e-4 secondes.

Le résultat est tout à fait convaincant. Il est donc possible d'assembler assez simplement dans un code C les moyens de calcul de la librairie GSL et les moyens d'affichage de Python.

Contenu et design par Dominique Lefebvre - www.tangenteX.com juillet 2017   Licence Creative Commons    Contact :

Cette œuvre est mise à disposition selon les termes de la Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Pas de Modification 3.0 France.