Add sigma v estimator functionality for dark matter use case #2075
Conversation
There was a problem hiding this comment.
Thanks @Bultako! This is a very nice contribution and feature for the dark-matter community. I've left a few in-line comments.
Another comment I have is, whether you could split the .run() method, because currently it's very long. At least introduce a method that compute the likelihood profiles and stores in a Table and a second method that starts from the table and compute the best fit values and upper limits form the likelihood profile. This way the the code becomes easier to debug and reuse e.g. when the brentq fails. The .run() method could stay and just call both of the steps one after another. Again you might want to take a look how the FluxPointEstimator is implemented.
gammapy/astro/darkmatter/utils.py
Outdated
| jfact = 3.41e19 * u.Unit("GeV2 cm-5") | ||
| channels = ["b", "t", "Z"] | ||
| masses = [70, 200, 500, 5000, 10000, 50000, 100000]*u.GeV | ||
| estimator = SigmaVEstimator(simulated_dataset, masses, channels, jfact=JFAC) |
There was a problem hiding this comment.
I think this example does not work, as the simulated_dataset is missing...I guess now that #2064 is merged it can be added?
gammapy/astro/darkmatter/utils.py
Outdated
| xsection = DMAnnihilation.THERMAL_RELIC_CROSS_SECTION | ||
| self.xsection = xsection | ||
|
|
||
| self._spatial_model = dataset.model.spatial_model |
There was a problem hiding this comment.
Is it really useful to introduce these "private" variables? Why not create and use local variables later?
gammapy/astro/darkmatter/utils.py
Outdated
| self._psf = dataset.psf | ||
| self._edisp = dataset.edisp | ||
|
|
||
| self._counts_map = WcsNDMap(self._geom, np.random.poisson(dataset.npred().data)) |
There was a problem hiding this comment.
I think the simulated dataset already contains a simulated counts map, so maybe remove?
There was a problem hiding this comment.
I had imagined calling this SigmaVEstimator.run() inside an external loop of n runs (see cell[13] in notebook https://github.com/peroju/dm-gammapy/blob/master/notebooks/DarkMatterUseCaseUpgraded.ipynb and the TODO comment in the code cell) Every call in the loop should create a different random poisson counts map from our simulated dataset, and the results of each individual fitting of the run would be averaged at the end.
I had thought of providing a different tool later (if needed) building this loop of n runs, using SigmaVEstimator.run() internally in the loop and throwing the averaged results.
Does this make sense?
gammapy/astro/darkmatter/utils.py
Outdated
| Dict with results as `~astropy.table.Table` objects for each channel. | ||
| """ | ||
|
|
||
| result = {} |
There was a problem hiding this comment.
Would it make sense to use a flattened data structure here? I.e. use an astropy.table.Table with a mass and channel column? For every combination of channel and mass one row is added to the table. I think this simplifies the code as well as the data structure that is returned.
There was a problem hiding this comment.
There is a certain hierarchy to be kept in the results (at least from the user-access point of view) where every channel has a list of (mass, sigma-v, err, etc.) values. This is the main reason why tables are linked to dictionary items (channels).
gammapy/astro/darkmatter/utils.py
Outdated
| ) | ||
| try: | ||
| fit = Fit(dataset_loop) | ||
| fit.datasets.parameters.apply_autoscale = False |
There was a problem hiding this comment.
As you use the fit.likelihood_profile() method now, it's not necessary to switch off the autoscaling. Maybe just remove...
gammapy/astro/darkmatter/utils.py
Outdated
| fit_result = fit.run(optimize_opts, covariance_opts) | ||
| likemin = fit_result.total_stat | ||
| profile = fit.likelihood_profile( | ||
| parameter="scale", bounds=5, nvalues=50 |
There was a problem hiding this comment.
I would definitely make the bounds and nvalues configurable. Maybe take a look at what we did for the FluxPointEstimator.
gammapy/astro/darkmatter/utils.py
Outdated
| parameter="scale", bounds=5, nvalues=50 | ||
| ) | ||
| xvals = profile["values"] | ||
| yvals = profile["likelihood"] - likemin - 2.71 |
There was a problem hiding this comment.
Again I would introduce a ts_diff variable to make it configurable (with using 2.71 as a default value). If you don't want to make it configurable at all please add as global variable TS_DIFF and don't hard-code the number in the method.
gammapy/astro/darkmatter/utils.py
Outdated
| scale_max = max(xvals) | ||
|
|
||
| scale_found = brentq( | ||
| interp1d(xvals, yvals, kind="cubic"), |
There was a problem hiding this comment.
I think cubic interpolation is probably not suitable for likelihood profile. I presume quadratic should be more stable (not sure whether it makes a difference though...). There is also a helper function _interp_likelihood_profile, that you could use.
Bultako
force-pushed
the
sigma_estimator
branch
2 times, most recently
from
208e452
to
a7fae18
Compare
|
Thanks @adonath for the review. |
Bultako
force-pushed
the
sigma_estimator
branch
3 times, most recently
from
19e0ec1
to
b68f711
Compare
Bultako
force-pushed
the
sigma_estimator
branch
from
b68f711
to
4e654bc
Compare
Bultako
force-pushed
the
sigma_estimator
branch
from
3bc2509
to
ae4cb6a
Compare
Bultako
force-pushed
the
sigma_estimator
branch
6 times, most recently
from
683b972
to
75d63f9
Compare
Bultako
force-pushed
the
sigma_estimator
branch
from
12a927f
to
1e2e0c9
Compare
Bultako
force-pushed
the
sigma_estimator
branch
from
8dfb336
to
a24c869
Compare
This PR adds a new class
SigmaVEstimatorinastro.darkmatter.utilsas a first attempt to address a general use case in the dark matter community. See a preliminary notebook usingSigmaVEstimatorhere.It still needs to be improved, mainly with a clear method on how to estimate upper limits when fitting extremely low signals compared to the background model. Once this point will be solved, tests could be added accordingly.
For the moment I leave it as an open/work in progress PR, in case it needs discussion at this stage.