[1] | 1 | /* ----------------------------------------------------------------------------*/ |
---|
| 2 | /* ----------------------------------------------------------------------------*/ |
---|
| 3 | //char Ydxalea[STRSIZE20+1] = "R0.0000001"; pour un dx (perturbation de Xo) aleatoire ... |
---|
| 4 | char Ydxalea[STRSIZE20+1] = "\0"; // mais initialement (par defaut) , cet alea est neutralised |
---|
| 5 | /* ----------------------------------------------------------------------------*/ |
---|
| 6 | /* ----------------------------------------------------------------------------*/ |
---|
| 7 | void 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 */ |
---|
| 31 | printf("\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 | /* ----------------------------------------------------------------------------*/ |
---|
| 103 | void 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 */ |
---|
| 139 | printf("\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 | /* ----------------------------------------------------------------------------*/ |
---|
| 235 | void 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 | /* ----------------------------------------------------------------------------*/ |
---|
| 529 | void 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 : |
---|
| 555 | printf("\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 | /* ----------------------------------------------------------------------------*/ |
---|
| 664 | int 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 | /* --------------------------------------------------------------------------*/ |
---|
| 682 | int 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 | /* --------------------------------------------------------------------------*/ |
---|
| 693 | int 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 | /* |
---|
| 834 | void 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 | /* |
---|
| 936 | void 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 | /* |
---|
| 1066 | void 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 | |
---|