SUBROUTINE clcunc(ear, esavar, fluxar, bflxar, respon, & igroup, ichanb, ichane, ftstat, & effar, iepar, parlmt, adelstt, ierflg, & qpeg, iflag, toler, nertry, delstt) IMPLICIT NONE REAL esavar(*), fluxar(*), bflxar(*) REAL effar(*), respon(*), ear(0:*) REAL parlmt(2), adelstt(2) REAL ftstat, toler, delstt INTEGER ierflg(2), nertry, iepar INTEGER igroup(*), ichanb(*), ichane(*) LOGICAL qpeg(2) c Routine to calculate the uncertainty for parameter IEPAR. c Arguments : c ear r i: Energies c esavar r r: Model component spectra c fluxar r r: Total model spectrum c bflxar r r: Total model background spectrum c respon r i: Response matrix elements c igroup i i: Number of response groups for each energy c ichanb i i: Start channel of response group c ichane i i: End channel of response group c ftstat i i/r: Fit statistic c effar r i: Effective areas c iepar i i: Model parameter under consideration c parlmt r r: Parameter deltas found c adelstt r r: Delta statistic for parameter deltas c ierflg i r: Error flag for each direction c qpeg l r: True if parameter pegged c iflag i r: Return flag - 0 == OK c -1 == EOF entered by user c 1 == parameter is frozen c 2 == new minimum found c toler r i: Accuracy to which we need to hit delta stat. c nertry i i: Number of tries before prompting user c delstt r i: Target delta statistic c kaa 4/26/94 INCLUDE '../../inc/xspec.inc' DOUBLE PRECISION parasf(MAXFPR), pval REAL asgn, hrdval, hrdlmt, curdel, xtrial, cpar REAL cdel, csig, sttsav, xlimit(2), ylimit(2) REAL random INTEGER ieftpr, ieigen, ieval INTEGER ifpar, iflag, nfree INTEGER isgn, lenact, icerr, neigen INTEGER iertry, ierror, ierr, seed(2) LOGICAL qdone, qeof, qfintp, qbrack, qfirst, qcont CHARACTER wrtstr*255, ctime*8 EQUIVALENCE(ctime, seed(1)) DOUBLE PRECISION fgpval, xgcrdc, mgpval INTEGER fgnfpr, mgfpar, dgndst, dgnbne, dgdpnt, fgneig CHARACTER xgstat*11 LOGICAL mgqfrz, xgprft, xgaskq, xgansw, fgqfrz EXTERNAL fgpval, mgpval, xgstat, xgcrdc, xgprft, xgansw EXTERNAL fgnfpr, mgfpar, mgqfrz, xgaskq, dgndst, dgnbne EXTERNAL fgqfrz, dgdpnt, fgneig iflag = 0 c Initialize the random number generator in case we need to use it later CALL gettim(ctime) seed(1) = seed(1) + seed(2) CALL rluxgo(3, seed(1), 0, 0) c Save the old fit parameter values (these are the UNADJUSTED values) DO ifpar = 1, fgnfpr() parasf(ifpar) = fgpval(ifpar,'a') ENDDO sttsav = ftstat c Get the memory to save the eigenvector data then save it neigen = fgneig() CALL udmget(neigen, 4, ieftpr, ierr) IF ( ierr .NE. 0 ) CALL xwrite( & ' WARNING: Failed to allocate memory for eftpr in CLCUNC.', 2) CALL udmget(neigen**2, 7, ieigen, ierr) IF ( ierr .NE. 0 ) CALL xwrite( & ' WARNING: Failed to allocate memory for eigen in CLCUNC.', 2) CALL udmget(neigen, 7, ieval, ierr) IF ( ierr .NE. 0 ) CALL xwrite( & ' WARNING: Failed to allocate memory for eval in CLCUNC.', 2) CALL fgeftp(neigen, MEMI(ieftpr), ierr) CALL fgevec(neigen, MEMD(ieigen), ierr) CALL fgeval(neigen, MEMD(ieval), ierr) IF ( ierr .NE. 0 ) CALL xwrite( & ' WARNING: Failed to load eigenvalue data in CLCUNC.', 2) c Get the best fit value and sigma for the chosen parameter cpar = SNGL(mgpval(iepar, 'a')) IF ( .NOT. mgqfrz(iepar) ) THEN cdel = SNGL(mgpval(iepar,'d')) CALL mpqfrz(iepar, ierr) IF ( ierr .NE. 0 ) & CALL xwrite(' CLCUNC : error detected from MPQFRZ', 15) ELSE WRITE (wrtstr, '(a,i4,a)') & ' *ERROR*:CLCUNC: Error for a frozen parameter,' & , iepar, ', is requested.' CALL xwrite(wrtstr, 2) iflag = 1 RETURN ENDIF csig = SNGL(mgpval(iepar,'s')) IF (csig.LE.0.) THEN WRITE (wrtstr, '(a,f12.4)') & ' *WARNING*: Parameter sigma indicates possible error:' & , csig CALL xwrite(wrtstr, 5) csig = cdel ENDIF c Calculate the number of free parameters nfree = 0 DO ifpar = 1, fgnfpr() IF ( .NOT.fgqfrz(ifpar) ) nfree = nfree + 1 ENDDO c Calculate uncertainty in both directions DO isgn = 1, 2 qpeg(isgn) = .FALSE. asgn = isgn*2 - 3 c Note that we can only check against the hard limit if the soft and hard c limits are identical. If they are not then the adjusted parameter value c (which is used in this routine) will lie between +/- infinity. IF (isgn.EQ.2) THEN hrdval = SNGL(mgpval(iepar,'h')) hrdlmt = 1.e32 IF ( hrdval .EQ. mgpval(iepar, 't') ) hrdlmt = hrdval - cpar ELSE hrdval = SNGL(mgpval(iepar,'l')) hrdlmt = 1.e32 IF ( hrdval .EQ. mgpval(iepar, 'b') ) hrdlmt = cpar - hrdval ENDIF qdone = hrdlmt .LE. 0. IF ( qdone ) THEN WRITE (wrtstr, '(a,1pg13.4,a)') & ' Parameter pegged at hard limit ', hrdval, & ' with delta ftstat= 0.' CALL xwrite(wrtstr, 5) xtrial = 0. IF ( isgn .EQ. 1 ) CALL mpecnf(iepar, 4, ierr) IF ( isgn .EQ. 2 ) CALL mpecnf(iepar, 5, ierr) ENDIF qfirst = .TRUE. qbrack = .FALSE. iertry = 0 c loop round iterations DO WHILE ( .NOT.qdone ) iertry = iertry + 1 c calculate the next offset. IF ( qfirst ) THEN c For the first trial set assume that the fit statistic is parabolic c through the minimum and the curvature given by the csig trial value. c If this takes it outside the hardlimit then use min of parameter c delta and hardlimit/2 (use hardlimit/2 and not hardlimit to avoid the c case when the parameter is a normalization - setting it to the hardlimit c would likely be zero which is not a good place to try and start). xtrial = csig*sqrt(delstt) IF( xtrial .GT. hrdlmt ) xtrial = min(cdel, hrdlmt/2.) qfirst = .FALSE. qfintp = .TRUE. c If not yet bracketed then just step out ELSEIF ( .NOT.qbrack ) THEN xtrial = MIN(2*xtrial, hrdlmt) WRITE (wrtstr, '(a,1pg13.6)') & ' Bracketing... trial = ', cpar+asgn*xtrial CALL xwrite(wrtstr, 15) c otherwise do an interpolation to find the next point ELSE CALL clcxtr(qfintp, delstt, curdel, csig, & hrdlmt, xtrial, xlimit(1), xlimit(2), & ylimit(1), ylimit(2), icerr) c error 1 in clcxtr is if the parameter is pegged at the hard limit. IF (icerr .EQ. 1) THEN qdone = .TRUE. qpeg(isgn) = .TRUE. WRITE (wrtstr, '(a,1pg13.6)') & ' Parameter pegged at hard limit ', hrdval CALL xwrite(wrtstr, 5) WRITE (wrtstr, '(a,1pg13.6)') & ' with delta ftstat=', curdel CALL xwrite(wrtstr, 5) IF ( isgn .EQ. 1 ) CALL mpecnf(iepar, 4, ierr) IF ( isgn .EQ. 2 ) CALL mpecnf(iepar, 5, ierr) GOTO 200 c error 2 is if the latest fit statistic value is no different from c the last. Try to recover by resetting the trial value to a random c point between the previous values. ELSEIF (icerr .EQ. 2) THEN CALL ranlux(random, 1) xtrial = xlimit(1) + (xlimit(2)-xlimit(1))*random ENDIF WRITE (wrtstr, '(a,1pg13.6)') & ' Interpolating... trial = ', cpar+asgn*xtrial CALL xwrite(wrtstr, 15) ENDIF c set new parameter value using the new offset. Other parameters are c reinitialized to the starting values - this should help avoid the c fit getting trapped in local minima. DO ifpar = 1, fgnfpr() CALL fppval(ifpar, parasf(ifpar), 'a', ierr) IF ( ierr .NE. 0 ) & CALL xwrite(' CLCUNC : error detected from FPPVAL', & 15) ENDDO c set the new parameter value - note the check that the value lies within the c hard limits, this avoids a potential problem with numerical accuracy if c the hard limit is very much smaller than the best fit value pval = DBLE(cpar+asgn*xtrial) IF ( mgpval(iepar, 'h') .EQ. mgpval(iepar, 't') ) & pval = MIN(pval, mgpval(iepar, 'h')) IF ( mgpval(iepar, 'l') .EQ. mgpval(iepar, 'b') ) & pval = MAX(pval, mgpval(iepar, 'l')) CALL mppval(iepar, pval, 'a', ierr) CALL evlmod(ear, .true., esavar, fluxar, bflxar, effar) CALL fldmod(fluxar, bflxar, igroup, ichanb, & ichane, respon, MEMR(dgdpnt('mod')), effar) CALL calstat(MEMR(dgdpnt('mod')), MEMR(dgdpnt('src')), & MEMR(dgdpnt('csv')), MEMR(dgdpnt('cor')), & MEMR(dgdpnt('bkg')), MEMR(dgdpnt('cbv')), & MEMR(dgdpnt('reg')), MEMR(dgdpnt('brg')), & MEMR(dgdpnt('are')), MEMR(dgdpnt('bre')), & ftstat) c if there is no free parameters other than iepar then we don't c need to do a fit for the other parameters IF ( nfree .EQ. 0 ) THEN ierror = 0 ELSE c renorm if we are automatically doing this before fits IF ( xgprft() ) THEN CALL renorm(MEMR(dgdpnt('src')), MEMR(dgdpnt('mod')), & MEMR(dgdpnt('csv')), MEMR(dgdpnt('cor')), & MEMR(dgdpnt('bkg')), MEMR(dgdpnt('cbv')), & MEMR(dgdpnt('reg')), MEMR(dgdpnt('brg')), & MEMR(dgdpnt('are')), MEMR(dgdpnt('bre')), & fluxar, bflxar, igroup, ichanb, & ichane, respon, esavar, ftstat, & ear, effar) ENDIF c minimize the other parameters for the model statistic CALL ftdata(ear, esavar, fluxar, bflxar, & respon, igroup, ichanb, ichane, ftstat, & .TRUE., ierror, effar) ENDIF c check whether we have found a new minimum IF ((sttsav-ftstat).GT.SNGL(xgcrdc())) THEN WRITE (wrtstr, '(a11,a,i5,a,1pg13.6)') & xgstat(), ' when model parameter ', & iepar, '=', mgpval(iepar,'v') CALL xwrite(wrtstr, 2) WRITE (wrtstr, '(a,f10.4,a,f10.4)') ' is ', & ftstat, ', which is < previous minimum', & sttsav CALL xwrite(wrtstr, 2) WRITE (wrtstr, '(a,f10.4,a)') & '(critical delta = ', xgcrdc(), ')' CALL xwrite(wrtstr, 2) ierror = -2 CALL mpecnf(iepar, 1, ierr) ENDIF c if new minimum or entered from ftdata then abort back to command level IF (ierror.LT.0) THEN GOTO 400 ENDIF c if the ftdata call returned an error 5 (number of flamda increasing c iterations exceeded the maximum) or an error 5 (flamda was about to c overflow) then we are probably OK so carry on but issue a warning. IF ( ierror .EQ. 5 .OR. ierror .EQ. 7 ) THEN CALL xwrite( & ' Minimization may have run into a problem, check your result', & 10) CALL mpecnf(iepar, 3, ierr) c on other ftdata errors then go onto other direction or next parameter ELSEIF ( ierror .GT. 0 ) THEN WRITE (wrtstr, '(a,i2)') & ' Minimization failed, crfitx error = ', ierror CALL xwrite(wrtstr, 10) CALL xwrite(' giving up on this direction or parameter', & 10) qdone = .TRUE. xtrial = 2.*hrdlmt IF ( isgn .EQ. 1 ) CALL mpecnf(iepar, 7, ierr) IF ( isgn .EQ. 2 ) CALL mpecnf(iepar, 8, ierr) GOTO 200 ENDIF c If curdel < 0 but by less than the critical amount (so not trapped above) c then set curdel = 0. Not sure that this is the right way of doing this c but trying to avoid the problem of error crashing out with new minimum c when the fit statistic is changed by a very small amount. curdel = max(0., (ftstat - sttsav)) c if we have not yet bracketed the target then see whether this attempt c managed it. If it did then set the brackets and go onto the step of c converging on the solution. Otherwise test whether we are at a hard limit c in which case give up. IF ( .NOT.qbrack ) THEN IF ( ftstat .GT. (sttsav + delstt) ) THEN qbrack = .TRUE. xlimit(1) = 0. xlimit(2) = xtrial ylimit(1) = 0. ylimit(2) = curdel WRITE(wrtstr,'(''Successful bracket : '', 1p2g13.6)') & xtrial*asgn, curdel CALL xwrite(wrtstr, 11) ELSE IF ( xtrial .GE. hrdlmt ) THEN qdone = .TRUE. qpeg(isgn) = .TRUE. WRITE (wrtstr, '(a,1pg13.6)') & ' Parameter pegged at hard limit ', hrdval CALL xwrite(wrtstr, 5) WRITE (wrtstr, '(a,1pg13.6)') & ' with delta ftstat=', curdel CALL xwrite(wrtstr, 5) IF ( isgn .EQ. 1 ) CALL mpecnf(iepar, 4, ierr) IF ( isgn .EQ. 2 ) CALL mpecnf(iepar, 5, ierr) GOTO 200 ENDIF ENDIF c if the target has already been bracketed then we are in the convergence c phase so check for convergence and too many iterations. ELSE IF (abs(curdel-delstt) .LE. toler ) THEN qdone = .TRUE. ELSE WRITE (wrtstr, '(i4,2x,1p2g13.6)') iertry, & asgn*xtrial, curdel CALL xwrite(wrtstr, 11) WRITE (wrtstr, '(6x,1p4g13.6)') & asgn*xlimit(1), ylimit(1), & asgn*xlimit(2), ylimit(2) CALL xwrite(wrtstr, 11) C If the calculated statistic for the trial value exceeds the upper limit, C ylimit(2) or is less than the lower limit, ylimit(1), then there is probably c something wrong and the user should use steppar to investigate the c situation. Return a value midway between the bracketing trials. IF ( curdel .GT. ylimit(2) .OR. & curdel .LT. ylimit(1) ) THEN WRITE(wrtstr, '(a,a11,a)') & ' Apparent non-monotonicity in ', xgstat(), ' space detected' CALL xwrite(wrtstr, 5) WRITE (wrtstr, '(a26,1p2g13.6)') & ' Current bracket values : ', asgn*xlimit(1)+cpar, & asgn*xlimit(2)+cpar CALL xwrite(wrtstr, 5) WRITE (wrtstr, '(a11,a11,a4,1p2g13.6)') & ' and delta ', xgstat(), ' : ', ylimit(1), ylimit(2) CALL xwrite(wrtstr, 5) WRITE (wrtstr, '(a18,1pg13.6,a7,a11,a3,1pg13.6)') & ' but latest trial ', asgn*xtrial+cpar, ' gives ', & xgstat(), ' = ', curdel CALL xwrite(wrtstr, 5) CALL xwrite( & ' Suggest that you check this result using the steppar command', & 5) xtrial = (xlimit(1)+xlimit(2))/2. curdel = (ylimit(1)+ylimit(2))/2. qdone = .TRUE. CALL mpecnf(iepar, 2, ierr) GOTO 200 ENDIF c If the tolerance has not been achieved but the top and bottom of the xlimit c range are equal then we cannot calculate the delta to the required tolerance c so have to give up. IF ( ABS(xlimit(1)-xlimit(2))/xlimit(2) .LT. 1e-5 ) & THEN WRITE(wrtstr, '(a)') & ' Warning: identical values of the parameter give different' CALL xwrite(wrtstr, 10) WRITE(wrtstr, '(a)') & ' values of the statistic. Please check your result' CALL xwrite(wrtstr, 10) IF ( isgn .EQ. 1 ) THEN WRITE(wrtstr, '(a)') & ' for the low end of the confidence range.' CALL xwrite(wrtstr, 10) CALL mpecnf(iepar, 7, ierr) ELSE WRITE(wrtstr, '(a)') & ' for the high end of the confidence range.' CALL xwrite(wrtstr, 10) CALL mpecnf(iepar, 8, ierr) ENDIF qdone = .TRUE. GOTO 200 ENDIF c if too many times round the loop give the user a chance to get out IF ( iertry .GE. nertry ) THEN CALL xwrtpr( & ' Number of trials exceeded before convergence.' & ) WRITE (wrtstr, '(a,1p2g13.6)') & ' Current trial values : ', asgn*xlimit(1)+cpar, & asgn*xlimit(2)+cpar CALL xwrtpr(wrtstr(:lenact(wrtstr))) WRITE (wrtstr, '(a11,a11,a3,1p2g13.6)') & ' and delta ', xgstat(), ' : ', ylimit(1), ylimit(2) CALL xwrtpr(wrtstr(:lenact(wrtstr))) IF ( xgaskq() ) THEN qcont = .false. CALL xquest( & 'Continue error search in this direction ?', 'n', & qcont, qeof) IF ( qeof ) THEN ierror = -1 IF ( isgn .EQ. 1 ) THEN CALL mpecnf(iepar, 7, ierr) ELSE CALL mpecnf(iepar, 8, ierr) ENDIF GOTO 400 ENDIF qdone = .NOT.qcont ELSE IF ( .NOT.xgansw() ) qdone = .TRUE. ENDIF iertry = 0 IF ( qdone ) THEN IF ( isgn .EQ. 1 ) THEN CALL mpecnf(iepar, 7, ierr) ELSE CALL mpecnf(iepar, 8, ierr) ENDIF ENDIF ENDIF ENDIF ENDIF c Come from (break condition) in do while loop 200 CONTINUE c End loop over iterations ENDDO c save the parameter limit - note the check that the value lies within the c hard limits, this avoids a potential problem with numerical accuracy if c the hard limit is very much smaller than the best fit value pval = DBLE(asgn*xtrial + cpar) CALL mppval(iepar, pval, 'a', ierr) pval = mgpval(iepar, 'v') pval = MIN(pval, mgpval(iepar, 'h')) pval = MAX(pval, mgpval(iepar, 'l')) parlmt(isgn) = SNGL(pval) adelstt(isgn) = curdel ierflg(isgn) = ierror c the answer has converged so save the value and reset the values of the c fit parameters to the values at invocation DO ifpar = 1, fgnfpr() CALL fppval(ifpar, parasf(ifpar), 'a', ierr) IF ( ierr .NE. 0 ) & CALL xwrite(' CLCUNC : error detected from FPPVAL', 15) ENDDO CALL evlmod(ear, .TRUE., esavar, fluxar, bflxar, effar) CALL fldmod(fluxar, bflxar, igroup, ichanb, ichane, & respon, MEMR(dgdpnt('mod')), effar) c end loop over directions ENDDO c abort point from new minimum or from user. 400 CONTINUE CALL mpqthw(iepar, ierr) IF ( ierr .NE. 0 ) & CALL xwrite(' CLCUNC : error detected from MPQTHW', 15) c if ierror=-2 then new minimum so leave parameters as they are otherwise c restore to condition at start of routine. IF ( ierror .NE. -2 ) THEN DO ifpar = 1, fgnfpr() CALL fppval(ifpar, parasf(ifpar), 'a', ierr) IF ( ierr .NE. 0 ) & CALL xwrite(' CLCUNC : error detected from FPPVAL', 15) ENDDO ftstat = sttsav ENDIF CALL evlmod(ear, .TRUE., esavar, fluxar, bflxar, effar) CALL fldmod(fluxar, bflxar, igroup, ichanb, ichane, respon, & MEMR(dgdpnt('mod')), effar) CALL fpeftp(neigen, MEMI(ieftpr)) CALL fpevec(neigen, MEMD(ieigen)) CALL fpeval(neigen, MEMD(ieval)) c release the memory CALL udmfre(ieftpr, 4, ierr) IF ( ierr .NE. 0 ) CALL xwrite( & ' WARNING: Failed to free memory for eftpr in CLCUNC.', 2) CALL udmfre(ieigen, 7, ierr) IF ( ierr .NE. 0 ) CALL xwrite( & ' WARNING: Failed to free memory for eigen in CLCUNC.', 2) CALL udmfre(ieval, 7, ierr) IF ( ierr .NE. 0 ) CALL xwrite( & ' WARNING: Failed to free memory for eval in CLCUNC.', 2) c Set the returned flag to -1 for an EOF and 2 for a new minimum IF ( ierror .EQ. -1 ) iflag = -1 IF ( ierror .EQ. -2 ) iflag = 2 RETURN END