20 !   celor mai mici patrate  folosind metoda Gauss-Newton
30 !             cu algoritmul Levenberg-Marquart
40 !                     D. CIURCHEA 1993
50 !**************************************************************
60 !*             P R O G R A M U L  P R I N C I P A L                       *
70 !**************************************************************
80 RANDOMIZE(.75)
90 CALL CLEAR
100 !CREEN 2
110 !
120 !-------------    Specificarea variabilelor    ---------------
130 !
140 AA=1!Incepe evaluarea preciziei masinii
150 FOR I=1 TO 100
160 REPS=2^(-I)
170 IF AA=AA+REPS THEN 190
180 NEXT I
190 REPS=REPS*2 :: PRINT I-1,REPS :: ! Lungimea mantisei masinii.
200 FOR I=1 TO 70 :: A$=A$&" " :: NEXT I
210 NPCT=121! Nr. punctelor experimentale.
220 NPAR=5! nr. parametrilor din modelul teoretic.
230 DIM D(128)! Vector care contin spectrul experimental.
240 DIM DC(128)! Vector care contine spectrul calculat.
250 DIM W(5),V(5),VT(5),C$(8)! Vect. pt. param. modelului.
260 DIM A(5,5),X(5),B(5),G(5)! Rez. sist. liniar.
270 DIM AINV(5,5)! Calculul matricii inverse.
280 NSIG=5! Nr. cifrelor semn. in raportarea rez. final.
290 PREC=10^(-NSIG)
300 MAXFN=5000! Nr. max. de evaluari ale func. permis de util.
310 PRINT "FITAREA UNUI SPECTRU MOSSBAUER EXPERIMENTAL"
320 PRINT "Programul este in executie"
330 DATA 585,8.4,85,39,69
341 ! Datele experimentale
342 DATA 577,578,575,574,576,574,575,573,574
343 DATA 574,571,578,575,573,575,571,573,575,573
345 DATA 570,573,569,572,569,573,569,575,565,568
346 DATA 569,564,569,560,555,554,545,538,535,536
347 DATA 534,542,551,556,558,564,566,566,572,569
348 DATA 569,571,566,571,575,568,571,572,569,570
349 DATA 565,569,568,571,572,577,577,572,573,575
350 DATA 572,570,568,563,565,567,562,557,551,542
351 DATA 537,538,535,539,544,549,557,556,567,564
352 DATA 567,568,566,571,566,570,565,569,570,572
353 DATA 569,573,569,574,570,570,572,572,569,574
354 DATA 574,574,568,570,574,570,577,571,575,573,568,574
359 HEX$="0123456789ABCDEF"
360 RESTORE 330 :: FOR I=1 TO NPAR!  Citeste param. adevarati ai spectrului
370 READ V(I)!  experimental ai nitroprusiatului de Na.
375 W(I)=V(I)
380 NEXT I
390 RESTORE 342 :: FOR I=1 TO NPCT :: READ BB :: D(I)=BB :: NEXT I
400 !GOSUB 2480 :: GOSUB 3300 :: GOTO 2304
590 ALM=10!100 ! Parametrul Levenberg-Marquart.
630 EPS=SQR(REPS)!*32! Cst. pt. evaluarea derivatei numerice.
640 !
650 !------------------       INITIALIZARE       -----------------
660 !
670 FOR I=1 TO NPAR! Initializarea vectorului de lucru.
680 W(I)=V(I)
690 NEXT I
700 GOSUB 2480!Valorarea functiei cu parametrii initiali.
710 F0=F
720 !
730 !--------------   INCEPUTUL PROCESULUI ITERATIV  -------------
740 !
750 !
760 !=======  CONSTRUIREA SISTEMULUI LINIAR DE ECUATII  ==========
770 !
780 FOR I=1 TO NPAR! Init. mat. Hessian si a term. liberi.
790 FOR J=I TO NPAR
800 A(I,J)=0
810 A(J,I)=0
820 NEXT J
830 W(I)=V(I)! Actualizarea vectorului de lucru.
840 B(I)=0
850 NEXT I
860 !
870 !------ Calculul derivatelor numerice  in fiecare punct ------
880 !
890 GAM=V(2)*V(2)/4
900 FOR IP=1 TO NPCT
910 F1=V(1)-V(2)*V(3)/((IP-V(4))^2+GAM)
920 F1=F1-V(3)*V(2)/((IP-V(5)*2+V(4))^2+GAM)
930 ! F1 este o referinta pentru calculul derivatelor numerice
940 DIF=D(IP)-F1
950 FOR I=1 TO NPAR!Calc. comp.gradientului in punctul IP.
960 AA=ABS(V(I)) :: IF AA<.1 THEN AA=.1! Scala pas.
970 PAS=AA*EPS :: W(I)=V(I)+PAS! Pas  deriv.num.
980 GAMW=W(2)*W(2)/4
990 F2=W(1)-W(2)*W(3)/((IP-W(4))^2+GAMW)
1000 F2=F2-W(3)*W(2)/((IP-W(5)*2+W(4))^2+GAMW)
1010 G(I)=(F2-F1)/PAS! Componenta gradientului.
1020 W(I)=V(I)! Refacerea vectorului de lucru.
1030 NEXT I
1040 !
1050 !--- Calculul matricii Hessian, A si al term. liberi, B  --
1060 !
1070 FOR I=1 TO NPAR
1080 FOR J=I+1 TO NPAR
1090 A(I,J)=A(I,J)+G(I)*G(J)! Mat. Hessian.
1100 A(J,I)=A(I,J)
1110 NEXT J
1120 A(I,I)=A(I,I)+G(I)*G(I)
1130 B(I)=B(I)+DIF*G(I)! Termenii liberi.
1140 NEXT I
1150 NEXT IP
1160 IFN=IFN+(NPAR*NPAR+NPAR)/2
1170 GOSUB 2640! Rezolvarea sist.liniar printr-o metoda adecvata.
1180 !
1190 !===================  AJUSTAREA PASULUI  ====================
1200 !
1210 IF INEF>6 OR IFN>MAXFN THEN 2060
1220 ALFA=1! Incercarea cu un pas cit mai mare
1230 GOSUB 3140
1240 F2=F
1250 IF F2>=F0 THEN 1290
1260 AL=ALFA
1270 F0=F2
1280 GOTO 1540
1290 ALFA=.5! Incercarea cu jumatatea pasului maxim
1300 GOSUB 3140
1310 F1=F
1320 IF F1>=F0 THEN 1360
1330 AL=ALFA
1340 F0=F1
1350 GOTO 1540
1360 ALFA=(3*F0-4*F1+F2)/4/(F0-F1-F1+F2)
1370 !       Ajust. pas. cu  o parabola prin abscisele  0, 0.5, 1.
1380 GOSUB 3140
1390 IF F>=F0 THEN 1440
1400 AL=ALFA
1410 F0=F
1420 GOTO 1540
1430 !
1440 !--------   Decizie pentru cautarea ineficienta   -----------
1450 !
1460 !IF F0 >= F1 AND F0 >= F2 THEN 1270
1470 ALM=ALM*10.1! Se amplifica param. L-M.
1480 INEF=INEF+1
1490 GOSUB 2410! Schimba semnul si marimea corectiilor.
1500 GOTO 1210!Reia ajustarea pasului.
1510 !
1520 ! ==========   P R O G R E S   I N  C A U T A R E   =========
1530 !
1540 ITN=ITN+1
1550 ALM=ALM/2
1560 INEF=0
1570 !CALL CLEAR
1580 FOR I=1 TO NPAR
1590 W(I)=V(I)+AL*X(I)
1600 !LOCATE 13+I,1
1610 DISPLAY AT(6+I,1):USING "v(#)=###.### ":I,V(I)
1620 !LOCATE 13 + I, 65
1630 !PRINT USING " w()=."; I; w(I)
1640 NEXT I
1650 DISPLAY AT(1,1):USING "alfa=###.###^^^^":AL
1660 DISPLAY AT(2,1):USING "Par.L-M=##.##^^^^":ALM
1670 DISPLAY AT(3,1):USING "itn=#####":ITN
1680 DISPLAY AT(4,1):USING "rez=#####.##":F0
1690 NQ=0! Testeaza  cifrele semnif. in parametrii calculati.
1700 FOR I=1 TO NPAR
1710 IF ABS(V(I)/W(I)-1)<PREC THEN NQ=NQ+1
1720 NEXT I
1730 IF NQ=NPAR THEN 2060
1740 !GOSUB 3270
1750 !
1760 !===========  CAUTAREA DE-A LUNGUL GRADIENTULUI  ============
1770 !
1780 !-------------------   O B S E R V A T I E ------------------
1790 !      In aceasta sectiune  urmarim  sa  obtinem  un  progres
1800 ! orict  de   mic  in   minimizarea   functionalei  c.m.m.p.
1810 !      Justificarea  procedeului  a  fost data de Kontorovici
1820 ! care a observat ca matricea Hessianului  nu se schimba prea
1830 ! mult in vecinatatea solutiei. Procedeul este in mod special
1840 ! util daca parametrul Levenberg-Marquart este nenul.
1850 !------------------------------------------------------------
1860 !FOR I=1 TO NPAR
1870 !V(I)=W(I)
1880 !W(I)=W(I)+X(I)*AL/2
1890 !NEXT I
1900 !GOSUB 2480
1910 !IF F<F0 THEN F0=F:GOTO 1650
1920 !
1930 !==========  ACTUALIZAREA PARAMETRILOR OPTIMI ===============
1940 !
1950 FOR I=1 TO NPAR
1960 V(I)=W(I)
1970 NEXT I
1980 KO=0
1990 INEF=0
2000 !===================  RELUAREA CAUTARII  ====================
2010 GOTO 780
2020 !
2030 !===================  SFIRSITUL CAUTARII  ===================
2040 !
2050 GOSUB 3200! Actualizarea param. optimi pt. raportul final.
2060 !FOR I=1 TO 7
2080 !NEXT I
2090 GOSUB 2930!INVERSARE
2100 !CALL CLEAR
2110 !LOCATE 1, 30
2120 DISPLAY AT(21,1):"TERMINAT"
2130 IF NQ<>NPAR THEN 2170
2140 A$="TERM.NORMALA."
2150 A$=A$&" ## CIFRE SEMNIF. IN ### ITERATII"
2160 DISPLAY AT(22,1):USING A$:NSIG,ITN :: GOTO 2230
2170 A$="TERM.ANORMALA "
2180 IF INEF=0 THEN 2210
2190 DISPLAY AT(22,1):A$&"datorata cautarii ineficiente"
2200 GOTO 2230
2210 IF IFN<MAXFN THEN 2230
2220 DISPLAY AT(22,1):A$&". S-a depasit nr. maxim de eval. ale fct."
2230 FOR I=1 TO NPAR
2240 Z$="V(#)=##.####^^^^+/-#.###^^^^"
2250 DISPLAY AT(6+I,1):USING Z$:I,W(I),SQR(ABS(AINV(I,I)))
2260 NEXT I
2270 DISPLAY AT(24,1):"Orice tasta pt viz.finala"
2280 CALL KEY(0,K,S) :: IF S=0 THEN 2280
2290 CALL CLEAR
2300 GOSUB 3270
2304 CALL SOUND(400,110,0) :: CALL SOUND(-400,1000,12) :: GOTO 2304
2330 END
2340 !************************************************************
2350 !*          SFIRSITUL  PROGRAMULUI  PRINCIPAL               *
2360 !************************************************************
2370 !
2380 !************************************************************
2390 !   SUBRUTINA PT. SCH. SEMN. IN CAZUL CAUTARII INEFICIENTE
2400 !************************************************************
2410 FOR I=1 TO NPAR
2420 X(I)=-X(I)/4
2430 NEXT I
2440 RETURN
2450 !************************************************************
2460 !      SUBRUTINA PENTRU CALCULUL FUNCTIONALEI C.M.M.P.
2470 !************************************************************
2480 F=0
2490 IFN=IFN+1
2500 GAMW=W(2)^2/4
2510 FOR I=1 TO NPCT
2520 AUX=W(1)-W(2)*W(3)/((I-W(4))^2+GAMW)
2530 AUX=AUX-W(3)*W(2)/((I-W(5)*2+W(4))^2+GAMW)
2540 DC(I)=AUX
2550 AUX=D(I)-AUX
2560 F=F+AUX*AUX
2570 NEXT I
2580 RETURN
2590 !************************************************************
2600 !       SUBRUTINA PENTRU REZOLVAREA SISTEMULUI LINIAR
2610 !************************************************************
2620 !--Procedeul ales este elim. Gaussiana cu pivotare partiala--
2630 !
2640 FOR I=1 TO NPAR! Initializare
2650 A(I,I)=A(I,I)+ALM
2660 !          Parametrul L.-M. este adaugat pe diag. principala.
2670 NEXT I
2680 FOR K=1 TO NPAR-1! Contorizarea eliminarilor.
2690 !
2700 !------------  Triunghiularizarea matricii A  --------------
2710 !
2720 FOR I=K+1 TO NPAR
2730 T=A(I,K)/A(K,K)
2740 FOR J=K+1 TO NPAR
2750 A(I,J)=A(I,J)-T*A(K,J)
2760 NEXT J
2770 B(I)=B(I)-T*B(K)
2780 NEXT I
2790 NEXT K
2800 !
2810 !--------  Construirea solutiilor sistemului liniar   -------
2820 !
2830 X(NPAR)=B(NPAR)/A(NPAR,NPAR)
2840 FOR I=NPAR-1 TO 1 STEP-1
2850 S=0
2860 FOR J=I+1 TO NPAR
2870 S=S+A(I,J)*X(J)
2880 NEXT J
2890 X(I)=(B(I)-S)/A(I,I)
2900 NEXT I
2910 RETURN
2920 !
2930 !************************************************************
2940 !     SUBRUTINA PENTRU CALCULUL INVERSEI MAT. HESSIAN
2950 !************************************************************
2960 FOR II=1 TO NPAR
2970 FOR J=1 TO NPAR
2980 B(J)=0
2990 NEXT J
3000 B(II)=1
3010 AINV(II,NPAR)=B(NPAR)/A(NPAR,NPAR)
3020 FOR I=NPAR-1 TO 1 STEP-1
3030 S=0
3040 FOR J=I+1 TO NPAR
3050 S=S+A(I,J)*AINV(II,J)
3060 NEXT J
3070 AINV(II,I)=(B(I)-S)/A(I,I)
3080 NEXT I
3090 NEXT II
3100 RETURN
3110 !************************************************************
3120 !      SUBRUTINA PT. INCERC. DIF. PASI IN OPTIM. PARAM.
3130 !************************************************************
3140 FOR I=1 TO NPAR
3150 W(I)=V(I)+ALFA*X(I)
3160 NEXT I
3170 GOSUB 2480
3180 RETURN
3190 !
3200 !************************************************************
3210 ! SUB. PT.CONSTRUIREA VECT. DE LUCRU CU VAL.FINALA A PASULUI
3220 !************************************************************
3230 FOR I=1 TO NPAR
3240 W(I)=V(I)+AL*X(I)
3250 NEXT I
3260 RETURN
3270 !************************************************************
3280 !  SUBRUTINA PENTRU REPREZENTAREA GRAFICA A SPECTRULUI
3290 !************************************************************
3300 ! FOR I = 1 TO NPCT
3310 !   PP2 = (D(I) - MI) / SCL * 150
3320 !   PSET (I * 5, 180 - PP2 + 1)
3330 !   PSET (I * 5, 180 - PP2 - 1)
3340 !   PSET (I * 5, 180 - PP2)
3350 !   PSET (I * 5 - 1, 180 - PP2)
3360 !   PSET (I * 5 + 1, 180 - PP2)
3370 ! NEXT I
3380 !PSET (1, 180 - (DC(1) - MI) / SCL * 150)
3390 !FOR I = 2 TO NPCT
3400 !  LINE -(I * 5, 1 * (180 - (DC(I) - MI) / SCL * 150))
3410 !CALL GRAF(D(),DC(),B,1,121,NPCT)!NEXT I
3420 !RETURN
9000 !SUB GRAF(PLOT(),XP(),B,XMIN,XMAX,NN)
9004 !DIM A$(8)
9005 FOR I=1 TO 14 :: CALL COLOR(I,2,16) :: NEXT I! negru pe alb
9006 CALL CLEAR :: CALL SCREEN(16)
9010 CALL CHAR(33,"AAAAAAAAAAAAAAAA") :: CALL CHAR(143,"")
9020 MAXI=-1E20 :: MINI=1E20 :: FOR I=1 TO NPCT :: MAXI=MAX(MAXI,D(I))
9021 MINI=MIN(MINI,D(I)) :: MAXI=MAX(MAXI,DC(I)) :: MINI=MIN(MINI,DC(I))
9025 NEXT I
9030 FOR I=1 TO NPCT :: D(I)=128-(D(I)-MINI)/(MAXI-MINI)*127+0.5
9031 DC(I)=128-(DC(I)-MINI)/(MAXI-MINI)*127+0.5
9035 B$="0"&B$ :: NEXT I
9037 BB=143
9040 FOR I=1 TO 224 STEP 16!2*121=242 !!
9041 MI=256 :: MA=0 :: C$(1),C$(2),C$(3),C$(4),C$(5),C$(6),C$(7),C$(8)=B$
9044 FOR J=1 TO 8 STEP 2 :: K=D(INT(I/2)+INT(J/2)+1) :: C$(J)=RPT$("1",K-1)&"1"&RPT$("0",128-K)
9045 K1=DC(INT(I/2)+INT(J/2)+1) :: C$(J+1)=RPT$("0",K1-1)&"1"&RPT$("0",128-K1)
9046 MI=MIN(MI,K) :: MA=MAX(MA,K) :: MI=MIN(MI,K1) :: MA=MAX(MA,K1)!CATE SUNT PE VERTICALA
9047 NEXT J
9058 FOR J=INT(MI/8)*8+1 TO MA STEP 8 :: M$="" :: FOR R=1 TO 8 :: K=J+R-1
9060 LOW=VAL(SEG$(C$(5),K,1))*8+VAL(SEG$(C$(6),K,1))*4+VAL(SEG$(C$(7),K,1))*2+VAL(SEG$(C$(8),K,1))+1
9070 HIGH=VAL(SEG$(C$(1),K,1))*8+VAL(SEG$(C$(2),K,1))*4+VAL(SEG$(C$(3),K,1))*2+VAL(SEG$(C$(4),K,1))+1
9072 M$=M$&SEG$(HEX$,HIGH,1)&SEG$(HEX$,LOW,1)
9073 NEXT R
9080 BB=BB-1 :: CALL CHAR(BB,M$) :: CALL VCHAR(J/8+1,INT(I/8)+1,BB) :: NEXT J
9085 IF MI>8 THEN CALL VCHAR(1,INT(I/8)+1,33,INT(MI/8))
9141 MI=256 :: MA=0 :: C$(1),C$(2),C$(3),C$(4),C$(5),C$(6),C$(7),C$(8)=B$
9144 FOR J=9 TO 16 STEP 2 :: K=D(INT(I/2)+INT(J/2)) :: C$(J-8)=RPT$("1",K-1)&"1"&RPT$("0",128-K)
9145 K1=DC(INT(I/2)+INT(J/2)) :: C$(J-7)=RPT$("0",K1-1)&"1"&RPT$("0",128-K1)
9146 MI=MIN(MI,K) :: MA=MAX(MA,K) :: MI=MIN(MI,K1) :: MA=MAX(MA,K1)!CATE SUNT PE VERTICALA
9147 NEXT J
9158 FOR J=INT(MI/8)*8+1 TO MA STEP 8 :: M$="" :: FOR R=1 TO 8 :: K=J+R-1
9160 LOW=VAL(SEG$(C$(5),K,1))*8+VAL(SEG$(C$(6),K,1))*4+VAL(SEG$(C$(7),K,1))*2+VAL(SEG$(C$(8),K,1))+1
9170 HIGH=VAL(SEG$(C$(1),K,1))*8+VAL(SEG$(C$(2),K,1))*4+VAL(SEG$(C$(3),K,1))*2+VAL(SEG$(C$(4),K,1))+1
9172 M$=M$&SEG$(HEX$,HIGH,1)&SEG$(HEX$,LOW,1)
9173 NEXT R
9180 BB=BB-1 :: CALL CHAR(BB,M$) :: CALL VCHAR(J/8+1,INT(I/8)+2,BB) :: NEXT J
9185 IF MI>8 THEN CALL VCHAR(1,INT(I/8)+2,33,INT(MI/8))
9900 NEXT I!:: FOR I=9 TO 14 :: CALL COLOR(I,2,4) :: NEXT I
9910 RETURN