-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathmncont.F
280 lines (279 loc) · 8.82 KB
/
mncont.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
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
*
* $Id: mncont.F,v 1.1.1.1 2000/06/08 11:19:19 andras Exp $
*
* $Log: mncont.F,v $
* Revision 1.1.1.1 2000/06/08 11:19:19 andras
* import of MINUIT from CERNlib 2000
*
* Revision 1.1.1.1 1996/03/07 14:31:29 mclareni
* Minuit
*
*
#include "minuit/pilot.h"
SUBROUTINE MNCONT(FCN,KE1,KE2,NPTU,XPTU,YPTU,IERRF,FUTIL)
#include "minuit/d506dp.inc"
CC Find NPTU points along a contour where the function
CC FMIN (X(KE1),X(KE2)) = AMIN+UP
CC where FMIN is the minimum of FCN with respect to all
CC the other NPAR-2 variable parameters (if any).
CC IERRF on return will be equal to the number of points found:
CC NPTU if normal termination with NPTU points found
CC -1 if errors in the calling sequence (KE1, KE2 not variable)
CC 0 if less than four points can be found (using MNMNOT)
CC n>3 if only n points can be found (n < NPTU)
CC
#include "minuit/d506cm.inc"
DIMENSION XPTU(NPTU), YPTU(NPTU), W(MNI),GCC(MNI)
CHARACTER CHERE*10
PARAMETER (CHERE='MNContour ')
LOGICAL LDEBUG
EXTERNAL FCN,FUTIL
C input arguments: parx, pary, devs, ngrid
LDEBUG = (IDBG(6) .GE. 1)
IF (KE1.LE.0 .OR. KE2.LE.0) GO TO 1350
IF (KE1.GT.NU .OR. KE2.GT.NU) GO TO 1350
KI1 = NIOFEX(KE1)
KI2 = NIOFEX(KE2)
IF (KI1.LE.0 .OR. KI2.LE.0) GO TO 1350
IF (KI1 .EQ. KI2) GO TO 1350
IF (NPTU .LT. 4) GO TO 1400
C
NFCNCO = NFCN
NFCNMX = 100*(NPTU+5)*(NPAR+1)
C The minimum
CALL MNCUVE(FCN,FUTIL)
U1MIN = U(KE1)
U2MIN = U(KE2)
IERRF = 0
CFROM = CHERE
NFCNFR = NFCNCO
IF (ISW(5) .GE. 0) THEN
WRITE (ISYSWR,'(1X,A,I4,A)')
+ 'START MNCONTOUR CALCULATION OF',NPTU,' POINTS ON CONTOUR.'
IF (NPAR .GT. 2) THEN
IF (NPAR .EQ. 3) THEN
KI3 = 6 - KI1 - KI2
KE3 = NEXOFI(KI3)
WRITE (ISYSWR,'(1X,A,I3,2X,A)')
+ 'EACH POINT IS A MINIMUM WITH RESPECT TO PARAMETER ',
+ KE3, CPNAM(KE3)
ELSE
WRITE (ISYSWR,'(1X,A,I3,A)')
+ 'EACH POINT IS A MINIMUM WITH RESPECT TO THE OTHER',
+ NPAR-2, ' VARIABLE PARAMETERS.'
ENDIF
ENDIF
ENDIF
C
C Find the first four points using MNMNOT
C ........................ first two points
CALL MNMNOT(FCN,KE1,KE2,VAL2PL,VAL2MI,FUTIL)
IF (ERN(KI1) .EQ. UNDEFI) THEN
XPTU(1) = ALIM(KE1)
CALL MNWARN('W',CHERE,'Contour squeezed by parameter limits.')
ELSE
IF (ERN(KI1) .GE. ZERO) GO TO 1500
XPTU(1) = U1MIN+ERN(KI1)
ENDIF
YPTU(1) = VAL2MI
C
IF (ERP(KI1) .EQ. UNDEFI) THEN
XPTU(3) = BLIM(KE1)
CALL MNWARN('W',CHERE,'Contour squeezed by parameter limits.')
ELSE
IF (ERP(KI1) .LE. ZERO) GO TO 1500
XPTU(3) = U1MIN+ERP(KI1)
ENDIF
YPTU(3) = VAL2PL
SCALX = 1.0/(XPTU(3) - XPTU(1))
C ........................... next two points
CALL MNMNOT(FCN,KE2,KE1,VAL2PL,VAL2MI,FUTIL)
IF (ERN(KI2) .EQ. UNDEFI) THEN
YPTU(2) = ALIM(KE2)
CALL MNWARN('W',CHERE,'Contour squeezed by parameter limits.')
ELSE
IF (ERN(KI2) .GE. ZERO) GO TO 1500
YPTU(2) = U2MIN+ERN(KI2)
ENDIF
XPTU(2) = VAL2MI
IF (ERP(KI2) .EQ. UNDEFI) THEN
YPTU(4) = BLIM(KE2)
CALL MNWARN('W',CHERE,'Contour squeezed by parameter limits.')
ELSE
IF (ERP(KI2) .LE. ZERO) GO TO 1500
YPTU(4) = U2MIN+ERP(KI2)
ENDIF
XPTU(4) = VAL2PL
SCALY = 1.0/(YPTU(4) - YPTU(2))
NOWPTS = 4
NEXT = 5
IF (LDEBUG) THEN
WRITE (ISYSWR,'(A)') ' Plot of four points found by MINOS'
XPT(1) = U1MIN
YPT(1) = U2MIN
CHPT(1) = ' '
NALL = MIN(NOWPTS+1,MAXCPT)
DO 85 I= 2, NALL
XPT(I) = XPTU(I-1)
YPT(I) = YPTU(I-1)
85 CONTINUE
CHPT(2)= 'A'
CHPT(3)= 'B'
CHPT(4)= 'C'
CHPT(5)= 'D'
CALL MNPLOT(XPT,YPT,CHPT,NALL,ISYSWR,NPAGWD,NPAGLN)
ENDIF
C
C ..................... save some values before fixing
ISW2 = ISW(2)
ISW4 = ISW(4)
SIGSAV = EDM
ISTRAV = ISTRAT
DC = DCOVAR
APSI = EPSI*0.5
ABEST=AMIN
MPAR=NPAR
NFMXIN = NFCNMX
DO 125 I= 1, MPAR
125 XT(I) = X(I)
DO 130 J= 1, MPAR*(MPAR+1)/2
130 VTHMAT(J) = VHMAT(J)
DO 135 I= 1, MPAR
GCC(I) = GLOBCC(I)
135 W(I) = WERR(I)
C fix the two parameters in question
KINTS = NIOFEX(KE1)
CALL MNFIXP (KINTS,IERR)
KINTS = NIOFEX(KE2)
CALL MNFIXP (KINTS,IERR)
C ......................Fill in the rest of the points
DO 900 INEW= NEXT, NPTU
C find the two neighbouring points with largest separation
BIGDIS = 0.
DO 200 IOLD = 1, INEW-1
I2 = IOLD + 1
IF (I2 .EQ. INEW) I2 = 1
DIST = (SCALX*(XPTU(IOLD)-XPTU(I2)))**2 +
+ (SCALY*(YPTU(IOLD)-YPTU(I2)))**2
IF (DIST .GT. BIGDIS) THEN
BIGDIS = DIST
IDIST = IOLD
ENDIF
200 CONTINUE
I1 = IDIST
I2 = I1 + 1
IF (I2 .EQ. INEW) I2 = 1
C next point goes between I1 and I2
A1 = HALF
A2 = HALF
300 XMIDCR = A1*XPTU(I1) + A2*XPTU(I2)
YMIDCR = A1*YPTU(I1) + A2*YPTU(I2)
XDIR = YPTU(I2) - YPTU(I1)
YDIR = XPTU(I1) - XPTU(I2)
SCLFAC = MAX(ABS(XDIR*SCALX), ABS(YDIR*SCALY))
XDIRCR = XDIR/SCLFAC
YDIRCR = YDIR/SCLFAC
KE1CR = KE1
KE2CR = KE2
C Find the contour crossing point along DIR
AMIN = ABEST
CALL MNCROS(FCN,AOPT,IERCR,FUTIL)
IF (IERCR .GT. 1) THEN
C If cannot find mid-point, try closer to point 1
IF (A1 .GT. HALF) THEN
IF (ISW(5) .GE. 0)
+ WRITE (ISYSWR,'(A,A,I3,A)') ' MNCONT CANNOT FIND NEXT',
+ ' POINT ON CONTOUR. ONLY ',NOWPTS,' POINTS FOUND.'
GO TO 950
ENDIF
CALL MNWARN('W',CHERE,'Cannot find midpoint, try closer.')
A1 = 0.75
A2 = 0.25
GO TO 300
ENDIF
C Contour has been located, insert new point in list
DO 830 MOVE= NOWPTS,I1+1,-1
XPTU(MOVE+1) = XPTU(MOVE)
YPTU(MOVE+1) = YPTU(MOVE)
830 CONTINUE
NOWPTS = NOWPTS + 1
XPTU(I1+1) = XMIDCR + XDIRCR*AOPT
YPTU(I1+1) = YMIDCR + YDIRCR*AOPT
900 CONTINUE
950 CONTINUE
C
IERRF = NOWPTS
CSTATU = 'SUCCESSFUL'
IF (NOWPTS .LT. NPTU) CSTATU = 'INCOMPLETE'
C make a lineprinter plot of the contour
IF (ISW(5) .GE. 0) THEN
XPT(1) = U1MIN
YPT(1) = U2MIN
CHPT(1) = ' '
NALL = MIN(NOWPTS+1,MAXCPT)
DO 1000 I= 2, NALL
XPT(I) = XPTU(I-1)
YPT(I) = YPTU(I-1)
CHPT(I)= 'X'
1000 CONTINUE
WRITE (ISYSWR,'(A,I3,2X,A)') ' Y-AXIS: PARAMETER ',KE2,
+ CPNAM(KE2)
CALL MNPLOT(XPT,YPT,CHPT,NALL,ISYSWR,NPAGWD,NPAGLN)
WRITE (ISYSWR,'(25X,A,I3,2X,A)') 'X-AXIS: PARAMETER ',
+ KE1,CPNAM(KE1)
ENDIF
C print out the coordinates around the contour
IF (ISW(5) .GE. 1) THEN
NPCOL = (NOWPTS+1)/2
NFCOL = NOWPTS/2
WRITE (ISYSWR,'(/I5,A,G13.5,A,G11.3)') NOWPTS,
+ ' POINTS ON CONTOUR. FMIN=',ABEST,' ERRDEF=',UP
WRITE (ISYSWR,'(9X,A,3X,A,18X,A,3X,A)')
+ CPNAM(KE1),CPNAM(KE2),CPNAM(KE1),CPNAM(KE2)
DO 1050 LINE = 1, NFCOL
LR = LINE + NPCOL
WRITE (ISYSWR,'(1X,I5,2G13.5,10X,I5,2G13.5)')
+ LINE,XPTU(LINE),YPTU(LINE),LR,XPTU(LR),YPTU(LR)
1050 CONTINUE
IF (NFCOL .LT. NPCOL) WRITE (ISYSWR,'(1X,I5,2G13.5)')
+ NPCOL,XPTU(NPCOL),YPTU(NPCOL)
ENDIF
C . . contour finished. reset v
ITAUR = 1
CALL MNFREE(1)
CALL MNFREE(1)
DO 1100 J= 1, MPAR*(MPAR+1)/2
1100 VHMAT(J) = VTHMAT(J)
DO 1120 I= 1, MPAR
GLOBCC(I) = GCC(I)
WERR(I) = W(I)
1120 X(I) = XT(I)
CALL MNINEX (X)
EDM = SIGSAV
AMIN = ABEST
ISW(2) = ISW2
ISW(4) = ISW4
DCOVAR = DC
ITAUR = 0
NFCNMX = NFMXIN
ISTRAT = ISTRAV
U(KE1) = U1MIN
U(KE2) = U2MIN
GO TO 2000
C Error returns
1350 WRITE (ISYSWR,'(A)') ' INVALID PARAMETER NUMBERS.'
GO TO 1450
1400 WRITE (ISYSWR,'(A)') ' LESS THAN FOUR POINTS REQUESTED.'
1450 IERRF = -1
CSTATU = 'USER ERROR'
GO TO 2000
1500 WRITE (ISYSWR,'(A)') ' MNCONT UNABLE TO FIND FOUR POINTS.'
U(KE1) = U1MIN
U(KE2) = U2MIN
IERRF = 0
CSTATU = 'FAILED'
2000 CONTINUE
CFROM = CHERE
NFCNFR = NFCNCO
RETURN
END