C #MINIY V2B C LAST UPDATE 21-MAY-73. C SUBROUTINE MINIY(AFUN,PAM,ND,STP,ALF,MD,F1,ICALL,FLAM,WRK) DIMENSION PAM(ND), STP(ND), ALF(MD,MD), WRK(ND,3) C C *PURPOSE. C TO SEARCH A PARAMETER SPACE OF SEVERAL DIMENSIONS FOR THE C MINIMUM OF A FUNCTION SUPPLIED BY THE USER. C C *PARAMETERS. C AFUN - EXTERNAL NAME OF FUNCTION SUPPLIED BY USER TO BE C MINIMIZED. IT IS CALLED AS "AFUN(PAM)" C PAM - REAL ARRAY OF PARAMETERS. IS SET INITIALLY BY THE USER C TO AN ESTIMATE OF THE POSITION OF THE MINIMUM AND IS C CHANGED BY "MINIY" C ND - THE NUMBER OF PARAMETERS. THE FIRST ND PARAMETERS OF C THE ARRAY "PAM" WILL BE VARIED IN THE SEARCH. C STP - AN ARRAY OF STEPSIZES FOR CORRESPONDING PARAMETERS. C THE STEPSIZES ARE USED TO ESTIMATE FIRST AND SECOND C DERIVATIVES OF THE FUNCTION "AFUN". AS AN OPTION C "MINIY" WILL ADJUST THE STEPS AS IT GOES. C ALF - A REAL MATRIX OF DIMENSION "MD"X"MD" WHICH IS RETURNED C WITH ERROR ESTIMATES. THE DIAGONAL ELEMENTS C ALF(I,I) RETURN THE ESTIMATED STANDARD DEVIATION IN C THE I'TH PARAMETER PAM(I), AND ALF(I,J) RETURNS THE C FRACTIONAL CORRELATION COEFFICIENT BETWEEN THE I'TH C AND J'TH PARAMETERS. THE STATISTICAL SIGNIFICANCE OF C THESE VALUES DEPENDS ON THE PROBLEM BEING SIMILAR TO C A LEAST SQUARES ONE AT LEAST IN THE VICINITY OF THE C MINIMUM. C MD - THE DIMENSION FOR "ALF". THIS IS INCLUDED AS A SEPARATE C PARAMETER FROM "ND" AS "ND" IS OFTEN VARIED. C F1 - THE RETURNED SMALLEST VALUE OF THE FUNCTION FOUND IN C THIS ITERATION. C ICALL - MUST BE SET TO ZERO BY THE USER AT THE START OF EACH C MINIMIZATION. IT IS INCREMENTED INTERNALLY TO COUNT THE C NUMBER OF ITERATIONS. C FLAM - IS INITIALLY SET BY THE USER TO INDICATE THE TYPE OF C MINIMIZATION STRATEGY AS: C <0 STEP SIZES ARE KEPT FIXED C =0 STEPS FIXED AND NO MARQUARDT CORRECTION C >0 STEPS VARIED AND MARQUARDT CORRECTION USED C WHEN THE STEPS ARE VARIED BY MINIY THEY ARE SET TO WHAT C MINIY THINKS IS REASONABLE FOR A LEAST SQUARES PROBLEM C MULTIPLIED BY THE USERS INITIAL "FLAM" VALUE. C WRK - A WORKING ARRAY OF AT LEAST 3*ND REAL ELEMENTS. C C INPUT PARAMETERS: C AFUN, PAM, ND, STP, MD, ICALL, FLAM C OUTPUT PARAMETERS: C PAM, STP, ALF, F1, ICALL, FLAM, WRK C C *METHOD. C THIS ROUTINE CARRIES OUT ONE ITERATION OF A GENERAL PURPOSE C MINIMIZATION STRATEGY. IT IS THE USERS RESPONSIBILITY TO C CHECK FOR CONVERGENCE AND EXCEPT IN SPECIAL CIRCUMSTANCES C IT WILL BE NECESSARY TO CALL "MINIY" SEVERAL TIMES TO GET A GOOD C ESTIMATE OF THE TRUE MINIMUM. AS A VERY ROUGH GUIDE ABOUT C (ND/2) ITERATIONS MAY BE TAKEN. CHECKS FOR CONVERGENCE SHOULD C CHECK THE MOVEMENT OF THE PARAMETERS AS WELL AS THE CHANGE OF C THE FUNCTION VALUE. C THIS IS A FULL NEWTON METHOD WITH A DIRECTION CHOSEN C BETWEEN THE GRADIENT AND THE DIRECTION TO THE BOTTOM OF C THE LOCAL QUADRIC SURFACE AS DONE FIRST BY MARQUARDT. THE C METHOD OF INTERPOLATING BETWEEN THESE DIRECTIONS IS UNIQUE IN C THAT THE VALUE LAMBDA USED BY MARQUARDT IS ESTIMATED FROM THE C SCALED MAGNITUDE OF THE GRADIENT SO THAT WHEN FAR FROM THE C MINIMUM LAMBDA IS AUTOMATICALLY LARGE AND THE SEARCH IS CLOSE TO C THE GRADIENT AND WHEN CLOSE TO THE MINIMUM THE SEARCH PROCEEDS TO C THE BOTTOM OF THE QUADRIC SURFACE. C AUTOMATIC SCALING IS DONE TO BRING THE SECOND DERIVATIVES C IN EACH DIRECTION TO UNITY. C THE ROUTINE WILL REDUCE TO A GRADIENT SEARCH IF THE LOCAL C QUADRIC SURFACE APPEARS TO GIVE A MAXIMUM INSTEAD OF A MINIMUM C AND CAN BE REDUCED TO THE PURE NEWTON METHOD BY SETTING "FLAM"=0 C INITIALLY. THERE ARE CERTAIN SAFEGUARDS SUCH AS CHECKING THAT C THE PREDICTED POINT GIVES A LOWER VALUE THAN THE INITIAL ONE ETC. C STEPS ARE TAKEN IN THE PARAMETER SPACE TO DETERMINE FIRST AND C SECOND DERIVATIVES. IF REQUIRED A MARQUARDT CORRECTION VALUE IS C OBTAINED FROM THE SCALED GRADIENT AND THE CURVATURE MATRIX IS C MODIFIED. IN ANY CASE THE CURVATURE MATRIX IS SCALED C AND INVERTED. STEPS ARE TAKEN ALONG THE INDICATED DIRECTION AND C THE BEST 3 POINTS USED FOR A PARABOLIC INTERPOLATION TO THE C MINIMUM. C C *ACCURACY. C DEPENDS ENTIRELY ON THE PROBLEM. C C *RESTRICTIONS. C THE STEPS ARE NOT ALLOWED TO BE ZERO. C C *ERROR CONDITIONS. C IF THE FUNCTION SURFACE IS NOT SUITABLE FOR THE ROUTINE SOME C SYMPTOMS MAY APPEAR. PARTICULARLY THE CORRELATION COEFFICIENTS IN C ALF MAY HAVE MAGNITUDE GREATER THAN 1. C C *NON STANDARD ROUTINES CALLED. C INVTB (MRMLIB) C C *TYPICAL TIMES. C ABOUT 2*ND + 0.5*ND*(ND-1) POINTS ARE CALCULATED PER PASS C THIS DETERMINES THE TIME TAKEN. C C *ORIGIN. M.R.MANNING. C THIS ROUTINE HAS BEEN DEVELOPED FROM EARLIER VERSIONS C "MINIZE" AND "MINIZ" AND NOW REPLACES THEM. THE ONLY DIFFERENCE C FROM "MINIZ" IS THAT THE INVERSION OF THE MODIFIED CURVATURE C MATRIX IS NOW DONE USING A COLUMN PIVOTING CROUT INVERSION, INSTEAD C OF A NO PIVOTING INVERSION. C C *COMMENTS. C OTHER ROUTINES ARE BETTER FOR LINEAR LEAST SQUARES C FITTING BUT "MINIY" CAN SOLVE SUCH A PROBLEM AND "FLAM" SHOULD C BE SET ZERO IN THIS CASE SO THE SOLUTION IS ACHIEVED IN THE C FIRST ITERATION. C IT SHOULD BE NOTED CAREFULLY THAT THE "ALF" MATRIX SHOULD C NOT BE REGARDED AS AN INDICATION OF THE ERRORS UNLESS "PAM" C IS CLOSE TO THE TRUE MINIMUM! C C #END. C LOGICAL QUAD, FIXSTP INTEGER RETURN C DATA NB, NT, EPS/ 2, 1, 1.E-4/ C GRADIENT ARRAY BET IS IN WRK(*,NB) AND TEMPORARY ARRAY TEM C IN WRK(*,NT) C C #START. IF (ICALL) 10,10,11 C C FIRST TIME ONLY CODE. C 10 F1=AFUN(PAM) QUAD= FLAM.EQ.0. C QUAD TRUE IF SURFACE TO BE TAKEN AS QUADRIC FIXSTP= FLAM.LE.0. C FIXSTP TRUE IS STEPS TO BE LEFT FIXED. FAC= FLAM C STEPS ARE SET BY MINIY AND THEN MULTIPLIED BY FAC. FN= ND C C CALCULATE CURVATURE AND GRADIENT ARRAYS C UPPER TRIANGLE OF ALF ONLY AT THIS STAGE. C 11 ICALL= ICALL + 1 FL= 0. DO 30 J= 1, ND PJ= PAM(J) SJ= ABS( STP(J) ) PAM(J)= PJ + SJ F2= AFUN(PAM) DO 20 K= 1, ND SK= STP(K) IF (K-J) 15, 20, 16 15 ALF(K,J)= (ALF(K,J)-F2)/SJ/SK GO TO 20 16 PK= PAM(K) PAM(K)= PK + SK F3=AFUN(PAM) ALF(J,K)=F1-F2+F3 PAM(K)=PK 20 CONTINUE PAM(J)=PJ-SJ F3=AFUN(PAM) PAM(J)=PJ C SCALE FACTOR FOR J'TH DIRECTION Y=SQRT(ABS(F2-F1-F1+F3)) IF(Y.LT.EPS*SJ) Y=EPS*SJ WRK(J,NT)=SJ/Y Z=(F3-F2)/(Y+Y) WRK(J,NB)=Z 30 FL=FL+Z*Z C LENGTH OF NORMALIZED GRADIENT AND LAMBDA VALUE Y=1. IF(Y.LT.F1) Y=F1 FLAM=SQRT(FL/FN/Y)/10. FL=FLAM*25. IF(QUAD) FLAM=0. 42 IF(FLAM-4.) 45,45,70 C C INVERT CURVATURE MATRIX C 45 DO 49 J=1,ND J1=J+1 IF(J1.GT.ND) GO TO 49 Y=WRK(J,NT) DO 48 K=J1,ND ALF(J,K)=ALF(J,K)*Y*WRK(K,NT) 48 ALF(K,J)=ALF(J,K) 49 ALF(J,J)=1. + FLAM C CALL INVTB(ALF,ND,MD,DET, WRK(1,3) ) C ** INVTB IS A LIBRARY MATRIX INVERSION WITH DETERMINANT DET ** IF(DET) 70,70,55 C C CALCULATE SEARCH DIRECTION AND STEP IN TEM. C AND SET 95% CONFIDENCE LIMITS AND CORRELATION COEFFICIENTS C IN THE ALF ARRAY. 55 DO 65 J=1,ND Z=0. Y=SQRT(ABS(ALF(J,J))) DO 60 K=1,ND Z=Z+ALF(J,K)*WRK(K,NB) IF(K-J) 57,60,58 57 ALF(J,K)=ALF(K,J)/Y ALF(K,J)=ALF(J,K) GO TO 60 58 ALF(J,K)=ALF(J,K)/Y 60 CONTINUE ALF(J,J)=Y*WRK(J,NT) 65 WRK(J,NT)=Z*WRK(J,NT) GO TO 80 C C IN CASE OF GRADIENT SEARCH ONLY C 70 DO 75 J=1,ND DO 74 K=1,ND 74 ALF(J,K)=0. ALF(J,J)=WRK(J,NT) 75 WRK(J,NT)=WRK(J,NB)*WRK(J,NT)/FL C C TAKE STEPS ALONG NEW DIRECTION. C 80 BIAS=-1. NSTP=0 FL=1. SGN=1. 81 DEL=SGN ASSIGN 504 TO RETURN GO TO 500 C 504 F2=AFUN(PAM) C IF QUADRIC SURFACE THIS IS THE BOTTOM SO GET OUT IF(QUAD) GO TO 503 C IF AFUN INCREASED REDUCE STEP SIZE BY HALF AND CYCLE BACK TO HERE. IF(F1-F2) 85,90,90 85 BIAS=0. SGN=-1. F3=F2 DO 86 J=1,ND 86 WRK(J,NT)=0.5*WRK(J,NT) FL=0.5*FL NSTP=NSTP+1 IF(NSTP-5) 81,87,87 C AFTER 5 REDUCTION STEPS GET OUT. 87 DEL=-2. ASSIGN 505 TO RETURN GO TO 500 505 GO TO 106 C 90 IF(BIAS) 91,100,91 C IF AFUN DECREASED LEAVE STEP SIZE AND KEEP GOING C TILL AN INCREASE FOUND. 91 ASSIGN 501 TO RETURN DEL=1. GO TO 500 C 501 F3=AFUN(PAM) IF(F3-F2) 93,100,100 93 F1=F2 F2=F3 FL=FL+1. NSTP=NSTP+1 IF(NSTP-10) 91,100,100 C C OK FIT PARABOLA TO LAST THREE POINTS ON SEARCH DIRECTION. 100 DEL=0.5*(F1-F3)/(F1-F2-F2+F3) + BIAS ASSIGN 502 TO RETURN GO TO 500 C 502 F3=AFUN(PAM) IF(F3-F2) 105,105,103 103 DEL=BIAS-DEL ASSIGN 503 TO RETURN GO TO 500 C 503 F1=F2 GO TO 106 105 F1=F3 106 IF(FIXSTP) RETURN C FL IS NOW APPROX THE NUMBER OF STEPS TAKEN ALONG SEARCH DIRECTION. C C STEPS ARE: 38% CONFIDENCE LIMITS (I.E. 0.5*VARIANCE) C BY 1/2 NO OF STEPS TAKEN ABOVE. C BY USER FACTOR FAC, SET BY ORIGINAL FLAM. Y=FAC*FL/4. DO 120 J=1,ND Z=Y*ALF(J,J) C BUILD FOLLOWING INERTIA INTO STEP SIZES FOR SAFETY. IF(Z.LT.STP(J)/5.) Z=STP(J)/5. IF(Z.GT.STP(J)*5.) Z=STP(J)*5. 120 STP(J)=Z RETURN C C GENERAL ROUTINE USED FIVE TIMES TO STEP PAM ARRAY 500 DO 510 J=1,ND 510 PAM(J)=PAM(J)+DEL*WRK(J,NT) GO TO RETURN C END