Télécharger ajuste.procedur

Retour à la liste

Numérotation des lignes :

  1. * AJUSTE PROCEDUR PV 14/03/13 21:15:01 7997
  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.99999*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.  
  556.  
  557.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales