1 : * AJUSTE PROCEDUR GOUNAND 12/07/27 21:15:00 7444 2 : * 3 : * =============================================================== 4 : * AJUSTEMENT D'UNE FONCTION DE LA FORME 5 : * 6 : * F = q1*f1(x,p) + q2*f2(x,p) + ... + qL*fL(x,p) + g(x,p) 7 : * 8 : * on a L parametres lineaires : q 9 : * et K parametres non lineaires: p 10 : * 11 : * definir les procedures FCT et DERI 12 : * (dont on peut préciser des noms différents) 13 : * 14 : * nombre maximum de points: mxn=100000 15 : * de parametres lineaires: mxl=5000 16 : * de parametres non lineaires: mxk=50 17 : * 18 : * =============================================================== 19 : * ENTREES 20 : * tab1 : TABLE d'indices : 21 : * 'K' : nb de paramètres internes 'non linéaires' aux fonctions 22 : * 'L' : nb de fonctions précédées de paramètres 23 : * 'linéaires' à évaluer 24 : * 'X' : TABLE des coordonnées (x,y,z,...) 25 : * 'F' : LISTREEL des N valeurs à caler 26 : * 'MESSAGES' : niveau de bavardage (défaut 0) 27 : * 'IMPRESSION' : fréquence avec laquelle impression les résultats 28 : * 'NOM_FCT' : nom de la proc. de calcul de la fonction 29 : * 'NOM_DERI' : nom de la proc. de calcul de sa dérivée par 30 : * rapport aux paramètres 31 : * 32 : * 'PMIN' : LISTREEL des K val. minimales des K paramètres 33 : * internes 'non linéaires' 34 : * 'PMAX' : LISTREEL des K val. maximales " " " 35 : * 'PRECISION' : LISTREEL des K val. de précision souhaitée pour eux 36 : * Ces 3 indices doivent être donnés si K > 0 37 : * (défaut 1d-7) 38 : * 'POIDS' : LISTREEL des N poids à attribuer à chaque point de valeur 39 : * 'MXTER' : nb max. d'itérations avant sortie (si K > 1) 40 : * 41 : * SORTIES 42 : * q : LISTREEL des valeurs des paramètres internes optimaux 43 : * p : idem (non linéaires) 44 : * 45 : * =============================================================== 46 : * 47 : * MODIFICATIONS : 48 : * 20/04/2006, P. Maugis : 49 : * mise en paramètre des noms de procédures FCT et DERI 50 : * 10/10/2006, P. Maugis : 51 : * niveau de messages 52 : * 53 : * =============================================================== 54 : 55 : 'DEBPROC' AJUSTE tab1*'TABLE' ; 56 : 57 : * Lecture des entrées 58 : * ------------------- 59 : k = tab1.'K'; 60 : l = tab1.'L'; 61 : xtab = tab1.'X'; 62 : y = tab1.'F'; 63 : 64 : 'SI' ('NON' ('EXISTE' tab1 'NOM_FCT')) ; 65 : * on conserve le nom par défaut 66 : NOMFCT = 'MOT' 'FCT' ; 67 : 'SINON' ; 68 : NOMFCT = 'MOT' ('TEXTE' ('CHAINE' tab1.'NOM_FCT')) ; 69 : 'FINSI' ; 70 : 71 : 'SI' ('NON' ('EXISTE' tab1 'NOM_DERI')) ; 72 : * on conserve le nom par défaut 73 : NOMDERI = 'MOT' 'DERI' ; 74 : 'SINON' ; 75 : NOMDERI = 'MOT' ('TEXTE' ('CHAINE' tab1.'NOM_DERI')) ; 76 : 'FINSI' ; 77 : 78 : 'SI' ('EXISTE' tab1 'IMPRESSION'); 79 : iimpr = tab1.'IMPRESSION'; 80 : 'SINON' ; 81 : iimpr = 20; 82 : 'FINSI' ; 83 : 84 : * Coefficients de pondération 85 : 'SI' ('EXISTE' TAB1 'POIDS'); 86 : POI = TAB1.'POIDS'; 87 : 'SINON'; 88 : JG = 'DIME' Y; 89 : POI = 'PROG' JG*1.d0 ; 90 : 'FINSI'; 91 : 92 : * Niveau de message 93 : 'SI' ('EXISTE' tab1 'MESSAGES') ; 94 : IMESS = 'ENTIER' (tab1.'MESSAGES') ; 95 : 'SINON' ; 96 : IMESS = 0 ; 97 : 'FINSI' ; 98 : 99 : * TROIS cas de figure selon le nb de paramètres non linéaires 100 : * ----------------------------------------------------------- 101 : 102 : * 1) Cas où PAS DE paramètres non linéaires 103 : 'SI' (k '<EG' 0); 104 : tbfonc = ('TEXTE' NOMFCT) xtab ; 105 : 106 : phi q = 'AJU1' l tbfonc y POI; 107 : 'SI' (IMESS '>EG' 1) ; 108 : 'MESS' 'AJUSTE:' ; 109 : 'MESS' ' Les ' (@ARR l 0) 110 : ' paramètres linéaires valent respectivement'; 111 : 'REPETER' bcl1 l ; 112 : mess ' ' ('EXTR' q &bcl1) ; 113 : 'FIN' bcl1 ; 114 : 'MESS' ' Il n y a pas de paramètre non linéaire'; 115 : 'MESS' ' Erreur d estimation : ' wphi ; 116 : 'FINSI' ; 117 : 118 : 'FINSI'; 119 : 120 : * 2) Cas avec UN SEUL paramètre non linéaire 121 : 'SI' (k 'EGA' 1); 122 : 123 : * on lit les options pour les paramètres non linéaires 124 : pmin = tab1.'PMIN'; 125 : pmax = tab1.'PMAX'; 126 : 'SI' ( 'EXISTE' tab1 'PRECISION'); 127 : dp = tab1.'PRECISION'; 128 : 'SINON' ; 129 : dp = 'PROG' K * 1D-7 ; 130 : 'FINSI'; 131 : 132 : p1 = 0.6180340051651; 133 : q1 = 0.3819659948349; 134 : ZERO = 2.E-56; 135 : rela = 2.E-14; 136 : 137 : wa = 'EXTR' pmin 1; 138 : wb = 'EXTR' pmax 1; 139 : wdp = 'EXTR' dp 1; 140 : 141 : prec = 'MAXI' ('PROG' 0 wdp); 142 : eaq = 'MAXI' ('PROG' zero 143 : (rela*( ('ABS' wb)+('ABS' (wa+prec)))) ); 144 : raq = wb - wa - prec; 145 : 'SI' (('ABS' raq) '<EG' eaq); 146 : raq = 0.; 147 : 'FINSI'; 148 : 'SI' (raq '<EG' 0); 149 : pp = 0.5 * (wa+wb); 150 : 'SINON'; 151 : wu = (p1*wa) + (q1*wb); 152 : wv = (q1*wa) + (p1*wb); 153 : 154 : wul = 'PROG' wu; 155 : wvl = 'PROG' wv; 156 : 157 : wtbfoncu = ('TEXTE' NOMFCT) xtab wul; 158 : wtbfoncv = ('TEXTE' NOMFCT) xtab wvl; 159 : 160 : wf wq = 'AJU1' l wtbfoncu y POI; 161 : wg wq = 'AJU1' l wtbfoncv y POI; 162 : 163 : 'REPETER' wbloc ; 164 : eaq = 'MAXI' ('PROG' zero (rela*( ('ABS' wf)+('ABS' wg) ) ) ); 165 : raq = wf - wg; 166 : 'SI' (('ABS' raq) '<EG' eaq); 167 : raq = 0.; 168 : 'FINSI'; 169 : 'SI' (raq '<EG' 0.); 170 : wb = wv; 171 : eaq = 'MAXI' ('PROG' zero 172 : (rela*( ('ABS' wb)+('ABS' (wa+prec)))) ); 173 : raq = wb - wa - prec; 174 : 'SI' (('ABS' raq) '<EG' eaq); 175 : raq = 0. ; 176 : 'FINSI'; 177 : 'SI' (raq '<EG' 0); 178 : wpp = 0.5 * (wa+wb); 179 : 'QUITTER' wbloc; 180 : 'SINON'; 181 : wv = wu; 182 : wg = wf; 183 : wu = (p1*wa) + (q1*wb); 184 : wul = 'PROG' wu; 185 : wtbfoncu = ('TEXTE' NOMFCT) xtab wul; 186 : wf wq = 'AJU1' l wtbfoncu y POI; 187 : 'ITERER' wbloc; 188 : 'FINSI'; 189 : 190 : 'SINON'; 191 : wa = wu; 192 : eaq = 'MAXI' ('PROG' zero 193 : (rela*( ('ABS' wb)+('ABS' (wa+prec)))) ); 194 : raq = wb - wa - prec; 195 : 'SI' (('ABS' raq) '<EG' eaq); 196 : raq = 0. ; 197 : 'FINSI'; 198 : 'SI' (raq '<EG' 0.); 199 : wpp = 0.5 * (wa+wb); 200 : 'QUITTER' wbloc; 201 : 'SINON'; 202 : wu = wv; 203 : wf = wg; 204 : wv = (q1*wa) + (p1*wb); 205 : wvl = 'PROG' wv; 206 : wtbfoncv = ('TEXTE' NOMFCT) xtab wvl; 207 : wg wq ='AJU1' l wtbfoncv y POI; 208 : 'ITERER' wbloc; 209 : 'FINSI'; 210 : 'FINSI'; 211 : 'FIN' wbloc; 212 : 'FINSI'; 213 : p = 'PROG' wpp; 214 : wtbfp = ('TEXTE' NOMFCT) xtab p; 215 : wphi q= 'AJU1' l wtbfp y POI; 216 : 217 : 'SI' (IMESS '>EG' 1) ; 218 : 'MESS' 'AJUSTE:' ; 219 : 'MESS' ' Les ' (@ARR l 0) 220 : ' paramètres linéaires valent respectivement'; 221 : 'REPETER' bcl1 l ; 222 : mess ' ' ('EXTR' q &bcl1) ; 223 : 'FIN' bcl1 ; 224 : 'MESS' ' Le paramètre non linéaire vaut' ('EXTR' p 1) ; 225 : 'MESS' ' Erreur d estimation : ' wphi ; 226 : 'FINSI' ; 227 : 228 : 'FINSI'; 229 : 230 : * 3) Cas avec PLUSIEURS paramètres non linéaire 231 : 'SI' (k > 1); 232 : 233 : * on lit les options pour les paramètres non linéaires 234 : pmin = tab1.'PMIN'; 235 : pmax = tab1.'PMAX'; 236 : 'SI' ( 'EXISTE' tab1 'PRECISION'); 237 : dp = tab1.'PRECISION'; 238 : 'SINON' ; 239 : dp = 'PROG' k * 1d-5 ; 240 : 'FINSI'; 241 : 242 : p = pmin + ( (pmax - pmin) * 0.555 ); 243 : 'SI' ( 'EXISTE' tab1 'MXTER'); 244 : mxter = tab1.'MXTER'; 245 : 'SINON'; 246 : mxter = 100; 247 : 'FINSI'; 248 : 249 : hzinf = 1.D150; 250 : qzero = 2.D-56; 251 : qrela = 2.D-14; 252 : hk = 'DIME' p; 253 : hiter = 0; 254 : 255 : htbfonc = ('TEXTE' NOMFCT) xtab p; 256 : htbderi = ('TEXTE' NOMDERI) xtab p; 257 : 258 : hphi1 hq = 'AJU1' l htbfonc y POI; 259 : 260 : hu hf = 'AJU2' hk hq htbfonc htbderi y POI; 261 : 262 : 'REPETER' hbloc2 ; 263 : h = hf; 264 : hkas = 0; 265 : 266 : 'REPETER' hbloc3 ; 267 : zq = 'PROG' hk*qzero ; 268 : qt = qrela *(('ABS' p)+('ABS' (pmin+dp))); 269 : aqp = 'MASQ' zq 'EGSUPE' qt; 270 : aqpc = 'MASQ' zq 'INFERIEUR' qt; 271 : eaq = aqp*zq + (aqpc*qt); 272 : raq = p - pmin - dp; 273 : aaqa = 'MASQ' ('ABS' raq) 'SUPE' eaq; 274 : raq = raq * aaqa; 275 : hs = 'MASQUE' raq 'EGINFE' 0; 276 : zq = 'PROG' hk*qzero; 277 : qt = qrela *(('ABS' p)+('ABS' (pmax-dp))); 278 : aqp = 'MASQ' zq 'EGSUPE' qt; 279 : aqpc = 'MASQ' zq 'INFERIEUR' qt; 280 : eaq = aqp*zq + aqpc*qt; 281 : raq = p - pmax + dp; 282 : aaqa = 'MASQ' ('ABS' raq) 'SUPE' eaq; 283 : raq = raq * aaqa; 284 : 285 : ht = 'MASQUE' raq 'EGSUPE' 0; 286 : htt = 'MASQUE' ht 'INFERIEUR' 0.5; 287 : hss = 'MASQUE' hs 'INFERIEUR' 0.5; 288 : hh = 'MASQUE' h 'SUPERIEUR' 0; 289 : hhh = 'MASQUE' h 'INFERIEUR' 0; 290 : hminh= h*hhh; 291 : hmaxh= h*hh; 292 : 293 : h = (hs*hmaxh) + ((hss*htt)*h) + (ht*hminh); 294 : 295 : hv = 0; 296 : hi = 0; 297 : 'REPETER' hblo hk; 298 : hi = hi+1 ; 299 : hv = ('EXTR' h hi)**2 + hv; 300 : 'FIN' hblo; 301 : eaq = 'MAXI' ('PROG' qzero (qrela*( ('ABS' hu)+('ABS' 0) ) ) ); 302 : raq = hu - 0; 303 : 'SI' (('ABS' raq) '<EG' eaq); 304 : raq = 0.; 305 : 'FINSI'; 306 : 'SI' (raq '<EG' 0); 307 : 'QUITTER' hbloc2; 308 : 'FINSI'; 309 : eaq = 'MAXI' ( 'PROG' qzero (qrela*( ('ABS' hv)+('ABS' 0) ) ) ); 310 : raq = hv - 0; 311 : 'SI' (('ABS' raq) '<EG' eaq); 312 : raq = 0.; 313 : 'FINSI'; 314 : 'SI' (raq '<EG' 0); 315 : 'SI' ( hkas '<EG' 0); 316 : 'QUITTER' hbloc2; 317 : 'SINON' ; 318 : 'ITERER' hbloc2; 319 : 'FINSI'; 320 : 'FINSI'; 321 : 322 : hrm = hzinf; 323 : hdr = hzinf; 324 : *gounand Masquage de h pour éviter les divisions par zéro 325 : zh = 'MASQUE' ('ABS' h) 'EGINFE' 2.D-56 ; 326 : uh = 'PROG' ('DIME' h) * 1.D0 ; 327 : mzh = '-' uh zh ; 328 : h2 = '+' ('*' h mzh) ('*' ('*' uh 2.D-56) zh) ; 329 : h = h2 ; 330 : * fin masquage 331 : 332 : hl1 = (pmin-p) / h; 333 : hl2 = (pmax-p) / h; 334 : 335 : hl11 = 'MASQUE' hl1 'EGSUPE' hl2; 336 : hl22 = 'MASQUE' hl11 'INFERIEUR' 0.5; 337 : hmaxl = (hl11*hl1) + (hl22*hl2); 338 : 339 : hlr1 = 'MASQUE' hmaxl 'EGINFE' ('PROG' hk*hrm); 340 : hlr2 = 'MASQUE' hlr1 'INFERIEUR' 0.5; 341 : hlr3 = (hlr1*hmaxl) + (hlr2*('PROG' hk*hrm)); 342 : hrm = 'MINI' hlr3; 343 : 344 : hldr = dp / ('ABS' h); 345 : hld1 = 'MASQUE' hldr 'EGINFE' ('PROG' hk*hdr); 346 : hld2 = 'MASQUE' hld1 'INFERIEUR' 0.5; 347 : hld3 = (hld1*hldr) + (hld2*('PROG' hk*hdr)); 348 : hdr = 'MINI' hld3; 349 : 350 : 'REPETER' hbloc1 ; 351 : * recherche du minimum d'une fonction (k variable) sur un segment 352 : * hg hq hr hphi2=seck hrm hdr l xtab y p h; 353 : p1 = 0.6180340051651; 354 : q1 = 0.3819659948349; 355 : ZERO = 2.E-56; 356 : rela = 2.E-14; 357 : ka = 0.; 358 : kb = hrm; 359 : 360 : kprec = 'MAXI' ('PROG' 0 (0.1*hdr)); 361 : 362 : eaq = 'MAXI' ('PROG' zero 363 : (rela*(('ABS' kb)+('ABS' (ka+kprec)))) ); 364 : raq = kb - ka - kprec; 365 : 'SI' (('ABS' raq) '<EG' eaq); 366 : raq = 0. ; 367 : 'FINSI'; 368 : 'SI' (raq '<EG' 0); 369 : kpp = 0.5 * (ka+kb); 370 : 'SINON'; 371 : ku = (p1*ka) + (q1*kb); 372 : kv = (q1*ka) + (p1*kb); 373 : 374 : kpu = p + (ku*h); 375 : kpv = p + (kv*h); 376 : 377 : ktbpu = ('TEXTE' NOMFCT) xtab kpu; 378 : ktbpv = ('TEXTE' NOMFCT) xtab kpv; 379 : 380 : kf kq = 'AJU1' l ktbpu y POI; 381 : kg kq = 'AJU1' l ktbpv y POI; 382 : 383 : 'REPETE' kbloc ; 384 : eaq = 'MAXI' ('PROG' zero (rela*( ('ABS' kf)+('ABS' kg))) ); 385 : raq = kf - kg; 386 : 'SI' (('ABS' raq) '<EG' eaq); 387 : raq = 0. ; 388 : 'FINSI'; 389 : 'SI' (raq '<EG' 0); 390 : kb = kv; 391 : 392 : eaq = 'MAXI' ('PROG' zero 393 : (rela*(('ABS' kb)+('ABS' (ka + kprec)))) ); 394 : raq = kb - ka - kprec; 395 : 'SI' (('ABS' raq) '<EG' eaq); 396 : raq = 0. ; 397 : 'FINSI'; 398 : 'SI' (raq '<EG' 0); 399 : kpp = 0.5 * (ka+kb); 400 : 'QUITTER' kbloc; 401 : 'SINON'; 402 : kv = ku; 403 : kg = kf; 404 : ku = (p1*ka) + (q1*kb); 405 : kpu = p + (ku*h); 406 : ktbpu = ('TEXTE' NOMFCT) xtab kpu; 407 : kf kq = 'AJU1' l ktbpu y POI; 408 : 'ITERER' kbloc; 409 : 'FINSI'; 410 : 'SINON'; 411 : ka = ku; 412 : eaq = 'MAXI' ('PROG' zero 413 : (rela*(('ABS' kb)+('ABS' (ka+kprec)))) ); 414 : raq = kb - ka - kprec; 415 : 'SI' (('ABS' raq) '<EG' eaq); 416 : raq = 0. ; 417 : 'FINSI'; 418 : 'SI' (raq '<EG' 0); 419 : kpp = 0.5 * (ka+kb); 420 : 'QUITTER' kbloc; 421 : 'SINON'; 422 : ku = kv; 423 : kf = kg; 424 : kv = (q1*ka) + (p1*kb); 425 : kpv = p + (kv*h); 426 : ktbpv = ('TEXTE' NOMFCT) xtab kpv; 427 : kg kq = 'AJU1' l ktbpv y POI; 428 : 'ITERER' kbloc; 429 : 'FINSI'; 430 : 'FINSI'; 431 : 'FIN' kbloc; 432 : 'FINSI'; 433 : 434 : kt = p + (kpp*h); 435 : ktbft = ('TEXTE' NOMFCT) xtab kt; 436 : kphi kq= 'AJU1' l ktbft y POI; 437 : 438 : hg = kt; 439 : hq = kq; 440 : hr = kpp; 441 : hphi2 = kphi; 442 : 443 : hiter = hiter + 1; 444 : 'SI' (('MULT' (hiter + IIMPR - 1 ) IIMPR) 'ET' (IMESS '>EG' 2)); 445 : 'MESS' ' itération numéro :' hiter; 446 : 'SI' (l > 0) ; 447 : 'MESS' ' valeurs des ' (@ARR l 0) ' paramètres ' 448 : 'linéaires '; 449 : 'LIST' hq; 450 : 'FINSI' ; 451 : 'SI' (k > 0) ; 452 : 'MESS' ' valeurs des ' (@ARR k 0) ' paramètres ' 453 : 'non linéaires '; 454 : 'LIST' hg; 455 : 'FINSI' ; 456 : 'FINSI'; 457 : 458 : 'SI' ((mxter '>EG' 1) 'ET' (hiter '>EG' mxter)); 459 : 'SI' (IMESS '>EG' 1) ; 460 : 'MESS' 'AJUSTE: ' 461 : 'LE NOMBRE MAXIMUM D ITERATIONS EST ATTEINT !'; 462 : 'FINSI' ; 463 : 'QUITTER' hbloc2; 464 : 'FINSI'; 465 : 466 : eaq = 'MAXI' ('PROG' qzero (qrela*( ('ABS' hr)+('ABS' hdr))) ); 467 : raq = hr - hdr; 468 : 'SI' (('ABS' raq) '<EG' eaq); 469 : raq = 0.; 470 : 'FINSI'; 471 : 'SI' (raq '<EG' 0); 472 : 'SI' (hkas '<EG' 0); 473 : 'QUITTER' hbloc2; 474 : 'SINON'; 475 : 'ITERER' hbloc2; 476 : 'FINSI'; 477 : 'FINSI'; 478 : eaq = 'MAXI' ('PROG' qzero 479 : (qrela*(('ABS' hphi2)+('ABS' hphi1)))); 480 : raq = hphi2 - hphi1; 481 : 'SI' (('ABS' raq) '<EG' eaq); 482 : raq = 0.; 483 : 'FINSI'; 484 : 'SI' (raq '>EG' 0); 485 : hrm = 0.5 * hrm; 486 : 'ITERER' hbloc1; 487 : 'SINON'; 488 : 'QUITTER' hbloc1; 489 : 'FINSI'; 490 : 491 : 'FIN' hbloc1; 492 : 493 : hphi1 = hphi2; 494 : p = hg; 495 : htbfonc = ('TEXTE' NOMFCT) xtab p; 496 : htbderi = ('TEXTE' NOMDERI) xtab p; 497 : hb hg = 'AJU2' hk hq htbfonc htbderi y POI; 498 : 499 : 'SI' (hr '>EG' (0.9999*hrm)); 500 : hf = hg; 501 : hu = hb; 502 : 'ITERER' hbloc2; 503 : 'FINSI'; 504 : 505 : hfg = 0 ; 506 : hi = 0 ; 507 : 'REPETER' hhblo hk; 508 : hi = hi + 1; 509 : hfg = hfg + (('EXTR' hf hi)*('EXTR' hg hi)); 510 : 'FIN' hhblo; 511 : 512 : hc = (hb-hfg) / hu; 513 : hu = hb; 514 : hf = hg; 515 : h = hg + (hc*h); 516 : hkas = 1; 517 : 518 : 'ITERER' hbloc3; 519 : 520 : 'FIN' hbloc3; 521 : 522 : 'FIN' hbloc2; 523 : 524 : htbfonc = ('TEXTE' NOMFCT) xtab p; 525 : hphi q = 'AJU1' l htbfonc y POI; 526 : 527 : 'SI' (IMESS '>EG' 1) ; 528 : 'MESS' 'AJUSTE:' ; 529 : 'SI' (l > 0) ; 530 : 'MESS' ' Les ' (@ARR l 0) 531 : ' paramètres linéaires valent respectivement'; 532 : 'REPETER' bcl1 l ; 533 : mess ' ' ('EXTR' q &bcl1) ; 534 : 'FIN' bcl1 ; 535 : 'SINON' ; 536 : 'MESS' ' Il n y a pas de paramètre linéaire'; 537 : 'FINSI'; 538 : 'SI' (k > 0) ; 539 : 'MESS' ' Les ' (@ARR k 0) 540 : ' paramètres non linéaires valent respectivement'; 541 : 'REPETER' bcl1 k ; 542 : mess ' ' ('EXTR' p &bcl1) ; 543 : 'FIN' bcl1 ; 544 : 'SINON' ; 545 : 'MESS' ' Il n y a pas de paramètre non linéaire'; 546 : 'FINSI' ; 547 : 'MESS' ' Erreur d estimation : ' hphi ; 548 : 'FINSI' ; 549 : 'FINSI'; 550 : 551 : 'FINPROC' q p; 552 : 553 : 554 : 555 :
© Cast3M 2003 - All rights reserved.
Disclaimer