source: Sophya/trunk/SophyaExt/CodeMinuit/code/mnsimp.F@ 2932

*
* $Id: mnsimp.F,v 1.1.1.1 2003-06-11 14:18:29 cmv Exp $
*
* $Log: not supported by cvs2svn $
* Revision 1.2  1996/03/15 18:02:54  james
*     Modified Files:
* mnderi.F eliminate possible division by zero
* mnexcm.F suppress print on STOP when print flag=-1
*          set FVAL3 to flag if FCN already called with IFLAG=3
* mninit.F set version 96.03
* mnlims.F remove arguments, not needed
* mnmigr.F VLEN -> LENV in debug print statement
* mnparm.F move call to MNRSET to after NPAR redefined, to zero all
* mnpsdf.F eliminate possible division by zero
* mnscan.F suppress printout when print flag =-1
* mnset.F  remove arguments in call to MNLIMS
* mnsimp.F fix CSTATU so status is PROGRESS only if new minimum
* mnvert.F eliminate possible division by zero
*
* Revision 1.1.1.1  1996/03/07 14:31:31  mclareni
* Minuit
*
*
#include "minuit/pilot.h"
      SUBROUTINE MNSIMP(FCN,FUTIL)
#include "minuit/d506dp.inc"
CC        Performs a minimization using the simplex method of Nelder
CC        and Mead (ref. -- Comp. J. 7,308 (1965)).
CC
#include "minuit/d506cm.inc"
      EXTERNAL FCN,FUTIL
      DIMENSION Y(MNI+1)
      DATA ALPHA,BETA,GAMMA,RHOMIN,RHOMAX / 1.0, 0.5, 2.0, 4.0, 8.0/
      IF (NPAR .LE. 0)  RETURN
      IF (AMIN .EQ. UNDEFI)  CALL MNAMIN(FCN,FUTIL)
      CFROM = 'SIMPLEX '
      NFCNFR = NFCN
      CSTATU= 'UNCHANGED '
      NPFN=NFCN
      NPARP1=NPAR+1
      NPARX = NPAR
      RHO1 = 1.0 + ALPHA
      RHO2 = RHO1 + ALPHA*GAMMA
      WG = 1.0/FLOAT(NPAR)
      IF (ISW(5) .GE. 0) WRITE(ISYSWR,100) EPSI
  100 FORMAT(' START SIMPLEX MINIMIZATION.    CONVERGENCE WHEN EDM .LT.'
     +,E10.2 )
         DO 2 I= 1, NPAR
         DIRIN(I) = WERR(I)
           CALL MNDXDI(X(I),I,DXDI)
           IF (DXDI .NE. ZERO) DIRIN(I)=WERR(I)/DXDI
         DMIN = EPSMA2*ABS(X(I))
         IF (DIRIN(I) .LT. DMIN)  DIRIN(I)=DMIN
    2    CONTINUE
C**       choose the initial simplex using single-parameter searches
    1 CONTINUE
      YNPP1 = AMIN
      JL = NPARP1
      Y(NPARP1) = AMIN
      ABSMIN = AMIN
      DO 10 I= 1, NPAR
      AMING = AMIN
      PBAR(I) = X(I)
      BESTX = X(I)
      KG = 0
      NS = 0
      NF = 0
    4 X(I) = BESTX + DIRIN(I)
      CALL MNINEX(X)
      CALL FCN(NPARX,GIN, F, U, 4, FUTIL)
      NFCN = NFCN + 1
      IF (F .LT. AMING)  GO TO 6
C         failure
      IF (KG .EQ. 1)  GO TO 8
      KG = -1
      NF = NF + 1
      DIRIN(I) = DIRIN(I) * (-0.4)
      IF (NF .LT. 3)  GO TO 4
C         stop after three failures
      BESTX = X(I)
      DIRIN(I) = DIRIN(I) * 3.0
      AMING = F
      GO TO 8
C
C         success
    6 BESTX = X(I)
      DIRIN(I) = DIRIN(I) * 3.0
      AMING = F
      CSTATU= 'PROGRESS  '
      KG = 1
      NS = NS + 1
      IF (NS .LT. 6)  GO TO 4
C
C         3 failures or 6 successes or
C         local minimum found in ith direction
    8 Y(I) = AMING
      IF (AMING .LT. ABSMIN)  JL = I
      IF (AMING .LT. ABSMIN)  ABSMIN = AMING
      X(I) = BESTX
      DO 9 K= 1, NPAR
    9 P(K,I) = X(K)
   10 CONTINUE
      JH = NPARP1
      AMIN=Y(JL)
      CALL MNRAZZ(YNPP1,PBAR,Y,JH,JL)
      DO 20 I= 1, NPAR
   20 X(I) = P(I,JL)
      CALL MNINEX(X)
      IF (ISW(5) .GE. 1)  CALL MNPRIN(5,AMIN)
      EDM = BIGEDM
      SIG2 = EDM
      NCYCL=0
C                                        . . . . .  start main loop
   50 CONTINUE
      IF (SIG2 .LT. EPSI .AND. EDM.LT.EPSI)     GO TO 76
      SIG2 = EDM
      IF ((NFCN-NPFN) .GT. NFCNMX)  GO TO 78
C         calculate new point * by reflection
      DO 60 I= 1, NPAR
      PB = 0.
      DO 59 J= 1, NPARP1
   59 PB = PB + WG * P(I,J)
      PBAR(I) = PB - WG * P(I,JH)
   60 PSTAR(I)=(1.+ALPHA)*PBAR(I)-ALPHA*P(I,JH)
      CALL MNINEX(PSTAR)
      CALL FCN(NPARX,GIN,YSTAR,U,4,FUTIL)
      NFCN=NFCN+1
      IF(YSTAR.GE.AMIN) GO TO 70
C         point * better than jl, calculate new point **
      CSTATU = 'PROGRESS  '
      DO 61 I=1,NPAR
   61 PSTST(I)=GAMMA*PSTAR(I)+(1.-GAMMA)*PBAR(I)
      CALL MNINEX(PSTST)
      CALL FCN(NPARX,GIN,YSTST,U,4,FUTIL)
      NFCN=NFCN+1
C         try a parabola through ph, pstar, pstst.  min = prho
      Y1 = (YSTAR-Y(JH)) * RHO2
      Y2 = (YSTST-Y(JH)) * RHO1
      RHO = 0.5 * (RHO2*Y1 -RHO1*Y2) / (Y1 -Y2)
      IF (RHO .LT. RHOMIN)  GO TO 66
      IF (RHO .GT. RHOMAX)  RHO = RHOMAX
      DO 64 I= 1, NPAR
   64 PRHO(I) = RHO*PBAR(I) + (1.0-RHO)*P(I,JH)
      CALL MNINEX(PRHO)
      CALL FCN(NPARX,GIN,YRHO, U,4,FUTIL)
      NFCN = NFCN + 1
      IF (YRHO .LT. AMIN)     CSTATU = 'PROGRESS  '
      IF (YRHO .LT. Y(JL) .AND. YRHO .LT. YSTST)  GO TO 65
      IF (YSTST .LT. Y(JL))  GO TO 67
      IF (YRHO .GT. Y(JL))  GO TO 66
C         accept minimum point of parabola, PRHO
   65 CALL MNRAZZ (YRHO,PRHO,Y,JH,JL)
      GO TO 68
   66 IF (YSTST .LT. Y(JL))  GO TO 67
      CALL MNRAZZ(YSTAR,PSTAR,Y,JH,JL)
      GO TO 68
   67 CALL MNRAZZ(YSTST,PSTST,Y,JH,JL)
   68 NCYCL=NCYCL+1
      IF (ISW(5) .LT. 2)  GO TO 50
      IF (ISW(5) .GE. 3 .OR. MOD(NCYCL, 10) .EQ. 0) CALL MNPRIN(5,AMIN)
      GO TO 50
C         point * is not as good as jl
   70 IF (YSTAR .GE. Y(JH))  GO TO 73
      JHOLD = JH
      CALL MNRAZZ(YSTAR,PSTAR,Y,JH,JL)
      IF (JHOLD .NE. JH)  GO TO 50
C         calculate new point **
   73 DO 74 I=1,NPAR
   74 PSTST(I)=BETA*P(I,JH)+(1.-BETA)*PBAR(I)
      CALL MNINEX (PSTST)
      CALL FCN(NPARX,GIN,YSTST,U,4,FUTIL)
      NFCN=NFCN+1
      IF(YSTST.GT.Y(JH)) GO TO 1
C     point ** is better than jh
      IF (YSTST .LT. AMIN)     CSTATU = 'PROGRESS  '
      IF (YSTST .LT. AMIN)  GO TO 67
      CALL MNRAZZ(YSTST,PSTST,Y,JH,JL)
      GO TO 50
C                                        . . . . . .  end main loop
   76 IF (ISW(5) .GE. 0)  WRITE(ISYSWR,'(A)')
     +                    ' SIMPLEX MINIMIZATION HAS CONVERGED.'
      ISW(4) = 1
      GO TO 80
   78 IF (ISW(5) .GE. 0)  WRITE(ISYSWR,'(A)')
     +                    ' SIMPLEX TERMINATES WITHOUT CONVERGENCE.'
      CSTATU= 'CALL LIMIT'
      ISW(4) = -1
      ISW(1) = 1
   80 DO 82 I=1,NPAR
      PB = 0.
      DO 81 J=1,NPARP1
   81 PB = PB + WG * P(I,J)
   82 PBAR(I) = PB - WG * P(I,JH)
      CALL MNINEX(PBAR)
      CALL FCN(NPARX,GIN,YPBAR,U,4,FUTIL)
      NFCN=NFCN+1
      IF (YPBAR .LT. AMIN)  CALL MNRAZZ(YPBAR,PBAR,Y,JH,JL)
      CALL MNINEX(X)
      IF (NFCNMX+NPFN-NFCN .LT. 3*NPAR)  GO TO 90
      IF (EDM .GT. 2.0*EPSI)  GO TO 1
   90 IF (ISW(5) .GE. 0)  CALL MNPRIN(5, AMIN)
      RETURN
      END