search

fermiPy / fermipy

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

Error with array slicing in irf.py possibly linked to numpy version #515

Closed grburgess closed 1 year ago

grburgess commented 1 year ago

Describe the error I am using the FermiBottle docker with the preinstalled fermipy. I need to install a few additional packages including multinest which results in bumping numpy from 1.22 -> 1.24. This then makes it incompatible with numba and is thus downgraded to 1.23.

Before all of this, my analysis runs fine. But after the above numpy changes (the only ones that I can see have an effect) I get the following error:

File ~/threeML/threeML/plugins/FermipyLike.py:221, in _get_fermipy_instance(configuration, likelihood_model)
    218 gta = GTAnalysis(configuration)  # noqa: F821
    220 # This will take a long time if it's the first time we run with this model
--> 221 gta.setup()
    223 # Substitute all spectra for point sources with FileSpectrum, so that we will be able to control
    224 # them from 3ML
    226 energies_keV = None

File /opt/anaconda/envs/fermi/lib/python3.9/site-packages/fermipy/gtanalysis.py:1085, in GTAnalysis.setup(self, init_sources, overwrite, **kwargs)
   1083 # Run setup for each component
   1084 for i, c in enumerate(self.components):
-> 1085     c.setup(overwrite=overwrite)
   1087 # Create likelihood
   1088 self._create_likelihood()

File /opt/anaconda/envs/fermi/lib/python3.9/site-packages/fermipy/gtanalysis.py:5131, in GTBinnedAnalysis.setup(self, overwrite, **kwargs)
   5128 self._tmax = self._ltc.tstop
   5130 self.logger.debug('Creating PSF model')
-> 5131 self._psf = irfs.PSFModel.create(self.roi.skydir, self._ltc,
   5132                                  self.config['gtlike']['irfs'],
   5133                                  self.config['selection']['evtype'],
   5134                                  self.energies)
   5136 # Bin data and create exposure cube
   5137 if not use_external_srcmap:

File /opt/anaconda/envs/fermi/lib/python3.9/site-packages/fermipy/irfs.py:563, in PSFModel.create(cls, skydir, ltc, event_class, event_types, energies, cth_bins, ndtheta, use_edisp, fn, nbin)
    560 else:
    561     psf = create_avg_psf(skydir, ltc, event_class, event_types,
    562                          dtheta, energies, cth_bins)
--> 563     wts = calc_counts(skydir, ltc, event_class, event_types,
    564                       egy_bins, cth_bins, fn)
    566 exp = calc_exp(skydir, ltc, event_class, event_types,
    567                energies, cth_bins)
    568 if sum(sum(np.abs(exp))) == 0:

File /opt/anaconda/envs/fermi/lib/python3.9/site-packages/fermipy/irfs.py:872, in calc_counts(skydir, ltc, event_class, event_types, egy_bins, cth_bins, fn, npts)
    869 exp = calc_exp(skydir, ltc, event_class, event_types,
    870                egy_bins, cth_bins)
    871 dnde = fn.dnde(egy_bins)
--> 872 cnts = loglog_quad(egy_bins, exp * dnde[:, None], 0)
    873 cnts = sum_bins(cnts, 0, npts)
    874 return cnts

File /opt/anaconda/envs/fermi/lib/python3.9/site-packages/fermipy/irfs.py:51, in loglog_quad(x, y, dim)
     49 xs0[dim] = slice(None, -1)
     50 xs1[dim] = slice(1, None)
---> 51 log_ratio = np.log(x[xs1] / x[xs0])
     52 return 0.5 * (y[ys0] * x[xs0] + y[ys1] * x[xs1]) * log_ratio

IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

In [2]: %debug
> /opt/anaconda/envs/fermi/lib/python3.9/site-packages/fermipy/irfs.py(51)loglog_quad()
     49     xs0[dim] = slice(None, -1)
     50     xs1[dim] = slice(1, None)
---> 51     log_ratio = np.log(x[xs1] / x[xs0])
     52     return 0.5 * (y[ys0] * x[xs0] + y[ys1] * x[xs1]) * log_ratio
     53

ipdb> xs1
[slice(1, None, None), None]
ipdb> x
array([8.65964323e+01, 1.15478198e+02, 1.53992653e+02, 2.05352503e+02,
       2.73841963e+02, 3.65174127e+02, 4.86967525e+02, 6.49381632e+02,
       8.65964323e+02, 1.15478198e+03, 1.53992653e+03, 2.05352503e+03,
       2.73841963e+03, 3.65174127e+03, 4.86967525e+03, 6.49381632e+03,
       8.65964323e+03, 1.15478198e+04, 1.53992653e+04, 2.05352503e+04,
       2.73841963e+04, 3.65174127e+04, 4.86967525e+04, 6.49381632e+04,
       8.65964323e+04, 1.15478198e+05])

It is odd that this would be related to the bumpy version... and the operation looks completely valid.

NOTE: This is for errors in fermipy and not for questions about data analysis or issues with fermitools.

To Reproduce

After pulling fermi bottle and activating the fermi condo env:

1) conda install multinest pymultinest numba. (this is required so that all versions are consistent... though it makes no sense from the requirements of the packages)

I am using Fermipy through 3ML, an example script would be:

from threeML import *

lat_catalog = FermiLATSourceCatalog()

ra, dec, table = lat_catalog.search_around_source("PKS 2155-304", radius=10.0)

model = lat_catalog.get_model()

model.free_point_sources_within_radius(3.0, normalization_only=True)

obs = Observation(nu_obs=nu_obs, xrt_obs=xrt_obs)

###############################

evfile = "/data/soprano/data/lat/L56b77462b260f299f049d8a2adf80b8c_FT1.fits"
scfile = "/data/soprano/data/lat/L2212240528347D48B77D46_SC00.fits"

config = FermipyLike.get_basic_config(
    evfile=evfile,
    scfile=scfile,
    ra=ra,
    dec=dec,
    fermipy_verbosity=1,
    fermitools_chatter=0,
)

config.display()

LAT = FermipyLike("LAT", config)

data = DataList(LAT)

ba = BayesianAnalysis(model, data)

Expected behavior

Everything works fine BEFORE the installation of other packages which results in a change in numpy 1.22 -> 1.23

Log files

Above

Environment (please complete the following information):

grburgess commented 1 year ago

update: if one forces numpy==1.22.4 in the environment, all works.

henrikef commented 1 year ago

@hayalaso could this be related to the weirdness you were seeing in the plots of the model maps (see attached image)?

In case it is related to the same underlying problem: @hayalaso was seeing this issue with numpy 1.23.5, while my model maps looked fine with numpy 1.24.

image
hayalaso commented 1 year ago

I'm not sure. I just re-run the program I have and I still got the same problem with numpy 1.24. I tried to print a model map before the fit and I still got the same map.

This is the code I used for that. 

gta = GTAnalysis('config.yaml', logging={'verbosity':3})
gta.setup()
gta.optimize()
gta.print_roi()
gta.write_roi("optimize_model",make_plots=True)
henrikef commented 1 year ago

I see. Sorry for the confusion @hayalaso , I guess you should open a new issue then.

htajima-gamma commented 1 year ago

I ran into the same issue. Here is an error message.

/Users/htajima/opt/anaconda3/envs/fermi-root/lib/python3.9/site-packages/fermipy/irfs.py:51: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use arr[tuple(seq)] instead of arr[seq]. In the future this will be interpreted as an array index, arr[np.array(seq)], which will result either in an error or a different result. log_ratio = np.log(x[xs1] / x[xs0]) /Users/htajima/opt/anaconda3/envs/fermi-root/lib/python3.9/site-packages/fermipy/irfs.py:52: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use arr[tuple(seq)] instead of arr[seq]. In the future this will be interpreted as an array index, arr[np.array(seq)], which will result either in an error or a different result. return 0.5 (y[ys0] x[xs0] + y[ys1] x[xs1]) log_ratio 2023-01-17 14:31:03 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 56874.155220 from DATE-END'. [astropy.wcs.wcs] 2023-01-17 14:31:03 INFO GTBinnedAnalysis._create_srcmaps(): Skipping gtsrcmaps. 2023-01-17 14:31:03 INFO GTBinnedAnalysis.setup(): Finished setup for component 00 2023-01-17 14:31:03 INFO GTBinnedAnalysis._create_binned_analysis(): Creating BinnedAnalysis for component 00. 2023-01-17 14:31:11 INFO GTAnalysis.setup(): Initializing source properties 2023-01-17 14:31:15 INFO GTAnalysis.setup(): Finished setup. 2023-01-17 14:31:15 INFO GTBinnedAnalysis.write_xml(): Writing /Users/htajima/pCloud/Work/GLAST/fermipy/test/mkn421/fit0_00.xml... 2023-01-17 14:31:15 INFO GTAnalysis.write_fits(): Writing /Users/htajima/pCloud/Work/GLAST/fermipy/test/mkn421/fit0.fits... WARNING: Format %s cannot be mapped to the accepted TDISPn keyword values. Format will not be moved into TDISPn keyword. [astropy.io.fits.column] WARNING: Format %f cannot be mapped to the accepted TDISPn keyword values. Format will not be moved into TDISPn keyword. [astropy.io.fits.column] WARNING: Format %s cannot be mapped to the accepted TDISPn keyword values. Format will not be moved into TDISPn keyword. [astropy.io.fits.column] 2023-01-17 14:31:28 INFO GTAnalysis.write_roi(): Writing /Users/htajima/pCloud/Work/GLAST/fermipy/test/mkn421/fit0.npy... 2023-01-17 14:31:28 INFO GTAnalysis.optimize(): Starting Joint fit ['isodiff', 'galdiff', '3FGL J1104.4+3812', '3FGL J1112.4+3449', '3FGL J1127.8+3618'] /Users/htajima/opt/anaconda3/envs/fermi-root/lib/python3.9/site-packages/scipy/interpolate/_fitpack2.py:298: UserWarning: The maximal number of iterations maxit (set to 20 by the program) allowed for finding a smoothing spline with fp=s has been reached: s too small. There is an approximation returned but the corresponding weighted sum of squared residuals does not satisfy the condition abs(fp-s)/s < tol. warnings.warn(message)

Areustle commented 1 year ago

@grburgess, @hayalaso, The numpy 1.24 issue should have been fixed in fermipy version 1.2.0. Please use pip to install that version and try again.

FermiBottle needs to be updated to fermipy version 1.2 as well. Please open an issue at https://github.com/fermi-lat/FermiBottle/issues to encourage this.

grburgess commented 1 year ago

thanks!

grburgess commented 1 year ago

Tested and this works. Thanks again.