* * $Id: mnline.fdoc,v 1.2 2002/06/07 22:17:50 sachs Exp $ * * $Log: mnline.fdoc,v $ * Revision 1.2 2002/06/07 22:17:50 sachs * Added comments to mnline.fdoc * * Revision 1.1 2002/05/14 15:32:56 mf * Added Minuit source, plus other stuff from a year ago. * * Revision 1.1.1.1 1996/03/07 14:31:30 mclareni * Minuit * * #include "minuit/pilot.h" SUBROUTINE MNLINE(FCN,START,FSTART,STEP,SLOPE,TOLER,FUTIL) #include "minuit/d506dp.inc" CC Perform a line search from position START CC along direction STEP, where the length of vector STEP CC gives the expected position of minimum. CC FSTART is value of function at START CC SLOPE (if non-zero) is df/dx along STEP at START CCC Actually is predicted change in value of FCN 1 step away CCC This function does NOT appear to provide for a zero value of SLOPE CCC However the othjer minuit functions do appear to provide the proper CCC value of slope, which is -GRAD*STEP, where GRAD is the Gradient of CCC FCN at the start point. FSTART+SLOPE is the first order predicted CCC value of FCN at the point START+STEP. CC TOLER is initial tolerance of minimum in direction STEP #include "minuit/d506cm.inc" EXTERNAL FCN,FUTIL DIMENSION START(*), STEP(*) XPQ(I) stores the SLAM value used to get the function output stored in YPQ(I). i.e. YPQ(I) = f(start + step * XPQ(I)) MAXPT is the maximum number of values stored in these arrays before the algorithm will exit. (see line 65) CHPQ is used to store the letters of the alphabet (see CHARAL below). Exactly why it is used remains a mystery. PARAMETER (MAXPT=12) DIMENSION XPQ(MAXPT),YPQ(MAXPT) CHARACTER*1 CHPQ(MAXPT) XCALS, FVALS, and COEFF are used when finding the parabolic fit for the three best points discovered sofar. DIMENSION XVALS(3),FVALS(3),COEFF(3) CHARAL stores the alphabet, CMESS is not understood CHARACTER*26 CHARAL CHARACTER*60 CMESS PARAMETER (SLAMBG=5.,ALPHA=2.) C SLAMBG and ALPHA control the maximum individual steps allowed. C The first step is always =1. The max length of second step is SLAMBG. C The max size of subsequent steps is the maximum previous successful C step multiplied by ALPHA + the size of most recent successful step, C but cannot be smaller than SLAMBG. LOGICAL LDEBUG DATA CHARAL / 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' / LDEBUG = (IDBG(1).GE.1) C starting values for overall limits on total step SLAM OVERAL = 1000. UNDRAL = -100. C debug check if start is ok IF (LDEBUG) THEN CALL MNINEX(START) CALL FCN(NPARX,GIN,F1,U,4,FUTIL) NFCN=NFCN+1 IF (F1 .NE. FSTART) THEN WRITE (ISYSWR,'(A/2E14.5/2X,10F10.5)') + ' MNLINE start point not consistent, F values, parameters=', + (X(KK),KK=1,NPAR) ENDIF ENDIF C . set up linear search along STEP FVMIN refers to the current best minimum XVMIN refers to the SLAM value used to get there FVMIN = f(start + step * XVMIN) FVMIN = FSTART XVMIN = ZERO NXYPT is the last index used in the *PQ arrays NXYPT = 1 CHPQ(I) stores the Ith letter of the alphabet if XPQ(I) and YPQ(I) are in use. Why? CHPQ(1) = CHARAL(1:1) We use start as out first data point XPQ(1) = 0. YPQ(1) = FSTART C SLAMIN = smallest possible value of ABS(SLAM) SLAMIN = ZERO This loop sets up SLAMIN based on the minimum ratio individual elements of START and STEP. DO 20 I= 1, NPAR IF (STEP(I) .EQ. ZERO) GO TO 20 RATIO = ABS(START(I)/STEP(I)) IF (SLAMIN .EQ. ZERO) SLAMIN = RATIO IF (RATIO .LT. SLAMIN) SLAMIN = RATIO We take out first step with SLAM = 1 20 X(I) = START(I) + STEP(I) if START is the zero vector (i.e. all ratios are zero) then the machine epsilon is used as the minimum ratio. IF (SLAMIN .EQ. ZERO) SLAMIN = EPSMAC SLAMIN = SLAMIN*EPSMA2 preserve NPAR NPARX = NPAR C CALL MNINEX(X) CALL FCN(NPARX,GIN,F1,U,4,FUTIL) NFCN=NFCN+1 Here we store the next value of FCN (with SLAM=1) NXYPT = NXYPT + 1 CHPQ(NXYPT) = CHARAL(NXYPT:NXYPT) XPQ(NXYPT) = 1. YPQ(NXYPT) = F1 Determine the best guess sofar IF (F1 .LT. FSTART) THEN FVMIN = F1 XVMIN = 1.0 ENDIF C . quadr interp using slope GDEL and two points Finally initialize SLAM to the value we've been assuming in our above calculations SLAM = 1. preserve TOLER TOLER8 = TOLER SLAMAX starts at SLAMBG, but can change SLAMAX = SLAMBG FLAST = F1 C can iterate on two-points (cut) if no imprvmnt CCC We might come back here based on our next result 25 CONTINUE CCC The following formula for computing the value of SLAM to try next CCC is EXACT if the function is purely quadratic along the line CCC being searched. However for more complex functions, it can fail CCC badly. The denominator (DENOM) of the expression for the new CCC value of SLAM will be zero if the predicted value of FCN CCC (FSTART + SLOPE*STEP) is exactly equal to the computed value. CCC This normally requires, that there must be at least one maximum CCC and one minimum in the interval. DENOM = 2.0*(FLAST-FSTART-SLOPE*SLAM)/SLAM**2 C IF (DENOM .EQ. ZERO) DENOM = -0.1*SLOPE SLAM = 1. IF (DENOM .NE. ZERO) SLAM = -SLOPE/DENOM IF (SLAM .LT. ZERO) SLAM = SLAMAX IF (SLAM .GT. SLAMAX) SLAM = SLAMAX IF (SLAM .LT. TOLER8) SLAM = TOLER8 IF (SLAM .LT. SLAMIN) GO TO 80 IF (ABS(SLAM-1.0).LT.TOLER8 .AND. F1.LT.FSTART) GO TO 70 IF (ABS(SLAM-1.0).LT.TOLER8) SLAM = 1.0+TOLER8 IF (NXYPT .GE. MAXPT) GO TO 65 Now we take the next step DO 30 I= 1, NPAR 30 X(I) = START(I) + SLAM*STEP(I) CALL MNINEX(X) CALL FCN(NPAR,GIN,F2,U,4,FUTIL) NFCN = NFCN + 1 and store it NXYPT = NXYPT + 1 CHPQ(NXYPT) = CHARAL(NXYPT:NXYPT) XPQ(NXYPT) = SLAM YPQ(NXYPT) = F2 and see if it's the best point sofar IF (F2 .LT. FVMIN) THEN FVMIN = F2 XVMIN = SLAM ENDIF if we've made no progress at all IF (FSTART .EQ. FVMIN) THEN FLAST = F2 TOLER8 = TOLER*SLAM OVERAL = SLAM-TOLER8 SLAMAX = OVERAL GO TO 25 ENDIF C . quadr interp using 3 points these are the first three points we use for the parabolic fit XVALS(1) = XPQ(1) FVALS(1) = YPQ(1) XVALS(2) = XPQ(NXYPT-1) FVALS(2) = YPQ(NXYPT-1) XVALS(3) = XPQ(NXYPT) FVALS(3) = YPQ(NXYPT) C begin iteration, calculate desired step 50 CONTINUE SLAMAX changes based upon ALPHA and XVMIN Note that this is not exactly how it is stated in the header comments. SLAMAX = MAX(SLAMAX,ALPHA*ABS(XVMIN)) This finds the parabolic fit (quadratic) y = COEFF(1) + COEFF(2) x + COEFF(3) x^3 SDEV is always 0, since only three points are given CALL MNPFIT(XVALS,FVALS,3,COEFF,SDEV) The sign of COEFF(3) determines the shape of the parabola. Negative means concave down. That means we're going to approach a maximum! This is a problem, so change SLAM by the maximum possible value in the direction away from the slope. IF (COEFF(3) .LE. ZERO) THEN SLOPEM = 2.0*COEFF(3)*XVMIN + COEFF(2) IF (SLOPEM .LE. ZERO) THEN SLAM = XVMIN + SLAMAX ELSE SLAM = XVMIN - SLAMAX ENDIF ELSE It's positive, so we're appreaching a minimum Change SLAM based on -b/2a (best guess for location of the minimum) SLAM = -COEFF(2)/(2.0*COEFF(3)) IF (SLAM .GT. XVMIN+SLAMAX) SLAM = XVMIN+SLAMAX IF (SLAM .LT. XVMIN-SLAMAX) SLAM = XVMIN-SLAMAX ENDIF Verify UNDRAL <= SLAM <= OVERAL IF (SLAM .GT. ZERO) THEN IF (SLAM .GT. OVERAL) SLAM = OVERAL ELSE IF (SLAM .LT. UNDRAL) SLAM = UNDRAL ENDIF A step is cut if it yields no benefit in the result i.e. the F3 value is greater than the three values we have already C come here if step was cut below 52 CONTINUE Toler9 is only used in the one line here There is no reason for it to be declared This whole section determines based on the currently three best X points whether or not the line search has already "attained tolerance". TOLER9 = MAX(TOLER8,ABS(TOLER8*SLAM)) DO 55 IPT= 1, 3 IF (ABS(SLAM-XVALS(IPT)) .LT. TOLER9) GO TO 70 55 CONTINUE C take the step First, check to see if we've exhausted the limit of function evaluations IF (NXYPT .GE. MAXPT) GO TO 65 DO 60 I= 1, NPAR 60 X(I) = START(I)+SLAM*STEP(I) CALL MNINEX(X) CALL FCN(NPARX,GIN,F3,U,4,FUTIL) NFCN = NFCN + 1 Store the value NXYPT = NXYPT + 1 CHPQ(NXYPT) = CHARAL(NXYPT:NXYPT) XPQ(NXYPT) = SLAM YPQ(NXYPT) = F3 C find worst previous point out of three We do this so that we can replace it FVMAX = FVALS(1) NVMAX = 1 IF (FVALS(2) .GT. FVMAX) THEN FVMAX = FVALS(2) NVMAX = 2 ENDIF IF (FVALS(3) .GT. FVMAX) THEN FVMAX = FVALS(3) NVMAX = 3 ENDIF C if latest point worse than all three previous, cut step IF (F3 .GE. FVMAX) THEN IF (NXYPT .GE. MAXPT) GO TO 65 IF (SLAM .GT. XVMIN) OVERAL = MIN(OVERAL,SLAM-TOLER8) IF (SLAM .LT. XVMIN) UNDRAL = MAX(UNDRAL,SLAM+TOLER8) SLAM = 0.5*(SLAM+XVMIN) GO TO 52 ENDIF C prepare another iteration, replace worst previous point XVALS(NVMAX) = SLAM FVALS(NVMAX) = F3 IF (F3 .LT. FVMIN) THEN FVMIN = F3 XVMIN = SLAM ELSE IF (SLAM .GT. XVMIN) OVERAL = MIN(OVERAL,SLAM-TOLER8) IF (SLAM .LT. XVMIN) UNDRAL = MAX(UNDRAL,SLAM+TOLER8) ENDIF Iterate until we reach our evaluation limit IF (NXYPT .LT. MAXPT) GO TO 50 C . . end of iteration . . . C stop because too many iterations 65 CMESS = ' LINE SEARCH HAS EXHAUSTED THE LIMIT OF FUNCTION CALLS ' IF (LDEBUG) THEN WRITE (ISYSWR,'(A/(2X,6G12.4))') ' MNLINE DEBUG: steps=', + (STEP(KK),KK=1,NPAR) ENDIF GO TO 100 C stop because within tolerance 70 CONTINUE CMESS = ' LINE SEARCH HAS ATTAINED TOLERANCE ' GO TO 100 80 CONTINUE CMESS = ' STEP SIZE AT ARITHMETICALLY ALLOWED MINIMUM' 100 CONTINUE Set our new minimum, and change DIRIN to match the STEP * SLAM that we used, as well as changing X to match the START + STEP * SLAM that we used. AMIN = FVMIN DO 120 I= 1, NPAR DIRIN(I) = STEP(I)*XVMIN 120 X(I) = START(I) + DIRIN(I) call MNINEX make certain that U is fixed CALL MNINEX(X) IF (XVMIN .LT. 0.) CALL MNWARN('D','MNLINE', + ' LINE MINIMUM IN BACKWARDS DIRECTION') IF (FVMIN .EQ. FSTART) CALL MNWARN('D','MNLINE', + ' LINE SEARCH FINDS NO IMPROVEMENT ') IF (LDEBUG) THEN WRITE (ISYSWR,'('' AFTER'',I3,'' POINTS,'',A)') NXYPT,CMESS CALL MNPLOT(XPQ,YPQ,CHPQ,NXYPT,ISYSWR,NPAGWD,NPAGLN) ENDIF RETURN END