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

Implement SpectralGaussian model class #2244

Merged
merged 6 commits into from Jun 21, 2019

Conversation

@JouvinLea
Copy link

@JouvinLea JouvinLea commented Jun 17, 2019

Implement Gaussian model in the model Class

@JouvinLea JouvinLea requested a review from registerrier Jun 17, 2019
val_at_2TeV=u.Quantity(0.3321457437841875, 'cm-2 s-1 TeV-1'),
integral_1_10TeV=u.Quantity(1.6371892282407274, 'cm-2 s-1'),
integral_infinity=u.Quantity(2, 'cm-2 s-1'),
eflux_1_10TeV=u.Quantity(7.676635460998642, 'TeV cm-2 s-1'),
Copy link
Author

@JouvinLea JouvinLea Jun 17, 2019

Maybe a better test for the energy flux?

energy_power=0,
n_points=100,
**kwargs
self,
Copy link
Member

@cdeil cdeil Jun 17, 2019

There's a lot of unrelated lines that show up in this diff. Could you please run black gammapy/spectrum/models.py gammapy/spectrum/tests/test_models.py to go back to the black coding style?

Copy link
Member

@cdeil cdeil left a comment

I put a few superficial suggestions.

I didn't check the equations and implementation yet.
Is there a test case we could choose with a "known good answer"?
I can't think of one, but maybe there is, e.g. on Wikipedia or somewhere?

r"""Integrate Gaussian analytically.
.. math::
F(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}} \phi(E)dE = \frac{N_0}{2} \left( erf(\frac{E_{max}-
Copy link
Member

@cdeil cdeil Jun 17, 2019

Please break the line earlier and make it a bit easier to read.

u_max = (emax - pars['mean'].quantity) / (np.sqrt(2) * pars['sigma'].quantity)
a = pars['norm'].quantity * pars['sigma'].quantity / np.sqrt(2 * np.pi)
b = pars['norm'].quantity * pars['mean'].quantity / 2
int = a * (np.exp(- u_min ** 2) - np.exp(-u_max ** 2)) + b * (erf(u_max) - erf(u_min))
Copy link
Member

@cdeil cdeil Jun 17, 2019

int is a keyword, shouldn't be used a local variable name. Just change int = to return here.

.. math::
G(E_{min}, E_{max}) = \int_{E_{min}}^{E_{max}}E \phi(E)dE =
\frac{N_0 \sigma}{\sqrt{2*\pi}}* \left( \exp(\frac{E_{min}-\bar{E}}{\sqrt{2} \sigma}) -
exp(\frac{E_{max} - \bar{E}}{\sqrt{2} \sigma}) \right) \\ + \frac{N_0 * \bar{E}}{2} \left( erf(\frac{E_{max}-\bar{E}}{\sqrt{2} \sigma})
Copy link
Member

@cdeil cdeil Jun 17, 2019

This is very hard to read. Do we need / want the equation, or could we link to Wikipedia or somewhere else instead?
If you keep it, suggest to remove \int_{E_{min}}^{E_{max}}E \phi(E)dE = and to break the other lines more and shorten them to make the LaTeX more readable (for review here, but also if someone looks at the help in ipython or Jupyter they will see it non-rendered)

Copy link
Author

@JouvinLea JouvinLea Jun 18, 2019

This is very hard to read. Do we need / want the equation, or could we link to Wikipedia or somewhere else instead?
If you keep it, suggest to remove \int_{E_{min}}^{E_{max}}E \phi(E)dE = and to break the other lines more and shorten them to make the LaTeX more readable (for review here, but also if someone looks at the help in ipython or Jupyter they will see it non-rendered)

ok, I clean a bit. I don't know on the doc I found it pretty clear to read no? I did the analytical calculation myself but I guess we could find it on wikipedia but I don't find it ugly as it is.

@cdeil cdeil added the feature label Jun 17, 2019
@cdeil cdeil added this to the 0.13 milestone Jun 17, 2019
@cdeil cdeil added this to To do in gammapy.makers via automation Jun 17, 2019
@cdeil
Copy link
Member

@cdeil cdeil commented Jun 17, 2019

@JouvinLea - Thanks!

There's a conflict in gammapy/spectrum/tests/test_models.py - please resolve it by rebasing against upstream master so that this will run in CI. I also suggest to squash the commits into one or two with a good commit message. And to merge this in as soon as it's clear the equations and implementation is good.

@registerrier - Could you please review this PR and merge it in when ready?

x-ref: #1411

(unsubscribing now)

@cdeil cdeil changed the title Gaussian spectral model Add Gaussian spectral model Jun 17, 2019
@JouvinLea JouvinLea force-pushed the gaussian_spectral_model branch from c6bf7e2 to 774159f Jun 18, 2019
@JouvinLea JouvinLea force-pushed the gaussian_spectral_model branch from 5d144a5 to 36334b2 Jun 18, 2019
@JouvinLea
Copy link
Author

@JouvinLea JouvinLea commented Jun 18, 2019

I put a few superficial suggestions.

I didn't check the equations and implementation yet.
Is there a test case we could choose with a "known good answer"?
I can't think of one, but maybe there is, e.g. on Wikipedia or somewhere?

@JouvinLea JouvinLea closed this Jun 18, 2019
gammapy.makers automation moved this from To do to Done Jun 18, 2019
@JouvinLea
Copy link
Author

@JouvinLea JouvinLea commented Jun 18, 2019

Hi,
Sorry for the merge conflict, I started on the morning I didn't think I could have create merge conflict in one day. Now it is fixed and I squashed also the commits!

@JouvinLea
Copy link
Author

@JouvinLea JouvinLea commented Jun 18, 2019

I put a few superficial suggestions.

I didn't check the equations and implementation yet.
Is there a test case we could choose with a "known good answer"?
I can't think of one, but maybe there is, e.g. on Wikipedia or somewhere?

I don't know, I recalculate evrything my self. I think for the integral it is ok since the result is quite obvious and in the test I add a check that the integral between -inf and inf is q=equal to the norm No.

For the energy flux, it is true that the test is not a super good test since I didn't what to test to know If analytical integration was correct.

@adonath adonath reopened this Jun 18, 2019
gammapy.makers automation moved this from Done to In progress Jun 18, 2019
Copy link
Contributor

@registerrier registerrier left a comment

Thanks @JouvinLea . Looks good to me.
Some comments inline. Also run black formatter as @cdeil suggested.

gammapy/spectrum/tests/test_models.py Outdated Show resolved Hide resolved
gammapy/spectrum/tests/test_models.py Outdated Show resolved Hide resolved
gammapy/spectrum/models.py Show resolved Hide resolved
gammapy/spectrum/models.py Outdated Show resolved Hide resolved
gammapy/spectrum/models.py Outdated Show resolved Hide resolved
@JouvinLea
Copy link
Author

@JouvinLea JouvinLea commented Jun 21, 2019

@registerrier
The test are cheked with http://m.wolframalpha.com equation are ok!

from scipy.special import erf

pars = self.parameters
u_min = (emin - pars["mean"].quantity) / (np.sqrt(2) * pars["sigma"].quantity)
Copy link
Contributor

@registerrier registerrier Jun 21, 2019

you should add a .to_value('') for u_min and u_max.
Apparently astropy does not always handle the erf correctly for units.

Copy link
Contributor

@registerrier registerrier left a comment

Thanks @JouvinLea !
The fail is CI is in fact due to the handling of quantities with erf. So if you can apply a .to_value('') for your u_min and u_max argument, this should work fine.
Then we can merge once CI passes.

Changed to have `u_min` and `u_max` as floats
With correct parenthesis...
@registerrier registerrier merged commit e388ccf into gammapy:master Jun 21, 2019
9 checks passed
gammapy.makers automation moved this from In progress to Done Jun 21, 2019
@registerrier
Copy link
Contributor

@registerrier registerrier commented Jun 21, 2019

Thanks @JouvinLea .

@adonath adonath changed the title Add Gaussian spectral model Implement SpectralGaussian model class Jul 25, 2019
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
Linked issues

Successfully merging this pull request may close these issues.

None yet

4 participants