| T O P I C R E V I E W |
| duci |
Posted - 07/26/2005 : 17:22:02
'Least squares minimization by using the gradient method. Derivatives are calculated numerically
0 '* Metoda gradientului, versiunea Liberty BASIC *
30 '************************************************************
40 NPCT = 5 '----Numarul punctelor experimentale
50 NPAR = 3 '----Numarul parametrilor modelului teoretic
60 NSIG = 3 '----Numarul cifrelor semnificative in param. calc.
70 PREC = 10 ^ (-1*NSIG) '---- Precizia in calculul parametrilor.
80 DIM D(NPCT):dim T(NPCT)' D(),T() Contin datele experimentale.
90 DIM W(NPAR) ' W() este un vector de lucru.
100 DIM X(NPAR) ' X() contine parametrii modelului.
110 DIM G(NPAR) ' G() contine componentele gradientului.
120 D(1)=12:D(2)=15.7:D(3)=18.8:D(4)=20.2:D(5)=19.2
130 T(1)=0:T(2)=10:T(3)=30:T(4)=100:T(5)=300
150 '***********************************************************
160 '* P R O G R A M U L P R I N C I P A L *
170 '***********************************************************
180 FOR I = 1 TO NPCT
200 D(I) = D(I) / 10 ' Operatie pt. scalarea comp. grad.
210 NEXT I
UpperLeftX=45
UpperleftY=45
WindowWidth=400
WindowHeight= 380
open "grafic final" for graphics as #1
print #1, "fill yellow;color darkblue"
print #1, "goto ";11;" ";21
print #1, "backcolor yellow;\D (dimensiunea de graunte)"
print #1, "down; goto ";11;" ";300
print #1, "goto ";311;" ";300
print #1, "\ T (timpul)"
print #1, "color red;backcolor darkred"
for I=1 to 5
print #1, "down; place "+str$(11+T(I))+" "+str$(300-D(I)*100)
print #1, "circlefilled 6 "
next I
print #1, " flush"
UpperLeftX=410
UpperleftY=100
WindowHeight= 360
open "Reprezentarea Experimental-Calculat" for graphics as #2
print #2, "fill blue"
220 X(1) = 2 '----Valorile initiale ale parametrilor
230 X(2) = 0.5
240 X(3) = 0.5
250 ALFA = 1
260 REPS = 1'------Evaluam lungimea mantisei masinii ----
270 while(REPS+1<>1): REPS = REPS / 2:wend
280 rem IF REPS + 1 <> 1 THEN 270
290 REPS = 2 * REPS
300 EPS = REPS^0.5 '---Constanta folosita la deriv. numerica.
360 ITN = 0 ' Initializam contorul iteratiilor
370 INEF = 0
380 FOR I = 1 TO NPAR '----Initializam vectorul de lucru
390 W(I) = X(I)
400 NEXT I
410 GOSUB 1290
420 F0 = F '----Evaluarea reziduului cu param. initiali
430 '
440 '================ INCEPUTUL PROCESULUI ITERATIV ============
450 '
460 GNORM = 0 '------ Initializam norma gradientului.
470 FOR I = 1 TO NPAR '-- Incepem calculul derivatelor numerice.
480 Z=W(I)
490 IF W(I)>0.1 THEN PAS=ABS(W(I))*EPS ELSE PAS=0.1*EPS
530 W(I) = X(I) + PAS
540 GOSUB 1290
550 G(I) = (F - F0) / PAS ' Eval. numerica a comp. grad.
560 GNORM = GNORM + G(I) * G(I)
570 W(I) = Z
580 NEXT I
590 FOR I = 1 TO NPAR '--Incercam un progres in dir. de cautare.
600 W(I) = X(I) - ALFA * G(I) ' Noua aprox. a coord. minim.
610 NEXT I
620 GOSUB 1290
670 IQUIT = 0 '-------Incercam criteriile de oprire.
680 FOR I = 1 TO NPAR
690 IF ABS(W(I) / X(I) - 1) < PREC THEN IQUIT=IQUIT+1
710 NEXT I
720 IF IQUIT = NPAR THEN 960 ' Terminat !
730 IF ALFA < 100 * REPS THEN 780
740 IF F < F0 THEN 830 '----In caz de progres.
750 INEF = INEF + 1'--- Daca cautarea este ineficienta.
760 ALFA = ALFA / 2'----In caz de regres, injumatatim pasul.
770 GOTO 590 '----Regres.
780 A$ = "Cautarea a devenit ineficienta."
790 GOTO 980 ' Terminatia anormala a cautarii.
800 '
810 '=================== PROGRES IN CAUTARE ====================
820 '
830 INEF = 0 '
840 A$ = "Cautarea a fost eficienta."
860 FOR I = 1 TO NPAR'------ Actualizarea parametrilor
870 X(I) = W(I)
880 NEXT I
890 GOSUB 1510
900 F0 = F
910 ALFA = 1
920 ITN = ITN + 1
921 PRINT " ITERATIA = ";ITN;" REZIDUU = ";F0
922 FOR I = 1 TO NPAR
923 PRINT "x(";I;")=";X(I); " G(";I;")=";G(I)
924 NEXT I
925 PRINT "---------------------------------------------"
930 GOTO 460' -------Reluarea cautarii
950 '
960 '=========== SFARSITUL PROCESULUI ITERATIV ===============
970 '
990 PRINT:PRINT " TERMINAT. "+A$
1000 PRINT "================================================="
1010 PRINT " PARAMETRII FINALI AI MODELULUI"
1020 PRINT "-------------------------------------------------"
1030 PRINT "I PARAMETRUL COMPONENTA GRADIENTULUI"
1040 PRINT "-------------------------------------------------"
notice "Procesul iterativ s-a incheiat:"+chr$(13)+A$
1045 GOSUB 1680
1050 X(1) = X(1) * 10 ' Refacem factorul de scala pt. X(1).
1060 V$ = " .^^^^ .^^^^"
1070 FOR I = 1 TO NPAR
1080 PRINT I;". X(";I;")="; X(I);" G(";I;")="; G(I)
1090 NEXT I
1100 PRINT "-------------------------------------------------"
1110 PRINT "NORMA GRADIENTULUI= "; GNORM^0.5
1120 PRINT "REZIDUUL LA ITERATIA FINALA="; F * 100
1130 PRINT "NUMARUL DE ITERATII="; ITN
1140 Z$ = "S-au efectuat cautari inef. ,ALFA=.^^^^"
1150 PRINT "S-au efectuat ";INEF;" cautari ineficiente. ,ALFA="; ALFA
notice "Mesaj de incheiere"+chr$(13)+"Apasati <ENTER> in fereastra principala pentru a termina executia programului"
rem Salveza o fereastrele grafica
print #1, "trapclose [quit]"
[quit] print #1, "getbmp drawing 1 1 380 380"
bmpsave "drawing", "d:\lb201w\duci\1.bmp"
CLOSE #1
print #2, "trapclose [quit1]"
[quit1] print #2, "getbmp drawing1 1 1 380 380"
bmpsave "drawing1", "d:\lb201w\duci\2.bmp"
CLOSE #2
input "Apasa <ENTER> pentru a termina executia programului";r$
1230 END
1240 '
1250 '**********************************************************
1260 '* SUBRUTINA PENTRU CALCULUL FUNCTIONALEI C.M.M.P. *
1270 '**********************************************************
1280 '
1290 F = 0 ' Initializarea sumei.
1300 FOR J = 1 TO NPCT
1310 A = W(1) * (1 - W(2) * EXP(-1*W(3) * T(J)))
1320 A = D(J) - A
1330 F = F + A * A
1340 NEXT J
1350 RETURN
1360 '
1510 '**********************************************************
1520 '* SUBRUTINA GRAFICA PENTRU REPREZENTAREA PUNCTELOR *
1530 '**********************************************************
print #2,"cls"
print #2, "fill green;up;goto ";18;" ";31
print #2, "color darkred;backcolor green;\Calculat"
print #2, "goto ";11;" ";11;
print #2,"down;box 311 311;up;goto 11 311;down; line 11 311 311 11"
print #2, "up; goto ";201;" ";295
print #2, "\ Experimental"
print #2, "flush"
for I=1 to NPCT
xx=X(1)*(1-X(2)*exp(-1*X(3)*T(I)))
y.p=11+D(I)*100
x.p=311-xx*100
print #2, "up;goto "+str$(x.p)+" "+str$(y.p)
print #2,"down;color blue;backcolor red;circlefilled 5"
next I
1550 return
1680 '**********************************************************
1690 '* SUBRUTINA PENTRU REPREZENTAREA GRAFICA FINALA *
1700 '**********************************************************
print #1, "color black; size 2"
xx=X(1)*(1-1*X(2))
print #1, "up; place "+str$(11)+" "+str$(300-xx*100)
for tt=0 to 300
xx=X(1)*(1-X(2)*exp(-1*X(3)*tt))
print #1, "down; goto "+str$(tt+11)+" "+str$(300-xx*100)
next tt
print #1, "flush"
1800 RETURN
|
| 3 L A T E S T R E P L I E S (Newest First) |
| duci |
Posted - 08/07/2005 : 10:27:17
Superexample
Least Squares Fitting of a Mossbauer Spectrum
130 REPS = 1 'Incepe evaluarea lungimii mantisei masinii.
140 WHILE REPS+1<>1:REPS=REPS/2:WEND
180 REPS = REPS * 2'Lungimea mantisei masinii
220 NPCT = 5 ' Nr. punctelor experimentale.
230 NPAR = 3 ' nr. parametrilor din modelul teoretic.
240 DIM D(NPCT):dim T(NPCT)' Vect. care contin datele experim.
250 DIM W(NPAR):dim V(NPAR):dim VT(NPAR)' Vect.de lucru pt. param.
260 DIM A(NPAR, NPAR):dim X(NPAR):dim B(NPAR)
270 rem DIM AINV(NPAR, NPAR)' Calculul matr. inverse.
271 D(1)=12:D(2)=15.7:D(3)=18.8:D(4)=20.2:D(5)=19.2
272 T(1)=0:T(2)=10:T(3)=30:T(4)=100:T(5)=300
280 NSIG = 5 ' Nr. cifre semnif. in raportarea rez. final.
290 PREC = 10 ^ (-1*NSIG)
300 MAXFN = 5000' Nr. maxim de evaluari ale functiei c.m.m.p.
UpperLeftX=4
UpperleftY=-100
WindowWidth=420
WindowHeight= 420
open "grafic final" for graphics as #1
print #1, "fill yellow;color darkblue"
print #1, "goto ";11;" ";21
print #1, "backcolor yellow;\D (dimensiunea de graunte)"
print #1, "down; goto ";11;" ";300
print #1, "goto ";311;" ";300
print #1, "\ T (timpul)"
print #1, "color red;backcolor darkred"
for I=1 to 5
print #1, "down; place "+str$(11+T(I))+" "+str$(300-D(I)*100)
print #1, "circlefilled 6 "
next I
print #1, " flush"
WindowHeight= 420
UpperLeftX=550
UpperleftY=-100
open "Reprezentarea Experimental-Calculat" for graphics as #2
print #2, "fill blue"
360 ALM = 1 ' Parametrul Levenberg-Marquart.
370 V(1) = 20' Valorile initiale ale parametrilor.
380 V(2) = 0.5
390 V(3) = 0.5
400 EPS = REPS^0.5 '---Constanta folosita la deriv. numerica.
480 '------------------ INITIALIZARE ---------------
490 '
500 FOR I = 1 TO NPAR ' Incepe initializarea vect. de lucru.
510 W(I) = V(I)
520 NEXT I
530 GOSUB 2180 ' Calculul valorii functiei cu param. initiali.
540 FF = F
550 '
560 '---------------- INCEPUTUL PROCESULUI ITERATIV ---------
570 '
580 FOR I = 1 TO NPAR 'Afisaj val. corectiilor la param. calc.
590 PRINT I;" ";V(I)
600 NEXT I
610 '
620 '======== CONSTRUIREA SISTEMULUI LINIAR DE ECUATII =======
630 '
640 FOR I = 1 TO NPAR' Init. mat. Hessian si a term. liberi.
650 FOR J = I TO NPAR
660 A(I, J) = 0
670 A(J, I) = 0
680 NEXT J
690 W(I) = V(I)' Actualizarea vectorului de lucru.
700 B(I) = 0
710 NEXT I
720 '
730 '------ Calculul derivatelor numerice in fiecare punct ----
740 '
750 FOR IP = 1 TO NPCT
760 F1 = V(1) * (1 - V(2) * EXP(-1*V(3)*T(IP)))' Referinta
770 ' pentru calculul derivatei numerice.
780 DIF = D(IP) - F1
790 FOR I = 1 TO NPAR ' Calculul comp. grad. in punctul IP.
800 AA = ABS(V(I))
810 IF AA>0.1 THEN PAS=EPS*AA ELSE PAS=0.1*EPS' Scalarea pasului in deriv.num.
830 W(I) = V(I) + PAS ' Pasul in derivata numerica.
840 F2 = W(1) * (1 - W(2) * EXP(-1*W(3)*T(IP)))
850 G(I) = (F2 -1*F1) / PAS ' Componenta gradientului.
860 W(I) = V(I) ' Refacerea vectorului de lucru.
870 NEXT I
880 '
890 '- Calculul matricii Hessian, A si al term. liberi, B ---
900 '
910 FOR I = 1 TO NPAR
920 FOR J = I + 1 TO NPAR
930 A(I, J) = A(I, J) + G(I) * G(J)'Elem. matr.
940 A(J, I) = A(I, J) ' Hessian, simetrica.
950 NEXT J
960 A(I, I) = A(I, I) + G(I) * G(I)
970 B(I) = B(I) + DIF * G(I) ' Termenii liberi.
980 NEXT I
990 NEXT IP
1000 IFN = IFN + (NPAR * NPAR + NPAR) / 2
1010 GOSUB 2310 ' Rez. sist. liniar de ecuatii A(,)*X()=B().
1020 '
1030 '================= AJUSTAREA PASULUI ====================
1040 '
1050 IF INEF > 6 OR IFN > MAXFN THEN 1840
1060 ALFA = 1 ' Incercarea cu un pas cit mai mare
1070 GOSUB 3070
1080 F2 = F
1090 IF F2 < FF THEN AL = ALFA: FF = F2: GOTO 1290
1100 ALFA = 0.5 ' Incercarea cu jumatatea pasului maxim
1110 GOSUB 3070
1120 F1 = F
1130 IF F1 < FF THEN AL = ALFA: FF = F1: GOTO 1290
1140 'Ajustarea pasului cu o parabola prin abscisele ALFA=0,0.5,1.
1150 ALFA = (3*FF-4*F1+F2)/(FF-F1-F1+F2)/4
1160 GOSUB 3070
1170 IF F < FF THEN AL = ALFA: FF = F: GOTO 1290
1180 '
1190 '---------- Decizie pentru cautarea ineficienta -------
1200 '
1210 IF FF > F1 AND FF > F2 THEN 1270
1220 ALM = ALM * 1.1 ' Marim val. param. Levenberg-Marquart.
1230 INEF = INEF + 1
1240 GOSUB 2110 ' Schimbam semnul si marimea corectiilor X().
1250 GOTO 1050 'Reluam cautarea unui progres cu aceste cor.
1260 '
1270 ' ========= P R O G R E S I N C A U T A R E ========
1280 '
1290 ITN = ITN + 1
1300 ALM = ALM / 5 ' Param. Levenberg-Marquart este micsorat.
1310 INEF = 0
1315 PRINT "I v vechi() v nou() "
1320 FOR I = 1 TO NPAR
1330 W(I) = V(I) + AL * X(I)
1360 print I;" "; V(I);" "; W(I)
1370 NEXT I
1380 PRINT " alfa= "; AL;
1390 PRINT " Parametrul Levenberg-Marquart= "; ALM
1400 PRINT " itn= "; ITN;
1410 PRINT " rez= "; FF
1420 NQ = 0' Test pentru cifrele semnif. in param. calculati.
1430 FOR I = 1 TO NPAR
1440 IF ABS(V(I) / W(I) - 1) < PREC THEN NQ = NQ + 1
1450 NEXT I
1460 IF NQ = NPAR AND ITN > 2 THEN 1840
1470 '
1480 '=========== CAUTAREA DE-A LUNGUL GRADIENTULUI ==========
1580 PRINT "Kantorovici "; KO
1590 FOR I = 1 TO NPAR
1600 VT(I) = W(I)
1610 W(I) = W(I) + X(I) * AL / 2
1620 NEXT I
1630 KO = KO + 1
1640 GOSUB 2180
1645
1650 IF F < FF AND KO<10 THEN FF = F: GOTO 1290
1660 '
1670 '============== ACTUALIZAREA PARAMETRILOR OPTIMI =========
1680 '
1690 GOSUB 3190
1700 FOR I = 1 TO NPAR
1710 V(I) = VT(I)
1720 NEXT I
1730 GOSUB 3310
1750 KO = 0
1760 INEF = 0
1780 GOTO 550'================ RELUAREA CAUTARII
1800 '================ SFIRSITUL CAUTARII ================
1820 GOSUB 3130 ' Actualizarea param. optimi pt. raportul final.
1830 IFIN = 1
1840 REM notice "TERMINAT"
1890 IF NQ <> NPAR THEN 1940
1900 A$ = "TERMINATIE NORMALA."
1910 A$ = A$ + " S-AU OBTINUT ";NSIG;" CIFRE SEMNIF. IN ";ITN;" ITERATII"
1920 GOSUB 3310
1930 print A$: GOTO 1990
1940 A$ = "TERMINATIE ANORMALA. -> "
1950 B$ = "Datorata cautarii ineficiente"
1960 C$ = " S-a depasit nr. maxim de eval. ale fct."
1970 IF INEF > 0 THEN PRINT A$ + B$: GOTO 1990
1980 IF IFN > MAXFN THEN PRINT A$ + C$
1990 FOR I = 1 TO NPAR
2010 PRINT "Parametrul "; I;"= ";V(I)
2020 NEXT I
2030 GOSUB 3470
'NOTICE "TERMINAT"
INPUT R$
CLOSE #2
CLOSE #1
2040 END
2060 '* SFIRSITUL PROGRAMULUI PRINCIPAL *
2080 '
2090 '== SUB. PT. SCHIMBAREA SEMNULUI IN CAZUL CAUTARII INEF. ==
2110 FOR I = 1 TO NPAR
2120 X(I) = -1*X(I) / 5
2130 NEXT I
2140 RETURN
2150 '
2160 '==== SUBRUTINA PENTRU CALCULUL FUNCTIONALEI C.M.M.P. ===
2170 '
2180 IFN = IFN + 1
2190 F = 0
2210 FOR I = 1 TO NPCT
2220 AUX = W(1) * (1 - W(2) * EXP(-1*W(3) * T(I)))
2230 AUX = D(I) - AUX
2240 F = F + AUX * AUX
2250 NEXT I
2260 RETURN
2270 '
2280 '=== SUBRUTINA PENTRU REZOLVAREA SISTEMULUI LINIAR ===
2290 '
2300 ' Procedeul este eliminarea Gaussiana cu pivotare partiala
2310 '
2320 FOR I = 1 TO NPAR' Initializare
2330 A(I, I) = A(I, I) + ALM' Parametrul L-M este adaugat
2340 ' pe diagonala principala.
2360 NEXT I
2370 FOR K = 1 TO NPAR - 1' Contorizarea eliminarilor.
2380 '
2390 '--------------- Cautarea liniei pivot ---------------
2400 '
2410 L = K
2420 I = K
2430 I = I + 1
2440 IF ABS(A(I, K)) > ABS(A(L, K)) THEN L = I
2450 IF I < NPAR - 1 THEN 2430
2460 IF L = K THEN 2590
2470 FOR J = K TO NPAR
2480 T = A(K, J)
2490 A(K, J) = A(L, J)
2500 A(L, J) = T
2510 NEXT J
2520 T = B(K): B(K) = B(L): B(L) = T
2570 '----------- Triunghiularizarea metricii A ---------
2590 FOR I = K + 1 TO NPAR
2600 T = A(I, K) / A(K, K)
2610 A(I, K) = T ' Detaliile transformarii sunt
2620 ' salvate sub diagonala mat. Hessian.
2630 FOR J = K + 1 TO NPAR
2640 A(I, J) = A(I, J) - T * A(K, J)
2650 NEXT J
2660 B(I) = B(I) - T * B(K)
2670 NEXT I
2680 NEXT K
2700 '------- Construirea solutiilor sistemului liniar ------
2720 X(NPAR) = B(NPAR) / A(NPAR, NPAR)
2730 FOR I = NPAR - 1 TO 1 STEP -1
2740 S = 0
2750 FOR J = I + 1 TO NPAR
2760 S = S + A(I, J) * X(J)
2770 NEXT J
2780 X(I) = (B(I) - S) / A(I, I)
2790 NEXT I
2800 RETURN
3050 '== SUB. PT. INCERCAREA PASILOR IN OPTIMIZAREA PARAM. ====
3070 FOR I = 1 TO NPAR
3080 W(I) = V(I) + ALFA * X(I)
3090 NEXT I
3100 GOSUB 2180
3110 RETURN
3130 '===== SUBRUTINA PENTRU CONSTRUIREA VECTORULUI DE LUCRU CU
3140 ' VALOAREA FINALA A PASULUI =================
3150 FOR I = 1 TO NPAR
3160 W(I) = V(I) + AL * X(I)
3170 NEXT I
3180 RETURN
3190 REM
3200 '====== PROCEDURA GRAFICA PENTRU STERGEREA PUNCTELOR ===
3220 return
3310 '====== SUBRUTINA PENTRU REPREZENTAREA GRAFICA EXPERIMENTAL-CALCULAT =====
print #2,"cls"
print #2, "fill green;up;goto ";18;" ";31
print #2, "color darkred;backcolor green;\Experimental"
print #2, "goto ";11;" ";11;
print #2,"down;box 311 311;up;goto 11 311;down; line 11 311 311 11"
print #2, "up; goto ";201;" ";295
print #2, "\ Calculat"
print #2, "up; goto ";50;" ";95
print #2, "backcolor yellow"
REM print #2,"\Iteratia "+str$(ITN)
print #2, "up; goto ";185;" ";200
REM print #2,"\Reziduu ="+str$(int(F0*10000)/100)
print #2, "flush"
for I=1 to NPCT
xx=V(1)*(1-V(2)*exp(-1*V(3)*T(I)))
y.p=11+D(I)*10
x.p=311-xx*10
print #2, "up;goto "+str$(x.p)+" "+str$(y.p)
print #2,"down;color blue;backcolor red;circlefilled 5"
next I
3460 RETURN
3470 '
3480 '====== SUBRUTINA PENTRU REPREZENTAREA GRAFICA FINALA =====
3490 '
print #1, "fill yellow;color darkblue"
print #1, "goto ";11;" ";21
print #1, "backcolor yellow;\D (dimensiunea de graunte)"
print #1, "down; goto ";11;" ";300
print #1, "goto ";311;" ";300
print #1, "\ T (timpul)"
print #1, "color red;backcolor blue"
for I=1 to 5
print #1, "down; place "+str$(11+T(I))+" "+str$(300-D(I)*10)
print #1, "circlefilled 6 "
next I
print #1, "color black; size 2"
xx=V(1)*(1-1*V(2))
print #1, "up; place "+str$(11)+" "+str$(300-xx*10)
for tt=0 to 300
xx=V(1)*(1-V(2)*exp(-1*V(3)*tt))
print #1, "down; goto "+str$(tt+11)+" "+str$(300-xx*10)
next tt
print #1, " flush"
3710 RETURN
|
| duci |
Posted - 07/26/2005 : 17:26:13
'This program will do a least squares minimization by using the Nelder and Mead algorithm10 '*************************************************************
20 ' Program pentru minimizarea functionalei celor mai mici
30 ' patrate folosind cautarea SIMPLEX (algoritmul Nelder & Mead)
40 ' D. Ciurchea 1992
50 '*************************************************************
60 NPAR = 3
70 NPCT = 5
90 DIM D(NPCT):DIM T(NPCT):DIM X(NPAR)
91 D(1)=12:D(2)=15.7:D(3)=18.8:D(4)=20.2:D(5)=19.2
92 T(1)=0:T(2)=10:T(3)=30:T(4)=100:T(5)=300
140 X(1) = 20: X(2) = .4: X(3) = .3 ' Parametrii initiali.
UpperLeftX=45
UpperleftY=45
WindowWidth=400
WindowHeight= 380
open "grafic final" for graphics as #1
print #1, "fill yellow;color darkblue"
print #1, "goto ";11;" ";21
print #1, "backcolor yellow;\D (dimensiunea de graunte)"
print #1, "down; goto ";11;" ";300
print #1, "goto ";311;" ";300
print #1, "\ T (timpul)"
print #1, "color red;backcolor darkred"
for I=1 to 5
print #1, "down; place "+str$(11+T(I))+" "+str$(300-D(I)*10)
print #1, "circlefilled 6 "
next I
print #1, " flush"
open "Reprezentarea Experimental-Calculat" for graphics as #2
200 GOSUB 230 ' ----------------- Cheama algoritmul SIMPLEX
210 GOSUB 1810 '----------------- Procedura grafica finala
211 GOSUB 1360
NOTICE "TERMINAT"
rem Salveza fereastrele grafice
print #1, "trapclose [quit]"
[quit] print #1, "getbmp drawing 1 1 380 380"
bmpsave "drawing", "d:\lb201w\duci\1.bmp"
CLOSE #1
print #2, "trapclose [quit1]"
[quit1] print #2, "getbmp drawing1 1 1 380 380"
bmpsave "drawing1", "d:\lb201w\duci\2.bmp"
CLOSE #2
input r$
220 END
230 DIM XX(NPAR, NPAR + 1)
240 '------ XX(,) contine virfurile simplexului;
250 DIM G(NPAR + 1):DIM H(NPAR):DIM W(NPAR):DIM V(NPAR)
260 '------ G() contine valorile F(W(i));
270 '------ H() contine W()mediu-W()i;
280 '------ W() contine W()mediu;
290 '------ V() este ajutator pt. calc. functionalei c.m.m.p.
300 FOR J = 1 TO NPAR + 1 ' Initializam cumva un simplex
310 FOR I = 1 TO NPAR
320 XX(I, J) = X(I)
330 IF J = I + 1 THEN XX(I, J) = X(I) * 1.01
340 NEXT I
350 NEXT J
360 IQUIT = 0'--------------- Initializam criteriul de oprire
370 ITN = ITN + 1
380 PRINT "TOT.ITERATII= "; ITN
390 PRINT "REZIDUU= "; F
400 IF IQUIT >= 5 THEN RETURN ' Terminarea cautarii
410 '***********************************************************
420 ' CALCULAM VALORILE FUNCTIONALEI IN VARFURILE SIMPLEXULUI
430 '***********************************************************
440 FOR J = 1 TO NPAR + 1
450 FOR I = 1 TO NPAR
460 V(I) = XX(I, J)
470 NEXT I
480 GOSUB 1420
490 G(J) = F
500 NEXT J
510 '***********************************************************
520 ' Incepem identificarea valorii celei mai defavorabile
530 ' pentru problema noastra, VAB, care conduce la valoarea cea
540 ' mai mare pentru functionala c.m.m.p., corespunzatoare
550 ' varfului XX(JVAB) al simplexului.
560 '***********************************************************
570 VAP = 1000000000000
580 VAB = -0.0000000000001
590 FOR J = 1 TO NPAR + 1
600 IF VAP > G(J) THEN VAP = G(J): JVAP = J
610 IF VAB < G(J) THEN VAB = G(J): JVAB = J
620 NEXT J
630 '------- Initializam vectorul centrului de greutate
640 FOR I = 1 TO NPAR
650 W(I) = 0
660 NEXT I
670 '--- Calculam coordonatele varfurilor simplexului, exceptand
680 '----- punctul cel mai defavorabil
690 FOR J = 1 TO NPAR + 1
700 FOR I = 1 TO NPAR
710 IF J <> JVAB THEN W(I) = W(I) + XX(I, J)
720 NEXT I
730 NEXT J
740 FOR I = 1 TO NPAR
750 AA = W(I) / NPAR '------------- W()mediu-W()i--------
760 H(I) = AA - XX(I, JVAB)
770 W(I) = AA '------------Coordonatele centrului de greutate
780 NEXT I
790 GG = 1 '-------Incercam un simplex cu gama=1
800 GOSUB 1170
810 GOSUB 1420
820 F0 = F
830 IF F0 <= VAB THEN 900
840 CONTOR = CONTOR + 1 '----Dir. de cautare este rau aleasa.
850 IQUIT = IQUIT + 1 ' Cautarea este ineficienta.
860 PRINT "RESTRINGERI="; CONTOR
870 GG = -.5 '----- Inapoi cu gama=-0.5.
880 GOSUB 1250
890 GOTO 370 '----------------Restart pt. gama=-0.5
900 IQUIT = 0'-----------------Am progresat
910 REM GG=1:GOSUB 1680:GOTO 1070 ' Fara expansiuni ale simplexului
920 GG = 2 '---------Incercam o expansiune.
930 GOSUB 1170
940 GOSUB 1420
950 FP = F
960 GG = .5 '--------Incercam o contractie
970 GOSUB 1170
980 GOSUB 1420
990 FM = F
1000 IF FP < F0 THEN 1100 ELSE IF FM < F0 THEN 1010 ELSE 1070
1010 ICON = ICON + 1
1020 PRINT "CONTRACTII="; ICON
1030 GG = .5
1040 GOSUB 1250 ' Calculeaza noul varf simplex prin contractie
1050 GOTO 1150 ' Afiseaza si continua
1060 REM Simplex normal
1070 GG = 1
1080 GOSUB 1250
1090 GOTO 1150
1100 IEXP = IEXP + 1
1110 PRINT "EXPANSIUNI="; IEXP
1120 GG = 2
1130 GOSUB 1250': GOSUB 3000
1140 'GOTO 1610
1150 GOSUB 1360 '----------Afiseaza parametrii V(I)
1160 GOTO 370 '---------------RESTART
1170 '*******************************************************
1180 ' Procedura care propune noul varf simplex urmand a fi
1190 ' incercat in functionala c.m.m.p.
1200 '*******************************************************
1210 FOR I = 1 TO NPAR
1220 V(I) = W(I) + GG * H(I)
1230 NEXT I
1240 RETURN
1250 '********************************************************
1260 ' Procedura de actualizare a varfului simplex.
1270 ' Inlocuieste coordonatele punctului cu valori
1280 ' neconvenabile cu noile valori obtinute prin cautare.
1290 '********************************************************
1300 GOSUB 1530 '-----Grafica de stergere a punctelor vechi
1310 FOR I = 1 TO NPAR
1320 XX(I, JVAB) = W(I) + GG * H(I)
1330 NEXT I
1340 GOSUB 1670 '-----Grafica calculat/experimental
1350 RETURN
1360 '*********************************************************
1370 ' Procedura pentru afisarea parametrilor V(I).
1380 '*********************************************************
1400 PRINT "PARAMETRII : "; V(1);" "; V(2);" "; V(3)
1410 RETURN
1420 '*********************************************************
1430 ' Procedura pentru calculul functionalei c.m.m.p.
1440 '*********************************************************
1450 '
1460 F = 0
1470 FOR II = 1 TO NPCT
1480 B = V(1) * (1 - V(2) * EXP(-1*V(3) * T(II)))
1490 B = D(II) - B
1500 F = F + B * B
1510 NEXT II
1520 RETURN
1530 '********************************************************
1540 ' PROCEDURA GRAFICA PENTRU STERGEREA PUNCTELOR
1550 '********************************************************
1560 '
1570 RETURN
1670 '**********************************************************
1680 ' PROCEDURA GRAFICA PENTRU REPREZENTAREA PUNCTELOR
1690 '**********************************************************
1700 '
print #2,"cls"
print #2, "fill green;up;goto ";18;" ";31
print #2, "color darkred;backcolor green;\Calculat"
print #2, "goto ";11;" ";11;
print #2,"down;box 311 311;up;goto 11 311;down; line 11 311 311 11"
print #2, "up; goto ";201;" ";295
print #2, "\ Experimental"
print #2, "flush"
for I=1 to NPCT
xa=X(1)*(1-V(2)*exp(-1*V(3)*T(I)))
y.p=11+D(I)*10
x.p=311-xa*10
print #2, "up;goto "+str$(x.p)+" "+str$(y.p)
print #2,"down;color blue;backcolor red;circlefilled 5"
next I
1710 RETURN
1810 '**********************************************************
1820 ' PROCEDURA GRAFICA PENTRU REPREZENTAREA GRAFICA FINALA
1830 '**********************************************************
1840 '
print #1, "color black; size 2"
xx=V(1)*(1-1*V(2))
print #1, "up; place "+str$(11)+" "+str$(300-xx*10)
for tt=0 to 300
xx=X(1)*(1-V(2)*exp(-1*V(3)*tt))
print #1, "down; goto "+str$(tt+11)+" "+str$(300-xx*10)
next I
print #1, "flush"
1850 RETURN
|
| duci |
Posted - 07/26/2005 : 17:25:02
'This code will fit a (simulated) Mossbauer spectrum by using least squares Gauss-Newton minimization and a Levenberg-Marquardt parameter
10 ' Program pentru fitarea unui spectru Mossbauer prin metoda
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 '**************************************************************
110 '
120 '------------- Specificarea variabilelor ---------------
130 '
140 AA = 1 'Incepe evaluarea preciziei masinii
150 FOR I = 1 TO 100
160 REPS = 1/(2 ^ (I))
170 IF AA = AA + REPS THEN 190
180 NEXT I
190 REPS = REPS * 2' Lungimea mantisei masinii.
210 NPCT = 121 ' Nr. punctelor experimentale.
220 NPAR = 5 ' nr. parametrilor din modelul teoretic.
230 DIM D(NPCT) ' Vector care contin spectrul experimental.
240 DIM DC(NPCT) ' Vector care contine spectrul calculat.
250 DIM w(NPAR):DIM V(NPAR):DIM VT(NPAR)' Vect. pt. param. modelului.
260 DIM A(NPAR, NPAR):DIM X(NPAR):DIM B(NPAR)' Rez. sist. liniar.
270 DIM AINV(NPAR, NPAR)' Calculul matricii inverse.
280 NSIG = 5 ' Nr. cifrelor semn. in raportarea rez. final.
290 PREC = 1/10 ^ (NSIG)
300 MAXFN = 500000' Nr. max. de evaluari ale func. permis de util.
310 PRINT "FITAREA UNUI SPECTRU MOSSBAUER SIMULAT"
330 DATA 575,8.4,85,39,60
340 DATA 480,14,60,52,73
350 SZ = 0
360 V(1)=575:V(2)=8.4 :V(3)=85 :V(4)= 39:V(5)= 60
390 G = V(2) ^ 2 / 4 ' Patratul semilargimii la semiinaltime
400 MI = 100000
410 MA =-1110
420 FOR I = 1 TO NPCT
430 ZG = 0.0 * (RND - .5)' Simularea erorilor experimentale.
440 BM = I - V(4) ' Pozitia picului de la viteze negative,
450 BP = I - V(5) * 2 + V(4)'respectiv pt.viteze pozitive
460 ' la centrul spectrului dat de V(5).
470 A = V(2) * V(3)
480 ' A= Aria picului,V(3) * semilargimea V(2) * 2.
490 D(I) = V(1) - A / (BM * BM + G) - A / (BP * BP + G)
500 ' Vectorul D() contine intensitatea adevarata.
510 ZG = ZG * SQR(D(I))
520 ' Erorile "experimentale" urmeaza o distributie Poisson.
530 SZ = SZ + ZG * ZG ' Suma erorilor introduse in spectru.
540 D(I) = D(I) + ZG' Intens. "masurata", cu er. "exper.".
550 IF MA < D(I) THEN MA = D(I) ' Scalarea vectorului D().
560 IF MI > D(I) THEN MI = D(I)
570 NEXT I
580 SCL = (MA - MI) ' Factorul de scala.
UpperLeftX=4
UpperleftY=-100
WindowWidth=690
WindowHeight= 550
open "Rafinarea parametrilor unui spectru Moessbauer simulat" for graphics as #1
print #1, "fill white;color darkblue"
print #1, "goto ";11;" ";21
print #1, "backcolor white;\Absorbtia (u.a.)"
print #1, "down; goto ";11;" ";400
print #1, "goto ";611;" ";400
print #1, "up; goto ";411;" ";421
print #1, "\ Viteza (mm/s)"
print #1, "color red;backcolor darkred"
for I=1 to NPCT
print #1, "down; place "+str$(15+I*5)+" "+str$(300 -(D(I) - MI)/SCL*250+10)
print #1, "circlefilled 3 "
next I
590 ALM = 1'10'100 ' Parametrul Levenberg-Marquart.
600 V(1)=565:V(2)=11 :V(3)=94 :V(4)= 22:V(5)= 63 ' Valorile initiale ale parametrilor.
630 EPS = SQR(REPS)' 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 * 1.01 ' 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
1580 FOR I = 1 TO NPAR
1590 w(I) = V(I) + AL * X(I)
1640 NEXT I
1650 PRINT "alfa = " ; AL; " Par.L-M = "; ALM;" itn = "; ITN; " 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 OR F0 < SZ THEN 2060
1740 GOSUB 3270'plotare fit
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 ' oricƒt 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 VT(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' ATENTIE ---------!!!!!!!
1920 '
1930 '========== ACTUALIZAREA PARAMETRILOR OPTIMI ===============
1940 '
1950 FOR I = 1 TO NPAR
1960 V(I) = VT(I)
1970 NEXT I
1980 KO = 0
1990 INEF = 0
2000 '=================== RELUAREA CAUTARII ====================
2010 GOTO 630
2020 '
2030 '=================== SFIRSITUL CAUTARII ===================
2040 '
2050 GOSUB 3200 ' Actualizarea param. optimi pt. raportul final.
2060 GOSUB 2940'CALCULUL INVERSEI
2110 PRINT
2120 PRINT "TERMINAT"
2130 IF NQ <> NPAR THEN 2170
2140 A$ = "TERMINATIE NORMALA."
2160 PRINT A$;" S-AU OBTINUT "; NSIG;" CIFRE SEMNIF. IN "; ITN;" ITERATII":GOTO 2230
2170 A$ = "TERMINATIE ANORMALA"':print A$
2180 IF INEF = 0 THEN 2210 else print A$:goto 2230
2190 PRINT A$ + " datorata cautarii ineficiente"
2200 GOTO 2230
2210 IF IFN < MAXFN THEN 2230
2220 PRINT A$ + ". S-a depasit nr. maxim de eval. ale fct."
2230 FOR I = 1 TO NPAR
2250 PRINT I;" "; w(I); " +/- ";SQR(ABS(AINV(I, I)))
2260 NEXT I
2300 'GOSUB 3270
rem print #1, "down"
print #1, "color blue; size 2"
xx=300 - (DC(1) - MI) / SCL * 250+10
print #1, "up; place "+str$(15)+" "+str$(xx)
for tt=1 to NPCT
xx=300 -(DC(tt) - MI)/SCL*250+10
print #1, "down; goto "+str$(tt*5+15)+" "+str$(xx)
next tt
print #1, "color red; size 1"
for I=1 to NPCT
print #1, "place "+str$(15+I*5)+" "+str$(300 -(D(I) - MI)/SCL*250+10)
print #1, "circle 3 "
next I
print #1, "flush"
input r$
close #1
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) = -1*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
'print DC(I)
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 print #1, "color LightGray; size 1"
xx=300 - (DC(1) - MI) / SCL * 250+10
print #1, "up; place "+str$(15)+" "+str$(xx)
for tt=1 to NPCT
xx=300 -(DC(tt) - MI)/SCL*250+10
print #1, "down; goto "+str$(tt*5+15)+" "+str$(xx)
next tt
RETURN
|
|
|
