cta-observatory / ctapipe

Low-level data processing pipeline software for CTAO or similar arrays of Imaging Atmospheric Cherenkov Telescopes
https://ctapipe.readthedocs.org
BSD 3-Clause "New" or "Revised" License
63 stars 267 forks source link

Overburden-to-height conversion functions #2422

Closed gschwefer closed 5 months ago

gschwefer commented 10 months ago

Implements the overburden to height conversion functions in the atmosphere module

maxnoe commented 10 months ago

I wanted to make a simple plot of height vs. grammage using this function to check for the edge behavior for the cubic interpolation, but it flat out failed:

from ctapipe.io import EventSource
import numpy as np
import astropy.units as u

source = EventSource("dataset://gamma_prod5.simtel.zst")
grammage = np.linspace(0, 1200, 100) * u.g/u.cm**2
source.atmosphere_density_profile.height_from_overburden(grammage)
/home/maxnoe/Uni/CTA/ctapipe/ctapipe/atmosphere.py:307: RuntimeWarning: divide by zero encountered in log10
  self._height_interp(np.log10(overburden.to("g cm-2").value)), u.km
Traceback (most recent call last):
  File "/home/maxnoe/Uni/CTA/ctapipe/invalid_values.py", line 7, in <module>
    source.atmosphere_density_profile.height_from_overburden(grammage)
  File "/home/maxnoe/.local/conda/envs/cta-dev/lib/python3.11/site-packages/astropy/units/decorators.py", line 313, in wrapper
    return_ = wrapped_function(*func_args, **func_kwargs)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/maxnoe/Uni/CTA/ctapipe/ctapipe/atmosphere.py", line 307, in height_from_overburden
    self._height_interp(np.log10(overburden.to("g cm-2").value)), u.km
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/maxnoe/.local/conda/envs/cta-dev/lib/python3.11/site-packages/scipy/interpolate/_polyint.py", line 80, in __call__
    y = self._evaluate(x)
        ^^^^^^^^^^^^^^^^^
  File "/home/maxnoe/.local/conda/envs/cta-dev/lib/python3.11/site-packages/scipy/interpolate/_interpolate.py", line 755, in _evaluate
    below_bounds, above_bounds = self._check_bounds(x_new)
                                 ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/maxnoe/.local/conda/envs/cta-dev/lib/python3.11/site-packages/scipy/interpolate/_interpolate.py", line 784, in _check_bounds
    raise ValueError("A value ({}) in x_new is below "
ValueError: A value (-inf) in x_new is below the interpolation range's minimum value (-4.864756005337933).
maxnoe commented 10 months ago

The error is mostly misleading, because it reports the transformed by log-values, so the user cannot really tell what is wrong easily

maxnoe commented 10 months ago

Please also rebase vs. current main to get the docs build working

maxnoe commented 10 months ago

@gschwefer

I played around a bit locally and wanted to open a PR against your branch but accidentally pushed here directly.

Could you have a look at 51dbae7 and see if this ok with you? If not I'll revert

sorry

gschwefer commented 10 months ago

The commit looks alright. Should there be a test similar to the one you added for the other profiles as well, i.e. should we explicitly guard against inputing a negative height to the exponential model or is this a potentially relevant case?

maxnoe commented 10 months ago

@gschwefer Yes, I think the three classes should behave the same regarding invalid inputs, otherwise they are not interchangeable

gschwefer commented 10 months ago

ImPACT reconstructs slant depth, but what that is converted to as an output is up to whatever we decide it should. Currently, we would need a function that is height_from_slant_depth.

gschwefer commented 10 months ago

Another question about the functions currently there: Is the line_of_sight_integral function really correct? Because the distance you put into it is a distance from the observer, but the integral functions uses heights a.s.l.

maxnoe commented 10 months ago

It seems you ran an older version of black that formats the powers slightly differently.

Running black via pre-commit should use the same version as in the CI, e.g. copy the exact command of the CI:

$ pre-commit run --show-diff-on-failure --color=always --files $(git diff origin/main --name-only)
gschwefer commented 10 months ago

Yep thanks for the tip, I committed from a different ctapipe installation than normal where the pre-commit hooks are not setup properly

maxnoe commented 6 months ago

Hi @gschwefer Could you fix the conflicts?

gschwefer commented 6 months ago

I'll have a look later today

codecov[bot] commented 6 months ago

Codecov Report

Attention: Patch coverage is 98.55072% with 2 lines in your changes are missing coverage. Please review.

Project coverage is 92.66%. Comparing base (f8d1049) to head (0799687). Report is 7 commits behind head on main.

Files Patch % Lines
src/ctapipe/atmosphere.py 97.22% 2 Missing :warning:
Additional details and impacted files ```diff @@ Coverage Diff @@ ## main #2422 +/- ## ========================================== + Coverage 92.59% 92.66% +0.06% ========================================== Files 232 232 Lines 20034 20220 +186 ========================================== + Hits 18551 18736 +185 - Misses 1483 1484 +1 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

maxnoe commented 6 months ago

docs build failure is unrelated, I am investigating

maxnoe commented 6 months ago

The docs failure is fixed in main, please rebase

maxnoe commented 6 months ago

@kosack Have your change requests been addressed?

kosack commented 6 months ago

@kosack Have your change requests been addressed?

Will review today. I want to check in detail that everything works as expected

maxnoe commented 5 months ago

@kosack I think we are at the point where we have accumulated enough significant changes to warrant a new release, it would be great to get this in, did you finish your review?

maxnoe commented 5 months ago

@gschwefer Looks good besides the couple of minor issues the linter has found:

https://github.com/cta-observatory/ctapipe/actions/runs/8566542976/job/23476496256?pr=2422

gschwefer commented 5 months ago

Did the pre-commit hooks change somehow? Or why is this not caught by them?

maxnoe commented 5 months ago

The CI is running the pre-commit hooks:

pre-commit run --show-diff-on-failure --color=always --files $(git diff origin/main --name-only)
gschwefer commented 5 months ago

Running the same thing, it doesn't complain for me...but I'll just fix it manually

maxnoe commented 5 months ago

The PR runs a merge request build: it runs on the merged result of your branch into main.

It could be that this branch is not up to date with main and the config changed in main.

maxnoe commented 5 months ago

This branch is 23 commits ahead of, 45 commits behind cta-observatory/ctapipe:main.

Yes, that seems to be the case. Update your branch and you should get the same errors as in the CI:

git fetch upstream
git rebase upstream/main
gschwefer commented 5 months ago

@maxnoe Looks like this tests fail is a github issue to me

maxnoe commented 5 months ago

Yes, seems to be something wrong / changed with the mac os base images, will check

maxnoe commented 5 months ago

I am a bit confused: https://github.com/actions/runner-images?tab=readme-ov-file#available-images

Says the macos-latest should be a macos 14 arm based image, but after I restarted the job, it says it's macos 12 on x86...

Maybe the rollout is not yet complete? And the failed run was running on arm / macos 14.

Edit: no, both were running on macos 12

ctao-dpps-sonarqube[bot] commented 5 months ago

Failed

Analysis Details

2 Issues

Coverage and Duplications

Project ID: cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

View in SonarQube

maxnoe commented 5 months ago

Ok, was flaky, the latest run passed

maxnoe commented 5 months ago

Sonarqube is complaining about duplicated code in the tests, a large chunk is exactly equal:

    fit_reference = np.array(

        [

            [0.00 * 100000, -140.508, 1178.05, 994186, 0],

            [9.75 * 100000, -18.4377, 1265.08, 708915, 0],

            [19.00 * 100000, 0.217565, 1349.22, 636143, 0],

            [46.00 * 100000, -0.000201796, 703.745, 721128, 0],

            [106.00 * 100000, 0.000763128, 1, 1.57247e10, 0],

        ]

    )

    profile_5 = atmo.FiveLayerAtmosphereDensityProfile.from_array(fit_reference)

Could be a module constant or a fixture

gschwefer commented 5 months ago

I'll have a look

ctao-dpps-sonarqube[bot] commented 5 months ago

Passed

Analysis Details

2 Issues

Coverage and Duplications

Project ID: cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

View in SonarQube