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 spectral model absorbed by EBL that can be fit #988

Merged
merged 9 commits into from Apr 27, 2017
Merged

Conversation

@jjlk
Copy link
Contributor

@jjlk jjlk commented Apr 25, 2017

Hi @cdeil ,
As we discussed here is the PR to add a spectral model absorbed by the EBL that can be fit with sherpa. The fit is running well.

The problem arises when I want to propagate errors (for a butterfly or for flux points). I use a NDDataArray in the SpectralModelAbsorbedbyEbl.evaluate function and it can't handle the error propagation. Do you think that we can try to fix this?

I added a dedicated test to reproduce the bug. The implementation of the model could be improve and we could discuss that further when the error propagation problem will be solved.

Thanks!

@cdeil cdeil self-assigned this Apr 25, 2017
@cdeil cdeil added the feature label Apr 25, 2017
@cdeil cdeil added this to the 0.6 milestone Apr 25, 2017
Copy link
Member

@cdeil cdeil left a comment

@jjlk - I had started a code review before we skyped. Here's the inline comments.

As I said in the call, my main suggestion here would be to move the EBLAbsorption code to a separate class, and then to have the EBL absorbed spectral model take an instance of that class.

And leave the part of getting spectral_model.evaluate_error to work to a future PR. (for now, you can leave a TODO comment and a test with a @pytest.mark.xfail marker.

"""
def __init__(self, spectral_model, redshift,
ebl_model='$GAMMAPY_EXTRA/datasets/ebl/ebl_franceschini.fits.gz'):
"""

This comment has been minimized.

@cdeil

cdeil Apr 25, 2017
Member

We always document parameters in the class-level docstring and have no docstrings for __init__.
So please move the content.

from gammapy.spectrum.models import SpectralModelAbsorbedbyEbl
from gammapy.scripts import CTAPerf
from gammapy.scripts.cta_utils import CTAObservationSimulation, Target, ObservationParameters
from gammapy.spectrum import SpectrumObservationList, SpectrumFit, SpectrumResult

This comment has been minimized.

@cdeil

cdeil Apr 25, 2017
Member

Move the imports to the top of the file.

redshift=0.1,
ebl_model='$GAMMAPY_EXTRA/datasets/ebl/ebl_dominguez11.fits.gz')

# should be changed, I hack to avoid double absorption by EBL

This comment has been minimized.

@cdeil

cdeil Apr 25, 2017
Member

If you have such comments as "should be changed" here, please always start with the string "TODO".
This allows us to find all such places in the future and (hopefully some day) clean everything up.

energy_bounds = EnergyBounds.from_lower_and_upper_bounds(lower=energy_lo,
upper=energy_hi,
unit=u.keV)
energy = energy_bounds.log_centers

This comment has been minimized.

@cdeil

cdeil Apr 25, 2017
Member

Local variable energy is not used. Remove?

@cdeil cdeil modified the milestones: 0.7, 0.6 Apr 25, 2017
jlk added 4 commits Apr 26, 2017
@jjlk
Copy link
Contributor Author

@jjlk jjlk commented Apr 26, 2017

Hi @cdeil ,
As we discussed, I introduced a new class to deal with absorption (not only EBL so I removed 'EBL' everywhere), Absorption. I adapted all the ~CTA-dedicated classes and all the tests. The AbsorbedSpectralModel class have been adapted also. As you suggested, one of the test has been marked as failed (error propagation for NDDataArray).

If you are ok, I'll adapt the CTA notebook once the PR is merged. ++

Copy link
Member

@cdeil cdeil left a comment

Left some inline comments.
Main goal for now is to make all tests and the docs build pass.

--------
Create and plot EBL absorption model for a redshift of 0.2
.. create::

This comment has been minimized.

@cdeil

cdeil Apr 26, 2017
Member

Does this exist / work? I've never seen the .. create:: Sphinx directive. Usually we use the .. plot:: Sphinx directive.

interpolation_mode='log', name='energy')
]

self.absorption = NDDataArray(axes=axes, data=data)

This comment has been minimized.

@cdeil

cdeil Apr 26, 2017
Member

Suggest to rename self.absorption to self.data here.
Otherwise you get absorption.absorption in code that uses this class, which is confusing.


# Target, PKS 2155-304 from 3FHL
name = 'test'
absorption = Absorption.read('$GAMMAPY_EXTRA/datasets/ebl/ebl_dominguez11.fits.gz')

This comment has been minimized.

@cdeil

cdeil Apr 26, 2017
Member

Can you make it convenient to access the "built in" EBL absorption models we provide, via

class Absorption:
    @classmethod
    def read_builtin(name):
        """"Read one of the built-in absorption models.
        Parameters
        -----------
        name : {'dominquez11', ...}

This comment has been minimized.

@jjlk

jjlk Apr 27, 2017
Author Contributor

Can we keep that? Since people my want to use their own model.

This comment has been minimized.

@cdeil

cdeil Apr 27, 2017
Member

Yes, we should keep the read method.

I was suggesting to add another method that calls read, where the user can just select the builtin absorption model via a string, instead of having to know where the file is.

This would also serve as documentation which absorption models we have in Gammapy.
At the moment even I don't know which files are available and are supposed to work with this class.
So if you add the extra read_builtin method and then document there the list of models we have, and add a test that shows that reading and evaluating works for those, we are in better shape, no?

This comment has been minimized.

@jjlk

jjlk Apr 27, 2017
Author Contributor

Oké done.

obs_param=obs_param)

# Model we want to fit
model = AbsorbedSpectralModel(spectral_model=PowerLaw(index=2.5 * u.Unit(''),

This comment has been minimized.

@cdeil

cdeil Apr 26, 2017
Member

Suggest to make the PowerLaw object on a line before and then pass it in here (less indented code)

fit.est_errors()
result = fit.result[0]

# here it explodes since NDDataArray can't handle error propagation

This comment has been minimized.

@cdeil

cdeil Apr 26, 2017
Member

You commented that you moved this failing part to an xfailed test?
I don't see that in the diff on Github!?

This comment has been minimized.

@jjlk

jjlk Apr 27, 2017
Author Contributor

Yess, I added the line @pytest.mark.xfail before the test_spectral_model_absorbed_by_ebl() test. Is it the way it's working?

This comment has been minimized.

@cdeil

cdeil Apr 27, 2017
Member

Just xfailing the whole test is bad. I have no idea if the code is working. I'll break as we refactor.

Please change the test so that the fit executes, either by commenting out this part at the end with a TODO marker, or by refactoring the test into several tests (one for fit, one fore ...)

This comment has been minimized.

@jjlk

jjlk Apr 27, 2017
Author Contributor

I added a TODO and I removed the @pytest.mark.xfailmarker

return ObservationParameters(alpha=alpha, livetime=livetime,
emin=emin, emax=emax)

@requires_data('gammapy-extra')

This comment has been minimized.

@cdeil

cdeil Apr 26, 2017
Member

@requires_data doesn't work on fixtures. Only use it on tests.

jlk
@jjlk
Copy link
Contributor Author

@jjlk jjlk commented Apr 27, 2017

HI @cdeil,
I did the the changes. Are we good?

@cdeil
cdeil approved these changes Apr 27, 2017
@cdeil cdeil modified the milestones: 0.6, 0.7 Apr 27, 2017
@cdeil
Copy link
Member

@cdeil cdeil commented Apr 27, 2017

I've re-started the CI builds

table = Absorption.read(filename='$GAMMAPY_EXTRA/datasets/ebl/ebl_dominguez11.fits.gz').table_model(parameter=0.2)
ax = plt.gca()
opts = dict(energy_range=[0.01, 100] * u.TeV, energy_unit='TeV', ax=ax)

This comment has been minimized.

@cdeil

cdeil Apr 27, 2017
Member

Does this example work?
You're missing an import import astropy.units as u, no?

What I do to test examples locally is to copying the code snipped and then typing %paste and enter in IPython.

This comment has been minimized.

@jjlk

jjlk Apr 27, 2017
Author Contributor

Fixed

@requires_data('gammapy-extra')
def test_cta_simulation():
text = str(cta_simu())
assert '*** Observation summary report ***' in text

This comment has been minimized.

@cdeil

cdeil Apr 27, 2017
Member

Shouldn't this test assert on something that was computed?
(the same holds for the other tests above that only assert on str, but especially this one because it actually does some coputations, no?)

This comment has been minimized.

@jjlk

jjlk Apr 27, 2017
Author Contributor

Yes it should. But there are some randomisation and I didn't add the possibility to fix the seed. My goal was to merge these tools with Johannes class, SpectrumSimulation. Maybe for the 0.7 release? I added a small test to check if we get more than 5 sigma for detection

jlk
@jjlk
Copy link
Contributor Author

@jjlk jjlk commented Apr 27, 2017

Hi @cdeil ,
Is it okay now?

Parameters
----------
param : `float`

This comment has been minimized.

@cdeil

cdeil Apr 27, 2017
Member

"param" -> "parameter"

def table_model(self, parameter, unit='TeV'):
"""
Returns `~gammapy.spectrum.models.TableModel` for a given parameter
from the input aborbed model

This comment has been minimized.

@cdeil

cdeil Apr 27, 2017
Member

aborbed -> absorbed.

Generally, try to fit the description on the first line. End it with a dot. Then if needed a blank line and more text after that. Fill words like "Returns" can always be omitted. The return type can be listed in a Returns Section after Parameters.

jlk
@jjlk
Copy link
Contributor Author

@jjlk jjlk commented Apr 27, 2017

HI @cdeil,
Done and done. You'll love the plot ++

@cdeil
Copy link
Member

@cdeil cdeil commented Apr 27, 2017

@jjlk - Merci pour les absorbcions de l'EBL! O la la.

@cdeil cdeil merged commit 70f7efd into gammapy:master Apr 27, 2017
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Linked issues

Successfully merging this pull request may close these issues.

None yet

2 participants