Attribute VB_Name = "Module1" '****************************************************************** '* MULTI-DIMENSIONAL CURVE FITTING BY SIMPLEX METHOD * '* -------------------------------------------------------------- * '* Reference: * '* * '* J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, * '* OKLAHOMA STATE UNIVERSITY * '* -------------------------------------------------------------- * '* Explanations: * '* * '* THIS PROGRAM FITS THE MODEL: * '* * '* A * '* Z = ---------------- * '* 1 + BX + CY * '* * '* WRITTEN AS: * '* * '* FIT(J) = X(1)/(1.0 + X(2)*T(J,1) + X(3)*T(J,2)) * '* * '* TO DATA POINTS (T(J,1),T(J,2),Y(J)). * '* * '* T(J,1) AND T(J,2) ARE THE TWO INDEPENDENT VARIABLES X AND Y * '* IN DATA SPACE. Y(J) IS THE MEASURED DEPENDENT VARIABLE Z AND * '* FIT(J) IS THE VALUE TO BE FITTED TO Y(J) BY MINIMIZING THE * '* WEIGHTED SUM OF SQUARES, * '* * '* J=NPTS * '* * '* PHI = SUM ((FIT(J)-Y(J))/YSIG(J))**2 * '* * '* J=1 * '* * '* THE VALUE OF PHI IS STORED IN THE VARIABLE FOBJ. * '* THE MODEL Z=F(X,Y) IS IMPLEMENTED IN FUNCTION FUNK. * '* THE RESULTS ARE BEST COEFFICIENTS A, B, C, STORED IN X(1), * '* X(2), X(3) TO FIT THE DATA PROVIDED IN FILE SIMPTEST.DAT. * '* -------------------------------------------------------------- * '* SAMPLE RUN: * '* File Simptest.dat contains: * '* 12 * '* 0.0 1.0 16.6 0.2 * '* 1.0 0.0 12.4 0.2 * '* 1.0 1.0 8.2 0.2 * '* 0.0 2.0 10.1 0.2 * '* 2.0 0.0 7.3 0.2 * '* 2.0 2.0 4.7 0.2 * '* 1.0 3.0 5.1 0.2 * '* 3.0 1.0 4.3 0.2 * '* 3.0 3.0 3.0 0.2 * '* 5.0 1.0 2.5 0.2 * '* 1.0 5.0 3.4 0.2 * '* 5.0 5.0 2.1 0.2 * '* * '* Output file Simpout.lst contains at the end: * '* --/-- * '* * '* 543 FUNCTION COMPUTATIONS. * '* * '* FINAL VALUE OF FOBJ = 6.50403378306642 * '* * '* FINAL VALUES OF X(J): * '* 48.2028883483963 2.87632824106749 1.90303605417914 * '* * '* * '* + 3.6498 * '* X( 1 ) = 48.20288 ( STDEVS = 1 ) * '* - 3.04467 * '* * '* * '* + 0.26377 * '* X( 2 ) = 2.87632 ( STDEVS = 1 ) * '* - 0.22392 * '* * '* * '* + 0.19635 * '* X( 3 ) = 1.90303 ( STDEVS = 1 ) * '* - 0.16812 * '* * '* -------------------------------------------------------------- * '* Visual Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '****************************************************************** 'ote: To change example, adjust function FUNK to your problem and check ' NV, NPTS and size of T(300,3). DefDbl A-H, O-Z DefInt I-N Option Base 1 ' global variables Public ERGSAV(20) Public FIDINT(2) Public X(20), XMAX(20), XMIN(20), DELTX(20), DELMIN(20) Public MASK(20) Public ERFRAC, ERGUES, FOBJ, FSAVE, STDEVS, TOLFID Public JXFID, MAXIND, MATRX, NV, NTRACE, NTRACF, NTRSET Public K, NFMAX, NFLAT, JVARY, NXTRA, KFLAG, NOREP, KERFL Public FLAMBD, FNU, RELDIF, RELMIN Public METHD, KALCP, KORDIF, MAXIT, LEQU, MAXSUB, MAXUPD, NPTS ' common variables for SMPLX, SIBEG... Public FZ(20) Public HUGE, ALPHA, BETA, GAMMA Public MAXPT, NVPLUS, NSSW, NF Public Y(300), YSIG(300) Public T(300, 2) Public Z(20, 21) Sub FUNK() '------------------------------------------------------------- ' J. P. CHANDLER, COMPUTER SCIENCE DEPT., ' OKLAHOMA STATE UNIVERSITY ' ' THIS VERSION OF SUBROUTINE FUNK COMPUTES THE FITTED VALUES ' AND THE WEIGHTED SUM OF SQUARES FOR THE PROBLEM OF FITTING ' THE MODEL ' ' FIT(J) = X(1) / (1# + X(2) * T(J, 1) + X(3) * T(J, 2)) ' ' TO DATA POINTS (T(J,1),T(J,2),Y(J)). ' '------------------------------------------------------------- Dim FITM(300) ' FOBJ WILL CONTAIN THE WEIGHTED SUM OF SQUARES, PHI. ' SIMPLEX REQUIRES THAT FUNK COMPUTE FOBJ. ' Here FOBJ is a global variable. FOBJ = 0# For J = 1 To NPTS FITM(J) = X(1) / (1# + X(2) * T(J, 1) + X(3) * T(J, 2)) FOBJ = FOBJ + ((FITM(J) - Y(J)) / YSIG(J)) ^ 2 Next J End Sub Function T1(ARG) TMP = ARG - FSAVE If TMP <= 0# Then TMP = 0# T1 = Abs(Sqr(TMP) - STDEVS) End Function Sub FIDO(JXFID, STDEVS, TOLFID, ERGUES, MAXIND, NTRSET, NTRACF, FIDINT()) '----------------------------------------------------------------- ' FIDO 4.5 DECEMBER 1991 ' A.N.S.I. STANDARD FORTRAN 77 ' ' J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, ' OKLAHOMA UNIVERSITY, STILLWATER, OKLAHOMA 74078 ' ' SUBROUTINE FIDO IS USED IN CONJUNCTION WITH SIMPLEX ' TO COMPUTE CONFIDENCE HALF-INTERVALS ("ERRORS") FOR ' THE PARAMETERS IN A FITTING PROBLEM, USING THE METHOD OF ' SUPPORT PLANES. ' THE QUANTITY BEING MINIMIZED (FOBJ) MUST EITHER BE ' CHI-SQUARE OR IT MUST BE TWICE THE NEGATIVE OF THE NATURAL ' LOGARITHM OF THE LIKELIHOOD FUNCTION. ' ' --------------------------------------------------------------- ' ' INPUT QUANTITIES..... JXFID,STDEVS,TOLFID,ERGUES, ' MAXIND,NTRSET,NTRACF,X(*), ' ' OUTPUT QUANTITY...... FIDINT( ) ' ' JXFID -- THE INDEX OF THE PARAMETER, X(JXFID), FOR ' WHICH CONFIDENCE HALF-INTERVALS ARE TO BE ' COMPUTED ' ' STDEVS -- THE NUMBER OF STANDARD DEVIATIONS TO WHICH ' THE HALF-INTERVALS ARE TO CORRESPOND ' ' TOLFID -- CONVERGENCE TOLERANCE (USE TOLFID=0.05) ' ' ERGUES -- ROUGH ESTIMATE OF THE LENGTH OF THE ' HALF-INTERVALS ' ' MAXIND -- =1 TO COMPUTE THE HALF-INTERVAL WITH THE ' SAME SIGN AS ERGUES, ' =2 TO COMPUTE BOTH HALF-INTERVALS ' ' NTRSET -- VALUE OF NTRACE TO BE USED IN THE ' FITTING SUBROUTINE (STEPIT OR MARQ, ' ETC.) WHEN CALLED BY FIDO (Not used here). ' ' NTRACF -- A SWITCH THAT CONTROLS PRINTING... ' ' =-1 TO PRINT NOTHING IN FIDO EXCEPT ANY ' WARNING OR ERROR MESSAGES, ' ' = 0 FOR INITIAL AND SUMMARY PRINTOUT, ' ' =+1 TO PRINT EVERY FIDO ITERATION ' ' X( ) -- THE POSITION OF THE MINIMUM OF FOBJ ' ' ' FIDINT(J) -- RETURN THE COMPUTED LENGTHS OF THE ' CONFIDENCE HALF-INTERVALS. ' ' FOR A LEAST SQUARES PROBLEM, FOBJ IS EQUAL TO CHI-SQUARE, ' THE WEIGHTED SUM OF SQUARES... ' ' NPTS ' FOBJ = Sum((FIT(JPT) - YDATA(JPT)) / YSIGMA(JPT))^2 ' JPT = 1 ' ' THE STANDARD ERRORS YSIGMA( ) MUST BE CORRECTLY SCALED. ' IF THE YSIGMA( ) ARE NOT KNow1N ACCURATELY, NORMALIZE THEM ' SO THAT THE VALUE OF FOBJ AT THE MINIMUM IS EQUAL TO THE ' NUMBER OF DEGREES OF FREEDOM, ' N.D.F. = (NO. DATA POINTS) - (NO. ADJUSTABLE PARAMETERS) ' ' FIDO IS ESSENTIALLY A ROOT FINDING ROUTINE. IT USES INVERSE ' QUADRATIC INTERPOLATION OR EXTRAPOLATION, INVERSE LINEAR ' INTERPOLATION, AND BISECTION (INTERVAL HALVING), AS ' SUCCESSIVELY MORE DIFFICULTIES ARE ENCOUNTERED IN FINDING ' THE ROOT. ' ' PRIMARY REFERENCES.... ' ' W. T. EADIE ET AL., "STATISTICAL METHODS IN EXPERIMENTAL ' EXPERIMENTAL PHYSICS" (AMERICAN ELSEVIER, 1971), ' CHAPTER 9 ' 'P.R.BEVINGTON , "DATA REDUCTION AND ERROR ANALYSIS IN ' THE PHYSICAL SCIENCES" (MCGRAW-HILL, 1969), ' PAGES 242-245 ' ' OTHER REFERENCES.... ' ' H. SCHEFFE, "THE ANALYSIS OF VARIANCE" (WILEY, 1959) ' ' H. STONE, DISCUSSION ON PAPER BY E. M. L. BEALE, ' J. ROY. STATIS. SOC. B., V. 22 (1960), P. 41, PP. 84-5 ' ' G. E. P. BOX AND G. A. COUTIE, PROC. INST. ELEC. ENGRS. ' V. 103, PART B, SUPPL. NO. 1 (1956). ' ' G. W. BOOTH, G. E. P. BOX, M. E. MULLER, AND ' T.I.PETERSON , "FORECASTING BY GENERALIZED REGRESSION ' METHODS NON-LINEAR ESTIMATION" (PRINCETON/IBM), ' FEBRUARY 1959, IBM MANUAL ' ' D.W.MARQUARDT , "LEAST-SQUARES ESTIMATION OF NONLINEAR ' PARAMETERS", SHARE DISTRIBUTION 3094 ' ' D. W. MARQUARDT, R. G. BENNETT, AND E. J. BURRELL, "LEAST ' SQUARES ANALYSIS OF ELECTRON PARAMAGNETI' RESONANCE ' SPECTRA", J. OF MOLECULAR SPECTROSCOPY 7 (1961) 269 ' '------------------------------------------------------------------------- 'Labels: 20,40,70,100,120,140,150,170,190, 210,240,250,270,290,310,330,350, ' 370,380,400,430 (inherited from Fortran 77). Dim XSAVE(20), XKA(20), XKB(20), XBEST(20) ' * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ' ' SET SOME PARAMETERS. ' ' RSMALL IS USED IN SETTING A PHONY VALUE IF A VALUE OF FOBJ ' IS FOUND THAT IS LESS THAN THE "GLOBAL MINIMUM" VALUE. ' RSMALL = 0.0001 ' CONVRG = SATISFACTORY RATE OF CONVERGENCE CONVRG = 0.5 ' BRACK = FACTOR FOR BRACKETTING SHOTS BRACK = 2# ' MAXTRY = MAXIMUM NUMBER OF ITERATIONS MAXTRY = 20 ' FACMAX = MAXIMUM FACTOR FOR QUADRATIC INTERPOLATION FACMAX = 4# ' FRMIN = MINIMUM FRACTION FOR LINEAR INTERPOLATION FRMIN = 0.1 RZERO = 0# RHALF = 0.5 UNITR = 1# ' SAVE SOME INPUT QUANTITIES, AND INITIALIZE EVERYTHING NTRSAV = NTRACE NTRACE = NTRSET MATSAV = MATRX MATRX = 0 MSKSAV = MASK(JXFID) MASK(JXFID) = 1 NACTIV = 0 For J = 1 To NV XBEST(J) = X(J) XSAVE(J) = X(J) If MASK(J) = 0 Then NACTIV = NACTIV + 1 Next J Call FUNK FSAVE = FOBJ FBEST = FOBJ FPSGSQ = FOBJ + STDEVS ^ 2 X(JXFID) = XSAVE(JXFID) + STDEVS * ERGUES ' LOOP OVER THE FIDINT(JJ) JJ = 1 20 If ERGUES <> RZERO Then GoTo 40 Print #2, " " Print #2, " ERROR IN INPUT TO FIDO... ERGUES IS EQUAL TO ZERO." End 40 SGN1 = UNITR If ERGUES < RZERO Then SGN1 = -UNITR If NTRACF < 0 Then GoTo 70 Print #2, " " Print #2, " SUBROUTINE FIDO:" Print #2, " JXFID ="; JXFID; " STDEVS = "; STDEVS; " TOLFID = "; TOLFID Print #2, " GLOBAL MINIMUM OF FOBJ ="; FSAVE; " IS AT X(JXFID) = "; XSAVE(JXFID) Print #2, " FIRST GUESS FOR LENGTH OF HALF-INTERVAL = "; ERGUES Print #2, " " Print #2, " NV ="; NV; " NACTIV ="; NACTIV 70 KLOSB = 0 MLIN = 0 For K = 1 To NV XKA(K) = XSAVE(K) Next K FB = FSAVE FKA = FSAVE ' BEGIN THE ITERATION LOOP FOR LOCATING THE PLANE ' PERPENDICULAR TO THE X(JXIFD) AXIS IN WHICH THE MINIMUM ' VALUE OF FOBJ IS EQUAL TO FPSGSQ = FSAVE + STDEVS^2, ' WHERE FSAVE IS THE VALUE OF FOBJ AT THE GLOBAL MINIMUM. For ITER = 1 To MAXTRY FC = FB TEMP = FSAVE If NACTIV = 0 Then Call FUNK Else Call SMPLX End If FSAVE = TEMP 'restore current FSAVE after SMPLX If FOBJ >= FBEST Then GoTo 100 FBEST = FOBJ For K = 1 To NV XBEST(K) = X(K) Next K 100 DX = X(JXFID) - XSAVE(JXFID) DC = FOBJ - FPSGSQ If NTRACF >= 1 Then Print #2, " " Print #2, " ITERATION "; ITER; " X("; JXFID; ") = "; X(JXFID) Print #2, " DISTANCE FROM GLOBAL MINIMUM = "; DX Print #2, " MINIMUM FOBJ IN THIS PLANE = "; FOBJ Print #2, " VALUE SOUGHT = "; FPSGSQ; " DIFFERENCE = "; DC End If If FOBJ = FSAVE Then GoTo 140 If FOBJ < FSAVE Then GoTo 120 ' TEST FOR CONVERGENCE ' IF(ABS(FOBJ-FPSGSQ).LE.TOL) ............................ DIFF = FOBJ - FPSGSQ If DIFF < RZERO Then DIFF = -DIFF If DIFF <= TOLFID Then GoTo 330 If FOBJ >= FPSGSQ Then GoTo 170 GoTo 150 120 DIFF = FSAVE - FOBJ Print #2, " " Print #2, " ***** FOBJ LESS THAN VALUE AT INPUT POINT BY "; DIFF Print #2, " X(J) AT THIS POINT: "; For K = 1 To NV Print #2, " "; X(K); Next K Print #2, " " 140 FOBJ = FSAVE + RSMALL * (FPSGSQ - FSAVE) ' CHECK FOR TIGHTER BRACKETTING OF FPSGSQ FROM BELOW 150 If (X(JXFID) - XKA(JXFID)) * SGN1 <= RZERO Then GoTo 190 For K = 1 To NV XKA(K) = X(K) Next K FKA = FOBJ GoTo 190 ' CHECK FOR BRACKETTING OF FPSGSQ FROM ABOVE 170 If KLOSB = 1 Then If (X(JXFID) - XKB(JXFID)) * SGN1 >= RZERO Then GoTo 190 End If KLOSB = 1 For K = 1 To NV XKB(K) = X(K) Next K FKB = FOBJ 190: If T1(FOBJ) < T1(FC) Then FB = FOBJ ' CHECK THE RATE OF CONVERGENCE. IF IT IS SATISFACTORY, AND ' IF LINEAR INTERPOLATION HAS BEEN USED, USE IT AGAIN. If (ITER >= 2) And (T1(FOBJ) > CONVRG * T1(FC)) Then GoTo 210 If MLIN = 1 Then GoTo 250 ' USE INVERSE QUADRATIC INTERPOLATION TEMP = FOBJ - FSAVE If TEMP <= 0# Then FAC = 1# Else FAC = STDEVS / Sqr(TEMP) End If ' FAC=AMIN1(FACMAX,AMAX1(FAC,1./FACMAX)) ................... If FAC < UNITR / FACMAX Then FAC = UNITR / FACMAX If FAC > FACMAX Then FAC = FACMAX XX = XSAVE(JXFID) + FAC * (X(JXFID) - XSAVE(JXFID)) ' CHECK THAT THE PROPOSED POINT IS INSIDE THE BRACKETTED INTERVAL If (XX - XKA(JXFID)) * SGN1 <= RZERO Then GoTo 210 If KLOSB = 1 Then If (XX - XKB(JXFID)) * SGN1 >= RZERO Then GoTo 240 End If For K = 1 To NV If (K = JXFID) Or (MASK(K) = 0) Then X(K) = XSAVE(K) + FAC * (X(K) - XSAVE(K)) End If Next K GoTo 310 210 If KLOSB = 1 Then GoTo 240 ' CONVERGENCE IS POOR, AND FPSGSQ HAS NOT YET BEEN BRACKETTED. ' TRY TO BRACKET IT. FRAC = BRACK If FOBJ >= FPSGSQ Then FRAC = UNITR / BRACK For K = 1 To NV If (K = JXFID) Or (MASK(K) = 0) Then X(K) = XSAVE(K) + FRAC * (X(K) - XSAVE(K)) Next K If NTRACF >= 1 Then Print #2, " " Print #2, " BRACKETTING SHOT...." End If GoTo 310 240 If MLIN = 1 Then GoTo 270 ' TRY LINEAR INTERPOLATION BETWEEN THE TWO BRACKETTING POINTS 250 MLIN = 1 FRAC = (FPSGSQ - FKA) / (FKB - FKA) ' FRAC=AMAX1(FRMIN,AMIN1(1.0-FRMIN,FRAC)) ............. If FRAC > UNITR - FRMIN Then FRAC = UNITR - FRMIN If FRAC < FRMIN Then FRAC = FRMIN If NTRACF >= 1 Then Print #2, " " Print #2, " LINEAR INTERPOLATION...." End If GoTo 290 ' CONVERGENCE IS POOR, AND LINEAR INTERPOLATION HAS BEEN USED. ' BISECT THE BRACKETTED INTERVAL. 270 FRAC = RHALF If NTRACF >= 1 Then Print #2, " " Print #2, " INTERVAL BISECTED...." End If 290 For K = 1 To NV If (K = JXFID) Or (MASK(K) = 0) Then X(K) = XKA(K) + FRAC * (XKB(K) - XKA(K)) End If Next K 310 Next ITER ' END OF ITERATION LOOP. THE ITERATION FAILED TO CONVERGE Print #2, " " Print #2, " CONVERGENCE FAILURE IN SUBROUTINE FIDO." FIDINT(JJ) = RZERO GoTo 350 ' CONVERGENCE ACHIEVED ... A SATISFACTORY PLANE HAS BEEN LOCATED 330 FIDINT(JJ) = X(JXFID) - XSAVE(JXFID) If NTRACF >= 0 Then Print #2, " " Print #2, " THE LENGTH OF THE CONFIDENCE HALF-INTERVAL FOR" Print #2, " X("; JXFID; ") IS "; FIDINT(JJ); " (STDEVS = "; STDEVS; ")" End If 350 JJ = JJ + 1 If JJ > MAXIND Then GoTo 370 ' REFLECT X THROUGH XSAVE, AND SEARCH FOR THE OTHER HALF-INTERVAL For K = 1 To NV If (K = JXFID) Or (MASK(K) = 0) Then X(K) = XSAVE(K) + (XSAVE(K) - X(K)) Next K ERGUES = X(JXFID) - XSAVE(JXFID) GoTo 20 ' END OF THE JJ LOOP. PRINT SUMMARY RESULTS 370 If (MAXIND < 2) Or (FIDINT(1) = RZERO) Then GoTo 400 If (FIDINT(1) < RZERO) And (FIDINT(2) <= RZERO) Then GoTo 400 If (FIDINT(1) > RZERO) And (FIDINT(2) >= RZERO) Then GoTo 400 If (FIDINT(1) > RZERO) And (FIDINT(2) < RZERO) Then GoTo 380 XX = FIDINT(1) FIDINT(1) = FIDINT(2) FIDINT(2) = XX 380 XX = -FIDINT(2) ' Print summary of result from FIDO (error calculation) Print #2, " " Print #2, " + "; Int(FIDINT(1) * 100000) / 100000 Print #2, " X("; JXFID; ") = "; Int(XSAVE(JXFID) * 100000) / 100000; " ( STDEVS ="; STDEVS; ")" Print #2, " - "; Int(XX * 100000) / 100000 ' RESTORE THE SAVED ENTRY VALUES 400 For J = 1 To NV X(J) = XBEST(J) Next J Call FUNK If FOBJ >= FSAVE Then GoTo 430 Print #2, " " Print #2, " SUBROUTINE FIDO FOUND A BETTER MINIMUM." Print #2, " OLD FOBJ = "; FSAVE; " NEW FOBJ = "; FOBJ Print #2, " X(J) = "; For J = 1 To NV Print #2, " "; X(J); Next J Print #2, " " 430 MASK(JXFID) = MSKSAV NTRACE = NTRSAV MATRX = MATSAV End Sub 'FIDO() Sub SMPLX() '--------------------------------------------------------------- ' SIMPLEX 2.12 DECEMBER 1991 ' ' A.N.S.I. STANDARD FORTRAN 77 ' ' COPYRIGHT (C) 1965, 1975, 1991 J. P. CHANDLER ' (PRESENT ADDRESS ... ' COMPUTER SCIENCE DEPARTMENT, ' OKLAHOMA STATE UNIVERSITY, ' STILLWATER, OKLAHOMA 74078 ' (405)-744-5676 ) ' ' SIMPLEX FINDS LOCAL MINIMA OF A SMOOTH FUNCTION OF SEVERAL ' PARAMETERS. IT WILL ALSO HANDLE SOME NON-SMOOTH FUNCTIONS. ' THE METHOD IS OFTEN VERY SLOW IF THE NUMBER OF PARAMETERS IS ' AT ALL LARGE (GREATER THAN ABOUT SIX TO EIGHT). ' ' "A SIMPLEX METHOD FOR FUNCTION MINIMIZATION", ' J. A. NELDER AND R. MEAD, ' THE COMPUTER JOURNAL 7 (1965) 308-313 ' ' FOR APPLICATIONS, SEE ' ' "'DIRECT SEARCH' SOLUTION OF ' NUMERICAL AND STATISTICAL PROBLEMS", ' ROBERT HOOKE AND T. A. JEEVES, JOURNAL OF THE ' ASSOCIATION FOR COMPUTING MACHINERY 8 (1961) 212-229 ' ' "THE NELDER-MEAD SIMPLEX PROCEDURE FOR FUNCTION ' MINIMIZATION", ' D. M. OLSSON AND L. S. NELSON, ' TECHNOMETRICS 17 (1975) 45-51 ' ' ' SIMPLEX 2.9 IS AVAILABLE FROM THE ' QUANTUM CHEMISTRY PROGRAM EXCHANGE ' DEPT. OF CHEMISTRY, INDIANA UNIVERSITY ' BLOOMINGTON, INDIANA 47401 ' ' * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ' ' INPUT QUANTITIES..... NV,NTRACE,MASK(*),X(*), ' XMAX(*),XMIN(*),DELTX(*), ' DELMIN(*),NFMAX,NFLAT,NXTRA,KW ' ' OUTPUT QUANTITIES.... X(*), FOBJ, KFLAG, NOREP ' ' NV -- THE NUMBER OF PARAMETERS, X(J) ' ' NTRACE -- =0 FOR NORMAL OUTPUT, ' =+1 FOR TRACE OUTPUT, ' =+2 FOR DEBUG OUTPUT, ' =-1 FOR NO OUTPUT EXCEPT ERROR MESSAGES ' ' MATRX -- USED BY STEPIT PROPER BUT NOT BY SIMPLEX ' ' FOBJ -- THE VALUE OF THE FUNCTION TO BE MINIMIZED ' ' MASK(J) -- NONZERO IF X(J) IS TO BE HELD FIXED ' ' X(J) -- THE J-TH PARAMETER ' ' XMAX(J) -- THE UPPER LIMIT ON X(J) ' ' XMIN(J) -- THE LOWER LIMIT ON X(J) ' ' DELTX(J) -- THE INITIAL STEP SIZE FOR X(J) ' ' DELMIN(J) -- THE LOWER LIMIT (CONVERGENCE TOLERANCE) ON ' THE STEP SIZE FOR X(J) ' ' ERR(J,K) -- SCRATCH STORAGE ' ' NFMAX -- THE MAXIMUM NUMBER OF FUNCTION ' COMPUTATIONS TO BE ALLOWED ' ' NFLAT -- NONZERO IF THE SEARCH IS TO TERMINATE WHEN ' ALL TRIAL STEPS GIVE IDENTICAL FUNCTION ' VALUES. ' THE RECOMMENDED VALUE OF NFLAT IS ' USUALLY NFLAT=1 . ' ' JVARY -- USED BY STEPIT AND STP BUT NOT BY SIMPLEX ' (SIMPLEX SETS JVARY TO ZERO ' PERMANENTLY) ' ' NXTRA -- NUMBER OF EXTRA POINTS TO BE ADDED TO ' THE SIMPLEX ' (NXTRA.GT.0 CAUSES A MORE THOROUGH ' SEARCH) ' ' KFLAG -- RETURNED .GT. ZERO FOR A NORMAL EXIT, ' RETURNED .LT. ZERO FOR AN ABNORMAL EXIT ' ' NOREP -- RETURNED .GT. ZERO IF THE FUNCTION WAS NOT ' REPRODUCIBLE ' ' KERFL -- NOT USED BY THIS ROUTINE ' ' --------------------------------------------------------------- ' SIMPLEX SHOULD USUALLY BE RUN USING A FLOATING POINT ' PRECISION OF AT LEAST TEN SIGNIFICANT DIGITS. ON MOST ' COMPUTERS, THIS REQUIRES THE USE OF DOUBLE PRECISION. 'Labels: 10, 40, 50, 70, 100, 130, 140, 150, 160, 170, 190, 220, 250, 280, 310, 340, 360, 380 'Labels: 390, 400, 420, 440, 450, 460, 470, 490, 510, 530, 550, 560, 570, 590, 610, 640, 660, ' 90, 720 (inherited from Fortran 77). Dim ZBAR(20), ZSTAR(20) ' SIMPLEX CALLS NO FUNCTIONS, EITHER EXTERNAL OR INTRINSIC. ' THE SUBROUTINES CALLED ARE FUNK, SIBEG, SIFUN, AND DATSW. ' SIMPLEX TERMINATES IF SENSE SWITCH NUMBER -NSSW- IS ON. ' THE STATEMENT CALL DATSW(NSSW,JUMP) RETURNS JUMP=1 IF ' SENSE SWITCH NUMBER -NSSW- IS ON, AND JUMP=2 IF IT IS OFF. ' IF NO SENSE SWITCH IS TO BE USED, THE USER SHOULD SUPPLY ' A DUMMY SUBROUTINE FOR DATSW. ' THIS SUBROUTINE CONTAINS NO REAL OR DOUBLE PRECISION ' CONSTANTS. Call SIBEG If KFLAG < 0 Then GoTo 720 FSAVE = FOBJ KERFL = 0 ' SET FIXED QUANTITIES.... ' METHD=1 TO USE THE METHOD OF NELDER AND MEAD, ' METHD=2 TO USE A MODIFIED METHOD. ' METHD=1 CAN CAUSE THE BEST KNow1N POINT TO BE DISCARDED, ' AND HENCE IS NOT RECOMMENDED. ' METHD=2 IS RECOMMENDED. METHD = 2 RZERO = 0# UNITR = 1# RTWO = 2# JL = NVPLUS LATER = 0 ' BEGIN THE NEXT ITERATION. ' FIND JH AND JL, THE INDICES OF THE POINTS WITH THE HIGHEST ' AND LOWEST FUNCTION VALUES, RESPECTIVELY. 10 JH = JL For J = 1 To NVPLUS If FZ(J) > FZ(JH) Then JH = J If FZ(J) < FZ(JL) Then JL = J Next J ' CHECK FOR POSSIBLE TERMINATION If KFLAG <> 0 Then GoTo 690 If (JH = JL) And (NFLAT > 0) Then GoTo 660 ' CHECK FOR TIES -- MORE THAN ONE POINT IN THE SIMPLEX HAVING ' THE HIGHEST FUNCTION VALUE. JHTIE = 0 For J = 1 To NVPLUS If (J <> JH) And (FZ(J) >= FZ(JH)) Then JHTIE = 7 Next J If (JH <> JL) And (JHTIE = 0) Then GoTo 70 ' THERE IS A TIE FOR HIGHEST FUNCTION VALUE, CHOOSE H FAR FROM L. DX = RZERO JHSAV = JH For J = 1 To NVPLUS If (J = JL) Or (FZ(J) < FZ(JH)) Then GoTo 50 For K = 1 To NV If MASK(K) <> 0 Then GoTo 40 SCALK = X(K) If SCALK = RZERO Then SCALK = DELTX(K) DZ = (Z(K, J) - Z(K, JL)) / SCALK If DZ < RZERO Then DZ = -DZ If DZ <= DX Then GoTo 40 JH = J DX = DZ 40 Next K 50 Next J If NTRACE >= 2 Then Print #2, " " Print #2, " TIE BREAKING.... JHSAV ="; JHSAV; " JH ="; JH End If If DX <= RZERO Then GoTo 640 70 If NTRACE < 2 Then GoTo 100 Print #2, " " Print #2, " FZ("; JL; ") = "; FZ(JL) Print #2, " Z(J,JL) = "; For J = 1 To NV Print #2, " "; Z(J, JL); Next J Print #2, " " Print #2, " " Print #2, " FZ("; JH; ") = "; FZ(JH); Print #2, " Z(J,JH) = "; For J = 1 To NV Print #2, " "; Z(J, JH); Next J Print #2, " " ' ADD EXTRA POINTS TO THE SIMPLEX, IF DESIRED. ' ANY EXTRA POINTS WILL BE SUPERIMPOSED ON THE HIGHEST POINT, ' P(JH), SO THAT THEY WILL SPLIT AS SOON AS POSSIBLE. 100 If LATER <> 0 Then GoTo 140 LATER = 7 If NVPLUS + NXTRA > MAXPT Then NXTRA = MAXPT - NVPLUS If NXTRA <= 0 Then GoTo 130 JLOW = NVPLUS + 1 NVPLUS = NVPLUS + NXTRA For J = JLOW To NVPLUS FZ(J) = FZ(JH) For K = 1 To NV Z(K, J) = Z(K, JH) Next K Next J 130: FSAVE = FZ(JL) ' CALCULATE PBAR, THE CENTROID OF ALL POINTS EXCEPT P(JH). ' THE MEAN VALUE IS COMPUTED USING A STABLE UPDATE FORMULA BY ' D. H. D. WEST, COMMUNICATIONS OF THE A.C.M. 22 (1979) 532-535. 140 For J = 1 To NV ZBARJ = X(J) If MASK(J) <> 0 Then GoTo 190 ZMAX = Z(J, JL) ZMIN = Z(J, JL) XS = RZERO For K = 1 To NVPLUS If K = JH Then GoTo 150 ZJK = Z(J, K) If ZJK > ZMAX Then ZMAX = ZJK If ZJK < ZMIN Then ZMIN = ZJK XS = XS + UNITR ZBARJ = ZBARJ + (ZJK - ZBARJ) / XS 150 Next K ' CHECK THE ROUNDING. ' ZBARJ MUST LIE IN THE CLOSED INTERVAL (ZMIN,ZMAX). If ZBARJ <= ZMAX Then GoTo 160 ZBARJ = ZMAX GoTo 170 160 If ZBARJ >= ZMIN Then GoTo 190 ZBARJ = ZMIN 170 Now1 = 1 If NTRACE >= 1 Then Print #2, " " Print #2, " ROUNDING OF X("; J; ") CORRECTED AT CHECKPOINT "; Now1; " WITH NF ="; NF End If 190 ZBAR(J) = ZBARJ Next J If NTRACE >= 2 Then Print #2, " " Print #2, " ZBAR(J) = "; For J = 1 To NV Print #2, " "; ZBAR(J); Next J Print #2, " " End If ' ATTEMPT A REFLECTION. FORM P* . ' THE FORMS USED BELOW FOR REFLECTION, EXPANSION, AND ' CONTRACTION ALL HAVE OPTIMAL ROUNDOFF PROPERTIES. JDIFF = 0 For J = 1 To NV XJ = X(J) If MASK(J) <> 0 Then GoTo 220 XJ = ZBAR(J) + ALPHA * (ZBAR(J) - Z(J, JH)) If XJ <> Z(J, JH) Then JDIFF = 7 X(J) = XJ 220 ZSTAR(J) = XJ Next J If JDIFF <> 0 Then GoTo 250 If NTRACE >= 1 Then Print #2, " ZSTAR = Z(JH), TREAT AS A CONTRACTION FAILURE" End If GoTo 440 250 Call SIFUN FSTAR = FOBJ KONTR = 0 If FSTAR < FZ(JL) Then GoTo 280 If NTRACE >= 1 Then Print #2, " " Print #2, " FOBJ ="; FSTAR; " REFLECTION FAILED, NF ="; NF End If For J = 1 To NVPLUS ' HERE WE DEVIATE FROM THE FLOWCHART IN THE ARTICLE BY NELDER ' AND MEAD. THEY USE FSTAR.LE.FZ(J) . ' THAT IS NOT CONSISTENT WITH THE TEXT OF THE ARTICLE, AND CAN ' CREATE AN INFINITE LOOP OF REFLECTIONS. If (J <> JH) And (FSTAR < FZ(J)) Then GoTo 340 Next J ' HERE WE DEVIATE FROM THE FLOWCHART IN THE ARTICLE BY NELDER ' AND MEAD. THEY USE FSTAR.LE.FZ(JH) (SEE COMMENT ABOVE). KONTR = 1 If FSTAR < FZ(JH) Then GoTo 340 GoTo 380 280 If NTRACE < 1 Then GoTo 310 Print #2, " " Print #2, " FOBJ ="; FSTAR; " REFLECTION SUCCEEDED, NF ="; NF Print #2, " X = "; For J = 1 To NV Print #2, " "; X(J); Next J Print #2, " " ' THE REFLECTED VALUE FSTAR IS LESS THAN THE LOWEST VALUE, ' FZ(JL), IN THE SIMPLEX. ' ATTEMPT AN EXPANSION. FORM P** . 310 For J = 1 To NV If MASK(J) = 0 Then X(J) = X(J) + (GAMMA - UNITR) * (X(J) - ZBAR(J)) Next J Call SIFUN ' CHOOSE ONE OF TWO STRATEGIES FOR ACCEPTING THE EXPANSION POINT P** . If ((METHD = 1) And (FOBJ < FZ(JL))) Or ((METHD = 2) And (FOBJ < FSTAR)) Then GoTo 490 If NTRACE >= 1 Then Print #2, " " Print #2, " FOBJ ="; FOBJ; " EXPANSION FAILED." End If ' THE REFLECTION FAILED BUT THE REFLECTED POINT WAS BETTER ' THAN SOME OTHER POINT, OR ELSE THE REFLECTION SUCCEEDED ' BUT THE EXPANSION FAILED. ' ACCEPT THE REFLECTED POINT. REPLACE P(JH) BY P* . 340 If NTRACE < 1 Then GoTo 360 Print #2, " " Print #2, " ACCEPT REFLECTED POINT:" For J = 1 To NV Print #2, " "; ZSTAR(J); Next J Print #2, " " 360 For J = 1 To NV Z(J, JH) = ZSTAR(J) Next J FZ(JH) = FSTAR ' IF THE REFLECTION FAILED AND THE REFLECTED POINT WAS BETTER ' THAN ONLY P(JH), CONTRACT THE SIMPLEX. If KONTR = 0 Then GoTo 530 ' ATTEMPT A CONTRACTION. FORM P** 380: JDIFF = 0 For J = 1 To NV If MASK(J) <> 0 Then GoTo 400 ZBARJ = ZBAR(J) ZJJH = Z(J, JH) XJ = ZBARJ + BETA * (ZJJH - ZBARJ) ' CHECK THE ROUNDING. ' XJ MIGHT HAVE BEEN ROUNDED UP TO ZJJH, WHICH COULD CREATE ' AN INFINITE LOOP. ' XJ MUST LIE IN THE CLOSED INTERVAL (ZJJH,ZBARJ), AND ' XJ MUST NOT BE EQUAL TO ZJJH UNLESS ZJJH.EQ.ZBARJ . ' EQUIVALENTLY, AN OPEN INTERVAL (ZJJH,ZBARJ) MUST EXIST AND ' XJ MUST LIE IN IT, OR XJ MUST BE EQUAL TO ZBARJ. If ((XJ > ZJJH) And (XJ < ZBARJ)) Or ((XJ > ZBARJ) And (XJ < ZJJH)) Then GoTo 390 If XJ = ZBARJ Then GoTo 390 Now1 = 2 If NTRACE >= 1 Then Print #2, " " Print #2, " ROUNDING OF X("; J; ") CORRECTED AT CHECKPOINT "; Now1; " WITH NF ="; NF End If XJ = ZBARJ 390 If XJ <> ZJJH Then JDIFF = 7 X(J) = XJ 400 Next J If JDIFF = 0 Then GoTo 420 Call SIFUN If FOBJ > FZ(JH) Then GoTo 420 If NTRACE < 1 Then GoTo 420 Print #2, " " Print #2, " FOBJ ="; FOBJ; " CONTRACTION SUCCEEDED, NF ="; NF GoTo 510 420 If NTRACE >= 1 Then Print #2, " " Print #2, " FOBJ ="; FOBJ; ", CONTRACTION FAILED." End If 440: JS = JL ' REPLACE ALL P(J) BY (P(J)+P(JL))/2.0 For J = 1 To NVPLUS If J = JL Then GoTo 470 For K = 1 To NV If MASK(K) <> 0 Then GoTo 460 ZKJ = Z(K, J) ZKJL = Z(K, JL) XK = ZKJL + (ZKJ - ZKJL) / RTWO ' CHECK THE ROUNDING. ' XK MIGHT HAVE BEEN ROUNDED UP TO ZKJ, WHICH COULD CREATE ' AN INFINITE LOOP. ' XK MUST LIE IN THE CLOSED INTERVAL (ZKJ,ZKJL), AND ' XK MUST NOT BE EQUAL TO ZKJ UNLESS ZKJ.EQ.ZKJL . ' EQUIVALENTLY, AN OPEN INTERVAL (ZKJ,ZKJL) MUST EXIST AND ' XK MUST LIE IN IT, OR XK MUST BE EQUAL TO ZKJL If ((XK > ZKJ) And (XK < ZKJL)) Or ((XK > ZKJL) And (XK < ZKJ)) Then GoTo 450 If XK = ZKJL Then GoTo 450 Now1 = 3 If NTRACE >= 1 Then Print #2, " " Print #2, " ROUNDING OF X("; J; ") CORRECTED AT CHECKPOINT "; Now1; " WITH NF ="; NF End If XK = ZKJL 450: Z(K, J) = XK X(K) = XK 460 Next K Call SIFUN If FOBJ < FZ(JS) Then JS = J FZ(J) = FOBJ 470 Next J If JS = JL Then GoTo 530 JL = JS If NTRACE >= 1 Then Print #2, " " Print #2, " FZ(JS) = "; FZ(JS) Print #2, " Z = "; For K = 1 To NV Print #2, " "; Z(K, JS); Next K Print #2, " " End If GoTo 530 490 If NTRACE < 1 Then GoTo 510 Print #2, " " Print #2, " FOBJ ="; FOBJ; ", EXPANSION SUCCEEDED." Print #2, " X = "; For J = 1 To NV Print #2, " "; X(J); Next J Print #2, " " ' REPLACE P(JH) BY P* OR P** 510 For J = 1 To NV Z(J, JH) = X(J) Next J FZ(JH) = FOBJ ' TEST FOR CONVERGENCE 530 For J = 1 To NV If MASK(J) <> 0 Then GoTo 550 ZMAX = Z(J, 1) ZMIN = Z(J, 1) For K = 2 To NVPLUS ZJK = Z(J, K) If ZJK > ZMAX Then ZMAX = ZJK If ZJK < ZMIN Then ZMIN = ZJK Next K DZ = ZMAX - ZMIN If DZ > DELMIN(J) Then GoTo 560 550 Next J GoTo 640 ' RETURN IF THE SENSE SWITCH IS ON 560 JUMP = 2 DATSW NSSW, JUMP If JUMP <= 1 Then GoTo 570 If NF > NFMAX Then GoTo 590 GoTo 10 570 KFLAG = -3 If NTRACE >= -1 Then Print #2, " " Print #2, " ABNORMAL TERMINATION... TERMINATED BY OPERATOR VIA SENSE SWITCH "; NSSW End If GoTo 610 590 KFLAG = -2 If NTRACE >= -1 Then Print #2, " " Print #2, " ABNORMAL TERMINATION... MORE THAN "; NFMAX; " CALLS TO THE FOBJ SUBROUTINE." End If 610 If NTRACE < -1 Then GoTo 690 Print #2, " " For J = 1 To NVPLUS Print #2, " SIMPLEX POINT "; J; " FOBJ="; FZ(J) Print #2, " X(J) = "; For K = 1 To NV Print #2, " "; Z(K, J); Next K Print #2, " " Next J GoTo 690 640 KFLAG = 1 If NTRACE >= 0 Then Print #2, " " Print #2, " TERMINATED WHEN THE DIMENSIONS OF THE SIMPLEX BECAME AS SMALL" Print #2, " AS THE DELMIN(J)." End If GoTo 690 660 KFLAG = 2 If NTRACE >= 0 Then Print #2, " " Print #2, " TERMINATED WHEN THE FUNCTION VALUES AT ALL POINTS OF THE SIMPLEX" Print #2, " WERE EXACTLY EQUAL." End If ' GO BACK AND COMPUTE JL AND JH ONE LAST TIME GoTo 10 ' FINISH UP AND EXIT 690 For J = 1 To NV X(J) = Z(J, JL) Next J Call SIFUN If (FOBJ <= FSAVE) And (FOBJ = FZ(JL)) Then GoTo 720 NOREP = NOREP + 2 If NTRACE >= -1 Then Print #2, " " Print #2, " WARNING.... FOBJ IS NOT A REPRODUCIBLE FUNCTION OF X(J)." Print #2, " NF="; NF; " FSAVE="; FSAVE; " FZ(JL)="; FZ(JL); " FOBJ="; FOBJ End If 720 If NTRACE >= 0 Then Print #2, " " Print #2, " " Print #2, " " Print #2, " "; NF; " FUNCTION COMPUTATIONS." Print #2, " " Print #2, " FINAL VALUE OF FOBJ = "; FOBJ Print #2, " " Print #2, " FINAL VALUES OF X(J):" For J = 1 To NV Print #2, " "; X(J); Next J Print #2, " " End If Print #2, " " End Sub 'SMPLX() Sub SIBEG() '-------------------------------------------------------------- ' SIBEG 1.3 OCTOBER 1991 ' ' A.N.S.I. STANDARD FORTRAN 77 ' ' COPYRIGHT (C) 1965, 1975, 1990 J. P. CHANDLER, ' COMPUTER SCIENCE DEPARTMENT, ' OKLAHOMA STATE UNIVERSITY ' ' SIBEG SETS DEFAULT VALUES AND PRINTS INITIAL OUTPUT ' FOR SIMPLEX. ' ' ------------------------------------------------------------ ' ' INPUT QUANTITIES..... FUNK,X(*),XMAX(*),XMIN(*),DELTX(*), ' DELMIN(*),NV,NTRACE,MASK(*), ' NFMAX,NFLAT,NXTRA,KW ' ' OUTPUT QUANTITIES.... FZ,HUGE,ALPHA,BETA,GAMMA,MAXPT, ' NVPLUS,NSSW,NF,KFLAG,NOREP, ' DELTX(*), AND DELMIN(*) '------------------------------------------------------------- 'Labels: 20, 30, 50, 60, 70, 80, 90, 100, 180, 200, 230, 260 ' * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ' ' SET FIXED QUANTITIES .... ' ' NSSW = SENSE SWITCH NUMBER FOR SEARCH TERMINATION ' (IRRELEVANT IF A DUMMY DATSW ROUTINE IS USED) NSSW = 6 ' HUGE = A VERY LARGE REAL NUMBER ' (DEFAULT VALUE FOR XMAX AND -XMIN, AND ' OUT-OF-BOUNDS VALUE FOR FOBJ) HUGE = 1E+35 ' NVMAX IS THE MAXIMUM VALUE OF NV. ' NVMAX IS ALSO THE DIMENSION OF X, XMAX, XMIN, DELTX, DELMIN, ' MASK, ZBAR, ZSTAR, AND THE FIRST DIMENSION OF ERR AND Z. NVMAX = 20 ' MAXPT IS THE MAXIMUM NUMBER OF POINTS IN THE SIMPLEX. ' MAXPT IS ALSO THE DIMENSION OF FZ, AND THE SECOND DIMENSION ' OF ERR AND Z. ' MAXPT MUST BE .GE. (NVMAX+1). MAXPT = 21 ' DELDF = DEFAULT VALUE FOR DELTX(J) DELDF = 0.01 ' ALPHA = REFLECTION PARAMETER ALPHA = 1# ' BETA = CONTRACTION PARAMETER BETA = 0.5 ' GAMMA = EXPANSION PARAMETER GAMMA = 2# RZERO = 0# ' NO REAL OR DOUBLE PRECISION CONSTANTS ARE USED BEYOND THIS ' POINT, IN THIS SUBROUTINE KFLAG = 0 NOREP = 0 KERFL = 0 ' CHECK SOME INPUT QUANTITIES, AND SET THEM TO DEFAULT VALUES ' IF DESIRED ' MAKE SURE THAT THE SENSE SWITCH IS OFF JUMP = 2 DATSW NSSW, JUMP If JUMP = 2 Then GoTo 30 Form1.Print " TURN OFF SENSE SWITCH "; NSSW 20 DATSW NSSW, JUMP If JUMP = 1 Then GoTo 20 30 JVARY = 0 If (NV >= 1) And (NV <= NVMAX) And (NVMAX <= MAXPT) Then GoTo 50 KFLAG = -1 Print #2, " " Print #2, " ERROR IN INPUT TO SUBROUTINE SIMPLEX, NV="; NV; " NVMAX="; NVMAX; " MAXPT="; MAXPT If NTRACE >= -1 Then Print #2, " NACTIV="; NACTIV End 50 For J = 1 To NV If MASK(J) <> 0 Then GoTo 100 If DELMIN(J) < RZERO Then DELMIN(J) = -DELMIN(J) ' CHECK THAT DELTX(J) IS NOT NEGLIGIBLE XPLUS = X(J) + DELTX(J) If XPLUS = X(J) Then GoTo 60 XPLUS = X(J) - DELTX(J) If XPLUS <> X(J) Then GoTo 80 60 If X(J) = RZERO Then GoTo 70 DELTX(J) = DELDF * X(J) GoTo 80 70 DELTX(J) = DELDF 80 If XMAX(J) > XMIN(J) Then GoTo 90 XMAX(J) = HUGE XMIN(J) = -HUGE 90 If X(J) > XMAX(J) Then X(J) = XMAX(J) If X(J) < XMIN(J) Then X(J) = XMIN(J) 100 Next J If NTRACE < 0 Then GoTo 180 Print #2, " " Print #2, " SIMPLEX MINIMIZATION, INITIAL VALUES...." Print #2, " MASK = "; For J = 1 To NV Print #2, " "; MASK(J); Next J Print #2, " " Print #2, " X = "; For J = 1 To NV Print #2, " "; X(J); Next J Print #2, " " Print #2, " XMAX = "; For J = 1 To NV Print #2, " "; XMAX(J); Next J Print #2, " " Print #2, " XMIN = "; For J = 1 To NV Print #2, " "; XMIN(J); Next J Print #2, " " Print #2, " DELTX = "; For J = 1 To NV Print #2, " "; DELTX(J); Next J Print #2, " " Print #2, " DELMIN = "; For J = 1 To NV Print #2, " "; DELMIN(J); Next J Print #2, " " ' CALCULATE THE INITIAL P(I) AND Y(I) (Z AND FZ, ' RESPECTIVELY). ' NF = NUMBER OF CALLS TO FUNK SO FAR 180 NF = 0 ' NACT = NUMBER OF ACTIVE X(J) NACTIV = 0 For J = 1 To NV If MASK(J) <> 0 Then GoTo 200 NACTIV = NACTIV + 1 For K = 1 To NV Z(K, NACTIV) = X(K) Next K Z(J, NACTIV) = Z(J, NACTIV) + DELTX(J) XS = X(J) X(J) = Z(J, NACTIV) Call SIFUN X(J) = XS FZ(NACTIV) = FOBJ 200 Next J If NTRACE >= 0 Then Print #2, " " Print #2, " "; NV; " VARIABLES, "; NACTIV; " ACTIVE." Print #2, " NFMAX="; NFMAX; " NFLAT="; NFLAT; " NXTRA="; NXTRA Print #2, " ALPHA="; ALPHA; " BETA="; BETA; " GAMMA="; GAMMA End If If NACTIV > 0 Then GoTo 230 KFLAG = -1 If NTRACE >= -1 Then Print #2, " " Print #2, " WARNING ... MASK(J).EQ.0 FOR ALL J IN PROCEDURE SMPLX." Print #2, " FOBJ WILL BE EVALUATED BUT NOT MINIMIZED." End If ' SET THE BASE POINT 230 NVPLUS = NACTIV + 1 For K = 1 To NV Z(K, NVPLUS) = X(K) Next K Call SIFUN FZ(NVPLUS) = FOBJ ' CHECK THE REPRODUCIBILITY OF FOBJ, AND PRINT THE REMAINING LINES FSAVE = FOBJ SIFUN If FOBJ = FSAVE Then GoTo 260 NOREP = 1 If NTRACE >= -1 Then Print #2, " " Print #2, " WARNING.... FOBJ IS NOT A REPRODUCIBLE FUNCTION OF X(J)." Print #2, " NF="; NF; " FSAVE="; FSAVE; " FOBJ="; FOBJ End If 260 If NTRACE >= 0 Then Print #2, " " Print #2, " FOBJ="; FOBJ Print #2, " " Print #2, " BEGIN MINIMIZATION...." End If End Sub 'SIBEG() Sub SIFUN() '----------------------------------------------------- ' SIFUN 1.2 DECEMBER 1991 ' ' A.N.S.I. STANDARD FORTRAN 77 ' ' COPYRIGHT (C) 1975, 1990 BY J. P. CHANDLER, ' COMPUTER SCIENCE DEPARTMENT, ' OKLAHOMA STATE UNIVERSITY ' ' SIFUN CALLS FUNK TO COMPUTE FOBJ, IF X IS IN BOUNDS. ' SIFUN IS CALLED BY SIMPLEX AND SIBEG. '----------------------------------------------------- 'Labels: 20, 30 For JF = 1 To NV If ((MASK(JF) = 0) And (X(JF) > XMAX(JF))) Or (X(JF) < XMIN(JF)) Then GoTo 20 End If Next JF Call FUNK NF = NF + 1 GoTo 30 20 FOBJ = HUGE 30 If NTRACE >= 2 Then Print #2, " " Print #2, " TRIAL FOBJ="; FOBJ Print #2, " X(J) = "; For JF = 1 To NV Print #2, " "; X(JF); Next JF Print #2, " " End If End Sub 'SIFUN() Sub DATSW(NSSW, JUMP) '----------------------------------------------------------- ' DUMMY VERSION OF SUBROUTINE DATSW -- ALL SWITCHES OFF. ' ' J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, ' OKLAHOMA STATE UNIVERSITY '---------------------------------------------------------- JUMP = 2 End Sub Sub STSET() '------------------------------------------------------------ ' STSET 3.2 DECEMBER 1991 ' ' STSET SETS SOME INPUT QUANTITIES TO DEFAULT VALUES, FOR ' SUBROUTINE SIMPLEX. ' ' J. P. CHANDLER, DEPARTMENT OF COMPUTER SCIENCE, ' OKLAHOMA STATE UNIVERSITY ' ' USAGE..... ' ' CALL STSET. ' THEN SET SOME INPUT QUANTITIES (NV, AT LEAST) AND RESET ANY ' OF THOSE SET IN STSET (BETTER VALUES OF X(J), ETC.) BEFORE ' CALLING SIMPLEX. '------------------------------------------------------------ HUGE = 1E+30 RZERO = 0# ' NVMAX IS THE MAXIMUM PERMISSIBLE VALUE OF NV, GIVEN THE ' PRESENT DIMENSIONS OF ARRAYS. ' NVMAX IS THE DIMENSION OF THE ARRAYS X(*), XMAX(*), XMIN(*), ' DELTX(*), DELMIN(*), AND MASK(*). NVMAX IS ALSO THE FIRST ' DIMENSION OF ERR(*,*). THE SECOND DIMENSION OF ERR(*,*) ' IS NVMAX+1 NVMAX = 20 ' THE USER MUST SET NV AFTER CALLING STSET NV = -1 NTRACE = 0 NFMAX = 32000 MAXIT = 50 MAXSUB = 30 METHD = 1 KALCP = 0 LEQU = 0 NFLAT = 1 MATRX = 105 NXTRA = 0 FLAMBD = 1# FNU = 10# KORDIF = 1 RELDIF = 0.00000001 RELMIN = 0.000001 For JX = 1 To NVMAX X(JX) = RZERO XMAX(JX) = HUGE XMIN(JX) = -HUGE DELTX(JX) = RZERO DELMIN(JX) = RZERO MASK(JX) = 0 Next JX End Sub 'STSET() 'main subroutine Sub Exec_Simplex() 'Open input/output files Open "SIMPTEST.DAT" For Input As #1 Open "SIMPOUT.LST" For Output As #2 Print #2, " *** SIMPLEX PROGRAM ***" ' READ IN THE VALUE OF NPTS, THEN READ IN THE DATA POINTS Input #1, NPTS Print #2, " NPTS = "; NPTS For K = 1 To NPTS Input #1, T(K, 1), T(K, 2), Y(K), YSIG(K) Print #2, " "; K; " "; T(K, 1); " "; T(K, 2); " "; Y(K); " "; YSIG(K) Next K Close #1 ' INITIALIZE FOR THE FIT Call STSET NV = 3 'Number of constants X(J) NTRACE = 0 'Normal output ' NTRACE = 1 For additional output NFMAX = 550 'Maximum number of FUNK evaluations X(1) = 10# 'Initial guess X(2) = 1# X(3) = 1# ' FIT THE MODEL TO THE DATA USING SIMPLEX METHOD Call SMPLX ' USE PROCEDURE FIDO TO COMPUTE ERRORS USING SIMPLEX ERFRAC = 0.1 For JXFID = 1 To NV ERGSAV(JXFID) = ERFRAC * Abs(X(JXFID)) If ERGSAV(JXFID) = 0# Then ERGSAV(JXFID) = ERFRAC Next JXFID STDEVS = 1# TOLFID = 0.05 MAXIND = 2 NTRACF = -1 '0 or 1 to print in FIDO NTRSET = -1 For JXFID = 1 To NV ERGUES = ERGSAV(JXFID) FIDO JXFID, STDEVS, TOLFID, ERGUES, MAXIND, NTRSET, NTRACF, FIDINT() Next JXFID Close #2 Form1.Print Form1.Print " Results in file Simpout.lst." Form1.Print " Program terminated." Form1.Print End Sub 'Exec_Simplex()