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

Change value clipping in LogScale class #2405

Merged
merged 8 commits into from Sep 28, 2019
Merged

Conversation

@QRemy
Copy link
Contributor

QRemy commented Sep 26, 2019

Change interpolation clipping to float64 precision
The clipping to float32 precision caused a bug when evaluating a table model containing the gamma-ray emissivity spectrum of the gas (that is <<10**-38/MeV/s/sr above 500GeV)

@QRemy QRemy requested a review from adonath Sep 26, 2019
@adonath adonath self-assigned this Sep 26, 2019
@adonath adonath modified the milestones: 0.15, 0.14 Sep 26, 2019
@cdeil

This comment has been minimized.

Copy link
Member

cdeil commented Sep 26, 2019

@QRemy - My thinking was that we should go to float32 for better performance (faster and less memory), for maps and IRFs and basically anything except the logL np.sum accumulator scalar (see #2374).

Can you please try to extract a minimal failing example (see here) from your current issue, and add that as a test case? I know this takes time, but without it there's a very good chance that in the coming weeks, months, years, someone would change something in Gammapy and this would break again. Naively I don't understand why float64 should be needed.

@cdeil cdeil added the bug label Sep 26, 2019
@QRemy

This comment has been minimized.

Copy link
Contributor Author

QRemy commented Sep 26, 2019

Here is an example:

energy=np.array([1.00000000e+05, 1.25892541e+05, 1.58489319e+05, 1.99526231e+05,
       2.51188643e+05, 3.16227766e+05, 3.98107171e+05, 5.01187234e+05,
       6.30957344e+05, 7.94328235e+05, 1.00000000e+06,1.25892541e+06,
       1.58489319e+06, 1.99526231e+06])

values=np.array([7.26683943e-35, 3.95773531e-35, 2.15549951e-35, 1.17394868e-35,
       5.44841625e-36, 2.43964992e-36, 1.09240767e-36, 4.89149905e-37,
       2.19027782e-37, 9.80745754e-38, 4.39150790e-38,1.96639562e-38,
       8.80497507e-39, 3.94262401e-39])

model=TableModel(energy=energy,values=values*u.Unit("MeV-1 s-1 sr-1"))
result=model.evaluate(energy,norm=1.,tilt=0.,reference=1*u.TeV)
print(result)
print(np.finfo(np.float32).tiny)

[7.26683943e-35 3.95773531e-35 2.15549951e-35 1.17394868e-35
 5.44841625e-36 2.43964992e-36 1.09240767e-36 4.89149905e-37
 2.19027782e-37 9.80745754e-38 4.39150790e-38 1.96639562e-38
 1.17549435e-38 1.17549435e-38] 1 / (MeV s sr)
1.1754944e-38

The two last values are below float32 precision so they are clipped to a higher value than expected.
Changing this to float64 here does not prevent to have the flux or npred of the final model (spatial+spectral) returned in float32 if you want.

I added this test:

np.testing.assert_allclose(values/values.max(),result.value/values.max())

Traceback (most recent call last):

  File "<ipython-input-532-185cc20d386d>", line 15, in <module>
    np.testing.assert_allclose(values/values.max(),result.value/values.max())

  File "/anaconda3/envs/gpy-dev/lib/python3.7/site-packages/numpy/testing/_private/utils.py", line 1501, in assert_allclose
    verbose=verbose, header=header, equal_nan=equal_nan)

  File "/anaconda3/envs/gpy-dev/lib/python3.7/site-packages/numpy/testing/_private/utils.py", line 827, in assert_array_compare
    raise AssertionError(msg)

AssertionError: 
Not equal to tolerance rtol=1e-07, atol=0

Mismatch: 14.3%
Max absolute difference: 0.00010751
Max relative difference: 0.66459864
 x: array([1.000000e+00, 5.446295e-01, 2.966213e-01, 1.615487e-01,
       7.497642e-02, 3.357237e-02, 1.503278e-02, 6.731261e-03,
       3.014072e-03, 1.349618e-03, 6.043216e-04, 2.705985e-04,
       1.211665e-04, 5.425500e-05])
 y: array([1.000000e+00, 5.446295e-01, 2.966213e-01, 1.615487e-01,
       7.497642e-02, 3.357237e-02, 1.503278e-02, 6.731261e-03,
       3.014072e-03, 1.349618e-03, 6.043216e-04, 2.705985e-04,
       1.617614e-04, 1.617614e-04])
@adonath

This comment has been minimized.

Copy link
Member

adonath commented Sep 26, 2019

Thanks @QRemy! I agree the clipping of the tiny values might be an issue, but I think the fix needs some more thought, because it re-introduces the problem for np.float32 values, that we tried to solve with the clipping:

from gammapy.utils.interpolation import LogScale
import numpy as np

LogScale.tiny = np.finfo(np.float64).tiny

values = np.array([1e-12, 0], dtype=np.float32)

scale = LogScale()
scale(values)

Gives:

array([-27.631021,       -inf], dtype=float32)

Where the -inf value breaks the linear interpolation, that is then done on the scaled values (I guess these are the test fails we see on Travis...). I guess there are two solutions:

  • If we really need high precision fluxes at high energy, we could change the default units for the Fermi isotropic diffuse model.
  • Somehow check the dtype of values (e.g. using np.dtype()) and set the clip values based on the dtype. I can't remember if there was a good motivation to set it to a fixed value, probably not. But then the issue well appear again if we maybe address #2374.
@adonath adonath modified the milestones: 0.14, 0.15 Sep 26, 2019
@QRemy QRemy force-pushed the QRemy:fix_interp_precision branch from ffd1034 to ed59b1e Sep 27, 2019
Copy link
Member

adonath left a comment

Thanks @QRemy, I have one specific comment to move the clipping directly to the LogScale class. Otherwise the looks good to me.

@@ -92,9 +94,21 @@ def __call__(self, points, method="linear", clip=True, **kwargs):
values = self._interpolate(points[0])
values = self.scale.inverse(values)

tiny = np.finfo(np.float32).tiny

This comment has been minimized.

Copy link
@adonath

adonath Sep 27, 2019

Member

As this is a specific problem of the LogScale, I think we should move the code for the clipping there. So my proposal would be to define self.is_tiny < (values < self.tiny) in LogScale._scale and the re-apply the mask in LogScale._inverse with values[self._is_tiny] = 0. This will also get rid of the try and except block for quantities, because the interpolation scale is called with the units stripped-of.

This comment has been minimized.

Copy link
@QRemy

QRemy Sep 27, 2019

Author Contributor

I moved the clipping code to LogScale.inverse, and I defined the mask only there as it need to be applied only to the output values. I added an other try and except block for the case where value is not an array.

@QRemy QRemy changed the title change interpolation clipping to float64 precision change interpolation clipping Sep 27, 2019
QRemy added 2 commits Sep 27, 2019
@adonath adonath changed the title change interpolation clipping Change value clipping in `LogScale` class Sep 27, 2019
Copy link
Member

adonath left a comment

Thanks @QRemy!

@adonath adonath merged commit ab57989 into gammapy:master Sep 28, 2019
9 checks passed
9 checks passed
Codacy/PR Quality Review Up to standards. A positive pull request.
Details
Scrutinizer Analysis: 1 updated code elements – Tests: passed
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
gammapy.gammapy Build #20190927.23 succeeded
Details
gammapy.gammapy (DevDocs) DevDocs succeeded
Details
gammapy.gammapy (Lint) Lint succeeded
Details
gammapy.gammapy (Test Python36) Test Python36 succeeded
Details
gammapy.gammapy (Test Windows36) Test Windows36 succeeded
Details
gammapy.gammapy (Test Windows37) Test Windows37 succeeded
Details
@adonath adonath changed the title Change value clipping in `LogScale` class Change value clipping in LogScale class Nov 18, 2019
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
3 participants
You can’t perform that action at this time.