star-gs

Literate program for a geometric sensitivity calculation
git clone git://git.meso-star.fr/star-gs.git
Log | Files | Refs | README | LICENSE

commit 97eaecb25d16f5f24182fc238fdbdd4d6cb53600
parent ee7d372819d3c73418d76b0c8cfbc063f58356b0
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Sun, 23 Apr 2023 10:04:12 +0200

Annexe résultats et intitulés des blocs de code

Relecture à la phrase de l'annexe qui présente les scripts résultats
pour alléger des tournures et corriger des fautes. Enfin, à l'exeption
des noms de fichiers, l'entièreté des blocs de code débutent désormais
par une majuscule.

Diffstat:
Msrc/sgs_compute_sensitivity_translation.nw | 273++++++++++++++++++++++++++++++++++++++++---------------------------------------
1 file changed, 138 insertions(+), 135 deletions(-)

diff --git a/src/sgs_compute_sensitivity_translation.nw b/src/sgs_compute_sensitivity_translation.nw @@ -390,7 +390,7 @@ Dans cette section nous écrivons un algorithme de Monte Carlo pour résoudre le problème décrit en section~\ref{sec:probleme} à partir du modèle de sensibilité géométrique décrit en section~\ref{sec:modele_sensib}. Notre algorithme est construit de manière analogue avec pour support la mise en {\oe}uvre pratique -de sa [[<<fonction de réalisation>>]] en langage C~\cite{C_language}. Pour +de sa [[<<Fonction de réalisation>>]] en langage C~\cite{C_language}. Pour cela, nous échantillonnerons des chemins qui démarrent directement de la source de sensibilité, à savoir la paroi du haut, et qui la propage jusqu'au récepteur positionné sur la paroi du bas. La mise en {\oe}uvre ici proposée est néanmoins @@ -406,24 +406,24 @@ s'éviter un travail d'abstraction et ainsi faciliter l'écriture analogue de notre algorithme. Une démarche adaptée en raison de la configuration simplifiée à laquelle nous nous intéressons. -<<fonction de réalisation>>= +<<Fonction de réalisation>>= static res_T realisation (struct ssp_rng* rng, const struct sgs_scene* scene, double* w) { - <<données locales à la fonction de réalisation>> + <<Données locales à la fonction de réalisation>> res_T res = RES_OK; - <<initialiser le poids>> + <<Initialiser le poids>> - <<échantillonner un chemin du problème couplé>> - <<calcul du poids>> + <<Échantillonner un chemin du problème couplé>> + <<Calcul du poids>> exit: - <<nettoyer les données locales de la fonction>> - <<renvoyer le poids>> + <<Nettoyer les données locales de la fonction>> + <<Renvoyer le poids>> return res; } @ @@ -441,10 +441,10 @@ sensibilité $A_h$. Ce chemin est alors complété par l'échantillonnage d'un chemin de dérivée spatiale dont la couplage avec le chemin de sensibilité forme un chemin du problème couplé. -<<échantillonner un chemin du problème couplé>>= -<<échantillonner une position sur la source de sensibilité>> -<<échantillonner un chemin de sensibilité>> -<<échantillonner un chemin de dérivée spatiale>> +<<Échantillonner un chemin du problème couplé>>= +<<Échantillonner une position sur la source de sensibilité>> +<<Échantillonner un chemin de sensibilité>> +<<Échantillonner un chemin de dérivée spatiale>> @ Comme point de départ du chemin du problème couplé, on commence donc par @@ -456,7 +456,7 @@ surface que l'on vient d'échantillonner, dans notre cas la surface supérieure $A_h$ identifié dans le code par la constante [[SGS_SURFACE_Z_MAX]] (voir figure~\ref{fig:configuration}). -<<échantillonner une position sur la source de sensibilité>>= +<<Échantillonner une position sur la source de sensibilité>>= /* Échantillonner uniformément une position sur la source de sensibilité */ sgs_geometry_sample_sensitivity_source(scene->geom, rng, &frag); d3_set(pos_h, frag.position); @@ -476,7 +476,7 @@ lambertienne [[dir_emit_h]] autour de la normale [[normal_h]] de la surface $A_h$, et de lancer un rayon dans cette direction. Nous stockons alors dans [[hit0]] l'intersection de ce rayon avec la géométrie de la scène. -<<échantillonner un chemin de sensibilité>>= +<<Échantillonner un chemin de sensibilité>>= /* Échantillonner une direction d'émission de sensibilité */ ssp_ran_hemisphere_cos(rng, normal_h, dir_emit_h, NULL); /* Lancer le rayon qui propage l'émission de sensibilité */ @@ -497,7 +497,7 @@ spéculaire de la direction d'émission $\vec{\omega}$ ([[dir_emit_h]]) avant de lancer un rayon dans cette direction jusqu'à l'intersection [[hit1]] avec une surface . -<<échantillonner un chemin de dérivée spatiale>>= +<<Échantillonner un chemin de dérivée spatiale>>= /* Calculer la direction spéculaire à dir_emit_h */ reflect(dir_spec_h, dir_emit_h, normal_h); /* Tracer le rayon spéculaire partant de surf_A_h */ @@ -515,15 +515,15 @@ contribution \emph{nulle} si le chemin de sensibilité n'atteint pas le récepteur ou si le chemin de dérivé spatiale n'atteint pas la source radiative, à savoir la paroi de droite. -<<initialiser le poids>>= +<<Initialiser le poids>>= sensib = 0; @ -<<échantillonner un chemin de sensibilité>>= +<<Échantillonner un chemin de sensibilité>>= if(!hit_receiver(scene, pos_h, dir_emit_h, &hit0)) { goto exit; } @ -<<échantillonner un chemin de dérivée spatiale>>= +<<Échantillonner un chemin de dérivée spatiale>>= if(!hit_source(scene, pos_h, dir_spec_h, &hit1)) { goto exit; } @@ -537,7 +537,7 @@ dont on aura besoin pour le calcul du poids. De même, on initialise la variable [[dir_emit_d]] à $-\vec{\omega}_s$, cette direction nous sera également utile pour évaluer la source de dérivée spatiale. -<<échantillonner un chemin de dérivée spatiale>>= +<<Échantillonner un chemin de dérivée spatiale>>= d3_normalize(normal_d, hit1.normal); d3_minus(dir_emit_d, dir_spec_h); @ @@ -546,20 +546,20 @@ Dans notre problème couplé, la contribution du chemin, {\ie} le poids Monte Carlo, va s'exprimer à travers la condition à la limite de sensibilité (équation \ref{eq:clsensib}) et des sources de chacun de ses couplages. -<<calcul du poids>>= -<<décomposition du vecteur de déformation $\vec{\chi}$>> +<<Calcul du poids>>= +<<Décomposition du vecteur de déformation $\vec{\chi}$>> -<<calcul de la dérivée surfacique de $\rho$>> -<<calcul des sources de dérivées spatiales>> +<<Calcul de la dérivée surfacique de $\rho$>> +<<Calcul des sources de dérivées spatiales>> -<<calculer le poids de sensibilité>> +<<Calculer le poids de sensibilité>> @ La décomposition du vecteur de déformation $\vec{\chi}$ permet d'obtenir le vecteur tangent $\vec{u}$ nécessaire dans l'expression de la source de sensibilité et de ses dérivées surfaciques (équation \ref{eq:clsensib}). -<<décomposition du vecteur de déformation $\vec{\chi}$>>= +<<Décomposition du vecteur de déformation $\vec{\chi}$>>= decomposition(chi, normal_h, dir_emit_h, &proj_chi_h); @ @@ -569,8 +569,8 @@ revient à travailler dans le plan. Et cette dérivée surfacique est le produit scalaire entre le gradient surfacique de $\rho$ et la direction de dérivation $\vec{u}$ transformée dans ce plan ([[u_2d]]). -<<calcul de la dérivée surfacique de $\rho$>>= -<<récupérer le gradient surfacique de $\rho$>> +<<Calcul de la dérivée surfacique de $\rho$>>= +<<Récupérer le gradient surfacique de $\rho$>> /* Transformer u dans le plan XY */ u_2d[0] = proj_chi_h.u[X]; @@ -584,7 +584,7 @@ Pour récupérer $\rho$ et son gradient, il nous suffit de transformer dans le plan la position d'émission [[pos_h]], et d'interroger les données associées à la position ainsi transformée ([[pos_h_2d]]). -<<récupérer le gradient surfacique de $\rho$>>= +<<Récupérer le gradient surfacique de $\rho$>>= /* Transformer pos_h dans le plan XY */ pos_h_2d[0] = pos_h[X]; pos_h_2d[1] = pos_h[Y]; @@ -598,11 +598,11 @@ $\vec{\chi}$ et $\vec{u}$) dont les sources sont données par les équations \ref{eq:cl_duL_haut}, \ref{eq:cl_duL_droite}, \ref{eq:cl_dchiL_haut} et \ref{eq:cl_dchiL_droite}. -<<calcul des sources de dérivées spatiales>>= -<<décomposition du vecteur $\vec{\chi}$>> -<<décomposition du vecteur $\vec{u}$>> -<<calcul de la dérivée surfacique de $Sb$ dans la direction $\vec{u}_e$>> -<<calcul de la dérivée surfacique de $Sb$ dans la direction $\vec{u}_{cs}$>> +<<Calcul des sources de dérivées spatiales>>= +<<Décomposition du vecteur $\vec{\chi}$>> +<<Décomposition du vecteur $\vec{u}$>> +<<Calcul de la dérivée surfacique de $Sb$ dans la direction $\vec{u}_e$>> +<<Calcul de la dérivée surfacique de $Sb$ dans la direction $\vec{u}_{cs}$>> @ Sur la paroi de droite, la décomposition du vecteur de déformation $\vec{\chi}$ @@ -610,7 +610,7 @@ permet d'obtenir le vecteur tangent $\vec{u}$ nécessaire dans l'expression de la source de de la dérivée spatiale dans la direction $\vec{\chi}$ et de sa dérivée surfacique (équation \ref{eq:cl_dchiL_droite}). -<<décomposition du vecteur $\vec{\chi}$>>= +<<Décomposition du vecteur $\vec{\chi}$>>= decomposition(chi, normal_d, dir_emit_d, &proj_chi_d); @ @@ -619,7 +619,7 @@ vecteur tangent $\vec{u}_e$ nécessaire dans l'expression de la source de la dérivée spatiale dans la direction $\vec{u}$ et de sa dérivée surfacique (equation \ref{eq:cl_duL_droite}). -<<décomposition du vecteur $\vec{u}$>>= +<<Décomposition du vecteur $\vec{u}$>>= decomposition(proj_chi_h.u, normal_d, dir_emit_d, &proj_uh_d); @ @@ -635,7 +635,7 @@ $S_b$ qui nous sera utile pour le calcul du poids. Pour rappel, dans notre configuration les deux dérivées spatiales $\partial_{\vec{\chi}} I$ et $\partial_{\vec{u}} I$ partagent une même position d'émission ([[pos_d]]). -<<calcul de la dérivée surfacique de $Sb$ dans la direction $\vec{u}_e$>>= +<<Calcul de la dérivée surfacique de $Sb$ dans la direction $\vec{u}_e$>>= /* Calculer la position sur l'émetteur */ pos_d[X] = pos_h[X] + dir_spec_h[X]*hit1.distance; pos_d[Y] = pos_h[Y] + dir_spec_h[Y]*hit1.distance; @@ -662,7 +662,7 @@ $\vec{u}_{cs}$ ([[dSb_uchid]]) comme étant le produit scalaire entre le gradient surfacique de $S_b$ et la direction $\vec{u}_{cs}$ transformée dans le plan de la paroi de droite ([[u_chid_2d]]). -<<calcul de la dérivée surfacique de $Sb$ dans la direction $\vec{u}_{cs}$>>= +<<Calcul de la dérivée surfacique de $Sb$ dans la direction $\vec{u}_{cs}$>>= /* Transformer u_cs dans le plan YZ */ u_chid_2d[0] = proj_chi_d.u[Y]; u_chid_2d[1] = proj_chi_d.u[Z]; @@ -679,7 +679,7 @@ ensuite multipliée par la surface $A_h$ et l'angle $\pi$ ([[PI]]) étant donné l'échantillonnage des densités de probabilités de la position sur la surface et de la direction dans l'hémisphère sortante de la paroi. -<<calculer le poids de sensibilité>>= +<<Calculer le poids de sensibilité>>= /* Calcul de la contribution du chemin couplé */ Sb_sensib = - proj_chi_h.beta * d_rho * Sb @@ -690,7 +690,7 @@ Sb_sensib = sensib = Sb_sensib * PI * get_Sr(scene); @ -<<renvoyer le poids>>= +<<Renvoyer le poids>>= w[0] = sensib; @ @@ -707,7 +707,7 @@ stocke le triangle sur lequel se trouve l'origine du rayon, une donnée d'entré utilisée pour éviter une auto-intersection, c'est à dire l'intersection du rayon avec le triangle dont il est issu. -<<données locales à la fonction de réalisation>>= +<<Données locales à la fonction de réalisation>>= /* Macro utilisée comme sucre syntaxique */ #define TRACE_RAY(Org, Dir, StartFrom, Hit) { \ double range[2]; \ @@ -718,15 +718,15 @@ rayon avec le triangle dont il est issu. } (void)0 @ -<<nettoyer les données locales de la fonction>>= +<<Nettoyer les données locales de la fonction>>= #undef TRACE_RAY @ \paragraph{} Enfin, nous déclarons ci-après l'ensemble des variables locales nécessaires à -notre [[<<fonction de réalisation>>]]: +notre [[<<Fonction de réalisation>>]]: -<<données locales à la fonction de réalisation>>= +<<Données locales à la fonction de réalisation>>= /* Variables de la source de sensibilité */ double dir_emit_h[3]; double pos_h[3]; @@ -998,7 +998,7 @@ La fonction [[decomposition]] réalise cette décomposition et les résultats ([[alpha]], [[beta]] et [[u]]) sont répertoriés dans la structure [[projection]]. -<<fonctions utilitaires>>= +<<Fonctions utilitaires>>= struct projection { double alpha; double beta; @@ -1060,7 +1060,7 @@ Le poids ([[weight_flux_part_spec]]) correspondant à la partie du flux venant de la réflection sur la paroi spéculaire est donc calculé au même moment que celui de la sensibilité. -<<données locales à la fonction de réalisation>>= +<<Données locales à la fonction de réalisation>>= double weight_flux_part_spec; @ @@ -1069,11 +1069,11 @@ Le poids de cette contribution correspond à la source radiative $S_b$ (equation $\pi$ ([[PI]]), la surface $A_h$ et le coefficient de réflection $\rho$ ([[rho]]). -<<calcul du poids>>= +<<Calcul du poids>>= weight_flux_part_spec = Sb * rho * PI * get_Sr(scene); @ -<<renvoyer le poids>>= +<<Renvoyer le poids>>= w[1] = weight_flux_part_spec; @ @@ -1081,7 +1081,7 @@ Ce poids est également initialisé en même temps que celui de la sensibilité sorte que les chemins réfléchis n'atteignant pas le récepteur aient une contribution nulle. -<<initialiser le poids>>= +<<Initialiser le poids>>= weight_flux_part_spec = 0; @ @@ -1099,18 +1099,18 @@ compléter le programme écrit jusqu'ici et ainsi proposer une mise en {\oe}uvre complète de l'algorithme de Monte Carlo objet du présent document. Ces étapes sont regroupées dans la fonction qui suit: -<<calculer la sensibilité à la translation>>= +<<Calculer la sensibilité à la translation>>= res_T compute_sensitivity_translation(struct sgs* sgs) { - <<variables locales au calcul de sensibilité>> + <<Variables locales au calcul de sensibilité>> res_T res = RES_OK; - <<exécuter l'intégration Monte Carlo>> - <<afficher les résultats de l'estimation>> + <<Exécuter l'intégration Monte Carlo>> + <<Afficher les résultats de l'estimation>> exit: - <<libérer les variables locales au calcul de sensibilité>> + <<Libérer les variables locales au calcul de sensibilité>> return res; error: goto exit; @@ -1118,7 +1118,7 @@ error: @ Dans cette fonction on commence par exécuter l'intégration Monte Carlo. Cette -étape consiste à invoquer notre [[<<fonction de réalisation>>]] autant de fois +étape consiste à invoquer notre [[<<Fonction de réalisation>>]] autant de fois que de nombre de réalisations demandées, et à accumuler les poids qu'elle retourne. D'apparence triviale, cette simple boucle s'avère plus compliquée pour qui souhaite paralléliser son calcul. Au delà des questions propres à une @@ -1139,7 +1139,7 @@ réalisation ([[run_realisation]]), le type de poids calculés ([[smc_doubleN]], {\ie} un vecteur de réels) et le nombre total de réalisations ([[nrealisations]]) défini comme paramètre d'entrée via la variable [[sgs]]. -<<exécuter l'intégration Monte Carlo>>= +<<Exécuter l'intégration Monte Carlo>>= /* Configurer et créer le système Star-MonteCarlo */ smc_args.logger = &sgs->logger; smc_args.allocator = sgs->allocator; @@ -1166,7 +1166,7 @@ la translation et la fraction de la luminance qui dépend du paramètre {$\PI$}. Nous pouvons dès lors récupérer l'état de notre estimateur pour afficher l'espérance et l'écart type de ces variables. -<<afficher les résultats de l'estimation>>= +<<Afficher les résultats de l'estimation>>= res = smc_estimator_get_status(estimator, &status); if(res != RES_OK) goto error; @@ -1183,7 +1183,7 @@ Ne reste plus qu'à déclarer les variables locales utilisées par notre fonctio de calcul et de libérer en sortie l'espace mémoire allouée dynamiquement pour ces variables. -<<variables locales au calcul de sensibilité>>= +<<Variables locales au calcul de sensibilité>>= /* Système */ struct smc_device_create_args smc_args = SMC_DEVICE_CREATE_ARGS_DEFAULT; struct smc_device* smc = NULL; @@ -1197,24 +1197,24 @@ struct smc_estimator* estimator = NULL; struct smc_estimator_status status = SMC_ESTIMATOR_STATUS_NULL; @ -<<libérer les variables locales au calcul de sensibilité>>= +<<Libérer les variables locales au calcul de sensibilité>>= if(estimator) smc_estimator_ref_put(estimator); if(smc) smc_device_ref_put(smc); @ Le lecteur attentif aura remarqué que l'intégrateur utilise la fonction -[[run_realisation]] et non directement la [[<<fonction de réalisation>>]] -développée dans ce document (voir [[<<exécuter l'intégration Monte Carlo>>]]). +[[run_realisation]] et non directement la [[<<Fonction de réalisation>>]] +développée dans ce document (voir [[<<Exécuter l'intégration Monte Carlo>>]]). [[run_realisation]] est une fonction intermédiaire qui ne fait qu'appeler la -[[<<fonction de réalisation>>]]. [[run_realisation]] est donc parfaitement +[[<<Fonction de réalisation>>]]. [[run_realisation]] est donc parfaitement dispensable sauf à la bibliothèque \texttt{Star-MonteCarlo} qui nous impose la signature de la fonction à utiliser, c'est à dire le type des paramètres d'entrées et de sorties. En d'autres termes, l'utilisation de cette fonction intermédiaire nous permet de faciliter l'écriture et la lecture de la -[[<<fonction de réalisation>>]] en la libérant de contraintes fonctionnelles +[[<<Fonction de réalisation>>]] en la libérant de contraintes fonctionnelles imposées par \texttt{Star-MonteCarlo}. -<<fonctions utilitaires>>= +<<Fonctions utilitaires>>= static res_T realisation (struct ssp_rng* rng, @@ -1247,13 +1247,13 @@ Cette partie décrit la structure du fichier C dans lequel l'algorithme de calcul de sensibilité est mis en oeuvre: <<sgs\_compute\_sensitivity\_translation.c>>= -<<liste des inclusions>> +<<Liste des inclusions>> -<<constantes>> +<<Constantes>> -<<fonctions utilitaires>> -<<fonction de réalisation>> -<<calculer la sensibilité à la translation>> +<<Fonctions utilitaires>> +<<Fonction de réalisation>> +<<Calculer la sensibilité à la translation>> @ En plus des fonctions de calcul écrites dans les sections précédents, notre @@ -1264,7 +1264,7 @@ l'émission thermique de la paroi de droite (équation~\ref{eq:S_b}). Le calcul de leur gradient surfacique, utilisé pour évaluer leur dérivée surfacique (section~\ref{subsec:poids}), est également ajouté comme fonction utilitaire: -<<fonctions utilitaires>>= +<<Fonctions utilitaires>>= static double get_rho (const struct sgs_scene* scene, @@ -1327,7 +1327,7 @@ pour déterminer si un rayon a intersecté le récepteur (fonction suffit de s'assurer que le rayon intersecte le côté droit de la boîte, ici identifié par la constante [[SGS_SURFACE_X_MAX]]. -<<fonctions utilitaires>>= +<<Fonctions utilitaires>>= static int hit_source (const struct sgs_scene* scene, @@ -1354,7 +1354,7 @@ l'intersection se situe sur le récepteur lui même. Pour cela il nous suffit de tester si la position d'intersection appartient au sous domaine de la surface sur lequel le récepteur est défini. -<<fonctions utilitaires>>= +<<Fonctions utilitaires>>= static int hit_receiver (const struct sgs_scene* scene, @@ -1390,7 +1390,7 @@ hit_receiver Dernière fonction utilitaire à écrire, la fonction qui calcule la surface $A_h$ de la paroi spéculaire: -<<fonctions utilitaires>>= +<<Fonctions utilitaires>>= static double get_Sr(const struct sgs_scene* scene) { @@ -1403,14 +1403,14 @@ On complète notre fichier en définissant les constantes [[X]], [[Y]] et [[Z]] utilisées tout du long pour simplifier la lecture du code lors des accès aux éléments d'un vecteur. -<<constantes>>= +<<Constantes>>= enum {X, Y, Z}; @ Enfin, on liste en début des sources les fichiers d'en-tête utilisés par notre code. -<<liste des inclusions>>= +<<Liste des inclusions>>= #include "sgs_c.h" #include "sgs_geometry.h" #include "sgs_log.h" @@ -1429,13 +1429,13 @@ code. \section{Script d'exécution du calcul et résultats} \label{annex:scripts} -Nous écrivons ici les scripts \textit{shell} qui exécutent notre programme de -calcul de sensibilité en vue de sa validation croisée par différences finies. Le -premier script [[<<mc.sh>>]] lance plusieurs fois le programme dont le présent -document décrit la fonction de réalisation; l'objet étant de calculer à -différentes valeur de $\PI$ la sensibilité du flux $\varphi$ reçu par le capteur -(équation~\ref{eq:flux}) ainsi que la fraction de ce même flux qui dépend de -$\PI$ (annexe~\ref{annex:flux}). +Nous écrivons ici les scripts \textit{shell} qui exécutent le programme de +calcul de sensibilité en vue de sa validation croisée par différences finies. +Le premier script [[<<mc.sh>>]] lance plusieurs fois le programme dont le +présent document décrit la fonction de réalisation; l'objet étant de calculer à +différentes valeurs de $\PI$ la sensibilité du flux $\varphi$ reçu par le +capteur (équation~\ref{eq:flux}) ainsi que la fraction de ce même flux qui +dépend de $\PI$ (annexe~\ref{annex:flux}). <<mc.sh>>= #!/bin/sh -e @@ -1443,13 +1443,14 @@ $\PI$ (annexe~\ref{annex:flux}). # Estimation par Monte Carlo de la sensibilité et de la fraction du # flux qui dépend de pi -<<configuration géométrique>> -<<paramètres des calculs Monte Carlo>> +<<Configuration géométrique>> +<<Paramètres des calculs Monte Carlo>> -<<pour différentes valeurs de pi>> do - <<lancer un calcul Monte Carlo>> - <<post traiter le résultat>> - <<passer à la valeur de pi suivante>> +<<Pour différentes valeurs de pi>> +do + <<Lancer un calcul Monte Carlo>> + <<Post traiter le résultat>> + <<Passer à la valeur de pi suivante>> done @ @@ -1462,17 +1463,17 @@ coin inférieur ([[lower]]) et le coin supérieur ([[upper]]) d'un parallélépi aligné aux axes, deux variables passées en arguments du programme via l'option [[-b]]. -<<lancer un calcul Monte Carlo>>= +<<Lancer un calcul Monte Carlo>>= out=$(./sgs \ -n "${nrealisations}" \ -b low="${lower}":upp="${upper}":pi="${pi}") @ avec [[nrealisations]] le nombre de réalisations utilisées par le calcul Monte -Carlo. Si ce document décrit en détail les sources C de sa [[<<fonction de -réalisation>>]], nous renvoyons le lecteur vers les autres fichiers C et l'aide -du programme [[sgs]] (affichée via l'option [[-h]]) pour plus d'informations -quant aux fonctionnement et options du programme. +Carlo. Si ce document décrit en détail les sources C de sa +[[<<Fonction de réalisation>>]], nous renvoyons le lecteur vers les autres +fichiers C et l'aide du programme [[sgs]] (affichée via l'option [[-h]]) pour +plus d'informations quant aux fonctionnement et options du programme. Une fois le calcul terminé, en post-traiter le résultat se résume à afficher sur la sortie standard la valeur $\PI$ courante suivie de l'estimation et de l'écart @@ -1480,7 +1481,7 @@ type de la sensibilité et du flux que nous venons d'estimer. Ci-après nous utilisons la commande [[sed]] pour extraire et afficher ces valeurs stockées dans la variable [[out]] à l'issu de notre calcul Monte Carlo. -<<post traiter le résultat>>= +<<Post traiter le résultat>>= prompt="[^~]\{1,\}~ " estimation="[^[:blank:]]\{1,\}" error="[^\n$]\{1,\}" @@ -1495,12 +1496,12 @@ lancées [[nsteps]] fois pour différentes valeurs de $\PI$, la valeur de $\PI$ chaque itération étant simplement sa valeur précédente incrémentée d'un pas constant égal à [[pi_step]]. -<<pour différentes valeurs de pi>>= +<<Pour différentes valeurs de pi>>= i=0 pi=0 -while [ "${i}" -lt "${nsteps}" ]; +while [ "${i}" -lt "${nsteps}" ] @ -<<passer à la valeur de pi suivante>>= +<<Passer à la valeur de pi suivante>>= i=$((i + 1)) pi=$(printf "%s + %s\n" "${pi}" "${pi_step}" | bc) @ @@ -1510,12 +1511,12 @@ caractérisent notre configuration géométrique ainsi que les paramètres qui pilotent nos différents calculs, tels que le nombre de calculs à lancer ou encore le nombre de réalisations par calcul. -<<configuration géométrique>>= +<<Configuration géométrique>>= h=1 # Hauteur du parallélépipède lower="0,0,0" upper="1,1,${h}" @ -<<paramètres des calculs Monte Carlo>>= +<<Paramètres des calculs Monte Carlo>>= nrealisations=100000000 nsteps=14 pi_step=0.1 @@ -1526,11 +1527,11 @@ finies la sensibilité du flux à $\PI$ à partir des estimations Monte Carlo en sortie de [[<<mc.sh>>]], des résultats lus via l'entrée standard. Pour pouvoir comparer les résultats, ce second script recopie en sortie les sensibilités estimées par Monte Carlo ([[sen_mc]]) en plus d'écrire ces mêmes sensibilités -calculées par différences finies ([[sen_fd]]). Leur écarts types respectifs -([[err_mc]] et [[err_fd]]) sont également des données de sortie. On notera enfin -que l'ensemble des résultats seront adimentionnés aux dimensions de la -configuration géométrique et seront par conséquent donnés pour une valeur de -$\PI$ indépendante de la hauteur du parallélépipède ([[pi_over_h]]). +calculées cette fois par différences finies ([[sen_fd]]). Leur écarts types +respectif ([[err_mc]] et [[err_fd]]) sont également des données de sortie. On +notera enfin que l'ensemble des résultats seront adimentionnés aux dimensions +de la configuration géométrique et seront par conséquent donnés pour une valeur +de $\PI$ indépendante de la hauteur du parallélépipède ([[pi_over_h]]). <<fd.sh>>= #!/bin/sh -e @@ -1541,44 +1542,44 @@ $\PI$ indépendante de la hauteur du parallélépipède ([[pi_over_h]]). # Liste des données pour chaque ligne écrite en sortie printf "pi_over_h sen_mc err_mc sen_fd err_fd\n" -<<fonctions utilitaires à [[fd.sh]]>> +<<Fonctions utilitaires à [[fd.sh]]>> -<<lire les paramètres d'entrée>> +<<Lire les paramètres d'entrée>> -<<pour chaque valeur de pi considérée>> do - <<lire la valeur du flux autour de pi>> - <<calculer par différences finies la sensibilité du flux à pi>> - <<lire l'estimation Monte Carlo de la sensibilité du flux à pi>> - <<écrire les résultats adimensionner aux dimensions du système>> - <<passer au calcul suivant>> +<<Pour chaque valeur de pi considérée>> +do + <<Lire la valeur du flux autour de pi>> + <<Calculer par différences finies la sensibilité du flux à pi>> + <<Lire l'estimation Monte Carlo de la sensibilité du flux à pi>> + <<Écrire les résultats adimensionner aux dimensions du système>> + <<Passer au calcul suivant>> done @ Pour les calculs en différences finies nous avons besoin de connaître les paramètres utilisés par les estimations Monte Carlo, comme le $\delta$ entre chaque valeur de $\PI$ ([[pi_step]]) ou encore la hauteur du parallélépipède -nécessaire pour adimensionner nos résultats ([[h]]). Pour passer ces données -d'un script à l'autre, on modifie le script [[mc.sh]] pour les écrire en -en-tête de ses sorties, en-tête qui peut alors être lu par le script -[[fd.sh]]. +nécessaire pour adimensionner les résultats ([[h]]). Pour passer ces données +d'un script à l'autre, on ajoute aux sorties de [[mc.sh]] un en-tête contenant +ces paramètres, en-tête qui peut alors être lu par le script [[fd.sh]]. -<<paramètres des calculs Monte Carlo>>= +<<Paramètres des calculs Monte Carlo>>= printf "%s %s\n" "${h}" "${pi_step}" @ -<<lire les paramètres d'entrée>>= +<<Lire les paramètres d'entrée>>= read -r header h=$(echo "${header}" | cut -d' ' -f1) pi_step=$(echo "${header}" | cut -d' ' -f2) @ -Autre donnée nécessaire pour adimensionner nos résultats, la valeur maximale du +Autre donnée nécessaire pour adimensionner les résultats, la valeur maximale du flux. Celle-ci correspond à la valeur de $\varphi$ lorsque la boîte est fermée, -c'est à dire lorsque $\PI = 0$ soit la valeur de $\PI$ du premier résultat en -sortie de [[mc.sh]]. Après l'en-tête nous lisons donc ce résultat que l'on -stocke dans la variable [[p]] avant d'en extraire la valeur $\varphi$ -([[phi_max]]). +c'est à dire lorsque $\PI$ vaut $0$, soit la valeur de $\PI$ du premier calcul +Monte Carlo lancé par [[mc.sh]]. Après avoir lu l'en-tête de ses sorties nous +lisons donc ce premier résultat que l'on stocke dans la variable [[p]] avant +d'en extraire la valeur $\varphi$ ([[phi_max]]). -<<lire les paramètres d'entrée>>= +<<Lire les paramètres d'entrée>>= read -r p phi_max=$(echo "${p}" | cut -d' ' -f4 | float_to_bc) @ @@ -1589,11 +1590,11 @@ résultat Monte Carlo pour la valeur de $\PI$ courante est stocké dans la variable [[c]] tandis que le résultat précédent et suivant sont respectivement stockés dans les variables [[p]] et [[n]]. -<<pour chaque valeur de pi considérée>>= +<<Pour chaque valeur de pi considérée>>= read -r c -while read -r n; +while read -r n @ -<<passer au calcul suivant>>= +<<Passer au calcul suivant>>= p="${c}" c="${n}" @ @@ -1602,28 +1603,30 @@ Il suffit alors d'extraire le flux ([[phi]]) et son erreur ([[err]]) aux valeurs de $\PI$ précédente ([[p]]) et suivante ([[n]]) pour calculer par différences finies la sensibilité et l'erreur associée pour la valeur de $\PI$ courante: -<<lire la valeur du flux autour de pi>>= +<<Lire la valeur du flux autour de pi>>= phi_p=$(echo "${p}" | cut -d' ' -f4 | float_to_bc) err_p=$(echo "${p}" | cut -d' ' -f5 | float_to_bc) phi_n=$(echo "${n}" | cut -d' ' -f4 | float_to_bc) err_n=$(echo "${n}" | cut -d' ' -f5 | float_to_bc) @ -<<calculer par différences finies la sensibilité du flux à pi>>= +<<Calculer par différences finies la sensibilité du flux à pi>>= sen_fd=$(echo "(${phi_n}-${phi_p})/(2*${pi_step})" | bc_cmd) err_fd=$(echo "(${err_n}+${err_p})/(2*${pi_step})" | bc_cmd) @ -On utilise alors le résultat Monte Carlo au $\PI$ considéré pour en retrouver -non seulement sa valeur mais aussi pour extraire l'estimation Monte Carlo de la -sensibilité de $\varphi$ à $\PI$. Ne reste plus qu'à écrire l'ensemble des -résultats attendus. +On utilise alors le résultat Monte Carlo au $\PI$ considéré pour retrouver +non seulement la valeur de $\PI$ mais aussi pour extraire l'estimation Monte Carlo de la +sensibilité de $\varphi$ à $\PI$: -<<lire l'estimation Monte Carlo de la sensibilité du flux à pi>>= +<<Lire l'estimation Monte Carlo de la sensibilité du flux à pi>>= pi=$(echo "${c}" | cut -d' ' -f1 | float_to_bc) sen_mc=$(echo "${c}" | cut -d' ' -f2 | float_to_bc) err_mc=$(echo "${c}" | cut -d' ' -f3 | float_to_bc) @ -<<écrire les résultats adimensionner aux dimensions du système>>= + +Ne reste plus qu'à écrire l'ensemble des résultats attendus: + +<<Écrire les résultats adimensionner aux dimensions du système>>= pi_over_h=$(echo "${pi}/${h}" | bc_cmd) sen_mc=$(echo "${sen_mc}/${phi_max}*${h}" | bc_cmd) err_mc=$(echo "${err_mc}/${phi_max}*${h}" | bc_cmd) @@ -1639,7 +1642,7 @@ flottante dans le format attendu par le programme \texttt{bc}, et la fonction [[bc_cmd]] qui exécute le programme \texttt{bc} et en post-traite le résultat pour supprimer les zéros qui suivent le dernier chiffre significatif. -<<fonctions utilitaires à [[fd.sh]]>>= +<<Fonctions utilitaires à [[fd.sh]]>>= float_to_bc() { sed 's/\([+-]\{0,1\}[0-9]\{0,\}\.\{0,1\}[0-9]\{1,\}\)'\