Initial commit
This commit is contained in:
233
gen/simplex.f
Normal file
233
gen/simplex.f
Normal file
@ -0,0 +1,233 @@
|
||||
SUBROUTINE SIMPLEX
|
||||
C ------------------
|
||||
|
||||
c minimization subroutine using a simple method by nelder and mead *).
|
||||
c it is reasonably fast when far from minimum, and may also be used to
|
||||
c converge to the exact minimum, but is mostly slower than "migrad".
|
||||
c
|
||||
c *) - nelder, j.a., mead, r. (1965), comput. j. 7, 308
|
||||
c - james, f., roos, m. (1985), cern documentation d506
|
||||
|
||||
include 'fit.inc'
|
||||
|
||||
integer jh,jl
|
||||
real p(maxpar,maxpar+1),y(maxpar+1)
|
||||
common/comsim/jh,jl,y,p
|
||||
|
||||
integer iflag, i, kg, ns, nf, k, ignal, ncycl, j, npfn, nparp1
|
||||
integer jhold
|
||||
real rho1, rho2, alpha, beta, gamma, rhomin, rhomax, wg
|
||||
real ynpp1, absmin, aming, bestx, f, sig2, pb, ystar, ystst
|
||||
real y1, y2, rho, yrho, ypbar
|
||||
|
||||
REAL PSTAR(MAXPAR) ,PSTST(MAXPAR) ,PBAR(MAXPAR) ,PRHO(MAXPAR)
|
||||
DATA ALPHA,BETA,GAMMA,RHOMIN,RHOMAX / 1.0, 0.5, 2.0, 4.0, 8.0/
|
||||
|
||||
IF (NPAR .LE. 0) RETURN
|
||||
NPFN=NFCN
|
||||
NPARP1=NPAR+1
|
||||
RHO1 = 1.0 + ALPHA
|
||||
RHO2 = RHO1 + ALPHA*GAMMA
|
||||
WG = 1.0/FLOAT(NPAR)
|
||||
IFLAG=4
|
||||
if (isyswr .ne. 0) WRITE(ISYSWR,100) EPSI
|
||||
100 FORMAT(/' Start SIMPLEX minimization'
|
||||
1/4X,'Convergence criterion: estimated distance to minimum edm <'
|
||||
1,E10.2)
|
||||
|
||||
DO 2 I= 1, NPAR
|
||||
IF (ISW(2) .GE. 1) DIRIN(I) = SQRT(V(I,I)*UP)
|
||||
IF (ABS(DIRIN(I)) .LT. 1.0E-5*ABS(X(I))) DIRIN(I)=1.0E-5*X(I) ! was 1e-8 M.Z. 6.02
|
||||
IF(ITAUR.LT. 1) V(I,I) = DIRIN(I)**2/UP
|
||||
2 CONTINUE
|
||||
IF (ITAUR .LT. 1) ISW(2) = 1
|
||||
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 FNCTN(X,F)
|
||||
NFCN = NFCN + 1
|
||||
IF (F .LE. 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
|
||||
NS = 6
|
||||
C SUCCESS
|
||||
6 BESTX = X(I)
|
||||
DIRIN(I) = DIRIN(I) * 3.0
|
||||
AMING = F
|
||||
KG = 1
|
||||
NS = NS + 1
|
||||
IF (NS .LT. 6) GO TO 4
|
||||
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 RAZZIA(YNPP1,PBAR)
|
||||
DO 20 I= 1, NPAR
|
||||
20 X(I) = P(I,JL)
|
||||
CALL INTOEX(X)
|
||||
IF (ISW(5) .GE. 1) CALL FIT_PRINT(0)
|
||||
SIGMA = SIGMA * 10.
|
||||
SIG2 = SIGMA
|
||||
IGNAL = 0
|
||||
NCYCL=0
|
||||
C. . . . . START MAIN LOOP
|
||||
50 CONTINUE
|
||||
IF (IGNAL .GE. 10) GO TO 1
|
||||
IF (SIG2 .LT. EPSI .AND. SIGMA.LT.EPSI) GO TO 76
|
||||
SIG2 = SIGMA
|
||||
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 FNCTN(PSTAR,YSTAR)
|
||||
NFCN=NFCN+1
|
||||
IF(YSTAR.GE.AMIN) GO TO 70
|
||||
C POINT * BETTER THAN JL, CALCULATE NEW POINT **
|
||||
DO 61 I=1,NPAR
|
||||
61 PSTST(I)=GAMMA*PSTAR(I)+(1.-GAMMA)*PBAR(I)
|
||||
CALL FNCTN(PSTST,YSTST)
|
||||
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*PSTAR(I) + (1. -RHO)*P(I,JH)
|
||||
CALL FNCTN(PRHO,YRHO)
|
||||
NFCN = NFCN + 1
|
||||
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 RAZZIA (YRHO,PRHO)
|
||||
IGNAL = MAX0(IGNAL-2, 0)
|
||||
GO TO 68
|
||||
66 IF (YSTST .LT. Y(JL)) GO TO 67
|
||||
IGNAL = MAX0(IGNAL-1, 0)
|
||||
CALL RAZZIA(YSTAR,PSTAR)
|
||||
GO TO 68
|
||||
67 IGNAL = MAX0(IGNAL-2, 0)
|
||||
675 CALL RAZZIA(YSTST,PSTST)
|
||||
68 NCYCL=NCYCL+1
|
||||
IF (ISW(5) .LT. 2) GO TO 50
|
||||
IF (ISW(5) .GE. 3 .OR. MOD(NCYCL, 10) .EQ. 0) CALL FIT_PRINT(0)
|
||||
GO TO 50
|
||||
C POINT * IS NOT AS GOOD AS JL
|
||||
70 IF (YSTAR .GE. Y(JH)) GO TO 73
|
||||
JHOLD = JH
|
||||
CALL RAZZIA(YSTAR,PSTAR)
|
||||
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 FNCTN (PSTST,YSTST)
|
||||
NFCN=NFCN+1
|
||||
IF(YSTST.GT.Y(JH)) GO TO 1
|
||||
C POINT ** IS BETTER THAN JH
|
||||
IF (YSTST .LT. AMIN) GO TO 675
|
||||
IGNAL = IGNAL + 1
|
||||
CALL RAZZIA(YSTST,PSTST)
|
||||
GO TO 50
|
||||
C. . . . . . END MAIN LOOP
|
||||
76 if (isyswr .ne. 0) WRITE(ISYSWR,120)
|
||||
120 format(/' SIMPLEX minimization has converged')
|
||||
GO TO 80
|
||||
78 if (nfcnmx .eq. 0) then
|
||||
if (isyswr .ne. 0) write(isyswr,'(/a)') ' SIMPLEX aborted'
|
||||
else
|
||||
if (isyswr .ne. 0) WRITE(ISYSWR,130)
|
||||
130 format(/' SIMPLEX terminates without convergence')
|
||||
endif
|
||||
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 FNCTN(PBAR,YPBAR)
|
||||
NFCN=NFCN+1
|
||||
IF (YPBAR .LT. AMIN) CALL RAZZIA(YPBAR,PBAR)
|
||||
IF (NFCNMX+NPFN-NFCN .LT. 3*NPAR) GO TO 90
|
||||
IF (SIGMA .GT. 2.0*EPSI) GO TO 1
|
||||
90 CALL FNCTN(X,AMIN)
|
||||
NFCN = NFCN + 1
|
||||
CALL FIT_PRINT(1-ITAUR)
|
||||
RETURN
|
||||
END
|
||||
|
||||
SUBROUTINE RAZZIA(YNEW,PNEW)
|
||||
C ------------------------------
|
||||
|
||||
include 'fit.inc'
|
||||
|
||||
integer jh,jl
|
||||
real p(maxpar,maxpar+1),y(maxpar+1)
|
||||
common/comsim/jh,jl,y,p
|
||||
|
||||
integer i, nparp1, j
|
||||
real ynew, us, pbig, plit
|
||||
real PNEW(MAXPAR)
|
||||
|
||||
DO 10 I=1,NPAR
|
||||
10 P(I,JH)=PNEW(I)
|
||||
Y(JH)=YNEW
|
||||
IF(YNEW.GE.AMIN) GO TO 18
|
||||
DO 15 I=1,NPAR
|
||||
15 X(I)=PNEW(I)
|
||||
CALL INTOEX(X)
|
||||
AMIN=YNEW
|
||||
JL=JH
|
||||
18 CONTINUE
|
||||
JH=1
|
||||
NPARP1=NPAR+1
|
||||
20 DO 25 J=2,NPARP1
|
||||
IF (Y(J) .GT. Y(JH)) JH = J
|
||||
25 CONTINUE
|
||||
SIGMA = Y(JH) - Y(JL)
|
||||
IF (SIGMA .LE. 0.) GO TO 45
|
||||
US = 1.0/SIGMA
|
||||
DO 35 I= 1, NPAR
|
||||
PBIG = P(I,1)
|
||||
PLIT = PBIG
|
||||
DO 30 J= 2, NPARP1
|
||||
IF (P(I,J) .GT. PBIG) PBIG = P(I,J)
|
||||
IF (P(I,J) .LT. PLIT) PLIT = P(I,J)
|
||||
30 CONTINUE
|
||||
DIRIN(I) = PBIG - PLIT
|
||||
IF (ITAUR .LT. 1 ) V(I,I) = 0.5*(V(I,I) +US*DIRIN(I)**2)
|
||||
35 CONTINUE
|
||||
40 RETURN
|
||||
45 WRITE (ISYSWR, 1000) NPAR
|
||||
GO TO 40
|
||||
1000 FORMAT (
|
||||
1 ' ***** FUNCTION VALUE DOES NOT SEEM TO DEPEND ON ANY OF THE '
|
||||
1, I3,' VARIABLE PARAMETERS'/15X,'VERIFY THAT STEP SIZES ARE'
|
||||
1,' BIG ENOUGH AND CHECK FUN LOGIC.'/1X, 79('*')/1X,79('H')///)
|
||||
END
|
||||
Reference in New Issue
Block a user