Application du Monte-Carlo symbolique avec Starlyx
Important: avant de traiter cet exercice voir le TP1
Le "Monte-Carlo symbolique" permet de remplacer la simulation Monte-Carlo , calcul d'une observable en fonction d'un paramètre, par une expression polynomiale de cette observable, fonction du paramètre.
Avec cette expression polynomiale il sera plus facile de calculer des sensibilités de l'observable par rapport au paramètre (dérivée) ou tout autre analyse avec des temps de calcul bien plus courts.
Quelques mots sur la méthode MC symbolique par projection sur une base de polynômes
Nous avons une observable exprimée par une équation intégrale des chemins $x$ sur le domaine des chemins $D_x$:
$ O(\alpha) = \int_D P_x(x,\alpha') f(x, \alpha') dx = \mathbb{E}(f(x, \alpha)) $
$\alpha$ est un paramètre que nous voulons faire varier entre $\alpha_{min}$ et $\alpha_{max}$: le symbole.
La densité de probabilité d'un chemin est : $P_x(x,\alpha)$
La fonction $f(x,\alpha')$ est le poids, échantillon de l'observable $O$ pour chaque chemin $x$.
Méthode de Monte Carlo: nous tirons les chemins $x_i$ selon la loi de probabilité $P_x(x,\alpha')$ et nous faisons la moyenne de f sur N tirages:
$ O(\alpha) \approx \frac{1}{N} \sum_{i=1}^N f(x_i, \alpha) $
Nous désirons trouver une approximation de $O(\alpha)$ par un polynôme de variable $\alpha$.
L'astuce est de déterminer les coordonnées de $O(\alpha)$ dans une base de polynômes orthonormée bien connue.
Nous choisissons à ce stade une base parmi les polynômes de Legendre, Tchebytchev, Laguerre, Fourier, Laplace.
Dans ce document, nous prenons les polynômes de Legendre, notés $P_n$.
Polynômes de Legendre
Produit scalaire entre 2 polynômes:
pour tous $n, m \in \mathbb{N} $
$ \langle P_n, P_m \rangle = \int_{-1}^{1} P_n(\alpha') P_m(\alpha') , w(\alpha') , d\alpha' = \delta_{nm} |P_n|^2 $
où
$ \delta_{nm} = \begin{cases} 1 & \text{si } n = m, \ 0 & \text{si } n \neq m. \end{cases} $
et
$w(\alpha')=1/2$
La norme est donnée par : $de de Monte Carlo: nous tirons les chemins xi selon la loi de probabilité Px(x,α′) et nous faisons la moyenne de f sur N tirages:
$ O(\alpha) \approx \frac{1}{N} \sum_{i=1}^N f(x_i, \alpha) $
$ |P_n|^2 = \int_{-1}^{1} (P_n(\alpha'))^2 w(\alpha') d\alpha'= \frac{2}{2n + 1} $
Les 3 premiers polynômes de Legendre :
$ P_0(\alpha') = 1 $
$ P_1(\alpha') = x $
$ P_2(\alpha') = \frac{1}{2}(3\alpha'^2 - 1) $
On définit la base orthonormale $\Psi_n$ à partir des polynômes de Legendre: [ \Psi_n = \frac{P_n}{|P_n|} = \sqrt{\frac{2n + 1}{2}} , P_n ]
Monte-Carlo symbolique
Ainsi nous pouvons exprimer $O(\alpha)$ comme une décomposition dans cette base de polynôme orthogonaux:
$ O(\alpha) = \sum_{n=0}^{\infty} a_n \Psi_n(\alpha) $
Maintenant, il faut trouver le moyen de calculer les coordonnées $a_n$.
Calcul des coefficients (a_n)
Les coefficients $a_n$ sont les coordonnées de $O(x,\alpha)$ par rapport à polynôme de la base, que nous calculons grâce au produit scalaire défini précédemment:
$ a_n = \langle O(\alpha), \Psi_n \rangle = \int_S \omega(\alpha') \Psi_n(\alpha') O(x, \alpha) , d\alpha' $
$ a_n = \int_S \omega(\alpha') \Psi_n(\alpha') \left[ \int_D P_x(x,\alpha) f(x, \alpha) dx \right] d\alpha' $
$ a_n = \int_S \int_D \omega(\alpha') \Psi_n(\alpha') P_x(x,\alpha) f(x, \alpha) dx d\alpha' $
On peut calculer cette intégrale par "double randomisation" en tirant d'abord $\alpha'$ dans $S$ selon la probabilité $\omega(\alpha')$ et ensuite un chemin $x(\alpha')$ selon la probabilité $P_x(x,\alpha)$. Le poids est donc $\Psi_n(x_i) f(x_i, \alpha)$
Nous pouvons donc calculer $a_n$ avec la série:
$ a_n = \frac{1}{N} \sum_{i=1}^{N} \Psi_n(x_i) f(x_i, \alpha) $
Application dans starlyx
Configurer le MC symbolique dans le fichier scene.xml de Starlyx
L'observable considérée dans le MC symbolique est définie dans un noeud XML au même niveau hierarchique que le noeud "scene".
On y précise la position de l'observable dans la scène, le nombre de réalisations (100000) utilisée pour calculer les coefficients de la projection sur la base de polynôme de Legendre et l'ordre maximum des polynômes.
Le "symbole" à faire varier est un fils du noeud "observable". Il est défini par l'attribut LSTAR du noeud "hg1" pour la longueur d'onde 650 nm. On fixe ses limites à VMIN=3mm et VMAX=30mm avec une répartition sur une échelle log.
<observable
SENSOR_NAME="sphereTdir"
FILTER="650"
NBR="100000"
NPOLMAX="4"
NORMAL_MC="1"
PLOT_NORMAL_MC="0"
PLOT_RECONSTRUCTION="1"
>
<symbol
OBJECT="hg1"
ATTRIBUTE="LSTAR"
LAMBDA="650"
VMIN="3"
VMAX="30"
LOG="1"
>
</symbol>
</observable>
Exercice
Configurer l'environnement:
cd /home/lafrier/environnement_tp_starlyx
source profile
Se mettre dans le dossier de ce TP:
cd tp_starlyx/TP3
Lancer starlyx avec l'option -s :
starlyx2 -s scene.xml
gnuplot va se lancer pour afficher le résultat: la courbe de variation de l'observable Tdir (signal de la sphère intégrante qui mesure la transmission directe) en fonction du paramètre LSTAR. On a la courbe donnée par quelques calculs de Monte-Carlo, pour vérification, et celle donnée par le polynome obtenu avec le Monte-Carlo symbolique. Voir les figures ci-bas.
Modifier le nombre de réalisations du calcul MCS à 1000 et relancer le calcul. Constater l'impact sur le polynôme obtenu.
Comparaison entre le polynôme obtenu par MCS et les simulations MC.
Avec 100000 réalisations:
Avec 1000 réalisations, le polynôme calculé "colle" moins bien:
