source: tag/release-8/yao8/product/include/Dyntst.h @ 1

Last change on this file since 1 was 1, checked in by lnalod, 15 years ago

Initial import of YAO sources
  • Property svn:eol-style set to native
File size: 52.4 KB
Line 
1/* ----------------------------------------------------------------------------*/
2/* ----------------------------------------------------------------------------*/
3//char Ydxalea[STRSIZE20+1] = "R0.0000001"; pour un dx (perturbation de Xo) aleatoire ...
4char Ydxalea[STRSIZE20+1] = "\0"; // mais initialement (par defaut) , cet alea est neutralised
5/* ----------------------------------------------------------------------------*/
6/* ----------------------------------------------------------------------------*/
7void Ytestad_mod(double pdx, double errelmax, int maxko, double pzedi)
8{                /* PURPOSE : test de l'ADjoint par module : <m'(Xo)x,y> = <x, m*(Xo).y>
9                    pdx :       Pourcentage De perturbation de Xo pour obtenir dx
10                    errelmax : Erreur relative max : on test tous les modules en chaque point
11                    et maxko : de coordonnee, et on considere que le test est ko si l'erreur
12                   relative est > a errelmax. Si on trouve plus de maxko tests ko,
13                   le process est interrompu (exit).
14                          pzedi    : Parametres pour se donner un niveau de zero informatique
15                                           lorsque les 2 |termes| sont < pzedi alors -> pas de test => OK !?
16
17        nb: maintenir en meme temps la fonction testad (s'il y a lieu)
18
19        called fct : YgradEQPstate_all, YdeltaEQPstate_all,YgradTOtab_all, YtabTOgrad_all
20                     before_it, Yforward, Ylinward, Yrazgrad_all, Ybackward
21                 */
22
23    int   witraj;
24    int   NBALLTIME;
25                YREAL *MXdx = new YREAL[YSIZECO];               /* tableau de sauvegarde */
26                if (MXdx==NULL){printf("Yao=>ABORT: pb on new MXdx for testad (see Yao system administrator)\n"); exit(-9);}
27          memset(MXdx, 0, YSIZECO*sizeof(YREAL));
28                Ytop0 = time((time_t *)NULL);
29
30                /* 0) Presentation */
31printf("\ntestad : NEW : use of new function set_dxalea (%s)\n", Ydxalea);
32                printf ("\n module Adjoint TEST :        dx = Xo * pdx  (with pdx = %e)", pdx);
33                printf ("\n --------------------    and  dy = m'(Xo).dx\n\n");
34
35                /* nb: On suppose Xo initialised par l'utilisateur (setstate ou xuserfct -> YSo) */
36
37                /* 1) some init */
38                Ytestad_module=1;
39                Yerrelmax=errelmax;
40                Ynbko=0; Ymaxko=maxko;
41                Ypzedi = pzedi;
42
43    YgradEQPOstate_target(pdx, Ydxalea); //Y'avait tout un baratin ici,
44    YdeltaEQPgrad_target(1.0);           //il faut le refaire ... ?
45
46                /* 2) MD (Modele Direct): Calculer M(Xo) */
47                Yset_modeltime(0);
48                before_it(1);                                //(YItRun);
49                Yforward (-1, 0);                    //MD : M(Xo) -> YS
50    Yact_operator('*', ON);
51    Yforward_operator('H');                                                  //forward operateurs a la fin; d'abord H seulement,
52    for (witraj=0; witraj<YNBTRAJ; ++witraj)                                                                                                                       //puis il faut mettre YW = YS pour les modules de
53    {   NBALLTIME = YTabTraj[witraj].nbsteptime + YTabTraj[witraj].nbuptime; //(M ou H) connectes vers des modules operateurs
54                    YwishEQPstate_traj_tocov(witraj, NBALLTIME-1, 1.0);                                                                  //de covariance (R ou K) ... puisque ... sinon ...
55    }
56    Yforward_operator('R');             //et on termine par les
57    Yforward_operator('K');     //    operateurs R et K
58
59                /* 3) LT (Lineaire Tangent): m'(Xo)dXo */
60                Yset_modeltime(0);
61                before_it(1);                                //(YItRun);
62                Ylinward(0);                                 // TL : M'(Xo)dXo -> YG
63    Ylinward_operator('*');  //+linward des operateurs a la fin
64
65                /* 4) YGn-> tableau Mxdx: avant RAZ, sauvegarde des dx obtenus sur les modules cout en fin de traj */
66    Y3windice=0;
67    for (witraj=0; witraj<YNBTRAJ; ++witraj)
68    {   NBALLTIME = YTabTraj[witraj].nbsteptime + YTabTraj[witraj].nbuptime;
69        YgradTOtab_traj(witraj, "Y#C", NBALLTIME-1, NBALLTIME , MXdx);                  /* YGn -> [MXdx] */
70    }
71
72                /* 5) RAZ tous les gradients pour l'Adjoint */
73                Yrazgrad_all();         /* YG = 0 */
74
75                /* 6) et Reprise des dx de fin de trajectoire des modules cout (qui vont etre retro-propagees) */
76    Y3windice=0;
77    for (witraj=0; witraj<YNBTRAJ; ++witraj)
78    {   NBALLTIME = YTabTraj[witraj].nbsteptime + YTabTraj[witraj].nbuptime;
79        YtabTOgrad_traj(witraj, "Y#C", NBALLTIME-1, NBALLTIME , MXdx);                          /* [MXdx] -> YGn */
80    }
81
82                /* 7) on se lance dans l'Adjoint */
83    Ybackward_operator('*'); //d'abord backward des operateurs en fin de trajectoire
84                Ybackward(-1, 0);        //au cours duquel le test sur chaque module se fait (YNBSTEPTIME);
85                //quid after_it ?
86
87                printf(" -->-->--> %i case(s) found \n\n", Ynbko);
88                /* 8) FIN */
89                /*    on libere la memoire               */
90                delete[] MXdx;
91                /*    et on laise le temps propre en ordre */
92                Yset_modeltime(0);
93
94                /* Restaurer (~) l'etat init ? En principe il n'est pas modified !?
95                for (witraj=0; witraj<YNBTRAJ; ++witraj)
96        YstateEQPdelta_traj(witraj, "Y#A", 0, YTabTraj[witraj].nbuptime, 1/pdx);
97                */
98                Ytestad_module=0;
99}
100
101/* ----------------------------------------------------------------------------*/
102/* ----------------------------------------------------------------------------*/
103void Ytestad(double pdx)
104{                /* PURPOSE : test de l'ADjoint : <M'(Xo)x,y> = <x, M*(Xo).y>
105                                        ici, dx correspond a x, et dx* a y
106                                        dans Yao la zone utilised pour Xo est YS, et YG est utilised pour dx lors
107                                        du TL (passe avant) puis pour dx* lors de l'Adjoint (passe arriere) d'ou
108                                        la necessite de passer par des tableaux pour sauvegarder et manipuler les
109                                        resultats des calculs.
110                    pdx :       Pourcentage De perturbation de Xo pour obtenir dx
111                    errelmax : Erreur relative max : si cette valeur est differente de 0 alors
112                    et maxko : on test tous les modules en chaque point de coordonnee, et on
113                               considere que le test est ko si l'erreur relative est > a errelmax.
114                               Si on trouve plus de maxko tests ko, le process est interrompu (exit).
115
116                                nb: maintenir en meme temps la fonction testad_mod (s'il y a lieu)
117
118                                called fct : YgradEQPstate_all, YgradTOtab_all, YtabTOgrad_all
119                                                                           Yrazgrad_all, before_it, Yforward, Ylinward, Ybackward
120                 */
121
122                int   witraj;
123                int   NBALLTIME;
124                double LTRes=0.0, AdRes=0.0;
125                YREAL *DX   = new YREAL[YSIZEPB];               /* tableau de sauvegarde de dx */
126                if (DX==NULL)  {printf("Yao=>ABORT: pb on new DX for testad (see Yao system administrator)\n"); exit(-9);}
127                YREAL *MXdx = new YREAL[YSIZECO];               /* tableau resultat du Lineaire Tangent M'(Xo)dX */
128                if (MXdx==NULL){printf("Yao=>ABORT: pb on new MXdx for testad (see Yao system administrator)\n"); exit(-9);}
129                YREAL *MAdx = new YREAL[YSIZEPB];               /* tableau resultat de l'Adjoint M*(Xo)d*X */
130                if (MAdx==NULL){printf("Yao=>ABORT: pb on new MAdx for testad (see Yao system administrator)\n"); exit(-9);}
131
132          memset(DX,   0, YSIZEPB*sizeof(YREAL));
133          memset(MXdx, 0, YSIZECO*sizeof(YREAL));
134          memset(MAdx, 0, YSIZEPB*sizeof(YREAL));
135
136                Ytop0 = time((time_t *)NULL);
137
138                /* 0) Presentation */
139printf("\ntestad : NEW : use of new function set_dxalea (%s)\n", Ydxalea);
140                printf ("\n Adjoint TEST :         dx = Xo * pdx  (with pdx = %e)", pdx);
141                printf ("\n --------------    and  dy = M'(Xo).dx\n\n");
142
143                Ytestad_module=0;               /* because c'est la meme instruction qui sert dans les 2 cas, il faut pouvoir distinguer */
144
145                /* nb: On suppose Xo initialised par l'utilisateur (setstate ou xuserfct -> YSo) */
146
147                /* 1) Init de dXo -> YGo */
148                Yrazgrad_all();
149
150    /* dXo = Xo * ~pdx (([1->UPTIME]:)YGo_mod=YSo_mod*pdx) */
151    YgradEQPOstate_target(pdx, Ydxalea); //calcul de dx avec un alea ...
152
153                /* 2) MD (Modele Direct): Calculer M(Xo) */
154                Yset_modeltime(0);
155                before_it(1);                           //(YitRun);
156                Yforward(-1, 0);                //MD : M(Xo) -> YS
157    Yact_operator('*', ON);
158    Yforward_operator('H');                                                  //forward operateurs a la fin; d'abord H seulement,
159    for (witraj=0; witraj<YNBTRAJ; ++witraj)                                                                                                                       //puis il faut mettre YW = YS pour les modules de
160    {   NBALLTIME = YTabTraj[witraj].nbsteptime + YTabTraj[witraj].nbuptime; //(M ou H) connectes vers des modules operateurs
161                    YwishEQPstate_traj_tocov(witraj, NBALLTIME-1, 1.0);                                                                  //de covariance (R ou K) ... puisque ... sinon ...
162    }
163    Yforward_operator('R');  //et on termine par les
164    Yforward_operator('K');  //    operateurs R et K
165
166                /* 3) LT (Lineaire Tangent): M'(Xo)dXo */
167                Yset_modeltime(0);
168                before_it(1);                           //(YitRun);
169                Ylinward(0);                            //TL : M'(Xo)dXo -> YG
170    Ylinward_operator('*'); //linward opera la fin !!?? (car on ne s'interesse qu'aux modules cout en fin de traj)
171
172                /* 4) YGn-> tableau Mxdx et calcul de < M'(Xo)dX, dy > (sur module cout en fin de traj) */
173                Y3windice=0;
174    for (witraj=0; witraj<YNBTRAJ; ++witraj)
175    {   NBALLTIME = YTabTraj[witraj].nbsteptime + YTabTraj[witraj].nbuptime;
176        YgradTOtab_traj(witraj, "Y#C", NBALLTIME-1, NBALLTIME , MXdx);                                  /* YGn -> [MXdx] */
177    }
178
179                for(int i=0; i<YSIZECO; ++i)
180                {LTRes += MXdx[i]*MXdx[i];      /* < M'(Xo)dX, dy > */
181                }
182
183                /* 6) au passage, on a peut etre aussi besoin de conserver donc de sauvegarder YGo (isnt-it?) */
184                Y3windice=0;
185    YgradTOtab_target(DX);                      /* YGo -> [DX]  */
186
187                /* 7) RAZ tous les gradients pour l'Adjoint */
188                Yrazgrad_all();         /* YG = 0 */
189
190                /* 8) et Reprise des gradients de fin de trajectoire des modules cout */
191                Y3windice=0;
192    for (witraj=0; witraj<YNBTRAJ; ++witraj)
193    {   NBALLTIME = YTabTraj[witraj].nbsteptime + YTabTraj[witraj].nbuptime;
194        YtabTOgrad_traj(witraj, "Y#C", NBALLTIME-1, NBALLTIME , MXdx);                          /* [MXdx] -> YGn */
195    }
196//
197                /* 9) on devrait pouvoir se lancer dans l'adjoint */
198    Ybackward_operator('*'); //d'abord backward des operateurs en fin de trajectoire
199                Ybackward(-1, 0);  /* AD (adjoint):-> dx =M*(Xo)d*X : Yjac * YG -> Ytbeta */
200                //quid after_it ?
201
202                /* 10) YGo-> tableau MAdx et calcul de < dx, M*(Xo).dy > */
203                Y3windice=0;
204    YgradTOtab_target(MAdx);             /* YGo -> [MAdx] */
205
206                for(int i=0; i<YSIZEPB; ++i)
207                {       AdRes += DX[i] * MAdx[i];                                                                       /* < dx, M*(Xo).dy > */
208                }
209
210                /* 11) Afficher le resultat */
211                printf (" < M'(Xo).dx, dy >  = %- 16.15e \n", LTRes);
212                printf (" < dx, M*(Xo).dy >  = %- 16.15e \n", AdRes);
213                printf ("                      ----------------------\n");
214                printf ("   ecart :            %- 16.15e \t", LTRes-AdRes);
215                printf (" relative error = %e \n\n", (LTRes-AdRes)/LTRes);
216
217
218                /* 12) on informe d'une possible cause d'echec par exemple si des obs sont */
219                if (YioInsertObsCtr != -1) /* charger dans une des fonction d'intervention */
220                {        printf("testad warning: obs must not be used for global Adjoint test\n");
221                           //because: ca modifie le dy au cours de la redescente qui ne sera plus le meme
222                           //         que celui utilised lors du calcul de LTRes.
223                           //         Cette contrainte n'est pas utile pour testad_mod car le test est local.
224                }
225
226                /* 13) FIN */
227                /*    on libere la memoire               */
228                delete[] DX;  delete[] MXdx;  delete[] MAdx;
229                /*    et on laise le temps propre en ordre */
230                Yset_modeltime(0);
231}
232
233/* ----------------------------------------------------------------------------*/
234/* ----------------------------------------------------------------------------*/
235void Ytestlt(double pdx, double alpha, double fdec, int nbloop, int modop, int syclop)
236{                                               //PURPOSE : test du Lineaire Tangent :
237                                                // lorsque dx~>0:  M(Xo+dx) - M(Xo) ~> M'(Xo)dx             (Taylor)
238                                                //                                                               M(Xo+dx) - M(Xo) -  M'(Xo)dx ~> 0
239                                                // 2eme maniere equivalente, mais (+)difficile a implementer a l'ordre 2 (!) :
240                                                //                                 [ (M(Xo+dXo) - M(Xo)) / M'(Xo)dx ] ~> 1
241                                                //
242                //pdx   : Pourcentage De perturbation de Xo pour obtenir dx
243                //alpha : parametre de tendance
244                //fdec  : facteur de decroissance de alpha (fdec=#10)
245                //nbloop: nombre de boucles de decroissance
246                //modop : mode operatoire: si 0: ->0, sinon ->1
247                //syclop: Secret Yao Systeme Code Operation (0: algo officiel
248                //        1: algo (initial) theoriquement equivalent a l'algo officiel mais plus lent
249                //                        2: algo (initial) bugged (i.e sans re-init de Xo ...) mais qui parfois "marche mieux" !!??
250               
251                int     witraj;
252                int     NBALLTIME;
253                double  apdx, svalpha;                  //...
254                double  *TrendO1, *TrendO2;     //Tableau des rms-tendances par module
255                int                     indT;                                                           //indice des TrendO tableaux
256                int             n                  = 0;                                 //indice (supposed) correspondant dans les tableaux de sauvegarde
257                int             noco   = 0;                                     //numero d'ordre des cout modul (parmi les modules cout)
258                int             sizeco;                                                 // taille atempo d'un module cout
259                double  galphaO1=0.0;                           //...
260                double  galphaO2=0.0;                           //...   
261                double  sigO1, sigO2;                   //...
262                int                     nbnanf;                                           //nombre de NaN ou Inf found pour le LT si = 0 (pour la 2eme maniere) ..
263                int     *Tnbnanf;           //Tableau des nombres de nanf par module
264                int                     isnanf;                                                 //booleen pour tester que les resultats de calcul sont finite (ie: !Nan && !Inf)
265               
266                // Tableaux de staokage des valeurs et resultats de calcul ...
267                YREAL *MX = new YREAL[YSIZECO];
268                if (MX==NULL){printf("Yao=>ABORT: pb on new MX for testlt (see Yao system administrator)\n"); exit(-9);}
269                YREAL *TLX = new YREAL[YSIZECO];
270                if (TLX==NULL){printf("Yao=>ABORT: pb on new TLX for testlt (see Yao system administrator)\n"); exit(-9);}
271                YREAL *MXdx = new YREAL[YSIZECO];
272                if (MXdx==NULL){printf("Yao=>ABORT: pb on new MXdx for testlt (see Yao system administrator)\n"); exit(-9);}
273                YREAL *MpXdx = new YREAL[YSIZECO];
274                if (MpXdx==NULL){printf("Yao=>ABORT: pb on new MpXdx for testlt (see Yao system administrator)\n"); exit(-9);}
275
276          memset(MX,   0, YSIZECO*sizeof(YREAL));
277          memset(TLX,  0, YSIZECO*sizeof(YREAL));
278          memset(MXdx, 0, YSIZECO*sizeof(YREAL));
279          memset(MpXdx,  0, YSIZECO*sizeof(YREAL));
280
281                Ytop0 = time((time_t *)NULL);
282
283                /* A) INITIALISATION */
284                /*   0) presentation : */
285                printf ("\n Lineaire Tangent TEST : with: pdx = %e, (a)lpha = %e", pdx, alpha);
286                printf ("\n -----------------------       fdec= %e, (n)bloop= %i\n", fdec, nbloop);
287                printf (" nb 1) dx = Xo * pdx \n");
288                printf ("    2) a(n) = a(n-1) / fdec \n");
289               
290                /*   1) determination du nombre de module cout : */
291                int     nbcout=0;
292                for (int w1=0; w1<YNBMODUL; ++w1) nbcout += YTabMod[w1].is_cout;
293               
294                /* allocation des tableaux de tendance pour chaque module */
295                TrendO1 = new double[nbcout*nbloop];
296                TrendO2 = new double[nbcout*nbloop];           
297                Tnbnanf = new int[nbcout*nbloop];               
298                               
299                /*   2) On suppose Xo initialised par l'utilisateur (setstate ou xuserfct -> YS) */
300                //for (witraj=0; witraj<YNBTRAJ; ++witraj)      /* sauvegarde Xo (i.e: YSo) dans YDo: (([1->UPTIME]:)YDo_mod=YSo_mod) */       
301    //    YdeltaEQPstate_traj(witraj, "Y#T", 0,  YTabTraj[witraj].nbuptime, 1);
302    YdeltaEQPstate_target(1);
303               
304                /*   3) sauvegarde du parametre de tendance */
305                svalpha = alpha;                               
306               
307                /*   4) Calcul de M(Xo) et sauvegarde */
308                Yset_modeltime(0);
309                before_it(1);      //... (YItRun);
310                Yforward (-1, 0);        //MD : M(Xo) -> YS
311    Yact_operator('*', ON);
312    Yforward_operator('H');                                                  //forward operateurs a la fin; d'abord H seulement,
313    for (witraj=0; witraj<YNBTRAJ; ++witraj)                                                                                                                       //puis il faut mettre YW = YS pour les modules de
314    {   NBALLTIME = YTabTraj[witraj].nbsteptime + YTabTraj[witraj].nbuptime; //(M ou H) connectes vers des modules operateurs
315                    YwishEQPstate_traj_tocov(witraj, NBALLTIME-1, 1.0);                                                                  //de covariance (R ou K) ... puisque ... sinon ...
316    }
317    Yforward_operator('R');  //et on termine par les
318    Yforward_operator('K');  //    operateurs R et K
319
320                Y3windice=0;
321    for (witraj=0; witraj<YNBTRAJ; ++witraj)                                                                                                                                    //YSn -> MX     (nb: que module cout
322    {   NBALLTIME = YTabTraj[witraj].nbsteptime + YTabTraj[witraj].nbuptime;    // en fin de trajectoire, qu'il y
323        YstateTOtab_traj(witraj, "Y#C", NBALLTIME-1, NBALLTIME , MX);                           // ait ou pas d'observations)                   
324    }
325
326                /*   5) Calcul de M'(Xo)dx une seule fois et sauvegarde : (algo officiel) */
327                if (syclop==0)
328                {  Yrazgrad_all();//RAZ de YG : YG = 0
329                   // on va mettre dx comme petite perturbation a propager 1 fois par le LT
330       YgradEQPstate_target(pdx);
331               
332                   // propagation par le LT
333                   Yset_modeltime(0);
334                   before_it(1);            //... (YItRun);
335                   Ylinward(0);                   //M'(Xo)dX -> YG
336       Ylinward_operator('*');  //linward des operateurs en fin de trajectoire
337
338                         Y3windice=0;
339       for (witraj=0; witraj<YNBTRAJ; ++witraj)
340       {   NBALLTIME = YTabTraj[witraj].nbsteptime + YTabTraj[witraj].nbuptime;
341           YgradTOtab_traj(witraj, "Y#C", NBALLTIME-1, NBALLTIME , MpXdx);              //sauvegarde: YGn -> [MpXdx]                   
342       }
343                }
344
345                /* B) BOUCLE DE CALCUL DES TENDANCES PAR MODULE */
346                for (int iloop=0; iloop<nbloop; ++iloop)                //^: nbloop : nombre de boucle
347                {               printf("o");fflush(stdout); //evolution loop
348                               
349                                /* 1) evolution des p* */
350                                alpha = alpha / fdec;                           //decroissance de alpha
351                                apdx  = pdx * alpha;
352
353                                /* 2) LT (Lineaire Tangent): M'(Xo)a.dXo */
354                                if (syclop==0) //algo officiel : cette boucle suffit pour avoir M'(Xo)a.dXo -> [TLX]
355                                {  for(int wi=0; wi<YSIZECO; ++wi) TLX[wi] = MpXdx[wi]*alpha;
356                                }                                                               
357                                if (syclop>0)  //algo non officiel (plus lent, eventuellement bugged)
358                                {        if (syclop==1) // avec cet algo,  il faut retablir YS (qui est cassed en fin de boucle)
359                                   {      YstateEQPdelta_target(1); //on restaure Xo (i.e: YSo) (([1->UPTIME]:)YSo_mod=YDo_mod)
360                                                                Yset_modeltime(0);
361                                                                before_it(iloop+1);          //...(YItRun);
362                                                                Yforward (-1, 0);                    //MD : M(Xo) -> YS
363                                                          Yforward_operator('H');                                                  //forward operateurs a la fin; d'abord H seulement,
364                                                                for (witraj=0; witraj<YNBTRAJ; ++witraj)                                                                                                                           //puis il faut mettre YW = YS pour les modules de
365                                                                {   NBALLTIME = YTabTraj[witraj].nbsteptime + YTabTraj[witraj].nbuptime; //(M ou H) connectes vers des modules operateurs
366                                                                                YwishEQPstate_traj_tocov(witraj, NBALLTIME-1, 1.0);                                                                      //de covariance (R ou K) ... puisque ... sinon ...
367                                                                }
368                                                                Yforward_operator('R'); Yforward_operator('K');         //et on termine par les operateurs R et K
369                                         } // sinon (syclop=2) on est dans la version bugged qui parfois marche mieux
370
371                                         // calcul de a.dx et propagation par le LT
372                                         Yset_modeltime(0); 
373                                         Yrazgrad_all();                  //RAZ de YG : YG = 0
374                                       
375           //for (witraj=0; witraj<YNBTRAJ; ++witraj) //a.dXo = Xo * alpha.pdx (([1->UPTIME]:)YGo_mod=YSo_mod*apdx)
376           //     YgradEQPstate_traj(witraj, "Y#T", 0, YTabTraj[witraj].nbuptime, pdx);
377           YgradEQPstate_target(pdx);
378                                                                               
379                                         before_it(iloop+1);        //... (YItRun);
380                                         Ylinward(0);                                       //TL : M'(Xo)a.dXo -> YG
381           Ylinward_operator('*');  //linward operateurs la fin
382
383                                         Y3windice=0;
384           for (witraj=0; witraj<YNBTRAJ; ++witraj)
385           {   NBALLTIME = YTabTraj[witraj].nbsteptime + YTabTraj[witraj].nbuptime;
386               YgradTOtab_traj(witraj, "Y#C", NBALLTIME-1, NBALLTIME , TLX);        //YGn -> [TLX]                     
387           }
388                                }
389
390                                /* 3) calcul de M(Xo+alpha.dXo) */
391                                //   3.1) calcule de Xo+alpha.dXo
392        YstateEQPdelta_target(1+apdx); //Xo=Xo+alpha.dXo : YDo * (1+alpha.pdx) -> YSo
393
394                    //   3.2) Modele Direct: M(Xo+dXo) et stockage
395                                Yset_modeltime(0);
396                                before_it(iloop+1);      //...(YItRun);
397                                Yforward (-1, 0);    //MD : M(Xo+alpha.dXo) -> YS
398        Yforward_operator('H');                                                  //forward operateurs a la fin; d'abord H seulement,
399        for (witraj=0; witraj<YNBTRAJ; ++witraj)                                                                                                                           //puis il faut mettre YW = YS pour les modules de
400        {   NBALLTIME = YTabTraj[witraj].nbsteptime + YTabTraj[witraj].nbuptime; //(M ou H) connectes vers des modules operateurs
401                        YwishEQPstate_traj_tocov(witraj, NBALLTIME-1, 1.0);                                                                      //de covariance (R ou K) ... puisque ... sinon ...
402        }
403        Yforward_operator('R');  //et on termine par les
404                          Yforward_operator('K');  //    operateurs R et K
405
406                                Y3windice=0;
407        for (witraj=0; witraj<YNBTRAJ; ++witraj)
408        {   NBALLTIME = YTabTraj[witraj].nbsteptime + YTabTraj[witraj].nbuptime;
409            YstateTOtab_traj(witraj, "Y#C", NBALLTIME-1, NBALLTIME , MXdx);             //YSn -> MXdx
410        }
411
412                                /* 4) on calcule les rms-tendances par module, et on les stocke */
413                                //   (nb: pour un meme module, toutes les sorties sont supposed homogenes)                             
414                                n                = 0; //indice (supposed) correspondant dans les tableaux de sauvegarde
415                                noco = 0;       //numero d'ordre des cout modul (parmi les modules cout)                       
416                                for (int w1=0; w1<YNBMODUL; ++w1)
417                                {               if (YTabMod[w1].is_cout)
418                                        {       
419                                                        sizeco = YTabMod[w1].nb_stout * YTabMod[w1].axi;                                                //   taille
420                                                if (YTabMod[w1].axj>0) sizeco = sizeco * YTabMod[w1].axj;               //   atempo
421                                                if (YTabMod[w1].axk>0) sizeco = sizeco * YTabMod[w1].axk;               //  du module
422                                                nbnanf = 0; //pour gerer NANF (NaN, Inf /0)
423                                                sigO1=sigO2=0.0;
424                                                for (int sz=0; sz<sizeco; ++sz)
425                                                {
426                                                                isnanf = (!finite(MXdx[n]) || !finite(MX[n]) || !finite(TLX[n])
427// || fabs(MXdx[n])<1.0e-15 || fabs(MX[n])<1.0e-15 || fabs(TLX[n])<1.0e-15                                                                      //ARAR
428                                                                                         );
429                                                          if (!isnanf)
430                                                          {
431                                                                                        if (modop==0)
432                                                                                        {       galphaO1 = MXdx[n] - MX[n] - TLX[n];                            //(1)ok
433                                                                                        sigO1 += (galphaO1*galphaO1);           
434                                                                                        }
435                                                                                        else                                                                                                                                                                    //(2) not as easy
436                                                                                        {        if (TLX[n] != 0.0)
437
438                                                                                                 {   galphaO1 = (MXdx[n] - MX[n]) / TLX[n];                                                                     
439                                                                                                         galphaO2 = galphaO1 - 1;       
440                                                                                                         sigO1 += (galphaO1*galphaO1);                         
441                                                                                                         sigO2 += (galphaO2*galphaO2);
442                                                                                         }
443                                                                                         else
444                                                                                                         ++nbnanf;                                                                                                             
445                                                                                  }
446                                                                                }
447                                                                                else
448                                                                                        ++nbnanf;                                                                               
449                                                                               
450                                                                                ++n;                                           
451                                                                }
452
453                                                                indT = (nbloop * noco) + iloop;
454                                                                Tnbnanf[indT]   = nbnanf;                                                               
455
456                                                                if (nbnanf < sizeco)
457                                                                {                                               
458                                                                        if (modop==0)
459                                                                        {       TrendO1[indT]   = sqrt(sigO1/sizeco-nbnanf);
460                                                                                //TrendO1[indT] = galphaO1; //courcircuitage RMS ok if only one component (pour mise au point)
461                                                                        TrendO2[indT]   = TrendO1[indT] / (alpha*alpha);
462                                                                        }
463                                                                        else
464                                                                        {       TrendO1[indT]   = sqrt(sigO1/(sizeco-nbnanf));
465                                                                        TrendO2[indT]   = sqrt(sigO2/(sizeco-nbnanf)) / alpha;
466                                                                        }
467                                                                }
468                                                                else
469                                                                {       Tnbnanf[indT]   = -nbnanf; /* dans ce cas la, toutes les mailles du module sont
470                                                                        nanf; sizeco-nbnanf est sensed etre = 0 ce qui est interdit pour la division;
471                                                                        pour le signaler, on affichera ci-apres un nombre de nanf negatif ! d'ou l'inversion
472                                                                        de signe; et pour faire bonne mesure, on met les Trend a 0. */
473                                                                        TrendO1[indT]=TrendO2[indT]=0.0;
474                                                                }                                                               
475                                                                                                               
476                                                ++noco;  //numero d'ordre (parmi les cout) des modules cout                             
477                                        }//fin if module is_cout
478                                }//fin boucle sur les modules                                                                                                                   
479                }//fin boucle des calculs avec decroissance des tendances
480                               
481                /* C) AFFICHAGE DES RESULTATS PAR MODULE */
482                if (modop==0)
483                {  printf("\n                             g = M(Xo+a.dx) - M(Xo) - M'(Xo)a.dx\n");
484                   printf("         a decrease     ordre 1: rms(g)-->0    ordre 2: [rms(g)/(a.a)]-->K \n");
485                }
486                else
487                {  printf("\n                             g = (M(Xo+a.dx) - M(Xo)) / M'(Xo)a.dx\n");
488                   printf("         a decrease     ordre 1: rms(g)-->1    ordre 2: [rms(g-1)/a]-->K \n");
489                }
490                noco = 0;       //numero d'ordre des cout modul (parmi les modules cout)                       
491                for (int w1=0; w1<YNBMODUL; ++w1)
492                {               if (YTabMod[w1].is_cout)
493                                {               if (YTabMod[w1].ctrord > 0)
494                                                {               //si le module est dans l'order (ce qui devrait etre le cas pour un module cout !?)
495                                                                printf("\n module cout ->       %s : \n",YTabMod[w1].Name);
496                                                                alpha = svalpha;
497                                                                for (int iloop=0; iloop<nbloop; ++iloop)
498                                                                {               indT = (nbloop * noco) + iloop;
499                                                                                alpha = alpha / fdec;
500                                                                                //printf("   %3i :  %13.6e           %13.6e     \n", iloop+1, TrendO1[indT], TrendO2[indT]);
501                                                                                if (modop==0)
502                                                                                    printf("   :%3i: %- 15.6e : -0-> %- 15.6e    -K-> %- 15.6e",
503                                                                                            iloop+1, alpha, TrendO1[indT], TrendO2[indT]);
504                                                                                else
505                                                                                    printf("   :%3i: %- 15.6e : -1-> %- 15.6e    -K-> %- 15.6e",
506                                                                                            iloop+1, alpha, TrendO1[indT], TrendO2[indT]);
507                                                                                if (Tnbnanf[indT]!=0) printf("  %i NANF", Tnbnanf[indT]);
508                                                                                printf("\n");                                                                           
509                                                                }
510                                                }
511                                                else
512                                                {               printf("\n module ->               %s : is not in the ORDER run over !?\n\n",YTabMod[w1].Name);
513                                                }
514                                                ++noco;
515                                }
516                }
517                               
518                /* D) FIN */
519        //    on restaure, autant que faire ce peut, l'etat initial
520    YstateEQPdelta_target(1);
521                //    on libere la memoire                             
522                delete[]MX;  delete[]TLX;  delete[]MXdx;        delete[]MpXdx;  delete[]TrendO1; delete[]TrendO2;
523                //    et on laise le temps propre en ordre
524                Yset_modeltime(0);
525}
526
527/* ----------------------------------------------------------------------------*/
528/* ----------------------------------------------------------------------------*/
529void Ytestof(double pdx, double alpha, double fdec, int nbloop, int modop, int syclop)
530{               //PURPOSE : test de l'Objective Fonction :
531                //pdx   : pourcentage De perturbation de Xo pour obtenir dx
532                //alpha : parametre de tendance
533                //fdec  : facteur de decroissance de pdec alias alpha
534                //nbloop: nombre de boucles de decroissance
535                //modop : mode operatoire: si 0: ->0, sinon ->1
536                //syclop: Secret Yao Systeme Code Operation (0: algo officiel par defaut)
537                //        1: algo (initial) theoriquement equivalent a l'algo officiel mais plus lent
538
539                double  apdx;
540                double  JXdx;
541                double  trendO1, trendO2;
542                double  DJdx=0.0;
543                double  JX=0.0;
544                YREAL           *Jpx = new YREAL[YSIZEPB];              // tableau pour J'(Xo)dx (avant c'etait pour les a.dx)
545                if (Jpx==NULL){printf("Yao=>ABORT: pb on new Jpx for testof (see Yao system administrator)\n"); exit(-9);}
546                YREAL           *DJ  = new YREAL[YSIZEPB];              // tableau de gradinat YG
547                if (DJ==NULL) {printf("Yao=>ABORT: pb on new DJ for testof (see Yao system administrator)\n"); exit(-9);}
548       
549          memset(Jpx, 0, YSIZEPB*sizeof(YREAL));
550          memset(DJ,  0, YSIZEPB*sizeof(YREAL));
551
552                Ytop0 = time((time_t *)NULL);
553
554                //   0) presentation :
555printf("\ntestof : NEW : use of new function set_dxalea (%s)\n", Ydxalea);
556                printf ("\n Objective Function TEST : with: pdx = %e, (a)lpha = %e", pdx, alpha);
557                printf ("\n -------------------------       fdec= %e, (n)bloop= %i\n", fdec, nbloop);
558                printf (" nb 1) dx = Xo * pdx \n");
559                printf ("    2) a(n) = a(n-1) / fdec \n");
560                if (modop==0)
561                { printf ("                              g = (J(Xo+a.dx) - J(Xo) - (@J/@Xo).a.dx\n");
562                  printf ("         a decrease        ordre 1: (g)-->0     ordre 2: [(g)/(a.a)]-->K \n");
563                }
564                else
565                { printf ("                              g = (J(Xo+a.dx) - J(Xo)) / (@J/@Xo).a.dx\n");
566                  printf ("         a decrease        ordre 1: (g)-->1     ordre 2: [(g-1)/a]-->K \n");
567                }
568               
569                YAL1Run=RUNL1_M1QN3; //ruse pour ne pas faire l'Yadjust inutilement dans basic_it
570                                                                                                 //de+, rappel: basic_it prevoit un razgrad_all apres la passe avant
571
572    //A) On suppose que les obs (y°) et les ebauches ebx (Xb) ont ete charges (comme Nal et d'hab)
573    //   avec loadobs, outoobs, outoebx ou une xuserfct.
574                //B) On suppose Xo initialised par l'utilisateur (setstate, loadstate ou xuserfct -> YS);
575                //   on sauvegarde Xo (i.e: YSo) dans YDo:
576    YdeltaEQPstate_target(1); //(([1->UPTIME]:)YDo_mod=YSo_mod) -> YDo
577
578                //C) algo officiel : calcul de J(Xo) et @J/@Xo.dx (une seule fois)
579                if (syclop==0)
580                {       // init de J(Xo) et (@J/@Xo).dX
581                        Ybasic_it();    //(MD+AD): J(Xo)->YTotalCost, et @J/@Xo->YGo
582                        JX = YTotalCost;
583      printf("                                 Cost function: J(Xo) = %e   J(Xo+a.dx)  (@J/@Xo).a.dx\n",JX);
584      YstateEQPOdelta_target(pdx, Ydxalea);
585
586                  Y3windice=0;
587      YstateTOtab_target(Jpx); //dX -> Jpx
588
589                  Y3windice=0;
590      YgradTOtab_target(DJ);     //@J/@Xo -> DJ
591
592      //for (int wi = 0; wi<YSIZEPB; ++wi) Jpx[wi] = Jpx[wi] * DJ[wi];  //@J/@Xo.dX
593      //attention, pour garder le meme random de dx par la suite, on intervertis Jpx et DJ
594      for (int wi = 0; wi<YSIZEPB; ++wi) DJ[wi] = Jpx[wi] * DJ[wi];     //@J/@Xo.dX
595                }
596
597                //D) boucle de test
598                for (int iloop=1; iloop<=nbloop; ++iloop)
599                {
600                                //1) evolution des p*
601                                alpha = alpha / fdec;                   //decroissance de alpha
602                                apdx  = pdx * alpha;
603
604                                //2 Xo + a.dX
605        Y3windice=0;
606        YstateEQAPTdelta_target(alpha,Jpx); //: YDo + alpha*Jpx -> YSo
607
608                                //3) iteration de base (MD+AD) -> J(Xo+a.dX)  (et @J/@(Xo+dX) inutile)
609                                Ybasic_it();    //(MD+AD): J(Xo+a.dX)->YTotalCost, et @J/@(Xo+.dX)->YGo (inutile)
610                                JXdx = YTotalCost;
611
612                                if (syclop==0) /*(nouvel algo officiel (default one))*/
613                                {        DJdx = 0.0;
614           //for (int wi = 0; wi<YSIZEPB; ++wi) DJdx += Jpx[wi] * alpha;
615           for (int wi = 0; wi<YSIZEPB; ++wi) DJdx += DJ[wi] * alpha; //nb: on a intervertis les roles de Jpx & DJ
616                                }
617                                if (syclop==1)  /*(ancien algo)*/
618                                {  printf ("testof: this old mode is no more maintened ; forget it\n"); break;
619                                   /*
620                                   //4)retauration de Xo :
621           YstateEQPdelta_target(1); //(([1->UPTIME]:)YSo_mod=YDo_mod)  ->YSo
622                                   //5) iteration de base (MD+AD) -> J(Xo)  et  @J/@Xo
623                                   Ybasic_it(); //(MD+AD): J(Xo)->YTotalCost, et @J/@Xo->YGo
624                                   JX = YTotalCost; //...
625                                   //6) calculer le denominateur (Djdx)
626                                   //   retauration de a.dx dans [Jpx]
627            YstateEQPdelta_target(apdx);  //(([1->UPTIME]:)YSo_mod=YDo_mod*ppd)->YSo    : a.dx -> YSo
628                       Y3windice=0;
629           YstateTOtab_target(Jpx); //YSo -> [Jpx]      :       a.dX -> [Jpx]
630                                   //    Mettre YGo ie @J/@Xo dans [DJ]
631                       Y3windice=0;
632           YgradTOtab_target(DJ);               //YGo -> [DJ]
633                                   //   Calcul de Djdx = @J/@Xo . a.dX
634                                   DJdx = 0.0; //@J/@Xo . p????
635                                   for (int wi = 0; wi<YSIZEPB; ++wi) DJdx += DJ[wi] * Jpx[wi];
636                                         */
637                                }
638
639                                //7) Resultat des tendences aux Ordre 1 et 2 et affichage du resultat
640                                //   1ere maniere
641                                if (modop==0)
642                                {  trendO1 = JXdx - JX - DJdx;
643                                   trendO2 = trendO1 / (alpha*alpha);
644                                   printf("   :%3i: %- 15.6e : -0-> %- 15.6e -K-> %- 15.6e %- e %- e\n", iloop, alpha, trendO1, trendO2, JXdx, DJdx);
645                                }
646                                else //   2eme maniere
647                                {  trendO1 = (JXdx-JX)/DJdx;
648                                   trendO2 = (trendO1 - 1)/alpha;
649                                   printf("   :%3i: %- 15.6e : -1-> %- 15.6e -K-> %- 15.6e %- e %- e\n", iloop, alpha, trendO1, trendO2, JXdx, DJdx);
650                                }
651                }
652
653                //E) FIN
654        //   on restaure, autant que faire ce peut, l'etat initial
655    YstateEQPdelta_target(1);
656                //   on libere la memoire
657                delete[] Jpx; delete[] DJ;
658                //   et on laise le temps propre en ordre
659                Yset_modeltime(0);
660}
661
662/* ----------------------------------------------------------------------------*/
663/* ----------------------------------------------------------------------------*/
664int Ydftestijkt(int imod)
665{   // verif que les cordonnees courantes (Yi, Yj, Yk, Yt) sont compatibles
666    // pour le module passed en parametre en se limitant aux coordonnées
667    // definies pour le module (fait pour testdf). Envoie un message d'erreur
668    // et un mauvais cr si ce n'est pas le cas
669    int codret=1;
670    if (Yi>=YTabMod[imod].axi) codret=0;
671    if (YTabMod[imod].dim>=2)
672                {        if (Yj>=YTabMod[imod].axj) codret=0;}
673    if (YTabMod[imod].dim>=3)
674                {        if (Yk>=YTabMod[imod].axk) codret=0;}
675                if (YTabMod[imod].nb_time>0)
676                {        if (YTemps>=YTabMod[imod].nb_time) codret=0;}  //ARAR
677                if (codret==0)
678                         printf("testdf: coordinates out for module %s; not proceeded\n", YTabMod[imod].Name);
679                return(codret);
680}
681/* --------------------------------------------------------------------------*/
682int Ydftesttt (int itraj)
683{   // test que le Temps courant est compatible avec une Trajectoire
684                //printf("itraj %i: %i %i %i \n", itraj, YTemps, YTabTraj[itraj].nbuptime,YTabTraj[itraj].nbsteptime);
685          if (YTemps >= YTabTraj[itraj].nbuptime + YTabTraj[itraj].nbsteptime
686          ||  YTemps < YTabTraj[itraj].nbuptime)
687                {  printf("testdf: time out for trajecyory %s; not proceeded\n", YTabTraj[itraj].name);
688                   return(0);
689                }
690                return(1);
691}
692/* --------------------------------------------------------------------------*/
693int Ytestdf(int argc, char *argv[])
694{               /*syntaxe: testdf i   j   k   t  codop [%]pdx ptol [nmmod-ko-max]
695                      [0]  [1] [2] [3] [4]  [5]    [6]   [7]      [8]            */
696
697                int   witraj, itrajelue;
698                int   modop=0;
699                int   yi, yj, yk, yt, wsyi, wsyj, wsyk, wsyt;
700                float pdx, ptol;
701                int     All, KeKo;
702                char  nmmod[LG_MAX_NAME+1];
703                int   nbko=0;
704                int       NbMaxKo=0;
705
706                //int           wt;
707
708                Ytop0 = time((time_t *)NULL);
709
710                yi=atoi(argv[1]);yj=atoi(argv[2]);yk=atoi(argv[3]);yt=atoi(argv[4]);
711
712    //if (yi<0 || yj<0 || yk<0 || yt<=0)
713    if (yi<=0 || yj<=0 || yk<=0 || yt<=0)
714    {  printf("testdf: coordinates (ijkt) must be positive \n");
715       return(0);
716    }
717
718                for (witraj=0; witraj<YNBTRAJ; ++witraj)    //sauvegarde dans YD de l'etat init
719        YdeltaEQPstate_traj(witraj, "Y#A", 0,  YTabTraj[witraj].nbuptime, 1);
720                               
721                --yi; --yj; --yk, --yt;
722                printf("\ntestdf: derivative calculus check starting : \n");
723
724                if (argv[6][0]=='%')
725                {        pdx=atof(&argv[6][1]);
726                         printf("perturbation percent  : pdx=%e \n", pdx);
727                }
728                else
729                {        modop=1; pdx=atof(argv[6]);
730                   printf("fixed perturbation    : pdx=%e \n", pdx);
731                }
732
733                ptol=atof(argv[7]);
734                printf("tolerence percent     : ptol=%e\n", ptol);
735
736                wsyi=Yi; wsyj=Yj; wsyk=Yk; wsyt=Yt; /* sauvegarde des coordonnees globales */
737                Yi=yi; Yj=yj; Yk=yk; YTemps=yt;           /* nb: dans le cas du codop = 'F' ca va etre ecrased */
738                       
739                All=1; KeKo=0; nmmod[0]='\0';
740                witraj = -1;
741                if (argc>=9)
742                {       if(Yimod(argv[8])>=0)
743                  {  All=0;
744                     strcpy(nmmod, argv[8]);
745                     witraj = Yitrajimod(Yimod(argv[8])); //trajectoire de ref = celle du module ou ...
746                  }             
747                        else
748                        {  KeKo=1;
749                           if (Yitraj(argv[8])>=0) witraj=Yitraj(argv[8]); //une trajectoire de ref it-self
750                           else //ce doit etre un maxFko > 0
751                           {  if (atoi(argv[8])<=0)
752                              {  printf("unknwon object (modul or trajectory) or bad value for MaxFko\n");
753                                 return(0);
754                              }
755                           }
756                        }
757                }
758
759                printf("\n Modulus |   grad   |  calculated  | approximeted |   |Ecart|    : test|     input    | (i j k t)\n");
760
761                if     (argv[5][0]=='f') /* test a la mode forward_order que sur une maille dont les coordonnees */
762                {       nbko = Ydfward_order_maille(modop, nmmod, All, KeKo, pdx, ptol); /* globales Yi, Yj, Yk, Yt sont deja positionnees */
763                  printf ("--> %i KO found at coordinate (%i %i %i %i)  \n", nbko, Yi+1, Yj+1, Yk+1, Yt+1);
764                }
765                else if(argv[5][0]=='F') /* test sur l'espace temps a partir des coordonnees passe en parametre...*/
766                {       /* Si le 9ième/10ième argument ne sont pas ou mal renseigned, on prend 100 ko par defaut avant de stoper */
767                        KeKo=1;
768                        if (argc==9)  NbMaxKo=atoi(argv[8]);
769                        if (argc==10) NbMaxKo=atoi(argv[9]);
770                        if (NbMaxKo==0) NbMaxKo=YNbItRun;       /* first default nbmaxko value  */
771                        if (NbMaxKo==0) NbMaxKo=100;                            /* second default nbmaxko value */
772
773                        /* une passe avant a la mode dfward */
774                        Yset_modeltime(0);
775                        before_it(1);
776                        if (witraj==-1) //trajectoire de ref par defaut
777                  {  witraj=0;
778                     printf("default trajectorie is %s \n", YTabTraj[witraj].name);
779                  }
780                  //il faut avancer jusqu'au temps voulu sans faire de test (avec trajectoire de ref)
781                        Yforward(witraj, yt);
782                        //puis continuer en faisant les tests
783                  while ((itrajelue=Yforward_elect()) >= 0)   // le while controle la progression
784                  {   YTemps  = YTabTraj[itrajelue].toptime;
785                                  YidTraj = itrajelue;
786          if (YDispTime)
787          {   printf (" >>> Traj %s, Current dfward time = %i \n", YTabTraj[itrajelue].name, YTemps+1);
788                                                 //Ydispcurstep(itrajelue);
789                                  }
790                                  forward_before(0); //arar
791                                  nbko += YTabTraj[itrajelue].dfward(modop, nmmod, All, KeKo, NbMaxKo-nbko, pdx, ptol, yi, yj, yk); //<-: dfward_order (pour LA trajectoire selectionnee)
792                                  forward_after(0); //arar
793          {++YTabTraj[itrajelue].toptime; //<-: ++YTemps
794           YTabTraj[itrajelue].curtime += YTabTraj[itrajelue].dtime; //<-: multi
795          }
796                                        if (nbko >= NbMaxKo) break;
797                  }
798                  /* retablissement coordonnees au bord si besoin pour affichage, ci-apres, en mode externe*/
799                        /*
800                        if (Yi>=YA1) Yi=YA1-1; if (Yi<0) Yi=0;
801                        #ifdef YA2
802                        if (Yj>=YA2) Yj=YA2-1; if (Yj<0) Yj=0;
803                        #endif
804                        #ifdef YA3
805                        if (Yk>=YA3) Yk=YA3-1; if (Yk<0) Yk=0;
806                        #endif
807                        if (Yt>=YNBALLTIME) Yt=YNBALLTIME-1;
808                  printf ("--> %i KO found [from] to reached coordinate (%i %i %i %i)  \n", nbko, Yi+1, Yj+1, Yk+1, Yt+1);
809                  */
810      /*ARAR*///printf ("--> %i KO found [from] to reached coordinate (%i %i %i %i)  \n", nbko, Yi+1, Yj+1, Yk+1, YTemps+1);
811      printf ("--> %i KO found to reached coordinate\n", nbko); /*ARAR*/
812                }
813                else
814                {       Ysetting(argv[5]);
815                  nbko = Ydfward_all(modop, nmmod, All, KeKo, pdx, ptol);
816                  printf ("--> %i KO found at coordinate (%i %i %i %i)  \n", nbko, Yi+1, Yj+1, Yk+1, Yt+1);
817                }
818               
819                Yi=wsyi; Yj=wsyj; Yk=wsyk; YTemps=wsyt; /* on laisse les globales comme on les a trouved */
820                printf("\n");
821
822                for (witraj=0; witraj<YNBTRAJ; ++witraj)  //restaur l'etat init autant que possible
823        YstateEQPdelta_traj(witraj, "Y#A", 0, YTabTraj[witraj].nbuptime, 1);                           
824
825                return(0);
826}
827/* ----------------------------------------------------------------------------*/
828/* ----------------------------------------------------------------------------*/
829
830/**OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO**/
831/**    OLD  VERSION      OLD  VERSION      OLD  VERSION      OLD  VERSION      OLD  VERSION      OLD  VERSION      OLD  VERSION **/
832/**OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO**/
833/*
834void Ytestad_mod(double pdx, double errelmax, int maxko, double pzedi)
835{                // PURPOSE : test de l'ADjoint par module : <m'(Xo)x,y> = <x, m*(Xo).y>
836                 // pdx :       Pourcentage De perturbation de Xo pour obtenir dx
837                 // errelmax : Erreur relative max : on test tous les modules en chaque point
838                 // et maxko : de coordonnee, et on considere que le test est ko si l'erreur
839     //            relative est > a errelmax. Si on trouve plus de maxko tests ko,
840     //            le process est interrompu (exit).
841                 // pzedi    : Parametres pour se donner un niveau de zero informatique
842                 //                      lorsque les 2 |termes| sont < pzedi alors -> pas de test => OK !?
843     //
844     // nb: maintenir en meme temps la fonction testad (s'il y a lieu)
845     //
846     // called fct : YgradEQPstate_all, YdeltaEQPstate_all,YgradTOtab_all, YtabTOgrad_all
847     //              before_it, Yforward, Ylinward, Yrazgrad_all, Ybackward
848                 //
849
850    int   witraj;
851    int   NBALLTIME;
852                YREAL *MXdx = new YREAL[YSIZECO];               // tableau de sauvegarde
853                if (MXdx==NULL){printf("Yao=>ABORT: pb on new MXdx for testad (see Yao system administrator)\n"); exit(-9);}
854          memset(MXdx, 0, YSIZECO*sizeof(YREAL));
855
856                Ytop0 = time((time_t *)NULL);
857
858                // 0) Presentation
859                printf ("\n module Adjoint TEST :        dx = Xo * pdx  (with pdx = %e)", pdx);
860                printf ("\n --------------------    and  dy = m'(Xo).dx\n\n");
861
862                // nb: On suppose Xo initialised par l'utilisateur (setstate ou xuserfct -> YSo)
863
864                // 1) some init
865                Ytestad_module=1;
866                Yerrelmax=errelmax;
867                Ynbko=0; Ymaxko=maxko;
868                Ypzedi = pzedi;
869
870                for (witraj=0; witraj<YNBTRAJ; ++witraj)
871                {               NBALLTIME = YTabTraj[witraj].nbsteptime + YTabTraj[witraj].nbuptime;
872        YgradEQPstate_traj (witraj, "Y#A", 0, NBALLTIME, pdx);  //il faut initialiser tous les dx et en
873        YdeltaEQPstate_traj(witraj, "Y#A", 0, NBALLTIME, pdx);  //avoir l'equivalent dans les delta
874    } //Avant (paradigme V7) on devait pouvoir se suffir d'initialiser jusqu'a uptime.
875      //Maintenant (paradigme V8) version multi-spatio-tempo + target sur trajectoire,
876      //il vaut mieux tout initialiser (i.e. toute la trajectoire) sachant que le LT
877      //ecrasera par les nelles valeurs qu'il aura calculees
878      //(c'est d'ailleurs aussi la raison pour laquelle un razgrad_all ne devrait pas etre utile !)
879                                                               
880                // 2) MD (Modele Direct): Calculer M(Xo)
881                Yset_modeltime(0);
882                before_it(1);                                //(YItRun);
883                Yforward (-1, 0);                    //MD : M(Xo) -> YS
884    Yact_operator('*', ON);
885    Yforward_operator('H');                                                  //forward operateurs a la fin; d'abord H seulement,
886    for (witraj=0; witraj<YNBTRAJ; ++witraj)                                                                                                                       //puis il faut mettre YW = YS pour les modules de
887    {   NBALLTIME = YTabTraj[witraj].nbsteptime + YTabTraj[witraj].nbuptime; //(M ou H) connectes vers des modules operateurs
888                    YwishEQPstate_traj_tocov(witraj, NBALLTIME-1, 1.0);                                                                  //de covariance (R ou K) ... puisque ... sinon ...
889    }
890    Yforward_operator('R');             //et on termine par les
891    Yforward_operator('K');     //    operateurs R et K
892
893                // 3) LT (Lineaire Tangent): m'(Xo)dXo
894                Yset_modeltime(0);
895                before_it(1);                                //(YItRun);
896                Ylinward(0);                                 // TL : M'(Xo)dXo -> YG
897    Ylinward_operator('*');  //+linward des operateurs a la fin
898
899                // 4) YGn-> tableau Mxdx: avant RAZ, sauvegarde des dx obtenus sur les modules cout en fin de traj
900    Y3windice=0;
901    for (witraj=0; witraj<YNBTRAJ; ++witraj)
902    {   NBALLTIME = YTabTraj[witraj].nbsteptime + YTabTraj[witraj].nbuptime;
903        YgradTOtab_traj(witraj, "Y#C", NBALLTIME-1, NBALLTIME , MXdx);                  // YGn -> [MXdx]
904    }
905
906                // 5) RAZ tous les gradients pour l'Adjoint
907                Yrazgrad_all();         // YG = 0
908
909                // 6) et Reprise des dx de fin de trajectoire des modules cout (qui vont etre retro-propagees)
910    Y3windice=0;
911    for (witraj=0; witraj<YNBTRAJ; ++witraj)
912    {   NBALLTIME = YTabTraj[witraj].nbsteptime + YTabTraj[witraj].nbuptime;
913        YtabTOgrad_traj(witraj, "Y#C", NBALLTIME-1, NBALLTIME , MXdx);                          // [MXdx] -> YGn
914    }
915
916                // 7) on se lance dans l'Adjoint
917    Ybackward_operator('*'); //d'abord backward des operateurs en fin de trajectoire
918                Ybackward(-1, 0);        //au cours duquel le test sur chaque module se fait (YNBSTEPTIME);
919                //quid after_it ?
920
921                printf(" -->-->--> %i case(s) found \n\n", Ynbko);
922                // 8) FIN
923                //    on libere la memoire
924                delete[] MXdx;
925                //    et on laise le temps propre en ordre
926                Yset_modeltime(0);
927
928                // Restaurer (~) l'etat init ? En principe il n'est pas modified !?
929                //for (witraj=0; witraj<YNBTRAJ; ++witraj)
930    //    YstateEQPdelta_traj(witraj, "Y#A", 0, YTabTraj[witraj].nbuptime, 1/pdx);
931                Ytestad_module=0;
932}
933*/
934/**OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO**/
935/*
936void Ytestad(double pdx)
937{                // PURPOSE : test de l'ADjoint : <M'(Xo)x,y> = <x, M*(Xo).y>
938                 //             ici, dx correspond a x, et dx* a y
939                 //             dans Yao la zone utilised pour Xo est YS, et YG est utilised pour dx lors
940                 //             du TL (passe avant) puis pour dx* lors de l'Adjoint (passe arriere) d'ou
941                 //             la necessite de passer par des tableaux pour sauvegarder et manipuler les
942                 //             resultats des calculs.
943                 //  pdx :      Pourcentage De perturbation de Xo pour obtenir dx
944                 //  errelmax : Erreur relative max : si cette valeur est differente de 0 alors
945                 //  et maxko : on test tous les modules en chaque point de coordonnee, et on
946                 //             considere que le test est ko si l'erreur relative est > a errelmax.
947                 //             Si on trouve plus de maxko tests ko, le process est interrompu (exit).
948     //
949                 //     nb: maintenir en meme temps la fonction testad_mod (s'il y a lieu)
950     //
951                 //     called fct : YgradEQPstate_all, YgradTOtab_all, YtabTOgrad_all
952                 //                                                Yrazgrad_all, before_it, Yforward, Ylinward, Ybackward
953                 //
954
955                int   witraj;
956                int   NBALLTIME;
957                double LTRes=0.0, AdRes=0.0;
958                YREAL *DX   = new YREAL[YSIZEPB];               // tableau de sauvegarde de dx
959                if (DX==NULL)  {printf("Yao=>ABORT: pb on new DX for testad (see Yao system administrator)\n"); exit(-9);}
960                YREAL *MXdx = new YREAL[YSIZECO];               // tableau resultat du Lineaire Tangent M'(Xo)dX
961                if (MXdx==NULL){printf("Yao=>ABORT: pb on new MXdx for testad (see Yao system administrator)\n"); exit(-9);}
962                YREAL *MAdx = new YREAL[YSIZEPB];               // tableau resultat de l'Adjoint M*(Xo)d*X
963                if (MAdx==NULL){printf("Yao=>ABORT: pb on new MAdx for testad (see Yao system administrator)\n"); exit(-9);}
964
965          memset(DX,   0, YSIZEPB*sizeof(YREAL));
966          memset(MXdx, 0, YSIZECO*sizeof(YREAL));
967          memset(MAdx, 0, YSIZEPB*sizeof(YREAL));
968
969                Ytop0 = time((time_t *)NULL);
970
971                // 0) Presentation
972                printf ("\n Adjoint TEST :         dx = Xo * pdx  (with pdx = %e)", pdx);
973                printf ("\n --------------    and  dy = M'(Xo).dx\n\n");
974
975                Ytestad_module=0;               // because c'est la meme instruction qui sert dans les 2 cas, il faut pouvoir distinguer
976
977                // nb: On suppose Xo initialised par l'utilisateur (setstate ou xuserfct -> YSo)
978
979                // 1) Init de dXo -> YGo
980                Yrazgrad_all();
981
982    // dXo = Xo * pdx (([1->UPTIME]:)YGo_mod=YSo_mod*pdx)
983    YgradEQPstate_target(pdx);
984
985                // 2) MD (Modele Direct): Calculer M(Xo)
986                Yset_modeltime(0);
987                before_it(1);                           //(YitRun);
988                Yforward(-1, 0);                //MD : M(Xo) -> YS
989    Yact_operator('*', ON);
990    Yforward_operator('H');                                                  //forward operateurs a la fin; d'abord H seulement,
991    for (witraj=0; witraj<YNBTRAJ; ++witraj)                                                                                                                       //puis il faut mettre YW = YS pour les modules de
992    {   NBALLTIME = YTabTraj[witraj].nbsteptime + YTabTraj[witraj].nbuptime; //(M ou H) connectes vers des modules operateurs
993                    YwishEQPstate_traj_tocov(witraj, NBALLTIME-1, 1.0);                                                                  //de covariance (R ou K) ... puisque ... sinon ...
994    }
995    Yforward_operator('R');  //et on termine par les
996    Yforward_operator('K');  //    operateurs R et K
997
998                // 3) LT (Lineaire Tangent): M'(Xo)dXo
999                Yset_modeltime(0);
1000                before_it(1);                           //(YitRun);
1001                Ylinward(0);                            //TL : M'(Xo)dXo -> YG
1002    Ylinward_operator('*'); //linward opera la fin !!?? (car on ne s'interesse qu'aux modules cout en fin de traj)
1003
1004                // 4) YGn-> tableau Mxdx et calcul de < M'(Xo)dX, dy > (sur module cout en fin de traj)
1005                Y3windice=0;
1006    for (witraj=0; witraj<YNBTRAJ; ++witraj)
1007    {   NBALLTIME = YTabTraj[witraj].nbsteptime + YTabTraj[witraj].nbuptime;
1008        YgradTOtab_traj(witraj, "Y#C", NBALLTIME-1, NBALLTIME , MXdx);                                  // YGn -> [MXdx]
1009    }
1010
1011                for(int i=0; i<YSIZECO; ++i)
1012                {LTRes += MXdx[i]*MXdx[i];      // < M'(Xo)dX, dy >
1013                }
1014
1015                // 6) au passage, on a peut etre aussi besoin de conserver donc de sauvegarder YGo (isnt-it?)
1016                Y3windice=0;
1017    YgradTOtab_target(DX);                      // YGo -> [DX]
1018
1019                // 7) RAZ tous les gradients pour l'Adjoint
1020                Yrazgrad_all();         // YG = 0
1021
1022                // 8) et Reprise des gradients de fin de trajectoire des modules cout
1023                Y3windice=0;
1024    for (witraj=0; witraj<YNBTRAJ; ++witraj)
1025    {   NBALLTIME = YTabTraj[witraj].nbsteptime + YTabTraj[witraj].nbuptime;
1026        YtabTOgrad_traj(witraj, "Y#C", NBALLTIME-1, NBALLTIME , MXdx);                          // [MXdx] -> YGn
1027    }
1028
1029                // 9) on devrait pouvoir se lancer dans l'adjoint
1030    Ybackward_operator('*'); //d'abord backward des operateurs en fin de trajectoire
1031                Ybackward(-1, 0);  // AD (adjoint):-> dx =M*(Xo)d*X : Yjac * YG -> Ytbeta
1032                //quid after_it ?
1033
1034                // 10) YGo-> tableau MAdx et calcul de < dx, M*(Xo).dy >
1035                Y3windice=0;
1036    YgradTOtab_target(MAdx);             // YGo -> [MAdx]
1037
1038                for(int i=0; i<YSIZEPB; ++i)
1039                {       AdRes += DX[i] * MAdx[i];                                                                       // < dx, M*(Xo).dy >
1040                }
1041
1042                // 11) Afficher le resultat
1043                printf (" < M'(Xo).dx, dy >  = %e \n", LTRes);
1044                printf (" < dx, M*(Xo).dy >  = %e \n", AdRes);
1045                printf ("                     --------------\n");
1046                printf ("   ecart :            %e \t", LTRes-AdRes);
1047                printf (" relative error = %e \n\n", (LTRes-AdRes)/LTRes);
1048
1049                // 12) on informe d'une possible cause d'echec par exemple si des obs sont
1050                if (YioInsertObsCtr != -1) // charger dans une des fonction d'intervention
1051                {        printf("testad warning: obs must not be used for global Adjoint test\n");
1052                           //because: ca modifie le dy au cours de la redescente qui ne sera plus le meme
1053                           //         que celui utilised lors du calcul de LTRes.
1054                           //         Cette contrainte n'est pas utile pour testad_mod car le test est local.
1055                }
1056
1057                // 13) FIN
1058                //    on libere la memoire
1059                delete[] DX;  delete[] MXdx;  delete[] MAdx;
1060                //    et on laise le temps propre en ordre
1061                Yset_modeltime(0);
1062}
1063*/
1064/**OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO**/
1065/*
1066void Ytestof(double pdx, double alpha, double fdec, int nbloop, int modop, int syclop)
1067{               //PURPOSE : test de l'Objective Fonction :
1068                //pdx   : pourcentage De perturbation de Xo pour obtenir dx
1069                //alpha : parametre de tendance
1070                //fdec  : facteur de decroissance de pdec alias alpha
1071                //nbloop: nombre de boucles de decroissance
1072                //modop : mode operatoire: si 0: ->0, sinon ->1
1073                //syclop: Secret Yao Systeme Code Operation (0: algo officiel par defaut)
1074                //        1: algo (initial) theoriquement equivalent a l'algo officiel mais plus lent
1075                               
1076                double  apdx;
1077                double  JXdx;
1078                double  trendO1, trendO2;
1079                double  DJdx=0.0;
1080                double  JX=0.0;
1081                YREAL           *Jpx = new YREAL[YSIZEPB];              // tableau pour J'(Xo)dx (avant c'etait pour les a.dx)
1082                if (Jpx==NULL){printf("Yao=>ABORT: pb on new Jpx for testof (see Yao system administrator)\n"); exit(-9);}
1083                YREAL           *DJ  = new YREAL[YSIZEPB];              // tableau de gradinat YG
1084                if (DJ==NULL) {printf("Yao=>ABORT: pb on new DJ for testof (see Yao system administrator)\n"); exit(-9);}
1085
1086          memset(Jpx, 0, YSIZEPB*sizeof(YREAL));
1087          memset(DJ,  0, YSIZEPB*sizeof(YREAL));
1088
1089                Ytop0 = time((time_t *)NULL);
1090
1091                //   0) presentation :
1092                printf ("\n Objective Function TEST : with: pdx = %e, (a)lpha = %e", pdx, alpha);
1093                printf ("\n -------------------------       fdec= %e, (n)bloop= %i\n", fdec, nbloop);
1094                printf (" nb 1) dx = Xo * pdx \n");
1095                printf ("    2) a(n) = a(n-1) / fdec \n");
1096                if (modop==0)
1097                { printf ("                             g = J(Xo+a.dx) - J(Xo) - (@J/@Xo).a.dx\n");
1098                  printf ("         a decrease        ordre 1: (g)-->0    ordre 2: [(g)/(a.a)]-->K \n\n");
1099                }
1100                else
1101                { printf ("                             g = (J(Xo+a.dx) - J(Xo)) / (@J/@Xo).a.dx\n");
1102                  printf ("         a decrease        ordre 1: (g)-->1    ordre 2: [(g-1)/a]-->K \n\n");
1103                }
1104
1105                YAL1Run=RUNL1_M1QN3; //ruse pour ne pas faire l'Yadjust inutilement dans basic_it
1106                                                                                                 //de+, rappel: basic_it prevoit un razgrad_all apres la passe avant
1107
1108    //A) On suppose que les obs (y°) et les ebauches ebx (Xb) ont ete charges (comme Nal et d'hab)
1109    //   avec loadobs, outoobs, outoebx ou une xuserfct.
1110                //B) On suppose Xo initialised par l'utilisateur (setstate, loadstate ou xuserfct -> YS);
1111                //   on sauvegarde Xo (i.e: YSo) dans YDo:
1112    YdeltaEQPstate_target(1); //(([1->UPTIME]:)YDo_mod=YSo_mod) -> YDo
1113
1114                //C) algo officiel : calcul de J(Xo) et @J/@Xo.dx (une seule fois)
1115                if (syclop==0)
1116                {       // init de J(Xo) et (@J/@Xo).dX
1117                        Ybasic_it();    //(MD+AD): J(Xo)->YTotalCost, et @J/@Xo->YGo
1118                        JX = YTotalCost;
1119      printf("                                 Cost Function: J(Xo)=%e     J(Xo+a.dx)  (@J/@Xo).a.dx\n",JX);
1120      YstateEQPdelta_target(pdx); //dX
1121
1122                  Y3windice=0;
1123      YstateTOtab_target(Jpx); //dX -> Jpx
1124
1125                  Y3windice=0;
1126      YgradTOtab_target(DJ);                    //@J/@Xo -> DJ
1127
1128                        for (int wi = 0; wi<YSIZEPB; ++wi) Jpx[wi] = Jpx[wi] * DJ[wi];  //@J/@Xo.dX
1129                }
1130
1131                //D) boucle de test
1132                for (int iloop=1; iloop<=nbloop; ++iloop)
1133                {
1134                                //1) evolution des p*
1135                                alpha = alpha / fdec;                   //decroissance de alpha
1136                                apdx  = pdx * alpha;
1137
1138                                //2 Xo + a.dX
1139        YstateEQPdelta_target(1+apdx); //: YDo * (1+apdx) -> YSo
1140
1141                                //3) iteration de base (MD+AD) -> J(Xo+a.dX)  (et @J/@(Xo+dX) inutile)
1142                                Ybasic_it();    //(MD+AD): J(Xo+a.dX)->YTotalCost, et @J/@(Xo+.dX)->YGo (inutile)
1143                                JXdx = YTotalCost;
1144
1145                                if (syclop==0) //(nouvel algo officiel (default one))
1146                                {        DJdx = 0.0;
1147                                         for (int wi = 0; wi<YSIZEPB; ++wi) DJdx += Jpx[wi] * alpha;
1148                                }
1149                                if (syclop==1)  //(ancien algo)
1150                                {        //4)retauration de Xo :
1151           YstateEQPdelta_target(1); //(([1->UPTIME]:)YSo_mod=YDo_mod)  ->YSo
1152
1153                                   //5) iteration de base (MD+AD) -> J(Xo)  et  @J/@Xo
1154                                   Ybasic_it(); //(MD+AD): J(Xo)->YTotalCost, et @J/@Xo->YGo
1155                                   JX = YTotalCost; //...
1156                                   //6) calculer le denominateur (Djdx)
1157                                   //   retauration de a.dx dans [Jpx]
1158            YstateEQPdelta_target(apdx);  //(([1->UPTIME]:)YSo_mod=YDo_mod*ppd)->YSo    : a.dx -> YSo
1159
1160                       Y3windice=0;
1161           YstateTOtab_target(Jpx); //YSo -> [Jpx]      :       a.dX -> [Jpx]
1162
1163                                   //    Mettre YGo ie @J/@Xo dans [DJ]
1164                       Y3windice=0;
1165           YgradTOtab_target(DJ);               //YGo -> [DJ]
1166
1167                                   //   Calcul de Djdx = @J/@Xo . a.dX
1168                                   DJdx = 0.0; //@J/@Xo . p????
1169                                   for (int wi = 0; wi<YSIZEPB; ++wi) DJdx += DJ[wi] * Jpx[wi];
1170                                }
1171
1172                                //7) Resultat des tendences aux Ordre 1 et 2 et affichage du resultat
1173                                //   1ere maniere
1174                                if (modop==0)
1175                                {  trendO1 = JXdx - JX - DJdx;
1176                                   trendO2 = trendO1 / (alpha*alpha);
1177                                   //printf("   :%3i: %- 15.6e : -0-> %- 15.6e    -K-> %- 15.6e\n", iloop, alpha, trendO1, trendO2);
1178                                   printf("   :%3i: %- 15.6e : -0-> %- 15.6e -K-> %- 15.6e %- e %- e\n", iloop, alpha, trendO1, trendO2, JXdx, DJdx);
1179                                }
1180                                else //   2eme maniere
1181                                {  trendO1 = (JXdx-JX)/DJdx;
1182                                   trendO2 = (trendO1 - 1)/alpha;
1183                                   //printf("   :%3i: %- 15.6e : -1-> %- 15.6e    -K-> %- 15.6e\n", iloop, alpha, trendO1, trendO2);
1184                                   printf("   :%3i: %- 15.6e : -1-> %- 15.6e -K-> %- 15.6e %- e %- e\n", iloop, alpha, trendO1, trendO2, JXdx, DJdx);
1185                                }
1186                }
1187
1188                //E) FIN
1189        //   on restaure, autant que faire ce peut, l'etat initial
1190    YstateEQPdelta_target(1);
1191                //   on libere la memoire
1192                delete[] Jpx; delete[] DJ;
1193                //   et on laise le temps propre en ordre
1194                Yset_modeltime(0);
1195}
1196*/
1197/**OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO**/
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
Note: See TracBrowser for help on using the repository browser.