Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add 3D fit example using gammapy.maps #1381

Merged
merged 4 commits into from Apr 18, 2018
Merged

Add 3D fit example using gammapy.maps #1381

merged 4 commits into from Apr 18, 2018

Conversation

@joleroi
Copy link
Contributor

@joleroi joleroi commented Apr 17, 2018

This PR modifies the example script for 3D fitting

  • Use gammapy.maps
  • Use the iminuit backend added in #1335

The fit converges but is quite slow and requires some tweaking of the input parameters. Most annoyingly I had to play a lot with the minimum of the amplitude parameter to keep the fit from setting amplitude=0. I guess adapting the minuit stepsize could already help, I will look into this.

UPDATE: After ee20503 the covariance matrix is used to determine the inital stepsize in the Minuit fit. This solves the problems described above.

INFO: Starting values
ParameterList
Parameter(name='lon_0', value=0.0, unit='deg', min=nan, max=nan, frozen=False)
Parameter(name='lat_0', value=0.0, unit='deg', min=nan, max=nan, frozen=False)
Parameter(name='sigma', value=1.0, unit='deg', min=nan, max=nan, frozen=False)
Parameter(name='index', value=2, unit='', min=nan, max=nan, frozen=False)
Parameter(name='amplitude', value=1e-10, unit='1 / (cm2 s TeV)', min=1e-12, max=nan, frozen=False)
Parameter(name='reference', value=1.0, unit='TeV', min=nan, max=nan, frozen=True)

Covariance: 
None [__main__]
INFO: Best fit values
ParameterList
Parameter(name='lon_0', value=0.19999767344948793, unit='deg', min=nan, max=nan, frozen=False)
Parameter(name='lat_0', value=0.0999974320886449, unit='deg', min=nan, max=nan, frozen=False)
Parameter(name='sigma', value=0.19999884967532772, unit='deg', min=nan, max=nan, frozen=False)
Parameter(name='index', value=3.000575227144423, unit='', min=nan, max=nan, frozen=False)
Parameter(name='amplitude', value=1.0008549544693324e-11, unit='1 / (cm2 s TeV)', min=1e-12, max=nan, frozen=False)
Parameter(name='reference', value=1.0, unit='TeV', min=nan, max=nan, frozen=True)

Covariance: 
None [__main__]
INFO: True values
ParameterList
Parameter(name='lon_0', value=0.2, unit='deg', min=nan, max=nan, frozen=False)
Parameter(name='lat_0', value=0.1, unit='deg', min=nan, max=nan, frozen=False)
Parameter(name='sigma', value=0.2, unit='deg', min=nan, max=nan, frozen=False)
Parameter(name='index', value=3, unit='', min=nan, max=nan, frozen=False)
Parameter(name='amplitude', value=1e-11, unit='1 / (cm2 s TeV)', min=nan, max=nan, frozen=False)
Parameter(name='reference', value=1.0, unit='TeV', min=nan, max=nan, frozen=True)

Covariance: 
None [__main__]

real    0m19.183s
user    0m18.992s
sys     0m0.185s
@joleroi joleroi added the cleanup label Apr 17, 2018
@joleroi joleroi added this to the 0.8 milestone Apr 17, 2018
@joleroi joleroi self-assigned this Apr 17, 2018
def load_cubes():
npred_cube = WcsNDMap.read('npred.fits')
exposure_cube = WcsNDMap.read('exposure.fits')
i_nan = np.where(np.isnan(exposure_cube.data))

This comment has been minimized.

@cdeil

cdeil Apr 17, 2018
Member

This shouldn't be necessary here, exposure shouldn't have nan.
Is this for pixels outside the max FOV offset cut?
Then we should just fill those with 0 instead of NaN.
Can you please either fix properly, or leave a TODO comment here to remember to make this change when filling exposure?

This comment has been minimized.

@joleroi

joleroi Apr 17, 2018
Author Contributor

I copied this from the old script, don't know what it does

This comment has been minimized.

@cdeil

cdeil Apr 17, 2018
Member

The exposure map is now filled properly, has no nan pixels:

>>> import numpy as np
>>> from gammapy.maps import Map
>>> m = Map.read("examples/example_3d_simulate/exposure.fits")
>>> np.isnan(m.data).sum()
0

So just remove these two old lines that do nan handling.

Parameters
----------
cubes : dict

This comment has been minimized.

@cdeil

cdeil Apr 17, 2018
Member

Why pass some inputs in a dict, but model separately?
Soon we'll pass other things in, e.g. PSF.
I'd suggest to just pass in all objects separately that we need here, at least for now as long as we have no "dataset" container for this.

This comment has been minimized.

@joleroi

joleroi Apr 17, 2018
Author Contributor

👍


def fit(self):
"""Run the fit"""
parameters, minuit = fit_minuit(parameters=self.model.parameters,

This comment has been minimized.

@cdeil

cdeil Apr 17, 2018
Member

Suggest to expose the Minuit object, at least for now for debugging having access to it for advanced operations (like hesse to compute parameter errors) or to check why an optimisaiton failed would be useful.

This comment has been minimized.

@joleroi

joleroi Apr 18, 2018
Author Contributor

Done

class CubeFit(object):
"""Perform 3D likelihood fit
This is my first go at such a class. It's geared to the SpectrumFit class

This comment has been minimized.

@cdeil

cdeil Apr 17, 2018
Member

Suggest to create a Sphinx link to the other class to make reading the docs more quick, i.e.

~gammapy.spectrum.SpectrumFit
@@ -76,5 +76,8 @@ def make_minuit_kwargs(parameters):
kwargs[par.name] = par.value
if par.frozen:
kwargs['fix_{}'.format(par.name)] = True
limits = par.parmin, par.parmax
limits = np.where(np.isnan(limits), None, limits)
kwargs['limit_{}'.format(par.name)] = limits

This comment has been minimized.

@cdeil

cdeil Apr 17, 2018
Member

Can you also pass the err_ options to Minuit?

In Minuit, these parameter errors are used to create the initial covariance matrix that the fitter then uses.
I.e. passing err_ is the usual way to set the "scale" of the parameters and make the fit converge or faster.
E.g. for fluxes it's usually 1e-11 but for other parameters it can be ~ 1.

This comment has been minimized.

@joleroi

joleroi Apr 17, 2018
Author Contributor

That's what I meant by initial stepsize

This comment has been minimized.

@joleroi

joleroi Apr 18, 2018
Author Contributor

Done in ee20503

@cdeil cdeil added feature and removed cleanup labels Apr 17, 2018
@cdeil cdeil changed the title Update 3D fitting example script Add 3D fit example using gammapy.maps Apr 17, 2018
@joleroi joleroi force-pushed the joleroi:3D_fit branch from ee20503 to 2fc1325 Apr 18, 2018
@joleroi joleroi merged commit 3697088 into gammapy:master Apr 18, 2018
1 of 2 checks passed
1 of 2 checks passed
continuous-integration/travis-ci/pr The Travis CI build is in progress
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
@joleroi joleroi deleted the joleroi:3D_fit branch Apr 18, 2018
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Linked issues

Successfully merging this pull request may close these issues.

None yet

2 participants