%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Al1 %%%%%
% LOG BOOK 5 -> Mise a jour 102 avec sensy_to_var, set_probe, free_parameters
% et aussi : fonction power, fonctions de grille.
% et corrections Kalman pour les obs.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass[a4paper,11pt]{article}
\usepackage{latexsym}
\usepackage{amsmath}
% amsfonts permet le \mathbb(C) attention, peu de lettres disponibles
\usepackage{amsfonts}
\usepackage{eepic,psfig,pstricks,multirow}
\usepackage{floatflt}
\usepackage{fancyheadings}
\usepackage{footnpag}
\usepackage{rotating}
\usepackage{euscript}
\usepackage{tabularx,xspace}
%\usepackage{minitoc}
\usepackage{makeidx}
\usepackage[frenchb]{babel}
\usepackage[latin1]{inputenc}

\topmargin -1cm
\addtolength{\textheight}{4cm}
\tolerance=8000
\pretolerance=4000
\textwidth 15cm
\oddsidemargin 0.5cm

%
\newcommand{\cad}{c'est-\`a-dire\xspace}
\newcommand{\Minik}{Mini\_ker\xspace}
\newcommand{\bra}[1]{<{#1}\mid}
\newcommand{\ket}[1]{\mid{#1}>}
\newcommand{\mat}[1]{\underline{\underline{\mathcal{#1}}}}
\newcommand{\diag}[1]{{#1}^{d}}
\newcommand{\mean}[1]{\overline{{#1}}}
\newcommand{\scal}[2]{{<\, #1 \, #2>}}
\newcommand{\scalm}[3]{{<\, #1 {#2}^{\dagger} \, #2 #3>}}
\newcommand{\qrayl}[3]{\frac{\scalm{#1}{#2}{#3}}{\scal{#1}{#3}}}
\newcommand{\trac}[1]{{\mathcal{T}r (#1)}}

\makeindex
\begin{document}

\def\todec{{\overline {\overline {\Upsilon  }}}_{,dec} }
\def\tox{{\overline {\overline {\Upsilon  }}}_{,x} }
\def\dt{{\delta t}}
\def\edtaux{{exp (-{\delta t  \over \taux } ) }}
\def\deta{{\delta \eta } }
\def\Deta{{\Delta \eta } }
\def\Dt{{\delta t } }
\def\dfi{{\delta \varphi } }
\def\Dfi{{\Delta \varphi } }
\def\dts{{\delta t\over 2 } }
\def\dtaf{{\overline {\overline {\partial _\eta f }}} }
\def\d(#1)G{{\overline {\overline {\partial_#1 G}}} }
\def\dA{{\overline {\overline {\partial A}}} }
\def\dtA{\partial_t A}
\def\VD{{\mathcal V \mathcal D}}


\title{Manuel pratique de Mini-ker\author{Al1, Janvier 2007\footnote{avec
des contributions de pat, stepH , stb, seb, jlj, jyg et ribs}}\date{Versions 1.02.xx}}
\maketitle 
{\it Ce manuel est une introduction à l'utilisation pratique de l'environnement
\Minik dans sa version 1.02. Le lecteur est renvoyé aux
différents documents concernant le TEF\index{TEF}, l'adjoint et les sensibilités
sur la page {\bf ZOOM}\index{ZOOM} en note de bas de page 2. Signalons également
l'existence d'un cahier de programmation --- sur demande}.

\vspace{0.2cm}
{\bf La version 1.02 comporte les facilités suivantes :}
\begin{itemize}
\item passage immédiat en double précision à la demande;
\item calcul des sensibilités en avant des variables à une liste de paramètres\footnote{On devrait alors dans
la foulée proposer en standard la méthode de Newton pour offrir une alternative
au fit de paramètres avec Minuit et ainsi réhabiliter LLSQ --- Least Linear Square fit par méthode de
Householder, routine du {\it Paris Sud Informatique} Université de Paris-Sud.}\index{quadratique}
\\ ({\bf free\_parameter}); plusieurs exemple de fit par cette méthode ont été effectuées (cf rapport
Céline de Mai 2007, et les exos Evex et ClimSI).
\item calcul des sensibilités de la fonction de coût à la même liste de paramètres\\ ({\bf free\_parameter, Zback});
\item un vecteur des ``observables'' étendant le vecteur des transferts \\ ({\bf set\_probe$<$});
\item le calcul de sensibilité en avant des variables à plusieurs types de perturbation\\ ({\bf Sensy\_to\_var$<$}).
Ceci permet notamment de calculer le propagateur et les exposants de Lyapunov le long de la trajectoire (cf exo de stepH),
mais aussi le gain de rétroaction non linéaire sous la forme réponse : $\displaystyle \frac{1}{1-g}(t,0)$;
\item calcul de la sensibilité de la matrice d'avance d'état {\bf dPi\_aspha} au {\bf fwd} paramètre
Ce calcul (symbolique) proposé par stb permet ensuite de calculer la sensibilité au paramètre choisi, soit
des valeurs propres, soit des valeurs singulières, de diverses matrices dérivées des Jacobiennes;
\item utilisation d'un filtre de Kalman sur une base de données (Pat); le filtrage a été
amélioré~(cf cahiers Propagateurs) pour assurer le mieux possible la positivité stricte de l'avance
de la var-covar;
\item recherche d'optimum paramétriques ou de loi de commande avec Minuit\footnote{\copyright CERN}.
\item calcul des fonctions SLT (gain, réponses) en Borel par balayage en dt ou par  la méthode des éléments
propres et inversion en variable temporelle (cf routines Bsweep et Boreleig).
\item dimensionnement automatique (ex DimEtaPhi générée);
\end{itemize}

\vspace{0.2cm}
{\bf Préannonce de la 2.00}

\hspace{0.5cm}Dans l'optique d'un développement avec composante industrielle, on n'envisage pas de version
1.03, mais directement celle d'une maquette (2.00) comportant un seul emboîtement de familles - à la ZOOM.

Ceci devrait être l'objet de discussions dans les mois à venir, avec ou sans l'accord
de l'ANR á laquelle nous avons soumis le projet {\bf Dyams}.

%
\newpage
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Principes et structure du Mini\_ker}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Un modèle formulé sous TEF sépare deux types d'équations : a) des équations
d'évolution référant aux modèles de {\bf cellules}\index{cellules} et : b) des modèles dits de
{\bf transferts}\index{transferts} exprimant des contraintes statiques. De façon générale, ces objets
associent à chaque variable une équation (un modèle) et ces objets sont rangés en
familles formant une structure
d'arbre\footnote{cf http:/www.lmd.jussieu.fr/documents.dir/tef-GB-partA5.ps.gz}.
Le traitement complet des modèles sous TEF renvoie à l'environnement logiciel
{\bf ZOOM}\footnote{cf http:/www.lmd.jussieu.fr/zoom}.

{\bf \Minik} constitue une
entrée plus facile à la modélisation sous TEF en simplifiant la structure
du modèle : une seule cellule, un seul transfert, une seule famille.
Alors que ZOOM peut traiter de l'ordre de 100 000 équations, \Minik
n'en peut traiter pratiquement que quelques centaines.

%Les vecteurs d'état\index{Etat} et de transferts sont deux tableaux Fortran {\bf eta(.)} et
%{\bf ff(.)}. Les modèles sont donnés pour chaque variable en utilisant
%l'instruction {\bf f\_set\index{f\_set}: name=expression;}.
%
\Minik effectue automatiquement les calculs de sensibilité aux conditions initiales
ou à un paramètre\index{paramètre} du modèle (en avant). Le modèle adjoint est automatiquement généré,
qui permet (en arrière) de calculer la sensibilité d'une fonctionnelle de coût
aux variables du modèle et à ses paramètres.

Toutes les dérivées partielles
nécessaires à la résolution du système (et des sensibilités et de l'adjoint) par le
TEF\index{TEF} sont effectuées symboliquement à la précompilationi.
L'utilisateur a donc terminé la programmation de son modèle
dès qu'il a entré les fonctions décrivant le modèle et donné
la valeur des paramètres ainsi que les conditions initiales.
%%% en{\bf eta(.)}\footnote{les transferts sont déterminés par les états.}.
La programmation pratique est décrite dans ce manuel en première section.

Le langage développé pour automatiser le développement du TEF à partir
des fonctions du modèle utilise {\bf Mortran}\index{Mortran}, un pré-compilateur Fortran développé
dans les années 1970 (cf manuel Mortran).
%%% au SLAC puis au CERN, et enfin par la présente collaboration.
La première annexe en explique le principe d'utilisation pour le \Minik

En section~\ref{analyse}, on indique comment se servir des différents outils
d'analyse des systèmes dynamiques traités, qui exploitent la formulation
par le TEF du modèle\index{TEF}. L'utilisation de Minuit en complèment de l'adjoint
est présentée en section~\ref{minuit}, et l'utilisation de données est étendue avec le
filtre de Kalman~(\ref{kalman}) et l'option {\bf ZObs}.
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{L'entrée d'un modèle dans \Minik}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Deux niveaux de langage étant disponibles, nous décrivons d'abord le
langage le plus proche de la résolution car il utilise une correspondance
immédiate entre les vecteurs $\eta, \varphi$ du TEF et les tableaux Fortran correspondants
{\bf eta(.)} et~{\bf ff(.)}.
Le deuxième niveau est alors introduit, qui permet une écriture plus
proche du modèle symbolique, et donne également accès aux facilités
pour mettre en \oe uvre un maillage uni-dimensionnel.
%---------------------------------------------
\subsection{Principe d'organisation du code}
%
Le programme principal de \Minik est subdivisé en morceaux, que l'on apelle séquences.
L'utilisateur peut se contenter d'un accès à un nombre réduit de séquences dans lesquelles
il code, en pseudo-Fortran, les instructions de dimentions du modèle et son expression
mathématique. Lors
de l'assemblage, ces instructions sont traduites par Mortran en Fortran à l'aide
de macros\index{macros (MTN)} qui remplacent, dans le corps des calculs, des instructions symboliques
par du code piloté par la description du modèle.

Les séquences principales sont rangées dans un fichier. Nous utilisons ici
implicitement le langage de {\bf cmz}\index{CMZ (CERN)}, logiciel de gestion de sources et de fabrication
d'exécutables du CERN\footnote{une version sans cmz, portable sur Windaube et Mac OSX
a été développée par Pat, se référer à la version anglaise du manuel.}.
Ainsi, une séquence est un élément de code formant
une {\bf +KEEP} de cmz, équivalent à un fichier file.h du préprocesseur Fortran.
Les séquences à programmer sont {\bf \$dimetaphi} pour entrer les dimentions\index{dimension}
des tableaux eta(.) et ff(.), et {\bf \$Zinit} pour programmer le modèle.
Deux autres séquences utiles ({\bf \$Zsteer, \$ZStep}) sont introduites à la fin de la boucle
d'avance temporelle
du programme principal\index{ZSteer}\index{ZStep} et utilisées pour effectuer des test et calculs numériques suuplémentaires,
gérer le niveau d'impression, interrompre la simulation etc.
%---------------------------------------------
\subsection{Description basique du modèle}
%
On démarre avec une cmz file du nom de son exercice (exemple donné : {\bf predator}
- en général, on part d'un des nombreux exemples fournis que l'on va modifier).
Sous CMZ, on va dans le patch {\bf ZinProc}\footnote{{\bf cd zinproc} puis {\bf ls}
pour lister le contenu du patch.}.
%%% {\bf ls} donne la liste des séquences.
Un modèle formulé sous le TEF
conduit au système matriciel :
\begin{eqnarray}
\partial_t \eta (t) &=&  g(\eta(t),\varphi(t)) \nonumber\\
\varphi(t) &=& f(\eta(t),\varphi(t))
\label{tef}
\end{eqnarray}
Le modèle prédateur-proie de Lotka-Volterra\index{Lotka-Volterra} peut s'écrire sous la forme suivante :
\begin{equation}
\left\{
\begin{array}{cc}
\partial_t \eta _{prey} (t) =&  a \eta _{prey} - a \varphi _{meet} \\
\partial_t \eta _{pred} (t) =&  -c \eta _{pred} + c \varphi _{meet} \nonumber\\
\end{array}
\right.
\end{equation}
\begin{equation}
\varphi _{meet} = \eta _{prey}\eta _{pred}
\label{pred}
\end{equation}
dans laquelle on a deux équations d'état\index{Etat}, \cad une équation d'évolution pour
chaque groupe proie et prédateur, et un seul transfert\index{transferts} représentant les rencontres
entre individus de groupes adverses.
%
\subsubsection*{séquence dimetaphi}
%
La séquence où on doit rentrer les dimensions\index{dimension} des vecteurs eta(.) et ff(.)
se trouve dans {\bf \$dimetaphi}, il s'agit respectivement des {\bf parameter np}\index{parameter FTN}
et {\bf mp}.
%* \footnote{maxstep ne concerne que l'adjoint, décrit plus loin. De
%* même les paramètres du maillage 1D.}.
\footnotesize
\begin{verbatim}
+KEEP,dimetaphi.
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!        nbre de variables d'etat       : np             +
!        nbre de variables de transfert : mp             +
!        nbre de parametres libres      : lp             +
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      parameter (np=2,mp=1);
      parameter (lp=0);
! _______________________________________________________
+KEND.
\end{verbatim}
\normalsize
Pour cet exemple, on a deux états et un seul transfert (et pas de paramètre libre)\index{paramètre}.

Pour que \Minik retienne ce calcul manuel, il faut à présent lever dans sa
selseq.kumac\footnote{fournie avec les fichiers \Minik, cf ftp anonyme pub/alain.} le
drapeau dimetaphi ({\bf sel dimetaphi}). Sinon, il est inutile de remplir cette séquence, le dimensionnement
sera automatiquement effectué par \Minik. 
%
\subsubsection*{séquence ZINIT}
%
Le modèle se programme dans {\bf \$Zinit}. On rentre d'abord les paramètres du
modèle\index{Lotka-Volterra} :
%\small
\footnotesize
\begin{verbatim}
!%%%%%%%%%%%%%%%%%%%%%%
! Parameters                               "! = carte commentaires         "
!%%%%%%%%%%%%%%%%%%%%%%

     dt=.01;                               "pas de temps initial           "
     nstep=10 000;                         "nbre d'itérations              "
     modzprint = 1000;                     "fréquence des printouts        " 
;
     time=0.;                              "initialisation du temps modèle "
;
! parametres du modele                    "le langage est du pseudo Fortran"
     apar = 1.5;                           "avec des fins d'intructions    "
     cpar = 0.7;                           "signalées par ';'              "

print*,'***************************************';
print*,'Lotka-Volterra model with parameters as:';
z_pr: apar,bpar;
print*,'***************************************';
\end{verbatim}
\normalsize
Une carte débutant par le signe ``!'' est un commentaire, les instructions Mortran\index{Mortran}
sont de format libre et terminées par un point-virgule. On a défini le pas
de temps initial $dt$, qui s'ajoutera à la variable {\bf time}\index{time, time-step} du temps de la
simulation à la fin du premier pas de temps\footnote{on fera varier ce pas de temps dans {\bf ZSTEER}}.
\index{ZSteer}Les paramètres\index{paramètre} du modèles sont aussi donnés\index{modzprint, Zprint}.

On remarque l'instruction d'impression {\bf z\_pr}
dont la macro Mortran\index{macros (MTN)} est en tête de la séquence ZINIT dans la cmz\index{CMZ (CERN)} file de chaque
exemple\footnote{ceci pour suggérer à l'utilisateur qu'il est très simple de se mettre
aussi a faire ses propres macros Mortran.}.
Les instructions ci-dessus produisent l'impression\index{paramètre}~:
\footnotesize
\begin{verbatim}
 ***************************************
 Lotka-Volterra model with parameters as:
    apar,cpar :
      1.5  0.699999988
 ***************************************
\end{verbatim}
\normalsize
On définit à présent le modèle, \cad ici une seule composante ff(1),
et deux d'états.
On programme directement les expressions à droite du signe égal des équations~(\ref{pred}),
composante par composante, en utilisant la macro {\bf Fun\_set}\index{Fun\_set} ({\it function set}):
\footnotesize
\begin{verbatim}
!%%%%%%%%%%%%%%%%%%%%%%
! Transfer definition
!%%%%%%%%%%%%%%%%%%%%%%
! ff(1) : rencontres
    fun_set Phi_tef(1) = eta(1)*eta(2);                  "modèle de transfert"
!%%%%%%%%%%%%%%%%%%%%%%
! Cell definition
!%%%%%%%%%%%%%%%%%%%%%%
! eta(1) : pray
! eta(2) : predator                                     "modèles de cellule"
    fun_set  deta_tef(1) =   apar*eta(1) - apar*ff(1);
    fun_set  deta_tef(2) = - cpar*eta(2) + cpar*ff(1);

!     Etat initial
!     ------------
     eta(1) = 1.;
     eta(2) = 1.;
;
    OPEN(50,FILE='title.tex',STATUS='UNKNOWN');   "utilisé par SLTC,SLTCIRC"
    write(50,5000) apar,cpar;                                "et Minuik    "
5000;format('Lotka-Volterra par:',2F4.1);
\end{verbatim}
\normalsize
Les états\index{cellules}\index{SLT, SLTC} sont aussi initialisés\index{Lotka-Volterra}, et on provoque l'écriture du titre choisi
dans le fichier {\bf title.tex} qui est obligatoire, car il sert à transmettre
de l'information aux sorties graphiques.

Les noms de fonctions {\bf dEta\_tef} et {\bf Phi\_tef}\footnote{la
convention majuscules-minuscules indistinctes est celle du Fortran.}
sont réservés, leur numérotation est en correspondance avec les
vecteurs eta(.) et ff(.).
La simulation sans calculs de sensibilité ou de l'adjoint est
prête\footnote{\emph{Il est également possible de programmer soi-même toutes
les matrices Jacobiennes\index{Jacobienne}
et vecteurs du TEF\index{TEF} en entrant dans son zinproc les séquences D\_Tr, D\_Mx, A\_Mx,
B\_Mx, Ct\_Mx, Omega\_vec et Gamma\_vec à la manière du \Minik originel
de stb. On attend d'avoir des exemples
pour décrire ce mode de programmation un peu près du filet.}}. En ce qui concerne le symbole des
noms de variables introduites par l'utilisateur, et pour éviter les symboles communs avec
ceux du reste du code, on a adopté la convention d'interdire le caractère ``o'' dans
les deux premiers caractères des noms de variables Fortran.
\vspace{0.5cm}

{\it {\bf Remarque sur l'initialisation des transferts :} dans le cas où le  modèle
comporte une matrice $D$ (\cad qu'il y a une relation implicite entre des
composantes de transferts) il peut être utile voire nécessaire d'initialiser
quelques transferts. Cet aspect ne regarde que la toute première itération de la résolution
du système des transferts. On a intéret à faciliter cette première passe en initialisant
les transferts-paramètres\index{paramètre} ainsi que les variables de transferts, qui prendront
sinon la valeur nulle comme première estimation. Cela peut donc conduire
à problème numérique si l'inverse d'une variable apparait dans les formules,
ou plus généralement rend une fonction singulière.

 On aura intérêt en cas de problème
à lever le drapeau {\bf sel debug} dans la selseq.kumac\index{kumac (macro CMZ)}\index{selseq (macro CMZ)}
 pour : a) recupérer en sortie
une trace de la D-Loop\index{D\_loop}, et b) activer un calcul d'estimation de conditionnement de
la matrice $D$ avant inversion (par LaPack)\index{LaPack}. On pourra également visualiser
la variation itérative des transferts dans le fichier {\bf tr\_in\_dloop.data}. L'option
{\bf debug} est donc très bavarde, ne pas oublier de réduire {\bf nstep}}.
%
\subsection{Lancement de la simulation et sorties}
%-------------------------------------------------
On lance la précompilation, compilation et fabrication de l'exécutable
en tapant {\bf mod} dans la fenêtre CMZ\index{CMZ (CERN)}. L'éxécution de cette macro\index{macros (MTN)} va assembler les
parties du code fragmenté en séquences, pour former un programme principal en langage Mortan,
que l'on trouve dans, avec l'exemple predator : predator/predator/cfs/predator\_o.tmp.
On retrouve dans ce fichier ses séquences dimetaphi et zinit, zsteer\index{ZSteer} en fin de la 
boucle de calcul, immergées dans le programme {\bf principal}.
L'assemblage est guidé par l'appel de mod à une macro cmz : {\bf selseq.kumac}, située
dans le répertoire de la cmz file. C'est dans cette kumac que l'on donne le path
d'accès au code complet de \Minik, qui impose les séquences de l'utilisateur, et
qui aiguille les options diverses décrites plus loin.

On a par ailleurs
le résultat du passage de Mortran\index{Mortran} dans le fichier : predator/predator/cfs/principal.f .
Après {\bf mod} sous cmz, un exécutable doit être présent sous le répertoire
{\bf predator/predator} par exemple. On lance cet éxécutable 

predator.exe $>$res.lis

On obtient en sortie des fichiers .data\index{data file} et le listing (res.lis).
Chaque fichier {\bf .data} comporte une première colonne donnant la variable {\bf time}\index{time, time-step}
\footnote{{\bf dres.data} comporte une deuxième colonne donnant {\bf dtime} qui a pu
varier en cours de simulation.}. Les colonnes suivantes  donnent
les valeurs de {\bf eta(.)} dans res.data, {\bf dEta(.)} dans dres.data,
la variation de eta(.) à chaque pas de temps, {\bf ff(.)} dans tr.data.
On a enfin sur sélection de l'option {\bf sel doteta} dans la selseq.kumac\index{Vel dot\_eta}
un fichier doteta.data donnant la vitesse $\partial_t \eta$.
D'autres fichiers {\bf .data} seront décrits plus loin. 

De manière analogue à l'option modzprint / Zprint, l'écriture des fichiers .data peut être
commandée par le couple modzout / Zout (cf {\bf ZSteer})\index{modzout, Zout}.
%
\newpage 
\subsection{Description symbolique du modèle}
%--------------------------------------------
Une couche de langage supplémentaire permet une description plus symbolique
d'un modèle\footnote{sur la suggestion de jlj.}, mais qui résultera exactement dans les mêmes vecteurs et fonctions
de base précédentes (et un programme Fortran identique).
L'exemple déjà donné peut en effet s'écrire de façon équivalente\index{transferts} :
\footnotesize
\begin{verbatim}
!%%%%%%%%%%%%%%%%%%%%%%
! Transfer definition
!%%%%%%%%%%%%%%%%%%%%%%
! rencontres
set_Phi
< var: ff_interact, fun: f_interact = eta_pray*eta_pred;
>;

!%%%%%%%%%%%%%%%%%%%%%%
! Cell definition
!%%%%%%%%%%%%%%%%%%%%%%

set_eta
< var: eta_pray, fun: deta_pray =   apar*eta_pray - apar*ff_interact;
  var: eta_pred, fun: deta_pred = - cpar*eta_pred + cpar*ff_interact;
>;
\end{verbatim}
\normalsize
Ces instructions génèrent en Fortran les instructions d'impression de la correspondance
entre les noms symboliques et les vecteurs et fonctions de base. Dans cet
exemple, on aura ainsi en sortie :
\footnotesize
\begin{verbatim}
  ---------------- Informing on Phi definition -----------------
    Var-name,           Function-name,       index in ff vector
             ff_interact              f_interact  1
  ----------------------------------------------------

  --------------- Informing on Eta definition ------------------
   Var-name,           Function-name,       index in eta vector
                eta_pray               deta_pray  1
                eta_pred               deta_pred  2
  ----------------------------------------------------
\end{verbatim}
\normalsize
Enfin une dernière simplification d'écriture est possible\index{set\_eta, set\_phi} :
\footnotesize
\begin{verbatim}
!%%%%%%%%%%%%%%%%%%%%%%
! Transfer definition
!%%%%%%%%%%%%%%%%%%%%%%
! rencontres
set_Phi
< eqn: ff_interact = eta_pray*eta_pred;
>;

!%%%%%%%%%%%%%%%%%%%%%%
! Cell definition
!%%%%%%%%%%%%%%%%%%%%%%

set_eta
< eqn: eta_pray =   apar*eta_pray - apar*ff_interact;
  eqn: eta_pred = - cpar*eta_pred + cpar*ff_interact;
>;
\end{verbatim}
\normalsize
où les fonctions associées auront un nom implicite {\bf {\$}ff\_interact} etc.
\newpage
%
\subsection{L'Options Grid1D}
%---------------------------------- 
Des macros\index{macros (MTN)} spécifiques ont été définies pour rentrer simplement un modèle
maillé, limité à une seule dimension\index{dimension}. En plus de la définition
générique des modèles attachés aux mailles, il est nécessaire en général
de définir des conditions aux limites, ou des sources ou forçages de part et d'autre
du maillage.

Nous présentons l'exemple d'une chaîne de maselottes\footnote{l'exemple est fournit dans {\bf masselottes}.} reliées entre elles par
des ressorts et amortisseurs, limitée à gauche par un mur et libre à droite.
La formulation TEF\index{TEF} du problème s'écrit dans l'espace de phase
(position, velocité) du n\oe ud $k$ :
\begin{equation}
\left\{
\begin{array}{cc}
\partial_t \eta _{k} ^{pos}  =  \eta _{k} ^{vel}\qquad& \\
\partial_t \eta _{k} ^{vel}  = ( \varphi_k ^{spr} -\varphi _{k+1} ^{spr}&+\,\,\varphi _{k} ^{dmp}-\varphi _{k+1} ^{dmp})\,/m_k \nonumber\\
\end{array}
\right.
\end{equation}
\begin{equation}
\left\{
\begin{array}{cc}
\varphi_k ^{spr} = -k_k (\eta _{k} ^{pos}- \eta _{k-1} ^{pos})\\
\varphi_k ^{spr} = -d_k (\eta _{k} ^{vel}- \eta _{k-1} ^{vel})
\end{array}
\right.
\label{mass}
\end{equation}
où $m_k$ est la masse du n\oe ud $k$, $r_k$ et $d_k$ la rigidité des ressorts et les
coefficients d'amortissement.
%
\subsubsection*{séquence ZINIT}
%
La convention de numérotage des mailles (pour cet exemple) est décrite en commentaire
du code suivant\index{set\_node}\index{node}\index{modzprint, Zprint}~:
\footnotesize
\begin{verbatim}
!%%%%%%%%%%%%%%%%%%%%%%
! Parameters
!%%%%%%%%%%%%%%%%%%%%%%
dimension rk(n_node),rd(n_node),rmassm1(n_node);
data rk/n_node*1./;
data rd/n_node*0.1/;
data rmassm1/n_node*1./;
     dt=.01;       time=0.;
     nstep=5 000;  modzprint = 1000;
!%%%%%%%%%%%%%%%%%%%%%%%%%%================================================
! Maillage convention inode
!%%%%%%%%%%%%%%%%%%%%%%%%%%                                 Open ended
!(2 Lo      Phi    Eta                         (n_node)
! Eta)  \|       .-----.       .-----.          .-----.        /
! wall  \|-\/\/\-|     |-\/\/\-|     |  . . .  -|     |-\/\/\- |dummy
!  pos  \|--***--|  1  |--***--|  2  |  . . .  -|  n  |--***-- |Phis
! speed \|   1   |_____|   2   |_____|      n   |_____|  n+1   \(2 Hi Phi)
! Dwn                                                           Up
!%%%%%%%%%%%%%%%%%%%%%%====================================================
! Cell definition
!%%%%%%%%%%%%%%%%%%%%%%

! LO cells (wall)
! ----------------
set_dwn_eta
< var: eta_pos_wall,  fun: deta_pos_wall  = 0.;
  var: eta_speed_wall, fun: deta_speed_wall= 0.;
>;
! node cells
! ----------
set_[3]node_Eta
< var: eta_move(inode),  fun: deta_move(inode) = eta_speed(inode);
  var: eta_speed(inode),
  fun: deta_speed(inode) = rmassm1(inode)
                             *( - ff_spring(inode+1) + ff_spring(inode)
                               - ff_dump(inode+1)  + ff_dump(inode)    );
>;
!---------------
! Initialisation
!---------------
eta_pos_wall = 0; eta_speed_wall = 0.;
<inode=1,n_node; eta_move(inode)=0.01; eta_speed(inode)=0.0;>;
;
Z_Pr: eta;
!%%%%%%%%%%%%%%%%%%%%%%
! Transfer definition
!%%%%%%%%%%%%%%%%%%%%%%
! Hi limiting Transfers : dummy ones
! ----------------------------------
set_Up_Phi
< var:ff_dummy_1, fun: f_dummy_1=0.;
  var:ff_dummy_2, fun: f_dummy_2=0.;
>;
! node transfers
! --------------
! convention de signe spring : comprime:= +
set_[3]node_Phi
< var: ff_spring(.),
  fun:
    f_spring(inode)= - rk(inode)*(eta_move(inode) - eta_move(inode-1));
  var: ff_dump(.),
  fun:
    f_dump(inode)  = - rd(inode)*(eta_speed(inode) - eta_speed(inode-1));
>;
\end{verbatim}
\normalsize
Comme on le voit sur l'exemple, le modèle attaché à une maille
est générique et peut nécessiter la création de conditions limites
particulières (transferts bidons à droite). Dans l'exemple, on a choisi mêmes nombres de
mailles de type cellule\index{cellules} ou transfert (c'est le [3] du set\_node).
Les noms de variables et
fonctions sont libres mais doivent impérativement être différents.
On dénomme une variable du maillage en une maille $k$ par le symbole
{\bf (inode:k)}, où $k$ est un chiffre ou un paramètre Fortran, ou encore une
opération sur des variables définies.
Le symbole {\bf inode}\index{node} est
réservé. Des instruction Fortran habituelles peuvent être effectuées
dans le bloc Mortran\index{Mortran} {\bf $<\,\,\, >$} des set\_xx\_yy.

{\bf paramètres NGrid, MGrid}\hspace{1.2cm}
Pour paramètrer les nombres de mailles, on utilisera d'une part l'instruction
{\bf parameter (NGrid=3,MGrid=3);} avec les deux formes {\bf set\_xxx}~:
%
\begin{verbatim}
! node cells
! ----------
set_[NGrid]node_Eta
< 
>;
! node transfers
! --------------
set_[MGrid]node_Phi
< 
>;
\end{verbatim}
%
La forme de l'instruction, comme le nom des variables {\bf NGrid, MGrid} sont strictes.
%
%--------------------------------
\subsubsection*{écriture simplifiée}
%--------------------------------
Lorsque l'utilisateur ne souhaite pas associer de nom à la
fonction définissant le modèle correspondant à une ou plusieurs composantes, il peut
remplacer une instruction :
\footnotesize
\begin{verbatim}
  var: ff_dump(.),
  fun: f_dump(inode)  = - rd(inode)*(eta_speed(inode) - eta_speed(inode-1));
\end{verbatim}
\normalsize
par :
\footnotesize
\begin{verbatim}
   eqn: ff_dump(.) = - rd(inode)*(eta_speed(inode) - eta_speed(inode-1));
\end{verbatim}
\normalsize
Dans ce cas, la fonction prendra le nom de la variable précédée du caractère\index{node}
\$ : {\bf \$ff\_dump(inode)}\footnote{L'initialisation éventuelle des transferts\index{transferts}
se fera en analogie avec celle des états\index{Etat}, l'indice de boucle étant impérativement {\bf inode}.}.
%--------------------------------
\subsubsection*{vérifications}
%--------------------------------
%
Il est possible d'utiliser conjointement les instructions de base
{\bf f\_set}\index{f\_set}, en prenant garde toutefois aux composantes définies
par les {\bf set\_}.
La correspondance avec les vecteurs et fonctions de base est donnée
en impression d'exécution, où sont précisées aussi les mailles à droite (Dwn)
et à gauche (Up) du maillage pour vérification :
\footnotesize
\begin{verbatim}
  --------------- Informing on Dwn Eta definition ---------------
 Var-name,            Function-name, index in eta vector
         eta_pos_wall        deta_pos_wall  1 [
       eta_speed_wall      deta_speed_wall  2 [

  -------------- Informing on Eta Nodes definition --------------
 Var-name,     Function, k2index of (inode: 0 [ 1,...n_node ] n_node+1)
            eta_move           deta_move    1 [   3 ...   7 ]   9
           eta_speed          deta_speed    2 [   4 ...   8 ]  10

  ---------------- Informing on Up  Phi  definition -------------
 Var-name,             Function-name, index in ff vector
          ff_dummy_1           f_dummy_1 ]    7
          ff_dummy_2           f_dummy_2 ]    8
         ff_move_sum          f_move_sum ]    9
        ff_speed_sum         f_speed_sum ]   10
  ----------------------------------------------------

 -------------- Informing on Phi Nodes definition ---------------
 Var-name,     Function, k2index of (inode: 0 [ 1,...m_node ] m_node+1)
           ff_spring            f_spring   -1 [   1 ...   5 ]   7
             ff_dump              f_dump    0 [   2 ...   6 ]   8
  ----------------------------------------------------
\end{verbatim}
\normalsize
En bref, on peut vérifier que les indices inférieurs nuls ou négatifs ne
sont pas utilisés en limite d'équation de maille courante.
Il est cependant
conseillé de vérifier les cohérences avec les tableaux eta(.) et ff(.)
dans ``principal.f'', c'est distrayant et aussi rassurant.
%
% -------------------------------------------
\subsubsection*{l'instruction {\bf mult\_sum}}
% -------------------------------------------
L'impression précédente correspond en réalité à une modification de l'exemple
donné, avec ajout de deux composantes de transfert :
\footnotesize
\begin{verbatim}
set_up_Phi
< var:ff_dummy_1, fun: f_dummy_1=0.;
  var:ff_dummy_2, fun: f_dummy_2=0.;
  var: ff_move_sum, fun: f_move_sum = (mult_sum:eta_move(inode:[.i]),[i]=1,3);
  var: ff_speed_sum,
                fun: f_speed_sum =  ff_spring(inode:1)*eta_speed(inode:2);
>;
\end{verbatim}
\normalsize
pour montrer\index{mult\_sum} l'utilisation de la macro\index{macros (MTN)} {\bf (mult\_sum: ... ,[i]=ideb,ifin)}.
où les arguments ideb et ifin sont des entiers chiffrés
 (ni paramètre ni variable FTN). Dans le cas présent, on obtient la somme explicite
des trois termes. Cette macro s'utilise également dans tout contexte de définition
de fonction, soit dans une {\bf Fun\_set:}, soit dans les blocs {\bf set\_xxx}, \cad hors
contexte Grid1D (cf in~\ref{probe}).

L'expression d'une mult\_sum -- obligatoirement entre parenthèses -- peut être combinée à
une expression plus complète, mais la forme de l'opérande doit rester simple
du point de vue du nombre de parenthèses\footnote{On peut sur commande demander
plus, ou se reporter à une autre macro {\bf mult\_op} plus complète, cf paragraphe~\ref{fset}.}.

% -------------------------------------------
\subsubsection*{séquence DimEtaPhi}
% -------------------------------------------
%
Nous donnons d'abord les détails de l'option d'écriture manuelle de la séquence {\bf dimetaphi},
en rappelant qu'elle est optionnelle et choisie en fonction de {\bf sel dimetaphi}
dans la selseq.kumac. En mode de génération automatique, le [3] du set\_node dans
l'exemple précédent indique
alors le nombre de mailles; cette indication est ignorée en mode manuel.

Les dimentions\index{dimension} des six groupes de variables doivent impérativement être
données dans la séquence \$DimEtaPhi (exemple non modifié)\index{node} :
\begin{verbatim}
+SELF,IF=GRID1D.
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! nodes parameters, and Limiting Conditions (Low and High)
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      parameter (n_node=3,n_dwn=2,n_up=0,n_mult=2);
      parameter (m_node=3,m_ldwn=0,m_up=2,m_mult=2);
! ________________________________________________________
\end{verbatim}
\normalsize
On a $n\_node=m\_node=3\,\,$ mailles\index{parameter FTN} de type respectivement cellule\index{cellules} et
transfert, deux cellules de conditions limites à gauche (Dwn),
aucune à droite (Up), le contraire en transferts.
La multiplicité est le nombre de variables
attachées à chaque maille cellule transfert ($n\_mult, m\_mult$), qui correspondra à une
boucle {\bf do io=1,n\_node} en Fortran dans le code généré.

Pour que le code généré prenne en compte les spécificités du
maillage, on doit lever le drapeau {\bf Grid1D}\index{Grid 1D}. Ceci
se fait dans la séquence {\bf selseq.kumac}\index{kumac (macro CMZ)}\index{selseq (macro CMZ)} du répertoire oú se trouve
la cmz file\index{CMZ (CERN)}~: l'instruction {\bf sel Grid1D} lève ce drapeau (elle est commentée
dans la selseq livrée avec les exos).

{\bf Version 102} Pour continuer à utiliser la séquence de dimentionnement manuel, il faut
délectionner le drapeau dimetaphi dans la selseq.kumac (sel dimetaphi). Sinon, les instruction Fortran
parameter seront automatiquement calculés et générés par \Minik, en prenant les nombres de
mailles indiqués dans les set\_[3]node.
\newpage

%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{nouveaux transferts : les observables}
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\label{probe}
%
La version {\bf 102} voit apparaître un nouveau vecteur ($\mu$), déclaré ainsi~:
\begin{verbatim}
set_probe
< var: p_sum, fun: f_sum= eta_pray + eta_pred;
>;
\end{verbatim}
% 
Il s'agit d'un transfert\index{observables}\index{transferts} de type particulier, qui n'intervient pas dans les couplages avec
les autres variables. Ils représentent typiquement un modèle d'appareil de mesure.

Leur intervention dans le code de {\bf principal} est marquée par {\bf ff(mp+[.j]}. Le fichier
des sorties s'apelle {\bf obs.data}, et il faut donner une dimension à {\bf parameter (mobs=0)} dans
la {\bf dimetaphi} si nécessaire.

Bien sûr, la fonction de coût est sensible à ce vecteur dans l'adjoint.

Lorsque le filtre de Kalamn est utilisé (avec sel kalman dans la selseq), ces variables seront
toutes utilisées -- et elles seules -- pour calculer le ``vecteur innovation''
(cf~\ref{kalman})\index{Kalman}.

\smallskip
On peut utiliser la macro {\bf mult\_sum} dans tout contexte, comme ici\index{mult\_sum}~:
\footnotesize
\begin{verbatim}
Set_Probe
<
 eqn: mu_dist=(eta(1)+5.02832E+00)**2 + (eta(2)+7.43784E+00)**2
               + (eta(3)- 1.95888E+01)**2;
 var: sum_xyz, fun: somme_eta = (mult_sum: eta[i],[i]=1,3);
!* ou encore  var: sum_xyz, fun: somme_eta = mult_sum: eta[i],[i]=1,3;
!* ou bien    var: sum_xyz, fun: somme_eta = ftn_var*(mult_sum: eta[i],[i]=1,3);
>;
\end{verbatim}
\normalsize
qui génère le code Fortran suivant, par exemple dans le calcul de cette fonction~:
\begin{verbatim}
c ---*-----------------------------------------------------------------I
10108 ff(mp+1)=(((eta(1)+5.02832E+00)**2+(eta(2)+7.43784E+00)**2+(eta(3)
     *-1.95888E+01)**2))
      ff(mp+2)=((eta(1)+eta(2)+eta(3)))
c ---*-----------------------------------------------------------------I
\end{verbatim}
%
{\bf Il est impératif de placer entre parenthèses l'expression concernée par la somme implicite.}


{\bf 102} Il est à présent plus commode d'utiliser des fonctions du genre précédent
codées en Fortran. Tel est le cas de la fonction {\bf quad\_dist} qui permet dd calculer
la distance quadratique\index{quadratique} entre deux points de l'espace des états ou des transferts.
Exemple~:
%
\begin{verbatim}
set_probe
< eqn: mu_dist = quad_dist(np,eta(1),eta_ref(1));
>;
\end{verbatim}
%
l'observable {\bf mu\_dist} donne la distance quadratique entre l'état courant et
un point de l'espace de phase pris comme référent. On peut y remplacer en argument
{\bf eta(1)} par un nom de variable symbolique (ou un transfer), mais celà n'évite pas
d'avoir en tête le vecteur $\eta$ (ou $\varphi$ et sa traduction en tableau Fortran\footnote{le principe de
dérivation de ce type de fonction intégrale est donné dans la section~\ref{derivation}.}.

On utilise par exemple cette fonction
pour évaluer la période (apparente\footnote{cette question débattue au Zircon de
l'Alexandrienne a particulièrement intéressé Céline, à savoir~: un système apparemment
périodique est-il périodique ? C'est à ces hauteurs de questionnement qu'on reconnait
l'unsecte.}) d'un système.
%
\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Options graphiques}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\label{graph}
Les différents fichiers {\bf .data}\index{data file} sont accessibles aux grapheurs courants.
Ainsi, les fichiers peuvent être visualisés avec {\bf gnuplot}\index{gnuplot}
en entrant :
%
\begin{verbatim}
		plot "res.data" using 1:(n+1) for eta(n)
		plot "tr.data" using 1:(n+1) for ff(n)
\end{verbatim}
où $n$ est un numéro de composante des vecteurs du TEF. Une premiere ligne indique les noms
de variables associées aux colonnes dans le format gnuplot.

Pour les fans de {\bf PAW}, logiciel graphique\index{PAW}
du CERN, \Minik génère les kumacs permettant la lecture des .data et leur 
structure en {\bf n-tuples}. Se reporter au manuel PAW pour plus
d'information. Les n-tuples ne sont cependant prêts à l'emploi
que pour les vecteurs de dimension 10 au plus (y-compris la variable $time$)\index{time, time-step}.
On édite ses kumacs pour les adapter au format souhaité - et on change
leurs noms, car elles seront écrasées au prochain run.

\medskip
{\it Cette option est en réalité obsolete, la génération des kumacs n'ayant pas été remise
à jour, faute d'utilisateur. On constate que non seulement la version ntuple de PAW est
buggée, mais que gnuplot est d'un usage plus simple. Cependant, pour faire des histogrammes
et des surfaces, PAW reste nettement supérieur. On conseille alors l'écriture d'un code
Fortran pour générer les histogrammes nécessaires et ensuite les visualiser (surfaces incluses)
dans PAW. On peut s'inspirer utilement de l'exemple fourni avec les \Minik sources
dans {\bf gains.f} (avec go.xqt). La procédure à suivre dans PAW est automatiquement
générée à l'éxécution dans la {\bf borel.kumac}}.
\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Des sorties listing d'aide}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\label{aide}
On a déjà vu l'aide fournie pour le maillage. Pour les modèles non-maillés, on
a un résumé des vecteurs en listing d'exécution, comme par exemple pour Lorenz-63\footnote{fourni, cf~{\bf lorhcl}.} :
%
\begin{verbatim}
  ---------------- Informing on Phi definition -----------------
    Var-name,           Function-name,       index in ff vector
            ff_courant_L           $ff_courant_L  1
               ff_T_czcx              $ff_T_czcx  2
            ff_Psi_Tczcx           $ff_Psi_Tczcx  3
              ff_Psi_Tsz             $ff_Psi_Tsz  4
  ----------------------------------------------------

  --------------- Informing on Eta definition ------------------
   Var-name,           Function-name,       index in eta vector
           eta_courant_L               deta_conv  1
              eta_T_czcx          deta_energy_zx  2
                eta_T_sz           deta_energy_z  3
  ----------------------------------------------------
\end{verbatim}
%
%\newpage
Un résumé des équations est également disponible dans le fichier {\bf Model.hlp},
comme le suivant, pour le même Lorenz :
%
\begin{verbatim}
==================== set_Phi                                                                                                                                                         
 1 ff_courant_L  $ff_courant_L   eta_courant_L
 2 ff_T_czcx     $ff_T_czcx      eta_T_czcx
 3 ff_Psi_Tczcx  $ff_Psi_Tczcx   eta_T_czcx*eta_courant_L
 4 ff_Psi_Tsz    $ff_Psi_Tsz     eta_T_sz*eta_courant_L

==================== set_Eta                                                                                                                                                 
 1 eta_courant_L deta_conv       pi_prandtl*ff_T_czcx-pi_prandtl*eta_courant_L
 2 eta_T_czcx    deta_energy_zx  -ff_Psi_Tsz+pi_Rayleigh_ratio*ff_courant_L-eta_T_czcx
 3 eta_T_sz      deta_energy_z   ff_Psi_Tczcx-pi_wave_nb_ratio*eta_T_sz
\end{verbatim}
%
les fonctions définies par des {\bf F\_set} ne figurent dans ce fichier que si
elles sont déclarées par {\bf Fun\_set}. On peut encore abonder l'information
grâce à la macro {\bf Aux\_fun}, de même format,  qui ne génère rien pour \Minik mais est
écrite dans {\bf Model.hlp}\index{Fun\_set}\index{Aux\_fun}. Exemples~:
%
\begin{verbatim}
 Fun_set no_w_neg_0 = (1.-ostepin(-w(0)-dt*(adv_w(0)+buoy_w(0)+sol_w(0))+dm2,2.*dm2))
 Aux_fun dm2 = 0.01;
\end{verbatim}
%

\vspace{1.cm}
Dans la mesure du possible toutes les valeurs numériques sont repérées par un nom de variable.
Ainsi, en listing d'exec :
%
\begin{verbatim}
 Gamma :-8.19100E-02-1.42151E-01 3.87150E-02
         eta_courant eta_T_czcx  eta_T_sz   
       ------------------------------------------------
 Omega : 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00
         ff_courant_ ff_T_czcx   ff_Psi_Tczc ff_Psi_Tsz 
       ------------------------------------------------
\end{verbatim}
pour les vecteurs connus par exemple, ou pour la matrice de couplage~:
\begin{verbatim}
 >ker : Matrice de couplage       4 4 4 4
ff_courant_L Raw(  1,j=1,  4):   1.000     -9.9010E-03   0.000       0.000    
ff_T_czcx    Raw(  2,j=1,  4): -2.7972E-02   1.000       0.000      9.9900E-04
ff_Psi_Tczcx Raw(  3,j=1,  4):  0.1605      9.7359E-02   1.000     -5.7321E-03
ff_Psi_Tsz   Raw(  4,j=1,  4):   0.000     -0.1376      5.7225E-03   1.000    
             Var-Name          ff_courant_ ff_T_czcx   ff_Psi_Tczc ff_Psi_Tsz 
             ----------------------------------------------------------
\end{verbatim}
%
comme pour les Jacobiennes et autres matrices de \Minik,
de même encore comme déjà indiqué en tête des fichiers {\bf data}.

%\newpage
 Les grosses matrices\index{matrices} bénéficient
d'un traitement particulier : elles ne sont pas données en tableaux de chiffres
mais sous la forme de ``scatter plots''. Nous commentons ici un exemple extrait
d'UrbaGus~:\footnote{thèse de Gus)}
%\scriptsize
\begin{verbatim}
 ------ TEF Matrices and vectors in ker ---------
 ------------------------------------------------
  A_matrix in ker                 48 48 48 48       >>>   2229 zero values  <<<
 [MIN:-0.4252    (32,33) (@zyxwvutsrqpi, .1PQRSTUVWXYZ*)  214.0    (43,43):MAX]
    123456789012345678901234567890123456789012345678
  1|P                                               || hous(1st)               
  2|,1,                                             || rent(1st)               
  3|  1                                             || densite(1st)            
  4|   P                                            || hous(.)                 
  5|   ,1,                                          || rent(.)                 
  6|     1                                          || densite(.)              
  7|      P                                         || hous(.)                 
  8|      ,1,                                       || rent(.)                 
  9|        1                                       || densite(.)              
 10|         P                                      || hous(.)                 
 11|         ,1,                                    || rent(.)                 
 12|           1                                    || densite(.)              
tableau tronqué
 37|                                    P           || hous(.)                 
 38|                                     1,         || rent(.)                 
 39|                                      1         || densite(.)              
 40|                                       P        || hous(.)                 
 41|                                        1,      || rent(.)                 
 42|                                         1      || densite(.)              
 43|                                          *     || hous(last)              
 44|                                          ,1,   || rent(last)              
 45|                                            1   || densite(last)           
 46|                                             1  || hous_dum                
 47|                                              1 || rent_dum                
 48|                                               1|| densite_dum             
    123456789012345678901234567890123456789012345678
    .    .    .    .    .    .    .    .    .    . densite_dum             
    .    .    .    .    .    .    .    .    .    .rent_dum                
    .    .    .    .    .    .    .    .    .    hous_dum                
    .    .    .    .    .    .    .    .    .   densite(last)           
    .    .    .    .    .    .    .    .    .  rent(last)              
    .    .    .    .    .    .    .    .    . hous(last)              
sortie tronquée
    .    .  densite(.)              
    .    . rent(.)                 
    .    .hous(.)                 
    .    densite(.)              
    .   rent(.)                 
    .  hous(.)                 
    . densite(1st)            
    .rent(1st)               
    hous(1st)               

\end{verbatim}
\normalsize
Une première ligne donne la taille de la matrice (deux indices par dimention~: le premier
la partie imprimée, la deuxième la taille réelle de la matrice dans sa ``Dimension'' Fortran).
On donne aussi le nombre réréré de coefficients strictement nuls.

La deuxième ligne indique le code alphabétique corespondant aux valeurs des coefficients
de la matrice, encadrés par les valeurs min et max référées. Ici, la plus petite valeur
est de $-0.4252$ en position $(32,33)$ (@ dans le tableau), le max $=214.$ en position $(43,43)$
(* dans le tableau).

 Les minuscules
codent des valeurs négatives, les majuscules des positives. Enfin, pour respecter le TEF,
on code systématiquement la valeur $1$ par $1$ (et $-1$ par $i$), l'intervalle
max par $*$ (et $@$ le min si négatif).
 Enfin la valeur stricte $0$ est codée par un blanc, les valeurs inférieures
à $1.$ en valeur absolues sont codées par ``.'' lorsqu'elles sont positives, et par ``,''
si négatives.
Entre 1. et max (resp. -1. et min si nécessaire) on code 12 intervalles de P,Q,R,...,  à *.

On me dit ``c'est bien'', mais ``pourquoi les poils ?''. Les poils
servent à repérer les variables associées à chaque coefficient. Et comme ici dans le cas de modèle
maillé, la variable est indicée par  (1st) et (last), pour la première et dernière maille,
les mailles intermédiaires par (.).
%
\subsection*{Test de linéarité sur $\varphi$}
% ..........................................
Lors de la fin du pas de temps, qui fournit les transferts, un test numérique est effectué
par comparaison de $\varphi=f(\eta(t-dt)+d\varphi$ ({\bf ffl}), avec $d\varphi$
donné par la résolution de $d\eta, d\varphi$ (dans {\bf ker}) et
$\varphi=f(\eta,\varphi,t)$  ({\bf ff}) :
%
\begin{verbatim}
! ========================================================
! test linearite ffl - ff
! ========================================================
if (istep.gt.1)
< res=0.; <io=1,m; res = res +(ffl(io)-ff(io))**2/max(one,ff(io)*ff(io)); >;
  if (res .gt. TOL_FFL)
  < print*,'*** pb linearite : res > TOL_FFL a istep',istep,res,' > ',TOL_FFL;
    do io=1,m < z_pr: io,ff(io),ff(io)-ffl(io); >;
  >;
>;
\end{verbatim}
%
Comme on le voit, il s'agit d'un critère portant sur la somme des écarts quadratiques entre
valeurs relatives (ou non). Ce test n'englobe ainsi que les non-linéarités des modèles
de transfert. On peut utiliser la valeur de {\bf res} dans {\bf ZSTEER} pour
piloter la valeur de {\bf dt}, voire, revenir au pas de temps précédent\index{ffl (test de linéarité)}.
Les plus grands écarts de la liste sont signalés par $<<<$.

Il est aussi possible de définir une valeur du critère {\bf TOL\_FFL} dans {\bf ZINIT}
différente du défaut fixé à $10^{-3}$ indépendamment de la précision Fortran.
%
\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Principes utiles à la programmation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
\subsection{Comment rendre un modèle continu}
% ..........................................
\label{fset}
\footnote{cette section concerne donc aussi les instructions {\bf set\_eta< et semblables.}}On
traite d'un problème lié à la continuité des fonctions et de leurs dérivées partielles.
Comme exemple, on peut avoir une fonction step de la variable, comme
dans les modèles avec changement de phase, avec des codes en Fortran comprenant
un {\bf IF} sur la variable de température\index{modèle {\bf IF}}. Dans ce cas, on régularise le problème
en remplaçant la step par une ``smooth step''. On peut noter à ce sujet qu'il s'agit
d'une prise en compte implicite de l'échelle spatio-temporelle. A une échelle
macroscopique, la fonction step est juste un non-sens\index{Heaviside (fonction de)}.

On peut avantageusement utiliser
la fonction {\bf Heavyside}\footnote{je précise pour éviter les propagations
malencontreuses que ``Heavyside'' est un joke, pour dire une fonction de Heaviside ``lourde'',
\cad à la fois inerte et miraculeuse.}
de la //mini\_ker/mathlib :
\begin{verbatim}
        Delta = -1."K";
        A_Ice =  heavyside("in:" (T_K-Tf), Delta, "out:" dAIce_dT);
\end{verbatim}
dans l'exemple, Tf est la température de fusion de la glace, A\_ice donne la fraction
englacée d'un volume d'eau à la température T\_k. La fonction smooth-step est une quasi
tangente hyperbolique fonction de la variable adimentionnelle $x/\Delta$), normalisée de 0 à 1, avec
une pente maximum de 2.5, voire figure~\ref{heavy}.
%% FIG : HeavySide Function
\begin{figure}[h]
\psfig{figure=heavyside.ps,%
bbllx=60pt,bblly=180pt,bburx=526pt,bbury=650pt,height=9cm,width=11cm,clip=}
\caption{La fonction ``domptée'' de Heaviside et sa dérivée pour une variable adimentionnée}
\label{heavy}
\end{figure}
%
Pour avoir la dérivée, il suffit d'ajouter (en tête de zinit) la règle
de dérivation :
\begin{verbatim}
&'(HEAVYSIDE(#,#,#))(/#)' = '((#1)(/#4)*HEAVYDELTA(#1,#2,#3))'
\end{verbatim}
(cette règle fait à présent partie des macros\index{macros (MTN)} de \Minik, de même la fonction et
son point d'entrée (HeavyDelta).

Un autre type de problème consiste en la traduction de $var=min(f(x),g(x))$.
Dans ce cas, on ne veut pas dériver la smooth-step, et on écrira simplement
\footnotesize
\begin{verbatim}
var = HeavySide(f(x)-g(x),Delta,dum)*g(x) + (1.-HeavySide(f(x)-g(x),Delta,dum)*f(x);
\end{verbatim}
\normalsize
ou de manière équivalente :
\footnotesize
\begin{verbatim}
var = HeavySide(f(x)-g(x),Delta,dum)*g(x) + HeavySide(g(x)-f(x),-Delta,dum)*f(x);
\end{verbatim}
\normalsize
%
{\bf avertissement} La valeur de Delta en argument est importante et doit être soigneusement choise,
car de cette valeur dépend la valeur de la dérivée partielle qui interviend dans les Jacobiennes.

{\bf ATTENTION !} Il est fortement déconseillé de donner en argument des valeurs numériques.
En effet, avec passage en double précision\index{double précision}, les chiffres restent
en simple\footnote{sauf à préciser qu'ils sont en double, comme 1.D0, mais alors ils restent
en double si on repasse en simple.} alors qu'ils seront vus en double dans la
fonction~: erreur garantie.

\medskip
{\bf Autre fonction de transition}

La Heavyside est par construction symétrique autour de zéro. Il est de nombreux cas où on veut
une fonction permettant une transition strictement entre zéro et {\bf Delta}, comme
par exemple deux fonctions prenant chacune en compte le cas de l'argument positife
ou négatif. C'est le rôle de la fonction {\bf oStepin}, avec la définition~:

 {\bf ostepin(x,Delta)}$\displaystyle=\frac{x}{\Delta}-\sin(\frac{2\pi x}{\Delta})/2\pi$,

strictement nulle pour les valeurs négatives de $x$, égale à l'unité dès que $x/\Delta \ge 1)$, et suivant
la fonction $sin$ entre 0 et 1. Pente max de 2. pour $\frac{x}{\Delta}=0.5, \Delta=1.$\footnote{prendre garde
au fait que la fonction ne dépend que de $\frac{x}{\Delta}$, mais la dérivée change avec $\Delta$.}.
%% FIG : oStepin Function
\begin{figure}[h]
\psfig{figure=stepfun.eps,%
bbllx=30pt,bblly=175pt,bburx=540pt,bbury=690pt,height=9cm,width=11cm,clip=}
\caption{Une fonction step et sa dérivée pour une variable adimentionnée $x/\Delta$ (dérivée pour $\Delta=1$).}
\label{ostepin}
\end{figure}
%

La règle de dérivation utilise la fonction associée {\bf oDeltin}~:
%
\begin{verbatim}
&'(OSTEPIN(#,#))(/#)' = '((#1)(/#3)*ODELTIN(#1,#2))'
\end{verbatim}
%
avec la fonction dérivée $=(1-cos\frac{2\pi x}{\Delta})/{\Delta}$.

Ainsi en zéro la fonction et sa dérivée est nulle.
%
%---------------------------------
\subsection{Règles de dérivation}
\label{derivation}
%---------------------------------
Se souvenir que les règles de la dérivation formelle correspondent à un tableau de macros {\bf Mortran}
dans la séquence {\bf P=INCLUDE, D=\$Derive\_mac} de \Minik. En cas de doute, il faut vérifier
que les fonctions mathématiques utilisées dans les modèles y sont traitées. Il est par ailleurs
très simple d'en ajouter en analogie avec celles qui y sont\index{Dérivées partielles}.

Il est par exemple est assez courant d'avoir à utiliser la fonction intrinsèque Fortran {\bf abs()}.
Il est bon de vérifier que la règle de dérivation est\index{fonction abs}, soit dans la séquence,
soit l'ajouter dans votre ZINIT :
%
\begin{verbatim}
 &'(ABS(#))(/#)' = '((#1)(/#2)*SIGN(1.,#1))'
!                               "derivation de la fonction valeur absolue"
\end{verbatim}
%
De manière générale, il est souvent plus opportun de concevoir une macro supplémentaire que
de jouer au Fortran (du genre sqrt($f()^2$) pour avoir une abs). Dans ce cas, il est impératif
d'utiliser les noms génériques de fonctions implicites, ceci pour pouvoir passer en double
précision\index{double précision} sans souci (ex : {\bf dsin()} est interdit, c'est {\bf sin()}
qu'il faut utiliser).

On peut aussi entrer ses valeurs numériques systématiquement en double, ce  qui n'aura aucune
conséquence en simple. Exemple :
\begin{verbatim}
pi = acos(-1.D0);
\end{verbatim}
%
La variable déclarée {\bf real pi} passera de 7 à 16 chiffres significatifs automatiquement
en double\footnote{une macro transforme en effet {\bf real} en {\bf double precision} quand
on sélectionne le flag {\bf double} dans \Minik.}. Celà signifie aussi qu'il est conseillé
de déclarer {\bf real} -- et non pas {\bf dimension} ces variables là ou tableaux
(souvent des paramètres du modèle); ça se fait en tête du {\bf ZINIT}. {\it On utilisera utilement
l'option selseq : {\bf sel none} pour s'assurer que toutes les variables ont bien été déclarées
(cette option provoque la compilation de {\bf principal} avec la règle {\bf implicit none})}.

{\bf ATTENTION} Cette règle ne concerne que les fonctions implicites Fortran. Pour les autres fonctions,
le type d'un argument sous forme de chiffre conserve sa précision de déclaration~: 1. reste sur quatre octets
en double, 1.D0 reste en huit octets en simple précision. NE PAS utilisez d'arguments chiffrés en règle
générale (cf fonction {\bf Heavyside}).

{\bf IMPORTANT~: Ne jamais utiliser dans une équation de valeurs numériques sous la forme FTN ``1.2E-3'',
car la dérivation symbolique ne l'interprète pas correctement\footnote{``bug'' signalé par stb en Mai 2007}.
Ainsi, l'entrée de\index{Valeurs!numériques}
{\bf eqn: var\_eta1 = 1.2E-3*(var\_ff1-var\_eta2)} conduira à des dérivées égales à $\pm 3$ !} Un principe
général de bonne programmation est d'utiliser des variables fortran pour ranger ces valeurs. 
%
%-----------------------------------------------------------------
\subsubsection*{Dérivation de puissance fonctionnelle de fonction.}
%-----------------------------------------------------------------
La dérivation partielle par rapport à un exposant fonctionnel n'est pas fiable ($d_y g(x,y)^{f(y)}$).
On doit la remplacer par {\bf power(g,f)}, ou si l'on préfère, écrire\\ $exp(f(y)*log(g(x,y)))$.
La dérivation de la fonction {\bf power}\index{fonction power} est prise sour la forme suivante~:
%
%
\begin{equation*}
\partial_x f^g=g f^{g-1}\partial_x f +  f^g \log f\partial_x g = f^{g-1}(g\partial_x f + f\partial_x g)
\end{equation*}
%
et se trouve dans le tableau des règles de dérivation {\bf \$DERIVE\_MAC}.
%
%------------------------------------------
\subsection{Modèle de type face, interface}
%------------------------------------------
%
On entend par là un code fournissant des grandeurs indépendantes du modèle
de simulation. Par exemple, un code fournissant des grandeurs météo pour calculer
des sources aux parois extérieures, des flux solaires etc.

La meilleure manière est de coder une fonction fortran qui fasse l'interface informatique
entre le code externe (lecture de fichiers, calculs etc). Dans la mesure où
il fournit des variables indépendantes, les dérivées partielles de la détyermination
des Jacobiennes sont nulles, il n'ay a ainsi pas à fournir de règles de dérivation
de la fonction.

On aura donc un codage du genre~:
%
\begin{verbatim}
   set_phi
   <  eqn:temp_paroi = Code_temp(time, orientation);
   >;
\end{verbatim}
%
avec Code\_temp() qui va faire tout ce qui est nécessaire pour fournir
la température de l'air extérieur vue par une paroi en fonction de la date,
de son orientation, etc.

{\bf Attention} Il peut être prévu d'étudier les sensibilités à une oriention
de paroi, auquel cas la fonction devra calculer en plus sa dérivée par rapport
à cette grandeur. On écrira alors dans le tableau des macros {\bf Derive\_Mac}
une règle comme suit~:
%
\begin{verbatim}
&'(Code_temp(#,#))(/#)' = '((#1)(/#3)*Dcode_Temp(#1,#2))'
\end{verbatim}
%
avec un point d'entrée ajouté à la fonction Code\_temp pour fournir la valeur de
cette dérivée. La section suivante complète utilement ces principes.
%-----------------------------------
\subsection{Jacobiennes numériques}
%-----------------------------------
%
Certains modèles ne peuvent être analytiquement formulés et exigent une procédure
numérique\index{Dérivées partielles}.
C'est le cas des calculs radiatifs GLO basé sur la procéduree de Malkmus.
Dans de tels cas, les macros de dérivation ne peuvent opérer symboliquement
sur une expression analytique du modèle, et on leur substitue un calcul numérique., qui doit
être ajouté au modèle de détermination numérique des grandeurs -- ici des flux radiatifs.

Dans l'exemple suivant extrait de ``ClimSI''\footnote{thèse de StepH.}, il reste dans ZINIT une expression analytique
représentant le modèle numérique complet linéarisé localement par rapport à un état de
référence (les variables X\_ref) (état qui change en cours de simulation). Les coefficients
de linéraisation (flux\_xxx(.)) sont calculés par le modèle numérique.

Quand donc ? Il y a une subtilité~: l'équation du premier flux (LW\_oce) possède la
particularité par rapport aux deux autres de faire un appel à une fonction Fortran
(ici Flux\_Rad(...)) qui fait l'interface entre \Minik et le modèle numérique~:
%
\begin{verbatim}
set_phi
<
 * * * * *
! ===========================================================
! Long wave radiation flux between Ocean-Atmosphere and Space
!   (3 variables : Phi_Oce and Phi_Atm and Phi_Str)      TRAD
! ===========================================================
!*Equ:       Phi_Oce = -sigma(SST^4-T^4)
!*Equ:       Phi_Atm = sigma(SST^4-T^4-(T+Gammat*haut(Q))^4)
!*Equ:     and Malkmus narrow band model from Cherkaoui et al.

! -----------------------------------------------------------------
  eqn:  flux_LW_ocean = Flux_Rad(
              T_ocean_eta,T_tro_surf_phi,rLapse_tro,T_strato_eta,rLapse_str,
              HumRel_phi,xCO2,rN_B,rN_M,rN_H,IDO,time,istep,Zprint,ZnewFlux
                      ) + flux_oce_ref
                        + flux_oce(2)*(T_ocean_eta-Tocean_ref)
                        + flux_oce(3)*(T_tro_surf_phi-Tsurfair_ref)
                        + flux_oce(4)*(T_strato_eta-Tstrato_ref)
                        + flux_oce(5)*(HumRel_phi-HumRel_ref)
                        + flux_oce(6)*(xCO2-xCO2_ref)
                        + flux_oce(7)*(rN_B-rN_B_ref)
                        + flux_oce(8)*(rN_M-rN_M_ref)
                        + flux_oce(9)*(rN_H-rN_H_ref);



  eqn:  flux_LW_tropo = flux_trp_ref
                        + flux_trp(2)*(T_ocean_eta-Tocean_ref)
                        + flux_trp(3)*(T_tro_surf_phi-Tsurfair_ref)
                        + flux_trp(4)*(T_strato_eta-Tstrato_ref)
                        + flux_trp(5)*(HumRel_phi-HumRel_ref)
                        + flux_trp(6)*(xCO2-xCO2_ref)
                        + flux_trp(7)*(rN_B-rN_B_ref)
                        + flux_trp(8)*(rN_M-rN_M_ref)
                        + flux_trp(9)*(rN_H-rN_H_ref);

  eqn:  flux_LW_strato = flux_str_ref
                        + flux_str(2)*(T_ocean_eta-Tocean_ref)
                        + flux_str(3)*(T_tro_surf_phi-Tsurfair_ref)
                        + flux_str(4)*(T_strato_eta-Tstrato_ref)
                        + flux_str(5)*(HumRel_phi-HumRel_ref)
                        + flux_str(6)*(xCO2-xCO2_ref)
                        + flux_str(7)*(rN_B-rN_B_ref)
                        + flux_str(8)*(rN_M-rN_M_ref)
                        + flux_str(9)*(rN_H-rN_H_ref);
! -----------------------------------------------------------------

>;"end set_Phi"
\end{verbatim}
%
Les trois expressions ci-dessus seront toujours appelées dans le même ordre, et ainsi,
l'appel à la fonction Flux\_Rad aura pour conséquence un calcul complet des flux
par le modèle numérique, la détermination de l'état courant de référence, et les
coefficients du modèle linéarisé. Ces coefficients seront en général déterminés
par accroissements finis. Au modèle numérique encore d'optimiser les calculs, en
rappelant le modèle complet seulement lorsque l'état de référence aura significativement
changé. Dans l'exemple présenté, pour diminuer le nombre d'arguments, les variables
de référence ainsi que les coefficients sont mis en {\bf common} au début de ZINIT~:
%
\begin{verbatim}
! ===== Interface to Flux_rad subroutine ===============================
 Logical ZnewFlux;
 real FLUX(4,4,9), Flux_oce(9), Flux_trp(9), Flux_str(9);
! Referenced Values for linearisation
 real Tocean_ref, Tsurfair_ref, Tstrato_ref, HumRel_ref, xCO2_ref,
                flux_oce_ref, flux_trp_ref, flux_str_ref,
            rN_B_ref, rN_M_ref, rN_H_ref, rLapse_tro_ref, rLapse_str_ref;

 common/fluxrad/Flux_oce, Flux_trp, Flux_str,
                flux_oce_ref, flux_trp_ref, flux_str_ref,
           Tocean_ref, Tsurfair_ref, Tstrato_ref, HumRel_ref, xCO2_ref,
           rN_B_ref, rN_M_ref, rN_H_ref, rLapse_tro_ref, rLapse_str_ref;
! ======================================================================
\end{verbatim}
%
%-----------------------------------
\subsection{Fonctions de grille}
%-----------------------------------
%
Il est courant d'avoir à calculer une fonction globale portant sur une variable
de grille, comme une moyenne spatiale, une variance, etc. Une manière classique
de le faire est de creer une variable en chaque maille effectuant la somme du terme
de maille sur la maille et sa voisine (de droite ou de gauche). Alors, la derniere
maille contiend la fonction sommée sur les mailles et peut être traitée pour
fournir le résultat final avec une ultime fonction de bord. Cette procédure consomme
$n+1$ variables.
;
On a encore la fonction mult\_sum, mais celle ci genere des calculs de dérivation
en Fortran dont la longueur fait facilement exploser Mortran.

Il est préférable de faire usage d'une librairie de fonctions de grille sobres.
Deux premières fonctions\index{Grid 1D} effectuent la moyenne et la variance d'une variable sur $n$ mailles à partir
d'une variable explicitée, comme ça~:
%
\begin{verbatim}
var: util_moy, fun: avrg_phi(12,util2(inode:1));
\end{verbatim}
%
qui est transformée par une macro en~:
%
\begin{verbatim}
 ...    avrg_node(12,util2(inode:1),j,ff(1),m_mult);
\end{verbatim}
%
qui renvoie à la fonction Fortran. Son point d'entrée {\bf davrg\_node} fournit la dérivée
par rapport à ff(j) (l'argument {\bf j} est bidon). Le principe est de récupérer les indices
du tableau via la fonction {\bf locf} du CERN\footnote{Oui, ça marche aussi en double précision.}.

\medskip
La variance est fournie par {\bf vari\_eta()} et sa s\oe ur {\bf vari\_phi}; elle est donnée
par $\displaystyle\frac{1}{n}\sqrt{\sum_{i=1} ^n (\zeta_i-\overline{\zeta})^2}$, où $\overline{\,}$
indique la moyenne de la variable de maille $\zeta$~:
%
\begin{verbatim}
var: util_quad, fun: vari_phi(12,densite(inode:1));
\end{verbatim}
%
Le point d'entrée de la dérivée est {\bf dvari\_node}.
%

\bigskip
{\bf Fonctions de grille pondérées}
\medskip

Les mêmes fonctions que précédemment dont l'argument est un produit de variables du TEF,
pour la moyenne pondérée, dans cet exemple issu de {\bf Urbadyn} avec une utilité en variable
de transfert pondérée par une densité (d'habitants en l'occurence) qui est une variable
d'état~:
%
\begin{verbatim}
any = avrg_poids(12,Eta:densite(inode:1),Phi:util2(inode:1));
!                   ---                  ---
any = avrg_poids(12,Phi:util2(inode:1),Eta:densite(inode:1));
any = avrg_poids(12,Phi:util2(inode:1),Phi:util2(inode:1));
any = avrg_poids(12,Eta:densite(inode:1),Eta:densite(inode:1));
;
j=56; "choisi pour correspondre à un indice de ff(.) <-> util2(inode:2)"
any = F_D(avrg_poids(12,Phi:util2(inode:1),Phi:util2(inode:1)))(/ff(j));
;
i=13;
any = F_D(avrg_poids(12,Eta:densite(inode:1),Eta:densite(inode:1)))(/eta(i));
\end{verbatim}
%
Les points d'entrée sont {\bf avrg\_two, davrg\_two1, davrg\_two2} car il faut pouvoir
dériver suivant l'une et l'autre des deux variables de grille.
\newpage
Pour la variance pondérée
 $\displaystyle\frac{1}{n}\sqrt{\sum_{i=1} ^n (\xi_i\zeta_i-\overline{\xi\zeta})^2}$~:
%
\begin{verbatim}
any = vari_poids(12,Eta:densite(inode:1),Phi:util2(inode:1));
any = vari_poids(12,Phi:util2(inode:1),Phi:util2(inode:1));
any = vari_poids(12,Eta:densite(inode:1),Eta:densite(inode:1));
;
j=56;
any = F_D(vari_poids(12,Phi:util2(inode:1),Phi:util2(inode:1)))(/ff(j));
!*** Code Genere
!      any=(DAVRG_TWO1(12,ff(k2j_12(1)),ff(k2j_12(1)),j,FF(1),FF(1),M_MUL
!     *T,M_MULT)+DAVRG_TWO2(12,ff(k2j_12(1)),ff(k2j_12(1)),j,FF(1),FF(1),
!     *M_MULT,M_MULT))
;
i=13;
any = F_D(vari_poids(12,Eta:densite(inode:1),Eta:densite(inode:1)))(/eta(i));
\end{verbatim}
%
Les points d'entrée sont {\bf vari\_two, dvari\_two1, dvari\_two2}.
(cf les routines dans {\bf mathlib} et les règles de dérivation dans {\bf derive\_mac}.

\medskip
{\bf Attention~:} Il n'y a pas de diagnostic d'erreur concernant la nature {\bf Eta, Phi}
des variables.

\bigskip
{\bf Indice de Gini}
\medskip

La fonction {\bf Gini\_index} a été introduite dans \Minik par François Gusdorf. On considère
une population $N$ partageant une utilité commune ($u$), et la courbe de O.~Lorenz associée
donnant, sur $p$ groupes de population ordonnés en $u$ croissant, la part d'utilité
totale $N\overline{u}$ que partage la population cumulée $\int_0^uN(v)dv,\quad v\le u$.

Le coefficient de Gini donne un indice d'inégalité de partage du gâteau sous forme d'une
intégrale~: la surface comprise entre la première bissectrice (égalité parfaite) et
la courbe de Lorenz, affectée d'un facteur 2 (ainsi, le Gini\index{Gini (indice de)}
varie entre 0 et 1).
Exemple tiré de l'UrbaGus~:
%
\begin{verbatim}
set_Probe
<  eqn: u_quad = vari_phi(n_node,util2(inode:1));
   eqn: gini_dx = Gini_index(n_node,ETA:densite(inode:1),PHI:util2(inode:1));
>;
\end{verbatim}
%
dans lequel la population est en premier argument et l'utilité (ou un revenu)
en deuxième.

\bigskip
{\bf Nous ne saurions éviter de poser la question de l'intérêt de cet indice
issu de l'écononmie pour caractériser un état physique~: est-ce que les
diphasiques de l'unsecte pourraient brainstormer là-dessus ?}
\newpage
%-----------------------------------
\subsection{Facilités pour gnuplot}
%-----------------------------------
%
{\bf Surfaces dans gnuplot}
\medskip

Le format des fichiers {\bf .data} ne permet pas simplement de faire des graphes 2D
du genre profil x temps. On dispose d'un petit éxécutable\index{gnuplot} qui permet
de le faire aisément. Usage~:
%
\begin{verbatim}
 splot "<1dgrid res.data 4 12 5" u 1:2:3 w l
 ou encore
 splot "<grid1d res.data 4 12 5" u 1:2:3 w l
\end{verbatim}
%
en lieu et place de la lecture du fichier ici res.data, grid1d est un éxécutable dont on montre
l'installation plus bas, qui va extraire de ce fichier les valeurs indiquées par les
arguments qui suivent. Le premier, 4, est le numéro de la première colonne où se trouvent
les valeurs de la grandeur sélectionnée; 12 est le nombre de mailles --- qui pouura englober
bien sùr les conditions limites (up et/ou dwn); 5 est le nombre d'équations par mailles.
Ces informations se retrouvent dans le Model.hlp habituel, mais il s'agit du numéro de colonne
du fichier, comme à l'habitude pour gnuplot.
Le résultat est un profil en fonction de tous les instants de la simulation, en filaire.
Avec la version 4 de gnuplot, "set pm3d" trace des surfaces en couleur.

On peut aisément récupérer le nom de la grandeur en faisant~"
%
\begin{verbatim}
gnuplot> !1dgrid res.data 4 12 5|tail -1
 qui donne 
gnuplot> title "densite(1st)"
 on trace alors par 
gnuplot> splot "<1dgrid res.data 4 12 5" u 1:2:3 title "densite(.)" w l
\end{verbatim}
%
par exemple.

Arguments supplémentaires~:
%
\begin{verbatim}
 splot "<1dgrid res.data 4 12 5 100 1000 10" u 1:2:3 w l
\end{verbatim}
%
avec 100 le premier istep demandé, 1000 le nombre à tracer, et 10 de dix en dix.
Les défaults sont en principe corrects.

Exemple avec les masselottes~:
%
\begin{verbatim}
gnuplot> set xlabel "t"
gnuplot> set ylabel "maille" 
gnuplot> set zlabel "grandeur" 
gnuplot> !1dgrid res.data 2 4 2 1 1|tail -1  
gnuplot> set title "eta_pos_wall"
!
gnuplot> splot "<1dgrid res.data 2 4 2 1 500 40" u 1:2:3
                                                      title "eta_pos(.)" w l
gnuplot> replot "<1dgrid res.data 3 4 2 1 500 40" u 1:2:($3-0.02) 
                                                     title "eta_speed(.)" w l
\end{verbatim}
%
on obtiend en filaire~:
%
\begin{figure}[h]
\begin{minipage}[t]{.5\linewidth}
\psfig{figure=exemple1dgrid.eps,%
bbllx=60pt,bblly=60pt,bburx=400pt,bbury=300pt,width=7.5cm,height=6cm,clip=}
\end{minipage}
\begin{minipage}[b]{.5\linewidth}
\psfig{figure=exemple1dgrid2.eps,%
bbllx=60pt,bblly=60pt,bburx=410pt,bbury=301pt,width=7.5cm,height=6cm,clip=}
\end{minipage}\hfill
\caption{Exemples de sortie 1dgrid en filaire et avec pm3d.}
\label{gnu1d}
\end{figure}
%
\begin{figure}[h]
\begin{minipage}[t]{.6\linewidth}
\psfig{figure=profil1dgrid.eps,%
bbllx=60pt,bblly=60pt,bburx=400pt,bbury=300pt,width=7.7cm,height=4cm,clip=}
\end{minipage}\hfill
\begin{minipage}[b]{.4\linewidth}
{\bf Profil simple}

Obtenir un seul profil peut être extrêmement intéressant pour réinitialiser le run
après modification du nombre de mailles.

 Pour sortir le 50ième pas de temps, il suffit
de faire un simple plot pour le résultat ci-contre.
%
\end{minipage}
\caption{plot "$<$grid1d tr.data 42 24 22 50 1" u 2:3  title "invest\_tot(.)" w l}
\label{gnuprof1d}
\end{figure}
%

\bigskip
{\bf Installation~:}

 Récupérer timap.f dans /pub/alain, le compliler avec, par exemple,
%
\begin{verbatim}
 f77 -fno-backslash -o timap.exe timap.f
\end{verbatim}
%
il faut recopier ce .exe dans /usr/bin (ou local/bin suivant vos configs).
Il ne reste plus qu'à créer un ou plusieurs liens pour l'appeler autrement (1dgrid, grid1d...).


\bigskip
{\bf Version avec toutes les variables}
\medskip

La version {\bf timapx} donne en sortie toutes les variables, qui seront
alors accessibles dans gnuplot à partir de la colonne 3, les deux
premières étants réservées au temps et au numéro de maille. On peut ainsi
effectuer des opérations entre variables; exemple~:
%
\begin{verbatim}
 splot "<1dgridx res.data 5 12" u 1:2:3 w l
 (identique au précédent) ou encore
 splot "<grid1dx res.data 5 12" u 1:2:($3+$4) w l
\end{verbatim}
%
effectuant la somme des deux premières variables.

Comme on le voit, il y a un argument de moins, le premier étant directement
l'incrément (la première colonne sera toujours la deuxième du fichier
à lire). Les trois arguments de gestion du temps sont inchangés.
L'installation est semblable à celle de {\bf timap}.

\bigskip
{\bf Histogrammes dans gnuplot}
\medskip

Un autre éxécutable {\bf histo.f}, qui s'installe comme les précédents, permet
de faire des histogrammes à partir des fichiers colonne {\bf .data}.

 Usage~:
%
\begin{verbatim}
plot "<hist dres.data 12 100 1. 20. 1 100 5 " u 1:2 w hist
\end{verbatim}
%
effectue un histogramme\index{Histogramme} sur la 12ème colonne du fichier {\bf dres.data} avec 100 canaux
de 1. à 20. en prenant 100 données, du premier incrément et en allant de 5 en 5.
\newpage
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Principes de la résolution}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\label{prog}
Dans la boucle temporelle, {\bf principal} résout d'abord l'équation
matricielle de transfert\index{transferts}. Le calcul de la matrice $D$ est effectué en
utilisant les dérivées partielles $\delta_{\varphi} \varphi$ symboliquement
construites par le passage de Mortran\index{Mortran} (le lecteur intéressé par la technique
peut se reporter à la première annexe au présent document).

Le système TEF\index{TEF} se construit de manière analogue par calcul de toutes les
matrices Jacobiennes\index{Jacobienne} du système :
%
\begin{eqnarray}
H = \partial_{\eta} g &;& B_b = \partial_{\varphi} g; \nonumber\\
C^{\dagger} = \partial_{\eta} f &;& D = \partial_{\varphi}f ;
\label{tefmat}
\end{eqnarray}
%
Deux matrices sont alors construites.
Une première  correspondant au Système Linéaire Tangent (SLT)\index{SLT, SLTC} :
%
\begin{equation}
\begin{bmatrix} \partial_t \!\cdot\, -H  & -B_b \\ -C^{\dagger} & 1-
D\end{bmatrix} \quad 
\begin{bmatrix} \Deta  \\ \Dfi \end{bmatrix} \quad = \quad
\begin{bmatrix} u^{\eta}\\  u^{\varphi} \end{bmatrix}
\label{sltc1} 
\end{equation}
%
avec les vecteurs $u$ de sources ou forçages perturbants.
On élimine les transferts de ce système, pour obtenir la matrice d'avance
des systèmes linéaires d'état\index{Etat} classiques ($aspha$ dans le code)\index{aspha (matrice d'avance de phase)}.
Elle est sortie dans le fichier
{\bf aspha.data}\index{data file} pour analyse post-run (cf SLT, SLTC), et ne sert pas
à la résolution (comme ce fichier pet-être trés volumineux, on active sa sortie avec l'option
{\bf sel aspha} dans la selseq.kumac).

Le système est alors discrétisé en temps selon un schéma de Crank Nicolson,
pour la résolution numérique du système à chaque pas de temps. Cela
se fait par un appel à la routine {\bf oker}, qui élimine les cellules\index{cellules}
pour former la matrice de couplage (tableau {\bf couplage(.)} dans le code)
entre tous les transferts du système, résout le système linéaire réduit
obtenu, calcule les variations $\dfi$ puis de $\deta$. L'état est alors
avancé dans ``principal'' et on peut résoudre le système implicite des
transferts; fin du pas de temps.

La matrice de couplage n'est pas sortie. Cependant, il est possible de préciser
une composante de transfert\index{transferts} sur laquelle on analyse un
gain\index{gain de rétroaction} de rétroaction,
comme expliqué plus bas (option balayage de Borel\index{Borel}).
%
\subsection*{ZSTEER}
% ..................
Après incrémentation de la variable $\eta$, on passe par la séquence {\bf ZSteer} qui
doit au minimum (fourni avec le code) contenir :\index{ZSteer}
\begin{verbatim}
! ***************************************************************
! ZSTEER : gestion en fin de pas de temps.
! ***************************************************************
  ZPRINT = mod(istep+1,modzprint).eq.0;
  Zprint = zprint .or. mod(istep+1,modzprint).eq.1;
  Zprint = zprint .or. mod(istep+1,modzprint).eq.2;

  Zout = mod(istep,modzout).eq.0;

\end{verbatim}
%
qui ne fait ici que gérer la fréquence de sortie des impressions et des fichier {\bf .data}.
Dans la mesure où on a accès aux matrices\index{matrices} jacobiennes, à la matrice d'avance $A^{spha}$ et
éventullement à leur dérivées par rapport à un paramètres, de nombreux calculs spécifiques
peuvent ainsi y être effectués. Deux exemples sont donnés pour le calcul de sensibilité
des valeurs porpres à un paramètre (cf stb et seb, voir {\bf lorhcl}, ou celui des exposants
de Lyapunov dans l'exemple {\bf Predator}).

Bizarrement, on n'a pas d'exemple de variation distinguée du pas de temps {\bf dt}, ce qui devrait
interpeller l'unsecte.
%
\subsection*{ZSTEP}
% ..................
Une nouvelle séquence\footnote{102}. (vide) située avant les incrémentations temporelles permet sur critère numérique
de redémarrer un pas de temps\index{ZStep}, par l'instruction {\bf goto :RedoStep}\footnote{ceci explique pourquoi
le balayage de Borel et le calcul de {\bf aspha} ont été déplacés, le label :redoStep: est
situé juste après.}.
%
%-----------------------------------------------------------------
\subsubsection*{Le calcul en double précision.}
%-----------------------------------------------------------------
\Minik passe automatiquement en double précision en levant le drapeau {\bf sel double}.
Ceci se fait de deux façons :
\begin{itemize}
\item  à la place de {\bf mod}, on lance {\bf smod double}\index{double précision} dans
cmz\index{CMZ (CERN)}; on obtient le .exe
sous le répertoire : exo/double, où on effectue les opérations courantes.
\item on lève le drapeau double dans la selseq.kumac\index{selseq (macro CMZ)}, et on remplace
les calculs en simple par les doubles sous exo/exo classiquement après {\bf mod}.
\end{itemize}
 Il est recommandé de passer en double
dés qu'on attaque un problème d'optimisation, comme décrit dans la section~\ref{minuit}, mais aussi
plus généralement pour
vérifier la précision numérique obtenue, conjointement avec une diminution de {\bf dt}.

 Le code double est obtenu
en utilisant l'option {\bf implicit double precision(a-h,o-y)}, et par une macro Mortran transformant
{\bf real} en {\bf double precision}.
%

{\bf ATTENTION} les valeurs numériques\index{Valeurs!numériques} ne changent pas de codage ---
 ils restent dans la précision de la manière
dont ils sont écrits (1. rest toujours en simple, 1.D0 toujours en double. Il ne faut donc pas
utiliser de chiffres en argument de fonction, sauf pour les fonctions Fortran implicites ($sin, cos$)
qui changnet automatiquement tous les types.
\newpage
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Outils d'analyse dynamique des systèmes}\label{analyse}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\vspace{0.7cm}

\subsection{Balayage de Borel}
%------------------------------
%
On montre dans l'introduction au TEF-ZOOM\index{ZOOM} comment un schéma de Crank-Nicolson
donne symboliquement la solution du système transformé de Borel, en identifiant
la variable dTime\index{time, time-step} à la variable de Borel $\tau$ à un facteur $2$ près. Dès lors,
il est possible d'étudier numériquement la dépendance en $\tau$ de transformée des
divers coefficients de matrices réduites du système\index{Borel}.
Ansi, en un point de la trajectoire, on peut numériquement obtenir la transformée
de Borel en effectuant un balayage sur la variable {\bf dTime}.

Il est prévu de sortir simplement
cette dépendance en $\tau$ pour la matrice réduite à une variable scalaire ff(k).
On définit ainsi un gain de rétroaction $g(\tau)$\index{gain de rétroaction}. On explique dans le document
de la page web ZOOM (documents.dir/ClimSIre3.ps.gz) comment on définit et interprète
cette fonction de gain ou de facteur de rétroaction attaché à une composante
de transfert donnée. Un exemple d'analyse de système est donné dans le même document.

Pour les systèmes linéaires - dont le SLT\index{SLT, SLTC} est statique (ou dit aussi autonome), la
fonction de $\tau$ du gain ne dépend pas du temps. Mais en général, on analyse
la fonction $g(\tau;t)$ le long d'un segment de trajectoire, comme dans l'exemple
suivant :
\footnotesize
\begin{verbatim}
! -------------------
! Gain de retroaction
! Borel
! -------------------
ZBorel=.true.;
if ZBorel                                             " activation de l'option"
<  istep_B_deb=1000;                          " debut du balayage (increment #"
   istep_B_fin=1200;                          " fin                           "
   istep_B_inc=1;                             " à chaque pas de temps         "
;
   index_ff_gain=7;                     "composante de ff(.) soumise au calcul"
   tau_B_ini=0.2;                     "valeur initiale de la variable de Borel"
   tau_B_mult=sqrt(sqrt(2.));        "facteur multiplicatif de tau du balayage"
   itau_max=100;                                     "nombre de valeurs de tau"
   z_pr/Borel/:tau_B_mult,tau_B_ini*(tau_B_mult)**itau_max;
>;
\end{verbatim}
\normalsize
Dans cet exemple, le balayage\index{ZBorel} sera effectué 100 fois à chaque pas de temps,
entre la 1000ième et la 1200ième
itération sur la trajectoire. Le fichier {\bf  tau\_Borel.data}\index{data file} contient ces valeurs
de $\tau$, et les gains sont dans {\bf gains.data}, une ligne par point de la
trajectoire.

Lorsque l'on entre une valeur initiale de $\tau$ négative (tau\_B\_ini=-0.2;),
le balayage sera effectué sur 201 valeurs, de $-0.2$,tau\_B\_mult$*(-0.2),...,0.,$
et les valeurs positives symétriques.

Pour obtenir des surfaces $g(t,\tau)$, on peut utiliser le programme Fortran fourni
{\bf gains.f}\index{gain de rétroaction} et sa shell de lancement {\bf gains.xqt}, qui met les sorties au format
d'un histogramme de PAW, dans lequel on pourra utiliser la kumac\index{kumac (macro CMZ)} : {\bf borels.kumac}
générée par le code. Un exemple de résultat est dans le document
http://lmd.jussieu.fr/Zoom/doc/LorenzGains.ps.gz~.

{\bf autre écriture} Dans le cas d'utilisation de symboles, on peut remplacer la donnée par
 {\bf index\_ff\_gain=...} de la composante du transfert\index{transferts} par l'instruction :
\begin{verbatim}
   ZBorel for: symbole; 
\end{verbatim}
le symbole est le nom donné à la variable, avec en Grid1D la convention (inode $=$ .).
%
\subsubsection*{Autre calcul du gain dynamique de rétroaction}
% -----------------------------------------------------------
Avec la version 1.01 apparait un autre calcul automatique des pôles et résidus\index{résidu} suivant
la méthode décrite dans le texte Effluents (réservé à l'unsecte). On trouve ces résultats
dans les fichiers {\bf vpgains} qui donne $g(\tau)$ et {\bf vprepons.data}~($\frac{g}{1-g}(\tau))$,
voir la subroutine {\bf Boreleig}
\index{Borel} (version 101e). Le post-traitement est dans le patch {\bf gsrun}, qu'il faut
importer dans son exercice ({\bf smod gsrun}), et un exemple d'utilisation dans {\bf predator};
fichiers de sorties $g(t)$ et répons efficace $(t)$e dans {\bf vpgaint et vprept.data}\index{gsrun gains SLT}.
Exemple encore dans http://lmd.jussieu.fr/Zoom/doc/LorenzGains.ps.gz~et l'article refusé par la
{\it Phys.~Rev.~Lett.}
%
%\newpage
\subsection{Calcul automatique des sensibilités}
%------------------------------------------------
Une deuxième utilisation directe des matrices Jacobiennes\index{Jacobienne} du système
permet le calcul automatique de trois types de sensibilités\index{sensibilités!variables} :
\begin{itemize}
\item la sensibilité des variables d'état\index{Etat} à la perturbation
      des conditions initiales de l'une d'entre elles.
\item la sensibilité des mêmes à un pulse ou un step initial sur un transfert.
\item la sensibilité des mêmes à la modification d'un paramètre\index{paramètre} du modèle.
\end{itemize}
Ces calculs sont démarrés soit par défaut à $t=0$, soit à un autre instant.
Nous décrivons d'abord les deux premiers, pour lesquels on déclare chaque variable,
et le type de perturbation~:\index{Sensibilité} :
%
\begin{verbatim}
  Sensy_to_var
  <
!*     var: eta_courant_L, pert: init at 100;
!*     var: ff_T_czcx,     pert: pulse at 100 every 20;
!*     var: ff_Psi_Tczcx,  pert: step_eff;
!*     var: ff_Psi_Tczcx,  pert: step_Resp at 10 every 100;
! *** premiers tests identiques a lorhcl.ref
    var: ff_courant_L , pert: step_eff;
    var: ff_T_czcx    , pert: step_eff;
    var: ff_Psi_Tczcx , pert: step_eff;
    var: ff_Psi_Tsz   , pert: pulse at 100 every 50;
  >;
>;
\end{verbatim}
\normalsize
Les instructions se comprennent d'elles-même; la différence entre step {\bf eff} et {\bf resp}
renvoie aux analogues (non stationnaires) entre respectivement $\displaystyle \frac{g}{1-g}$
et  $\displaystyle \frac{1}{1-g}$. Une vérification est effectuée pour distinguer selon
l'applicabilité des modes de perturbation aux états (seul INIT est permis) et aux transferts
(tous les autres).

La mise-à-jour {\bf 102q} propose de plus une amplitude de step modulable~: pour une composante k
des transferts à laquelle on applique ce step, il suffit de rentrer
%
\begin{verbatim}
   step_psi(k)= 3.3d-04; "2xCO2 amplitude"
\end{verbatim}
%
\paragraph{Initialisation} Par défaut et sans précision de {\bf at III}, les
sensibilités de type état {\bf INIT} sont initialisees à la valeur unité ---
correspondant à la matrice unité pour un propagateur $\Phi(0,0$ complet.
On peut modifier ces valeurs en remplissant soit-même le tableau {\bf Phi\_t}
après  la fermeture du bloc {\bf sensy\_to\_var}. Exemple, la cinquième
variable d'état est une concentration de CO2 ($\approx 3.3E-04$~: une perturbation
de $1$ mène à des sensibs énormes; on met alors {\bf Phi\_t(5,5)$=3.3E-04$},
et on directement les sensibs à un doublement de CO2.

Dans le cas où on lance ce calcul à un pas de temps non initial ({\bf at 100}),
il faut effectuer la mêmeinitialisation de {\bf Phi\_t(.,.)}
dans {\bf \$ZSTEP}, mais aussi anuller la (les) colonnes correspondantes
de {\bf dPhi\_t}. On aura une séquence comme ça~:
%
\begin{verbatim}
if istep.eq.100
<     call veczero(dPhi_t(1,5),np); Phi_t(5,5)=3.3E-04;
      print*,' ---  resetting C.I. sensy to State #:',5,' at:'
                                 ,time,istep,' on:',5,'th Phi_t colon.';
>;
\end{verbatim}
%

\paragraph{Sorties} Les fichiers de sortie comprennent deux premières lignes de commentaire (gnuplot), et sont~:
%
\begin{itemize}
\item pour les états {\bf sens.data};
\item pour les transfers {\bf sigma.data}\index{transferts};
\item pour les probes {\bf sigmo.data};\index{observables}
\end{itemize}
%
Si l'option de génération de DimEtaPhi n'est pas retenue, il faut donner une valeur au paramètre {\bf NXP}
dans la séquence. Les paramètrfes nyp et nzp restent nuls - ils correspondent à des développements
possibles de l'adjoint.

Il est à souligner que cette analyse de sensibilité utilise les dérivées formelles
du modèle établies sous le TEF, et fournit des résultats plus rigoureux que par
une méthode consistant a comparer des trajectoires correspondant à des accroissements
finis des\index{paramètre} paramètres\footnote{Pour une explication succinte du calcul automatique des
sensibilités, voir le document http:lmd.jussieu.fr/zoom/documents.dir/sensibilite.ps.\index{ZOOM}
Le document plus complet Propagateurs est réservé à l'unsecte}.

%
\subsubsection*{sensibilité à une liste de  paramètres.}
Supposons que nous voulions étudier\index{paramètre} la sensibilité\index{sensibilités!paramètres} des variables à un changement initial du paramètre
{\bf apar} du modèle de Lotka-Volterra\index{Lotka-Volterra}. Il suffit de rentrer dans
son Zinit\index{fwd}\index{Free\_parameters} :
%
\begin{verbatim}
Free_parameter: apar;
\end{verbatim}
%
pour obtenir\index{Sensibilité} dans les fichiers {\bf sensp, sigmap} et {\bf sigmop.data} la sensibilité
des variables au cours de la simulation.
Il faut aussi donner le {\bf nombre de ces paramètres} dans \$dimetaphi : lp=1 ici, lorsque l'on
n'utilise pas la génération automatique de {\bf DimEtaPhi}.

Forme plus complète~:
\begin{verbatim}
Free_parameter: apar, cpar at 100; "[every 10 000]"
\end{verbatim}
%
Il peut en effet s'avérer qu'un modèle subisse une forte transition initial dont on souhaite
qu'elle n'intervienne pas dans la senbilité aux paramètres. Dans le cas présenté, on démarre
au 100ème pas de temps. Il peut de plus s'avérer qu'une sensibilité paramétrique ``explose'',
d'où la possibilité de la redémarrer (ici tous les 10 000 pas de temps, en commentaire).

{\it remarque étonnante~: supposez que dans vos équations, la valeur numérique 3.14159
apparaisse -- et si répétée, avec toujours exactement la même liste de caractères; eh bien,
on peut demander~: {\bf Free\_parameter: apar, 3.14159;} ! On aura en sortie les sensibilités
à une modification de CE patamètre. Qui trouvera mieux que Mortran ? Remarquons une
fois de plus que ça ne marche pas avec le type Fortran 0.31E01 !}

\newpage
%
\subsubsection*{sensibilités diverses à un paramètre.}
% ....................................................
Un drapeau {\bf dPi\_aspha} dans la selseq.kumac active le code calculant la sensibilité de la
matrice d'avance de phase $M$ ({\bf aspha}) par rapport au paramètre désigné comme
{\bf [mx:apar],cpar} dans la liste des {\bf free\_parameters}. La sortie est dans le\index{dPiAspha}
tableau Fortran {\bf d\_pi\_aspha(n,n)}. Le calcul effectué est le suivant :
%
\begin{eqnarray}
  d_{\varpi} M =&& d_{\varpi}H + d_{\varphi}H \, d_{\varpi} \varphi \\
                 &+& \left\{ d_{\varphi}B \, d_{\varpi} \varphi + d_{\varpi}B \right\} (I-D)^{-1} C^{\dagger}\nonumber\\
                 &+& B \left\{ (I-D)^{-1} d_{\varpi}C^{\dagger} - (I-D)^{-1} d_{\varpi}D\, (I-D)^{-1} C^{\dagger} \right\}\nonumber
\label{dpiaspha}\end{eqnarray}
On montre comme exemple d'utilisation\footnote{Une idée de stb, avec, ici, la version où
seul un transfert est explcitement dépendant de $\varpi$.} que la sensibilité de valeurs propres à ce paramètre est donnée par :
\begin{equation}
\scal{f_i}{e_i} d_{\varpi}\lambda_i = \bra{f_i}d_{\varpi} M \ket{e_i}
\end{equation}
avec $f_i$ et $e_i$ les couples de vecteurs propres à gauche et à droite de $M$.
On a une formule équivalente pour les valeurs singulières\index{Valeurs!propres}\index{Valeurs!singulières}.
Il faut déclarer les tableaux nécessaires dans {\bf ZINIT} et les
calculs sont effectués dans le {\bf \$ZSTEER}\index{ZSteer} (cf exo {\bf lorhcl}).
%
\newpage
\subsection{Utilisation du modèle adjoint}
%------------------------------------------
%
De façon générale, le modèle adjoint\index{adjoint} permet de calculer les gradients
d'une fonctionnelle de coût $J$ à des perturbations, le long du
SLTC\footnote{Système Linéaire Tangent Circulant le long d'une trajectoire}\index{SLT, SLTC} :
%
\begin{equation}
J = \psi[\eta(T),\varphi(T) ,h(T)] + \int_0 ^T {l[\eta(\tau),\varphi(T),h(\tau)]}\, d\tau
\end{equation}
%
Exemple, où on active l'option de calcul de l'adjoint {\bf ZBack}\index{ZBack}
\index{masselottes}\index{node} :
\begin{verbatim}
! ---------------------------------------------------
!   Adjoint  (correspond a l'exemple des masselottes)
! ---------------------------------------------------
Zback=.true.;
! Initialisation v(T) -> pour zinit-adjoint
  Fun_set cout_Psi = eta_move(inode:1);
! et fun_set cout_l integrande de la fontion de cout
  Fun_set cout_l = 0.;
\end{verbatim}
Il suffit de donner les deux fonctions correspondant au calcul de
la fonctionnelle $J$ : cout\_Psi (pour $\psi(\eta(T)$, le coût dit final), et
cout\_l correspondant à l'intégrande. Dans l'exemple particulier, la fonctionnelle est
réduite à la valeur finale d'une composante de l'état\index{Etat}. A l'instant $t=0$, le
vecteur adjoint à l'état sera égal à la sensibilité de l'état final
de cette deuxième composante (position de la première masselotte)
à la perturbation de l'état initial des positions et vitesses de tout
ces objets. De façon plus générale, on obtient le gradient de la fonction
de coût\footnote{Les explications détaillées du fonctionnement de
l'adjoint sont données dans le document : http:www.lmd.jussieu.fr/documents/Adjoint.pdf
ou .ps.gz}\index{adjoint}.
Lors du calcul, les variables v\_adj(.) et w\_adj(.) sont respectivement
les adjointes à eta(.) et ff(.). On les récupère dans les fichiers {\bf vadj.data}\index{data file} et
{\bf wadj.data}.

%
\subsubsection*{sensibilité à des paramètres.}
On obtient très simplement en plus la sensibilité de la fonction de coût à des paramètres\index{paramètre}
du modèle dans le fichier : {\bf gradpj.data} en donnant la liste des paramètres dîts alors
``libres'' dans le Zinit\index{Free\_parameters} :
%
\begin{verbatim}
    Free_parameters: apar, cpar ;    " pour predator"
    Free_parameters: rd(1), rk(2) ;   " pour masselottes"
\end{verbatim}
que l'on déclare après l'initialisation numérique des variables FTN définies comme paramètres
libres. (L'option [fwd: . ]\index{fwd} peut porter en plus sur un terme de la liste.
{\bf Il faut aussi donner le nombre de ces paramètres dans \$dimetaphi : lp=2 ici.}
%
\subsubsection*{ système avec loi de commande}

Si le problème comprend une loi de commande\index{loi de commande}, on doit activer
le drapeau {\bf ZCommand} dans ZINIT. Les deux types de commandes sont données par les
variables {\bf ux\_com(1: np)} et {\bf uy\_com(1: mp)} commandant les états\index{Etat}
et les transferts\index{transferts}, composante par composante.

Deux types de commandes sont
possibles : a) soit elles sont données
comme fonction des variables du système, et le drapeau {\bf ZLaw} doit être levé
(cf l'exemple {\bf pousse})\index{pousse (exo)},
 b) soit elles pilotent à chaque pas de temps le système par
valeurs numériques dans les fichiers {\bf uxcom} et {\bf uycom} (cf Adjoint et \Minik).

La commande par fonction résultera en l'écriture de ces deux derniers fichiers.
Pour entrer une loi de commande, lever le drapeau {\bf Zlaw} flag, et entrer
la loi dans la séquence \$ZCMD\_LAW du patch ZinProc (voir l'exemple
{\bf pousse}). Si l'option d'entrée par fichier est choisie, l'utilisateur est
censé effectuer lui-même la recherche d'optimisation\index{optimisation}
entre deux simulations. On a provisonnellement prévu un embryon de programme
dans le PATCH {\bf ucom}.
On pourrait par exemple choisir de piloter le paramètre {\bf apar} dans le
modèle de Lotka-Volterra\index{Lotka-Volterra}.
Dans ce cas, on remplace l'occurence du paramètre dans l'équation d'évolution des proies
par la variable réservée {\bf ux\_com(1)} et, soit définir une loi dans la séquence
\$ ZCMD\_LAW , soit écrire les valeurs instantannées du pilotage dans le
fichier uxcom.data\index{data file} avant calcul dans \Minik.

\emph{Remarque} Il est supposé que chaque composante de commande n'interviend que dans l'équation
de même indice. Si cette simplification s'avérait trop restrictive dans l'avenir, le
code devra être légèrement modifié. Il est d'ores et déjà possible d'utiliser des composantes
de transfert pour le faire (pendant qu'il est chaud).

\bigskip
A l'heure actuelle, n'ont été effectuées que des optimisations avec Minuit. Horws, depuis
l'avènement de la {\bf 102}, il devient possible d'utiliser la matrice complète de sensibilité
à la liste des paramètres devant être dèterminés, ce qui donne accès à la méthode de Newton (cf Cahiers),
en utilisant par exmple la subroutine du PSI~:{\bf OLLSQ}. Nous attendons un exercice
pratique pour décrire en détail cette méthode dand \Minik.
\newpage
%
\section{Optimisation avec l'adjoint et Minuit}\label{minuit}
%--------------------------------------------------
\vspace{0.7cm}

\emph{La description suivante reste succinte et le lecteur est
renvoyé au manuel Minuit du CERN pour une introduction plus complète\index{optimisation}. Il est
utile aussi de regarder les subroutines citées dans {\bf $//$pousse$/$minuik} et
{\bf lorhcl/minuik}, ce dernier effectuant un ajustement des conditions initiales
sur une trajectoire de référence d'un Lorenz 63}\index{adjoint}\index{Lorenz}.
\vspace{0.5cm}

Le CERN a développé un programme très élaboré permettant de maximiser une
fonctionnelle calculée dans la fonction générique dénommée FUN : {\bf Minuit}.
Minuit n'est pas une boite noire qui entreprendrait de faire un fit.
Il s'agit plutôt d'une boite à outil permettant sobrement d'enchainer des
algorithmes sans se préoccuper de les programmer, de mettre en forme des
sorties, et tout le reste.

On construit son menu d'appels à Minuit à partir d'un programme principal
du patch {\bf minuik (DSDQ)}. On a prévu une interface avec \Minik qui est
la subroutine {\bf FMINI} dont le rôle est de passer les arguments de l'un à
l'autre dans les deux sens (voir ces programmes dans ``pousse''\index{pousse (exo)}.
Il est dirigé par Minuit via le drapeau IFLAG\index{IFLAG}\index{paramètre}.
%
\footnotesize
\begin{verbatim}
IFLAG = 
 0 : premier run pour inititialiser la commande\index{loi de commande} (ou les paramètres).
 1 : calcul de la valeur FUN à minimiser (-> calcul direct seul).
 2 : demande de Fun et de ses gradients       (-> adjoint aussi).
 3 : calcul final avec meilleur fit. Mini_ker en profite pour
     donner un listing complet, calcule les .kumacs\index{kumac (macro CMZ)} etc.
 4 : Fun seult (jamais de calcul de gradients : pour simplex ou scan).
\end{verbatim}
\normalsize
Minuit dirige les appels à \Minik en suivant une stratégie, lourde par
défaut quand on minimise par MIGRAD avec gradients conjugués, calcul des
matrices de covariance\index{covariance} des paramètres etc. On a de nombreux enchaînements
de IFLAG 4 et 2 suivant l'algorithme : typiquement, on aura facilement
npar x 500 appels à la fonction (npar=nbre de paramètres) !
Mais en le pilotant correctement, 50 x npar peuvent suffire.

La fonction FUN est ici le ``Principal'' de \Minik, qui se transforme en subroutine
quand on exécute sous cmz\index{CMZ (CERN)} la commande 

{\bf $>$ smod minuik}

ou

{\bf $>$ mod minuik}

mais avec le sel flag (monitor levé dans la selseq.kumac\index{selseq (macro CMZ)}, cf option Monitor).

L'exemple donné dans "pousse" consiste en un problème de commande optimale,
la commande consistant en la force de poussée à déterminer à chaque instant
pour minimiser une fonction de coût (cf manuel Adjoint et ...).
On sélectionne arbitrairement 11 paramètres seulement parmi les 1001
instants de commande dans le tableau uxcom(2).
On demande à Minuit de minimiser - cout\_J\index{cout\_J} (cout\_J étant à maximiser),
et on interpole la commande entre ces valeurs.

{\bf FMINI} est standard à ceci près qu'il faut extraire du vecteur commande
sa seule composante active (ou plusieurs en général).
Par contre, {\bf DSDQ} est très dépendant du problème puisqu'il pilote la
maximisation. {\bf Attention} : il est impératif avant un fit\index{fit} de définir correctement
certains paramètres de Minuit. En particulier, les paramètres {\bf TOL} et {\bf UP}
sont essentiels. En gros, UP devrait etre un pourcentage du min de J attendu, et Minuit
considère que l'erreur sur les paramètres correspond au passage de $J=J_{min}$ à
$J=J_{min}+$ UP. Donc, typiquement, si on travaille en écarts quadratiques\index{quadratique} avec
des incertitudes bien choisies, on aura UP=1., c'est le défaut.
TOL (tolérance) est le paramètre critère de convergence : on a convergence dès que
la distance estimée au minimum (EDM, Estimated Distance to Minimum) est inférieure
à $0.001*$TOL$*$UP dans MIGRAD (à moins que l'on ai dépassé MAXCALL). Avec SIMPLEX,
on prend simplement comme critère TOL - ou $0.1*$UP par défaut.

La subroutine {\bf TINTER} se trouve dans la \Minik file (mathlib). Elle
fait l'interpolation entre les (par ex.) 11 commandes\index{loi de commande} fixées par Minuit, et de plus,
au point d'entrée {\bf TINTERG}, elle
propage les gradients de toute les commandes sur les paramètres\index{paramètre} sélectionnés
pour Minuit (ici, les 11). Ce problème d'interpolation entre les instants significatifs
et les instants numériques (fixés par un critère numérique sur le pas de temps)
est classique et devra toujours être considéré, faute de sous-estimer fortement la
sensibilité\index{Sensibilité} de la fonction de coût à la commande.

\vspace{0.5cm}
%\newpage
{\bf Stratégie de DSDQ pour cet exemple :}
\vspace{0.3cm}
Le dernier paramètre joue (comme c'est assez généralement le cas) un rôle
particulier : ici, il a peu d'effet direct sur cout\_J\index{cout\_J}. On le détermine
par un {\bf SCAN(n)}. Minuit calcule n fois cout\_J et prend la meilleure valeur
de ce paramètre, et produit une très jolie courbe ascii de convergence :
%\newpage
\scriptsize
\begin{verbatim}
 MINUIT WARNING IN SCAn    
 ============== VARIABLE  11 IS AT ITS LOWER ALLOWED LIMIT
 SCAN OF PARAMETER NO. 11,  Ucom11    
                  ...........................................................
                  .                          .                              .
  -75000.00     ...                          .                              .
                  .                          .                              .
  -85000.00     ...                          .                              .
                  .                          .                              .
  -95000.00     ...                          .                              .
                  .                          .                              .
  -105000.0     ...                          .                             *.
                  .                          .                      ******* .
  -115000.0     ...                          .                *******       .
                  .                          .        ********              .
  -125000.0     ...                          *********                      .
                  .             *************&                              .
  -135000.0     ... ************             .                              .
  -140000.0     .............................................................
                  /         /         /         /         /         /        
            -2.000     18.00     38.00     58.00     78.00     98.00    
                         ONE COLUMN= 2.000000       Overprint character is &   

 FCN=  -132574.8     FROM SCAn      STATUS=IMPROVED     100 CALLS      101 TOTAL
                     EDM= unknown      STRATEGY= 1      NO ERROR MATRIX       

  EXT PARAMETER               CURRENT GUESS      PHYSICAL LIMITS       
  NO.   NAME        VALUE          ERROR       NEGATIVE      POSITIVE  
   1    UCom1       50.000        1.0000        0.0000        110.00    
   2    UCom2       50.000        1.0000        0.0000        110.00    
   3    Ucom3       50.000        1.0000        0.0000        110.00    
   4    Ucom4       50.000        1.0000        0.0000        110.00    
   5    Ucom5       50.000        1.0000        0.0000        110.00    
   6    Ucom6       50.000        1.0000        0.0000        110.00    
   7    Ucom7       50.000        1.0000        0.0000        110.00    
   8    Ucom8       50.000        1.0000        0.0000        110.00    
   9    Ucom9       50.000        1.0000        0.0000        110.00    
  10    Ucom10      50.000        1.0000        0.0000        110.00    
  11    Ucom11     0.10490E-03    1.0000        0.0000        110.00    
                               ERR DEF= 0.100E+04
 **********
 **    6 **FIX   11.00    
 **********
\end{verbatim}
\normalsize
ceci fait, on a fixé ce paramètre\index{paramètre} qui ne sera plus libre dans Minuit.
%
On essaie ensuite un simplex\index{Minuit!simplex}, puis un migrad (très coûteux).
On termine le fit par un calcul des matrices de covariance\index{covariance} et d'erreur (MINOS\index{Minuit!minos}.
%
On finit par deux appels à contour de corrélation entre
couples de paramètres pour vérifier la quadrique de convergence :
\newpage
\scriptsize
\begin{verbatim}
 **********
 **   12 **CONTOUR   3.000       6.000       2.000       25.00    
 **********
 Y-AXIS: PARAMETER   6: Ucom6     
    84.97               2222222        
    82.05            2222  *  2222     
    79.13          222     *     222   
    76.21         22       *       22  
    73.29        22        *        22 
    70.37       22         *         22
    67.45       2        111111       2
    64.53      22      111 *  11       
    61.61      2      11   *   11      
    58.69     22     11    *    11     
    55.77     2      1     *     11    
    52.85     2      1     00     1    
    49.93     2******1*****00*****1****
    47.01     2      1     00     1    
    44.09     2      1     *     11    
    41.17     22     11    *     1     
    38.25      2      11   *    11     
    35.33      22      111 *  111      
    32.41       2        111111       2
    29.49       22         *         22
    26.57        22        *        22 
    23.65         22       *       22  
    20.73          222     *     222   
    17.81            2222  *   222     
    14.89               22222222       
              I            I          I
           44.01                    110.0    
                       80.25    
       X-AXIS: PARAMETER  3: Ucom3       ONE COLUMN=   2.640    
 FUNCTION VALUES: F(I)= -0.1672E+06 +   1000.     *I**2
 **********
 **   13 **CONTOUR   2.000       7.000       2.000       25.00    
 **********
 Y-AXIS: PARAMETER   7: Ucom7     
    75.00     3           22222222     
    72.08              2222   *  2222  
    69.17            222      *     222
    66.26          222        *       2
    63.35         22          *        
    60.44        22           *        
    57.53        2         111111      
    54.61       22       111  * 111    
    51.70       2       11    *   111  
    48.79      22      11     *     1  
    45.88      2       1      *     11 
    42.97      2      11     00      1 
    40.06     *2******1******00******1*
    37.15      2      11     00      1 
    34.23      2       1      *     11 
    31.32      22      11     *     1  
    28.41       2       11    *   111  
    25.50       22       111  *  11    
    22.59        2         1111111     
    19.68        22           *        
    16.76         22          *        
    13.85          222        *       2
    10.94            222      *      22
    8.029              2222   *  22222 
    5.117     3           22222222     
              I               I       I
           50.62                    110.0    
                       89.05    
       X-AXIS: PARAMETER  2: UCom2       ONE COLUMN=   2.375    
 FUNCTION VALUES: F(I)= -0.1672E+06 +   1000.     *I**2
\end{verbatim}
\normalsize
Les quadriques sont interrompues par la borne maximale donnée aux paramètres 2 et 3\index{paramètre}.
\subsection*{STRATEGIE de MINUIT}
%================================
Sans gradients fournis, il en effectue le calcul par différences finies, soit
4 appels par paramètre ($\pm$ delta puis $\pm$ $.1\times$ Delta, pour vérifier la linéarité).
En mode gradient, il calcule la direction de minimisation, mais effectue
les appels à la fonction en modifiant successivement chaque paramètre. On peut
utilement consulter res.lis de pousse\index{pousse (exo)}.

S'il n'arrive pas à la précision demandée, il entreprend un autre essai
en fixant successivement chacun des paramètres libres et en recalculant
un minimum avec les autres (recherche de pathologies).

On peut lui en demander moins en donnant la valeur 0 à STRATEGY
(défaut=1, max=2).

{\bf Remarque sur la précision}:
On est en général déçu par la précision obtenue.
Les calculs dans Minuit sont faits en double précision, et on pousse
l'utilisateur à utiliser l'option {\bf double}\index{double précision} ou une machine de 64 bits pour \Minik.

\subsection*{TOURS DE MAIN}
%=========================
Une fois lancée l'exécution, on peut à tout moment faire
-- soit un {\bf tail} sur uxcom.data\index{data file} pour suivre l'évolution du fit
-- soit plus complètement :

{\bf $>$ Crtl  z}  pour interrompre et regarder les divers fichiers
 on remet en 

{\bf $>$ bg}   \cad en background pour relancer.

Pour interrompre de nouveau, mettre en foreground puis interrompre :

{\bf $>$ fg}

{\bf $>$ Ctrl z}

Minuik demande à Principal une sortie complète des zprint\index{modzprint, Zprint} tous les
100 appels par Minuit (trois sorties successives).

Il est très pratique de ``tracer'' la progression de Minuit en suivant les exemples donnés
(subroutine FMINI) qui génèrent un fichier {\bf Minuik\_trace.data}\index{data file} (cf exos pousse et lorhcl).
\newpage
%
\section{Le filtre de Kalman}\label{kalman}
%------------------------------------------
%
{\it version encore non raccrochée à la 102. Le vecteur des obs est à présent
donné par les {set\_probe:<>}, mais à part ça les principes sont les mêmes. On a cependant
grandement amélioré les résultats avec un nouvel algorithme d'avance de la matrice
de var-covar. Attention, dans ce paragraphe, $\mu$ désogne un bruit alors que $\omega$
est le vecteur des observables (noté $\mu$ ailleurs.}

\footnote{traduite librement du texte anglais de pat}Le filtre de Kalman\index{Kalman} permet
d'assimiler des données en cours de simulation. 
On suppose ainsi l'existence d'un modèle réel avec perturbation stochastique
de ses états\index{Etat}, et qu'on dispose d'observations bruitées de sa réalité.
La situation prise en compte
dans \Minik correspond à des états continument perturbés et à des obervations
échantillonnées. On a donc un système sous TEF\index{TEF} du type :
%
\begin{eqnarray}
\partial_t \eta (t) &=&  g(\eta(t),\varphi(t)) + W(t) \mu\\
\varphi(t) &=& f(\eta(t),\varphi(t))\\
\omega(t) &=& h ( \eta(t) , \varphi(t)) + \nu
\end{eqnarray}
%
Les obsevations $\omega$ sont disponibles aux temps discrets $t=s_i$. La perturbation
stochastique $\mu$ est caractérisée par sa matrice de variance-covariance $Q$\index{covariance},
et $\nu$ par $R$. $W$ effectue l'impact sur $\eta$ des perturbations. A chaque
instant d'observation, le filtre met à jour la meilleure estimation des états et de leur
matrice de vcov au sens où l'erreur résiduelle\index{résidu} est orthogonale à l'accroissement
d'information déduite des obs bruitées\index{GetObs}.

L'exemple suivant est linéaire et perturbé-bruité. On a trois états, trois
transferts-copie\index{transferts},
mais on a deux sources de perturbations qui agissent sur les trois états.
%
\begin{equation}
\left\{
\begin{array}{cc}
\partial_t \eta_1 =& a_{11} \eta_1 + a_{12} \varphi_2 + a_{13} \varphi_3 + W_{11} \mu_1 + W_{12} \mu_2\\
\partial_t \eta_2 =& a_{21} \varphi_1 + a_{22} \eta_2 + a_{23} \varphi_3 + W_{21} \mu_1 + W_{22} \mu_2\\
\partial_t \eta_3 =& a_{31} \varphi_1 + a_{32} \varphi_2 + a_{33} \eta_3 + W_{31} \mu_1 + W_{32} \mu_2
\end{array}
\right.
\end{equation}
\begin{equation}
\left\{
\begin{array}{cc}
\omega _1 =& \varphi _1 + \nu_1\\
\omega _2 =& \eta _2 + \nu_2 \\
\omega _3 =& \eta _3 + \nu_3
\end{array}
\right.
\end{equation}
%
\subsection*{Codage du filtre de Kalman}
%
Il faut lever le drapeau {\bf sel kalman}\index{Kalman} dans la selseq\index{selseq (macro CMZ)}, et la variable
logique Zkalman=.true. dans le zinit. Il faut aussi donner la dimention\index{dimension} de la source
de perturbation ainsi que celle du vecteur des observations dans {\bf dimetaphi} :
\begin{verbatim}
parameter (nobsp=3,nerrp=2);
\end{verbatim}
On a ici trois obs et deux sources de perturbation.
%
\subsubsection*{L'observation et les perturbations}
Dans {\bf zinit}, on doit définir un ``modèle diagnostic d'observation'',
qui est un ensemble de
fonctions de $\eta$ et $\varphi$ qui devra être ajusté aux vraies obs provenant
de l'extérieur du run (censé interroger le modèle réel). Exemple\index{Fun\_set} :
\begin{verbatim}
set_probe
< eqn: Obs_tef(1) = ff(1) ;
  eqn: Obs_tef(2) = eta(2);
  eqn: Obs_tef(3) = eta(3);
>;
\end{verbatim}
%
Le nom {\bf Obs\_tef(.)} est bine sûr arbitraire.
La vcov du bruit est constante et entrée dans le tableau {\bf covobs(mobs,mobs)} :
\begin{verbatim}
covobs(1,1) = 0.3; covobs(1,2) = 0.;  covobs(1,3) = 0.;
covobs(2,1) = 0.;  covobs(2,2) = 0.1; covobs(2,3) = 0.;
covobs(3,1) = 0.;  covobs(3,2) = 0.;  covobs(3,3) = 0.2;
\end{verbatim}
%
En cours de simulation, il est nécesaire d'indiquer au moteur de l'assimilation que
des mesures sont disponibles. Cela se fait dans {\bf zsteer}\index{ZSteer}, en levant le drapeau ZGetObs\index{GetObs}
et en remplissant le vecteur {\bf vobs(mobsp)}, qui sera ainsi pris en compte
dans le traitement du pas de temps.

Enfin, la vcov des états\index{Etat} doit être donnée dans le tableau {\bf coveta(np,np)} :
\begin{verbatim}
coveta(1,1) = 0.2;   coveta(1,2) = 0.001;
coveta(2,1) = 0.001; coveta(2,2) = 0.1;
\end{verbatim}
avec en plus la matrice d'impact des perturbations sur les états ($W$) dans
le tableau {\bf mereta(np,nerrp)} :
\begin{verbatim}
mereta(1,1) = 1.;   mereta(1,2) = 0.;
mereta(2,1) = 0.;   mereta(2,2) = 1.;
mereta(3,1) = 0.5;  mereta(3,2) = 0.5;
\end{verbatim}
%
\subsubsection*{Initialisation}
Enfin, on doit initialiser la matcov des états qui sera ensuite calculée par le filtre :
tableau {\bf covfor(np,np)} :
\begin{verbatim}
covfor(1,1) = 1000.; covfor(1,2) = 10.;   covfor(1,3) = 10.;
covfor(2,1) = 10.;   covfor(2,2) = 5000.; covfor(2,3) = 5.;
covfor(3,1) = 10.;   covfor(3,2) = 5.;    covfor(3,3) = 2000.;
\end{verbatim}
%
\subsubsection*{Résultats}
Les résultats de la simulation ajustée sont dans les fichiers standard\index{data file} ({\bf res, tr, dres.data}),
on a en plus un fichier {\bf obs.data} rassemblant les observable calculées par le modèle ayant servi
à l'ajustement.
%
%\newpage
%
\section{L'option Monitor}
%--------------------------------------------------------
%
Il est possible de piloter les simulations à partir d'un programme principal appelant
\Minik (cette option transforme donc comme pour {\bf Minuik} le ``principal'' en
subroutine). Dans cas, l'appel se fait avec les arguments impératifs : 
\begin{verbatim}
call principal(Cost,ncall,IFlag,osuffix,info,idxerror);
\end{verbatim}
L'arguments {\bf osuffix} sert à distinguer les fichiers de sorties : il remplace
le {\bf .data} standard. Un tableau des principales causes de déroutement de \Minik
est donné par les arguments de sortie info et idxerror\index{LaPack} :
%
\begin{table}[ht]
\begin{tabular}{|c||c|c|c|}
\hline Source of error / warning                                   & index      &  idxerror      \\
 & &\\
%%--------------------------------------------------------------------------------------------
\hline state matrix inversion in ker                                & LaPack    &   1    \\
% & &\\
\hline system resolution error in ker                               &  LaPack    &   2    \\
% & &\\
\hline state matrix advance in phase space,                         & principal    &   3    \\
% & &\\
\hline analysis, state matrix advance (aspha)\index{aspha (matrice d'avance de phase)}, (I-D) inversion      &  LaPack &   21    \\
% & &\\
\hline analysis, Choleski factorization, non pos-def var-covariance mx\index{covariance} &  ---    &   22    \\
% & &\\
\hline analysis, error matrix inversion                             &  LaPack    &   23    \\
% & &\\
\hline variance covariance update system resolution (Kalman)        &   ---    &   24    \\
% & &\\
\hline \hline linearity problem for transfers\index{transferts}              &   principal    &   -1   \\
% & &\\
\hline Newton D\_loop\index{D\_loop} does not converge                             &  newt    &  -2    \\
% & &\\
%%----------------------------------------------------------------------------------------------
\hline 
\end{tabular}
\caption{Monitoring errors and warning from \Minik}\index{Kalman}
\label{idxerr}
\end{table}
%

La structure choisie est la suivante :
\begin{itemize}
\item a) on définit un {\bf patch} de nom quelconque (machin),
dans lequel on code son programme pilote;
\item b) On lève le flag {\bf sel monitor} dans
la selseq.kumac\index{selseq (macro CMZ)};
\item c) On assemble et éxécute par {\bf smod machin}.
Les sorties sont alors dans le répertoire machin, sous lequel on a un cfs etc
\end{itemize}
Les divers paramètres en plus des arguments nécessaires au pilotage de \Minik seront
rangés dans un {\bf common/interface/ }. On devrait aussi convenir de noms spècifiques
pour la liste des paramètres (genre {\bf pi(.)}).

Comme exemple d'utilisation, Pat pensait à une étude des bifurcations par variation
paramétrique. Mais de façon générale, pour l'utilisation de Minuit, on a intérêt a utiliser
cette option. Donnons l'exemple de StepH qui cherche à définir les paramètres\index{paramètre} d'un modèle
macro-éco ayant propention certaine a devenir chaotique, voire singulier pour certaines valeurs des
paramètres. Par principe, Minuit explore au maximum le domaine de définition donné et fait
ainsi régulièrement tomber le modèle dans des régions singulières. StepH nous a ainsi expliqué
au Zircon ``Plateforme'' qu'il était indispensable de gérer les erreurs pour empècher \Minik de
se vautrer. Les erreurs détectées (cf tableau~\ref{idxerr}) aiguillent vers {\bf ZSteer}\index{ZSteer}, où
l'utilisateur renvoie au programme père. Dans le cas cité, il faut ``signaler'' à Minuit que la
région explorée est mauvaise, ce que StepH fait en renvoyant un coût de $10^9$ et des gradients
de $10^{-3}$.
%
%\newpage
\section{L'option ZObs} ZGetObs\index{GetObs}\index{ZObs, ZGetObs}
%----------------------------------------------------------------
Voir Pat et STB, j'ai pas tout suivi l'utilisation de cette possibilité, mais il
semble s'agir d'un point de rencontre (dans \$ZSteer) de la simulation du
modèle avec des instants où des donnèes externes (mesure, référence) sont
disponibles. Par exemple, la fonction de coût (cout\_J) n'accumulerait des
écarts quads que pour ces instants. Il faut signaler au programme par {\bf Zobs=.TRue.;}
qu'il existe des données, et synchroniser, avec la variable logique {\bf ZGetObs}, le
tempo de la mesure avec celui de la simulation.
%\newpage
%
\section{Autres outils d'analyse}
%--------------------------------------------------------
%
\subsection*{Analyse de stabilité des modes rapides : option SLT}
%---------------------------------------------------------------
\label{sltc}
Les\index{SLT, SLTC} analyses précédentes sont réalisées en même temps qu'un
parcours de trajectoire. On peut également analyser de manière
plus classique la dynamique du {\bf modèle d'état} grâce à la
matrice $A_{st}$. La macro cmz\index{CMZ (CERN)} {\bf smod SLTC} permet de préparer
un programme principal ({\bf circul} dans le +PATCH SLTC), fourni
en standard mais prévu pour être modifié au gré de l'utilisateur.
En standard, on calcule une SVD (Décomposition en Valeurs Singulières)
à la fois de $A_{st}$ et $(A_{st} + A_{st} ^\dagger)$.
Les éléments singuliers\index{SVD} de ces deux matrices correspondent à l'analyse
des modes de croissance les plus rapides du système perturbé.
Voir directement le programme {\bf $/sltc/circul$} bien commenté.

Attention, avant de lancer un run sous ./..exo../sltc/, il est
impératif de créer deux liens, sur deux fichiers de /..exo../.
Sur {\bf title.tex} pour transmettre un titre au graphique, et
{\bf aspha.data}\index{data file}\index{aspha (matrice d'avance de phase)} pour l'accès à la matrice d'avance d'état\index{Etat}.

{\bf Sorties} : {\bf u,w,v.data} et  {\bf us,ws,vs.data} pour les éléments singuliers
des deux matrices.
leurs macros .kumac\index{kumac (macro CMZ)} correspondantes pour PAW)\footnote{Explications
dans le papier sur le SLTC (Al1 2003)}.
%
\subsection*{Analyse du propagateur sur le SLTC}
%----------------------------------------------- 
\label{sltcirc}
La matrice $A_{st}$ est utlisée pour calculer le propagateur sur le SLTC.
Le calcul est effectué par produits finis du développemnt au cinquième ordre
de $\Phi(t+\delta t,t)=\exp{A_{st} \delta t}$.
Divers éléments d'analyse sont fournis en standard, dont en particulier
les exposants de Lyapunov\index{Lyapunov (exposant de)} (cf programme {\bf circule} dans le patch SLTCIRC.
On prépare une exécution en lançant la macro smod SLTCIRC sous cmz\index{CMZ (CERN)}.
Exécutable et sorties (avec leurs kumacs) dans /..exo../sltc\index{data file} :
%
\footnotesize
\begin{verbatim}
              ntuple # var name
---------------------------------------------------------------------
phit.data       55       p[ij]
uphit.data      50       up       SVD de Phit 
wphit.data      51       wp        = U[w_diag]V'
vphit.data      52       vp
wr.data         53       wr[i]    éléments propres
wi.data         54       wi       partie réelle et imaginaire
---------------------------------------------------------------------
lphit.data      65       lp[ij]   1/t Ln Phit     + 1/LnDetPhit,<TrA>
ulphit.data     60       ulp      SVD de 1/t LnPhit
wlphit.data     61       wlp
vlphit.data     62       vlp
wrl.data        63       wrl[i]   VP de 1/t LnPhit
wil.data        64       wil
---------------------------------------------------------------------
lwphit.data     67       lwp[i]   1/t Ln w : Phit=UwV' (Lyapunov)
lwr.data        68       lwr      1/t VP de UV'
lwi.data        69       lwi
---------------------------------------------------------------------
\end{verbatim}
\normalsize
Le bloc 60-65 contient les éléments d'analyse du propagateur assimilé à
un système autonome $\Phi(t,0)=\exp{At}$ où $A$ serait la matrice d'état\index{Etat}
autonome d'un système donnant le même propagateur au temps $t$.
Le dernier bloc 67-69 donne les exposants de Lyapunov ordonnés du système
et les éléments propres associés.

Trois liens doivent être créés avant de lancer une exécution, les
deux du SLTC\index{SLT, SLTC}, avec un troisième sur {\bf dres.data}\index{data file}, ceci pour avoir
accès à l'incément de temps du calcul le long de la trajectoire\footnote{Pour
la recherche sur l'analyse du propagateur, consulter les documents
sur le SLTC et les Gains\index{gain de rétroaction} sur champs (Al1 2003-2004)}.
%
\subsection*{Analyses en ligne}
%-----------------------------------------------
Comme on l'a déjà remarqué, l'utilisateur a toute laltitude pour extraire les quelques
lignes des nombreux calculs de SLT et SLTC qui l'intéresse et les programmer dans la
séquence {\bf ZSteer}\index{ZSteer} en dimentionnant les tableaux nécessaires dans le {\bf ZINIT}.
Toujours penser aux règles pour passer sans problème en double précision (en particulier le nom
générique des fonctions FTN sin, cos, min, etc).
%
\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Annexe sur la programmation de \Minik}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
On donne ici des indications sur la manière dont les macros\index{macros (MTN)}, et en particulier
celles de la dérivation partielle, arrivent à faire le travail, sans
effets magiques.

Dans cette boucle extraite de {\bf principal}, on calcule la matrice
jacobienne\index{Jacobienne} $Bb(,j)=\partial_{\varphi(j)} \eta(i)$, ligne par ligne,
pour la partie maillée\index{TEF}\index{node} :
\begin{verbatim}
    do j=1,mp
    < do inode=1, n_node
      < ;mult[i]=1,n_mult : Bb(k2i_[.i](inode),j)
                           = F_D(nde_eta_fun_[.i](inode))(/ff(j));
      >
    >;
;
! Build TEF B mx from Bd
         call scamat(-dt/2.,Bb,B,n,n,m,m);
\end{verbatim}
la macro {\bf mult} génère $n_{mult}$ lignes de code calculant chacun
un élément de la matrice. La macro {\bf F\_D(fun)(/ff(j)} va effectuer
la dérivation symbolique de la fonction {\bf fun} par rapport à la variable $ff$.
Le nom générique de la fonction associée à la ième variable d'état d'une
maille (nde\_eta\_fun) est associée a celle définie par le bloc Mortran\index{Mortran}
set\_node\_eta$<>$ du zinit. Ainsi, à la précompilation par Mortran de
ces instruction, la partie fonction est remplacée strictement par
l'expression donnée dans le ZINIT, fonction de $\eta$ et $\varphi$ ou
de leurs symboles (var:) donnés dans le même bloc de ZINIT. L'appel à {\bf scamat}
est réalisé pour passer à la matrice $B$ du système discrétisé en temps du
TEF (multiplication par $\delta t/2$).

Fortran généré (cfs/principal.f) :
\begin{verbatim}
10125 DO 10129 j=1,mp
      DO 10131 inode=1,n_node
      Bb(k2i_1(inode),j)=(0.)
      Bb(k2i_2(inode),j)=( (((-F_delta(k2j_1(inode+1),j))+F_delta(k2j_1(
     *inode),j))+((-F_delta(k2j_2(inode+1),j))+F_delta(k2j_2(inode),j)))
     **(rmassm1(inode)))
10131 CONTINUE
10132 CONTINUE
10129 CONTINUE
\end{verbatim}
On reconnait les fonctions dérivées de deta\_move et deta\_speed
du ZINIT des masselottes\index{masselottes}. La dérivation est faite sans prise en compte
de l'indice de la variable ff, et la dérivée est multipliée par la fonction
déclarative (Fortran ``statement function'') {\bf F\_Delta(k,l)} qui vaut un pour $k= l$,
zéro sinon.
\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Structure des répertoire, installation et lancement\\ \Minik en résumé}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\emph{Partie reprise d'un {\bf Help} ancien et peut-être redondante}.\index{CMZ (CERN)}\index{Lorenz}
\index{Etat}\index{pousse (exo)}\index{adjoint}\index{SVD}\index{SLT, SLTC}\index{loi de commande}\index{optimisation}
%\footnotesize
\begin{verbatim}
La structure informatique de mini_ker est la suivante :

           .-----------------------------------------------.
           |  .../mini_ker/                                |
           |                                               |
           |   cmz file contiend la partie commune         |
           |   aux exos. En mode exo, elle reste inchangée.|
           `-----------------------------------------------'
                          dessous :
                  /------------------------\
                  |        |        |       |
                /exo1    /exo2  [/lorhcl /pousse]
\end{verbatim}
%\normalsize
       pour chacun : une cmz file avec :
\begin{itemize}
\item      le patch zinproc qui contient toute l'information du modèle
         [calcul direct trajectoire, sensibilités, adjoint
\item      patch sltc : SVD de la matrice d'évolution des états
                      (Singular Value Decomposition).
\item      patch sltcirc : calcul du propagateur et de ses éléments singuliers.
\item      command : gestion de recherche d'une commande optimale hors Minuit.
         [ce dernier est a l'état embryon correspondant à l'exo pousse]
\item      minuik : recherche de maximum de fonctionnelle par {\bf Minuit}
\end{itemize}
\subsection*{Démarrage après réception}
%--------------------------------------
Avant de lancer cmz\index{CMZ (CERN)}, on édite {\bf cmzlogon} pour spécifier un éditeur :
\begin{itemize}
\item par défaut à Jussieu, c'est {\bf ce cv} (Apollo like)
\item en décommentant les deux commandes xedit, on invoque cet éditeur
\end{itemize}
{\bf attention, les options qui sont éditeur-dépendantes sont capitales,
   il faut que cmz reste en attente lors de l'édition d'un fichier.}

on lance cmz (les helps sont en ligne).
se familiariser avec un des exemples fournis :
 pour éxécution de calcul trajectoire : {\bf exe mod} (ou mod)

  crée {\bf exo/exo, exo/exo/cfs}

    on lance l'exécution dans exo/exo/.  [exo.exe $>$res.lis]

    on obtient un ensemble de fichiers file.data\index{data file} avec des kumacs\index{kumac (macro CMZ)}
associées, qui donnent le contenu des fichiers file.

Si on utilise {\bf PAW} comme grapheur, il suffit de lancer une kumac
(ex : {\bf exe eta} [ou eta], liste des kumacs par macro/list)\index{kumac (macro CMZ)} pour obtenir les
trajectoires.
\emph{la mémoire sature pour les grandes trajectoires, il
peut être nécessaire de faire $>$hi/del * entre deux kumac}. Un message
donne un exemple pour sortir un ntuple, genre :
\begin{verbatim}
 > NT/PL 51.wp1%t !  1000    1 ! L !   [ou >help nt/pl]
\end{verbatim}
on peut faire au max du 4D :
\begin{verbatim}
 > NT/PL 51.x%y%z%t !  1000    1 ! * !
\end{verbatim}
pour obtenir un portrait de phase 3D avec le t $\rightarrow$ couleur.
Et (rappel) pour avoir une palette arc-en-ciel :
\begin{verbatim}
>set ncol 32
>palette 1
\end{verbatim}
\subsection*{Applications persos}
%--------------------------------
Faire son exo :(attention, le tout prend au moins 1$/$4h avant de sortir
les résultats...)
On commence par copier une cmz\index{CMZ (CERN)} ressemblante, changer le nom et modifier
la {\bf cmzlogon.kumac} (ou en passant par un .car, interroger les cmz fans).

dans le patch zinproc, on modifie les séquences pour :
\begin{enumerate}
\item donner la dimention\index{dimension} de eta phi  \$DIMETAPHI (*)
\item rentrer son modèle par les fun\_set ou équivalentes de \$ZINIT etc
   on a aussi à choisir des flags logiques : 
\begin{itemize}
   \item {\bf Zback}\index{ZBack} si on veut un aller-retour avec adjoint et dans ce cas,
    \begin{itemize}
    \item {\bf Zcommand} si on gère une loi de commande\index{loi de commande}
    \item {\bf Zlaw} si on explicite une loi (les valeurs sont alors
      inscrites dans les fichiers uxcom et uycom. Sinon, on lit ces fichiers
      en entrée; ils sont donc censés avoir été remplis par une procédure de
      recherche d'optimum (cf patch command).
    \end{itemize}
  \item Rq \$ZSTEER comme dans ZOOM\index{ZOOM}, se situe en fin de boucle (gestion du pas
   de temps par exemple ou des sorties / modzprint)\index{modzprint, Zprint}
\end{itemize}
\item le lancer comme auparavant.
\end{enumerate}
\emph{remarque sur le fonctionnement : avec {\bf mod}, on éxécute la kumac
selseq.kumac\index{selseq (macro CMZ)} qui précise qu'on va chercher les sources de ../ qui est
le noyau de calcul du \Minik}.

\subsection*{Sorties diverses}
%..............................
Dans {\bf exo/exo/cfs}, on a les codes Mortran\index{Mortran} et Fortran générés.

- le fichier exo\_o.tmp est le plus explicite, il contient le programme
"principal" en mtn. {\bf principal.f} est aussi a lire pour i) vérifier la valeur
donnée aux paramètres de dimentionnement losque l'on en sollicite le calcul automatique; 
ii) vérifier le calcul
symbolique des dérivées et opérateur du TEF. Les deux sont utiles pour
qui veut rapidement se familiariser avec le \Minik. Les symboles utilisés
sont très proches du papier ``Adjoint''.

Comme toujours, au lieu de se précipiter sur le graphique, on a intérêt
à consulter un peu les sorties {\bf res.lis}. Et en particulier voir si on
à un message de warning sur ffl (phi). Un test est en effet effectué à chaque incrément
pour vérifier que {\bf ffl(time-dt)+dPhi} reste proche de {\bf ffl(eta(time))}\index{time, time-step}.
Sinon, ça veut dire que pb de linéarité\index{Test de linéarité} (diminuer alors dt). {\bf Remarque : }Ce dernier
critère est nécessaire, malheureusement pas suffisant, cala dépend de la manière dont
les non-linéarités sont distribuées entre cellules\index{cellules} et transferts\index{transferts}. Il reste donc
conseillé d'effectuer le test de comparaison des résultats par variation du pas de temps avec passage
en double précision.
%
\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Validation numérique des algorithmes}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
La validation des méthodes de calcul de la simulation a été faite en comparant
le modèle de Lorenz 63\index{Lorenz} en simulation effectuée avec ZOOM\index{ZOOM} sur station Apollo
et \Minik sous Linux. On trouve six décimales identiques au bout de 3000
pas de temps (time=s) avant que les résultats ne divergent. Cela constitue une
excellente validation, restituant en particulier le même attracteur que celui
abondamment décrit dans la littérature.

Il est également possible de mener une intéressante comparaison entre les trois manières
de calculer une sensibilité dans \Minik :
\begin{itemize}
\item avec la matrice $A^{spha}$ et le calcul du propagateur sur le SLTC\index{SLT, SLTC}, on a la
matrice de sensibilité des états\index{Etat} $(t)$ aux valeurs initiales;
\item l'option Zsensib\index{ZSensib}\footnote{aujourd'hui obsolete.} permet de calculer cette matrice colonne par colonne;
\item et ZBack\index{ZBack}, avec une fonction finale de cout restreinte à un seul état donne
une ligne du propagateur.
\end{itemize}
Voici ces trois matrices (transposées)\index{data file}\index{masselottes} pour un autre exercice
avec option {\bf Grid1D} :
\scriptsize
\begin{verbatim}                                 
========================================================================================================
         Comparaisons entre trois sensibilites / Masselottes Thu Oct 21 18:11:15 CEST 2004 Al1
========================================================================================================
sens.data (5000 calculs)
1 0.10000E+01  0.00000E+00  0.13165E+01 -0.42256E-01  0.15780E+01 -0.69412E-01  0.17284E+01 -0.79530E-01
2 0.00000E+00  0.10000E+01  0.13165E+00 -0.42258E-02  0.15780E+00 -0.69412E-02  0.17284E+00 -0.79530E-02
3 0.00000E+00  0.00000E+00 -0.55054E-01  0.15101E-01 -0.11106E+00  0.17038E-01 -0.15040E+00  0.10118E-01
4 0.00000E+00  0.00000E+00 -0.42257E-01 -0.53544E-01 -0.69413E-01 -0.10936E+00 -0.79531E-01 -0.14939E+00
5 0.00000E+00  0.00000E+00 -0.11106E+00  0.17038E-01 -0.20546E+00  0.25219E-01 -0.26147E+00  0.27156E-01
6 0.00000E+00  0.00000E+00 -0.69413E-01 -0.10936E+00 -0.12179E+00 -0.20294E+00 -0.14894E+00 -0.25875E+00
7 0.00000E+00  0.00000E+00 -0.15041E+00  0.10119E-01 -0.26147E+00  0.27156E-01 -0.31652E+00  0.42257E-01
8 0.00000E+00  0.00000E+00 -0.79531E-01 -0.14939E+00 -0.14894E+00 -0.25875E+00 -0.19120E+00 -0.31229E+00
vadj (5000 calculs)
  1*           2*           3*           4*            5              6*            7             8*
  0.10000E+00  0.00000E+00  0.13166E+01 -0.42131E-01  0.15788E+01 -0.69155E-01  0.17294E+01 -0.79180E-01
  0.00000E+01  0.10000E+01  0.13166E+00 -0.42131E-02  0.15788E+00 -0.69156E-02  0.17294E+00 -0.79181E-02
  0.00000E+00  0.00000E+00 -0.55089E-01  0.15106E-01 -0.11127E+00  0.17001E-01 -0.15053E+00  0.10024E-01
  0.00000E+00  0.00000E+00 -0.42131E-01 -0.53579E-01 -0.68061E-01 -0.10940E+00 -0.77685E-01 -0.14943E+00
  0.00000E+00  0.00000E+00 -0.11110E+00  0.17001E-01 -0.20577E+00  0.25130E-01 -0.26180E+00  0.27025E-01
  0.00000E+00  0.00000E+00 -0.69155E-01 -0.10940E+00 -0.11928E+00 -0.20301E+00 -0.14575E+00 -0.25883E+00
  0.00000E+00  0.00000E+00 -0.15043E+00  0.10024E-01 -0.26180E+00  0.27025E-01 -0.31704E+00  0.42131E-01
  0.00000E+00  0.00000E+00 -0.79179E-01 -0.14943E+00 -0.14575E+00 -0.25883E+00 -0.18734E+00 -0.31241E+00

 phit (5000 calculs) sltcirc
  0.10000E+01  0.00000E+00  0.13165E+01 -0.42254E-01  0.15780E+01 -0.69421E-01  0.17284E+01 -0.79551E-01
  0.00000E+00  0.10000E+01  0.13165E+00 -0.42249E-02  0.15780E+00 -0.69421E-02  0.17284E+00 -0.79554E-02
  0.00000E+00  0.00000E+00 -0.55048E-01  0.15090E-01 -0.11106E+00  0.17035E-01 -0.15040E+00  0.10131E-01
  0.00000E+00  0.00000E+00 -0.42264E-01 -0.53540E-01 -0.69431E-01 -0.10936E+00 -0.79559E-01 -0.14941E+00
  0.00000E+00  0.00000E+00 -0.11106E+00  0.17035E-01 -0.20546E+00  0.25222E-01 -0.26146E+00  0.27166E-01
  0.00000E+00  0.00000E+00 -0.69419E-01 -0.10936E+00 -0.12180E+00 -0.20293E+00 -0.14897E+00 -0.25875E+00
  0.00000E+00  0.00000E+00 -0.15041E+00  0.10130E-01 -0.26146E+00  0.27166E-01 -0.31651E+00  0.42255E-01
  0.00000E+00  0.00000E+00 -0.79553E-01 -0.14939E+00 -0.14897E+00 -0.25875E+00 -0.19123E+00 -0.31229E+00

 calcul de Phi (avec Psi) en avant version 1.01 (Juin 2005)
  0.10000E+01  0.00000E+00  0.13165E+01 -0.42258E-01  0.15780E+01 -0.69416E-01  0.17284E+01 -0.79535E-01
  0.00000E+00  0.10000E+01  0.13165E+00 -0.42259E-02  0.15780E+00 -0.69416E-02  0.17284E+00 -0.79536E-02
  0.00000E+00  0.00000E+00 -0.55054E-01  0.15101E-01 -0.11106E+00  0.17039E-01 -0.15040E+00  0.10119E-01
  0.00000E+00  0.00000E+00 -0.42259E-01 -0.53544E-01 -0.69417E-01 -0.10936E+00 -0.79536E-01 -0.14939E+00
  0.00000E+00  0.00000E+00 -0.11106E+00  0.17039E-01 -0.20546E+00  0.25220E-01 -0.26147E+00  0.27158E-01
  0.00000E+00  0.00000E+00 -0.69417E-01 -0.10936E+00 -0.12179E+00 -0.20294E+00 -0.14895E+00 -0.25875E+00
  0.00000E+00  0.00000E+00 -0.15040E+00  0.10119E-01 -0.26147E+00  0.27158E-01 -0.31652E+00  0.42259E-01
  0.00000E+00  0.00000E+00 -0.79536E-01 -0.14939E+00 -0.14895E+00 -0.25875E+00 -0.19121E+00 -0.31229E+00
 idem, mais avec dt doublé, cad 2500 calculs
  0.10000E+01  0.00000E+00  0.13165E+01 -0.42263E-01  0.15780E+01 -0.69399E-01  0.17284E+01 -0.79485E-01
  0.00000E+00  0.10000E+01  0.13165E+00 -0.42263E-02  0.15780E+00 -0.69397E-02  0.17284E+00 -0.79486E-02
  0.00000E+00  0.00000E+00 -0.55067E-01  0.15128E-01 -0.11107E+00  0.17047E-01 -0.15040E+00  0.10087E-01
  0.00000E+00  0.00000E+00 -0.42263E-01 -0.53555E-01 -0.69398E-01 -0.10937E+00 -0.79485E-01 -0.14939E+00
  0.00000E+00  0.00000E+00 -0.11107E+00  0.17047E-01 -0.20547E+00  0.25215E-01 -0.26147E+00  0.27135E-01
  0.00000E+00  0.00000E+00 -0.69398E-01 -0.10937E+00 -0.12175E+00 -0.20294E+00 -0.14888E+00 -0.25876E+00
  0.00000E+00  0.00000E+00 -0.15040E+00  0.10088E-01 -0.26147E+00  0.27135E-01 -0.31654E+00  0.42263E-01
  0.00000E+00  0.00000E+00 -0.79486E-01 -0.14939E+00 -0.14888E+00 -0.25876E+00 -0.19115E+00 -0.31231E+00
========================================================================================================
\end{verbatim}
\normalsize
les colonnes de vadj$^{\dagger}$ sans signe $*$ correspondent à 4999 calculs, pour
calibrer les comparaisons.

Ces mêmes comparaisons ont été effectuées sur deux autres modèles avec le même
score.
On remarque en général que le calcul de sensibilité automatique s'approche le
mieux du propagateur - que l'on peut prendre comme référence du fait de l'ordre
élevé des dévelopements du calcul. Cela se conçoit si l'on remarque que les
calculs avec Zsensib\index{ZSensib} sont effectués en prenant en compte les dérivées des
matrices (par différences finies).

On a aussi comparé ces trois matrices pour l'exemple macro--éco de stepH (dénommé
gras\_double par Al1), exercice non-linéaire avec matrice $D$\footnote{cf document Propagateurs.}.

D'autres validations numériques concernant les différents modes de perturbation des transferts
et le tableau $\Psi(t,0)$ est disponible sur demande. On a pu noter en particulier que le passage
en double précision ne fait gagner en gros que deux ordres de grandeur, et que l'on bute
manifestement sur une complexité algorithmique.

\vspace{1.cm}
\hrule
%
%\begin{figure}[h]
%\centerline{\psfig{figure=gain2sltc.eps,%
%bbllx=2pt,bblly=2pt,bburx=386pt,bbury=294pt,width=10.5cm,clip=}}
%\caption{Couplage de $\psi$ sur $\varphi$ : deux coupures.}
%\label{couplfig}
%\end{figure}
%
\printindex
\tableofcontents
\end{document}
\end

