search

fermiPy / fermipy

Fermi-LAT Python Analysis Framework
http://fermipy.readthedocs.org/
BSD 3-Clause "New" or "Revised" License
51 stars 53 forks source link

ROI optimization and fit not working on Mac OS with scipy > 1.10.1 #566

Open ranieremenezes opened 7 months ago

ranieremenezes commented 7 months ago

When running fermipy 1.2.2 on a Mac OS, the functions optimize() and fit() from GTAnalysis crash with the error "ValueError: The function value at x=0.0 is NaN; solver cannot continue." (full log below). This problem does not happen if the analysis is done on Linux or if we install scipy v1.10.1 in the fermipy environment.

My analysis was the following: Mrk 421 time window: START, 14/10/2008 15:43:00 energy window: 100, 300000

Full log below:

fermipy version 1.2.2 ScienceTools version 2.2.0 2024-02-15 15:53:49 INFO GTAnalysis.setup(): Running setup. 2024-02-15 15:53:49 INFO GTBinnedAnalysis.setup(): Running setup for component 00 2024-02-15 15:53:49 INFO GTBinnedAnalysis._select_data(): Skipping data selection. 2024-02-15 15:53:49 INFO GTBinnedAnalysis._create_ltcube(): Skipping LT Cube. 2024-02-15 15:53:50 INFO GTBinnedAnalysis._create_expcube(): Skipping gtexpcube. WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs] WARNING: FITSFixedWarning: 'datfix' made the change 'Set DATEREF to '2001-01-01T00:01:04.184' from MJDREF. Set MJD-OBS to 54682.655283 from DATE-OBS. Set MJD-END to 54753.65486 from DATE-END'. [astropy.wcs.wcs] 2024-02-15 15:53:50 INFO GTBinnedAnalysis._create_srcmaps(): Skipping gtsrcmaps. 2024-02-15 15:53:50 INFO GTBinnedAnalysis.setup(): Finished setup for component 00 2024-02-15 15:53:50 INFO GTBinnedAnalysis._create_binned_analysis(): Creating BinnedAnalysis for component 00. Drm_Cache::update Measured counts < 0 4FGL J1204.4+3750 27 -7.52459e-18 1.38686e-17 0.129667 0.116145 0.0848453 0.0519819 0.0279117 0.0123597 0.00530023 0.00184387 0.000664404 0.000191539 6.01258e-05 1.55422e-05 4.04127e-06 8.71618e-07 2.15058e-07 5.05724e-08 1.29549e-08 3.12197e-09 7.25482e-10 1.58131e-10 3.10266e-11 5.43056e-12 8.23529e-13 1.09763e-13 1.31946e-14 1.36786e-15 1.43733e-16 1.38686e-17 2024-02-15 15:54:06 INFO GTAnalysis.setup(): Initializing source properties 2024-02-15 15:54:11 INFO GTAnalysis.setup(): Finished setup. 2024-02-15 15:54:11 INFO GTAnalysis.optimize(): Starting Joint fit ['isodiff', 'galdiff', '4FGL J1104.4+3812', '4FGL J1131.0+3815']

ValueError Traceback (most recent call last)

File ~/Documents/Mrk421/Mrk421/BinnedAnalysis_so_ate_o_fit.py:31 28 h[0].data = np.sum(counts,axis=0) 29 h.writeto('./cmap.fits',overwrite=True) ---> 31 gta.optimize(npred_threshold=500,shape_ts_threshold=200) 33 gta.write_roi('roi1') 35 srcs = gta.find_sources(sqrt_ts_threshold=8.0, min_separation=0.5,multithread=True)

File ~/miniforge3/envs/fermipy/lib/python3.9/site-packages/fermipy/gtanalysis.py:2339, in GTAnalysis.optimize(self, kwargs) 2336 break 2338 print ("Joint fit ", joint_norm_fit) -> 2339 self.fit(loglevel=logging.DEBUG, config['optimizer']) 2340 self.free_sources(free=False, loglevel=logging.DEBUG) 2342 # Step through remaining sources and re-fit normalizations 2343 # FIXME, EAC, use npred_wt here

File ~/miniforge3/envs/fermipy/lib/python3.9/site-packages/fermipy/gtanalysis.py:3074, in GTAnalysis.fit(self, update, **kwargs) 3072 if len(freePars) == 0: 3073 continue -> 3074 self.update_source(name, reoptimize=config['reoptimize']) 3076 # Update roi model counts 3077 self._update_roi()

File ~/miniforge3/envs/fermipy/lib/python3.9/site-packages/fermipy/gtanalysis.py:4106, in GTAnalysis.update_source(self, name, paramsonly, reoptimize, **kwargs) 4103 npts = self.config['gtlike']['llscan_npts'] 4104 optimizer = kwargs.get('optimizer', self.config['optimizer']) -> 4106 sd = self.get_src_model(name, paramsonly, reoptimize, npts, 4107 optimizer=optimizer) 4108 src = self.roi.get_source_by_name(name) 4109 src.update_data(sd)

File ~/miniforge3/envs/fermipy/lib/python3.9/site-packages/fermipy/gtanalysis.py:4264, in GTAnalysis.get_src_model(self, name, paramsonly, reoptimize, npts, **kwargs) 4261 pass 4262 # self.logger.error('Failed to update source parameters.', 4263 # exc_info=True) -> 4264 lnlp = self.profile_norm(name, savestate=True, 4265 reoptimize=reoptimize, npts=npts, 4266 optimizer=optimizer) 4268 src_dict['loglike_scan'] = lnlp['loglike'] 4269 src_dict['dloglike_scan'] = lnlp['dloglike']

File ~/miniforge3/envs/fermipy/lib/python3.9/site-packages/fermipy/gtanalysis.py:2439, in GTAnalysis.profile_norm(self, name, logemin, logemax, reoptimize, xvals, npts, fix_shape, savestate, **kwargs) 2436 xvals = self._find_scan_pts(name, npts=9) 2437 lnlp = self.profile(name, parName, 2438 reoptimize=False, xvals=xvals) -> 2439 lims = utils.get_parameter_limits(lnlp['xvals'], 2440 lnlp['dloglike'], 2441 cl_limit=0.99) 2443 if lnlp['npred'].sum() < 0.1: 2444 self.logger.warning("Source model for source %s has no counts in ROI" % name)

File ~/miniforge3/envs/fermipy/lib/python3.9/site-packages/fermipy/utils.py:804, in get_parameter_limits(xval, loglike, cl_limit, cl_err, tol, bounds) 802 # Refine the peak position 803 if np.sign(sd(xval[ilo])) != np.sign(sd(xval[ihi])): --> 804 x0 = find_function_root(sd, xval[ilo], xval[ihi]) 806 if np.isnan(x0): 807 x0 = xval[imax]

File ~/miniforge3/envs/fermipy/lib/python3.9/site-packages/fermipy/utils.py:700, in find_function_root(fn, x0, xb, delta, bounds) 697 xtol = 1e-10 * np.abs(xb + x0) 699 try: --> 700 return brentq(lambda t: fn(t) + delta, x0, xb, xtol=xtol) 701 except RuntimeError: 702 raise RuntimeError("Brentq failed ", x0, xb, fn(x0)+delta, fn(xb)+delta, xtol, delta)

File ~/miniforge3/envs/fermipy/lib/python3.9/site-packages/scipy/optimize/_zeros_py.py:802, in brentq(f, a, b, args, xtol, rtol, maxiter, full_output, disp) 800 raise ValueError(f"rtol too small ({rtol:g} < {_rtol:g})") 801 f = _wrap_nan_raise(f) --> 802 r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp) 803 return results_c(full_output, r)

File ~/miniforge3/envs/fermipy/lib/python3.9/site-packages/scipy/optimize/_zeros_py.py:98, in _wrap_nan_raise..f_raise(x, *args) 96 err._x = x 97 err._function_calls = f_raise._function_calls ---> 98 raise err 99 return fx

ValueError: The function value at x=0.0 is NaN; solver cannot continue.

clodoN1109 commented 7 months ago

Based on the log, NANs are being generated in fermipy/utils.py, declaration sd = spline.derivative(), which is the function passed to find_function_root, which calls brentq, that fails when finds NANs.

image

spline.derivative is a scipy method from the UnivariateSpline class.

Comparing versions, it's found that for Scipy > 1.10.1 the brentq method for finding the zeros of a function implements a test to determine whether the function evaluates to NAN (NotANumber). Although version 1.10.1 may be running without errors, it may be only because it is not testing for NANs.

image

I hope that this information helps.

omodei commented 6 months ago

Thanks! The error is raised by brentq, but the issue originates when we define the parameters range, in some cases adding "0" as the first element of the scan causes a crash. I am testing a solution and it should be fixed in the next merge.