-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathmnsimp.F
206 lines (206 loc) · 6.3 KB
/
mnsimp.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
*
* $Id: mnsimp.F,v 1.1.1.1 2000/06/08 11:19:20 andras Exp $
*
* $Log: mnsimp.F,v $
* Revision 1.1.1.1 2000/06/08 11:19:20 andras
* import of MINUIT from CERNlib 2000
*
* 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