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

Update interpolation code to pyirf v0.10 API #1184

Merged
merged 5 commits into from
Nov 30, 2023
Merged

Update interpolation code to pyirf v0.10 API #1184

merged 5 commits into from
Nov 30, 2023

Conversation

morcuended
Copy link
Member

@morcuended morcuended commented Nov 27, 2023

See #1160

Changed the interpolation code used in lstchain to support the new API of pyirf as explained in #1160. These changes keep the same functionality as before. They do not allow for new extrapolation functionality (that can be done elsewhere).

Copy link

codecov bot commented Nov 27, 2023

Codecov Report

Attention: 53 lines in your changes are missing coverage. Please review.

Comparison is base (045a5fd) 74.42% compared to head (674501f) 74.39%.
Report is 29 commits behind head on main.

Files Patch % Lines
lstchain/scripts/lstchain_create_run_summary.py 60.29% 27 Missing ⚠️
lstchain/scripts/lstchain_merge_run_summaries.py 75.75% 24 Missing ⚠️
lstchain/high_level/interpolate.py 83.33% 2 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##             main    #1184      +/-   ##
==========================================
- Coverage   74.42%   74.39%   -0.04%     
==========================================
  Files         129      129              
  Lines       13169    13199      +30     
==========================================
+ Hits         9801     9819      +18     
- Misses       3368     3380      +12     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

Comment on lines +365 to +368
interp = GridDataInterpolator(
grid_points=grid_points, params=gh_cuts, method=method
)
gh_cuts_interp = interp(target_point)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If extrapolation is needed I'd suggest to create a Estimator object (same for the AL_CUTS), just like it is done for RAD_MAX tables in pyirf. If you don't want extrapolation this should be fine.

@morcuended
Copy link
Member Author

@RuneDominik, ready for review. Please have a look at the changes.

Although you suggested leaving all the optional arguments by default, I have included the interpolation method in the effective area and RADMAX estimators (via interpolator_kwargs) to keep this possibility as it was before.

@morcuended morcuended added the ready for review Pull request is ready for review label Nov 27, 2023
@RuneDominik
Copy link
Member

RuneDominik commented Nov 27, 2023

LGTM, if you want this to work just as the old version. One could argue that this would be a good opportunity to support the whole API, including extrapolation and support for different inter-/extrapolator classes.

@morcuended
Copy link
Member Author

LGTM, if you want this to work just as the old version. One could argue that this would be a good opportunity to support the whole API, including extrapolation and support for different inter-/extrapolator classes.

Many thanks for the quick review. I agree with you. However, this change has been pending for quite some time and I'd say it's somewhat urgent to make a new minor release (with the interpolation code at least working as before). I suggest merging this PR as it is. Make the release. And extend the support for extrapolation in another PR. What do you think, @rlopezcoto?

@rlopezcoto
Copy link
Contributor

thanks @RuneDominik for the review and the comment. I support @morcuended's point of view, since this change is urgently needed to release a bugfix, I would just stick to solve this issue.

@morcuended
Copy link
Member Author

I tested these changes (using pyirf v0.10) on a Crab sample and compared them with results from lstchain v0.10.4 (pyirf v0.8 + applying the edisp fix script). Both use interpolation. The resulting SED is not the same.

pyirf v0.8 + applying the edisp fix script

flux_points_LP_8bin_gh_eff_0 7_th_cont_0 7_pyirf08

pyirf v0.10 (this PR)

flux_points_LP_8bin_gh_eff_0 7_th_cont_0 7_pyirf010

@RuneDominik was there any change in the interpolation code of pyirf (from v0.8 to v0.10) that can explain this change in the spectrum?

@RuneDominik
Copy link
Member

The was the edisp fix which was also applied to the interpolation code. As far as I see when going though the QuantileInterpolator history, this was the only one actually changing the code. I assume the IRFs had the applicable edisp corrections, so that

  • pyirf v0.8 + applying the edisp fix script used uncorrected IRFs (thus with edisp summed to 1) as input and the output was corrected to an integral of 1
  • pyirf v0.10 (this PR) used correct IRFs as input and no re-normalization was applied to the output?

@morcuended
Copy link
Member Author

pyirf v0.8 + applying the edisp fix script used uncorrected IRFs (thus with edisp summed to 1) as input and the output was corrected to an integral of 1

Correct. I applied the script correcting the edisp normalization to the DL3 files produced with pyirf v0.8

pyirf v0.10 (this PR) used correct IRFs as input and no re-normalization was applied to the output?

No re-normalization was applied to the output in this case.

@RuneDominik
Copy link
Member

RuneDominik commented Nov 29, 2023

I tested a bit and found that for data with a logarithmic binning (hence edisps) the results have changed when switching between pyirf 0.9 (pre normalization fix) and pyirf 0.10 (post normalization fix), as you can see here in the residuals between truth and estimated (and re-normalized if needed) Gaussian histograms.

Residuals

For linear binning I found no difference. Thinking of it this seems plausible, as the bin-width quite drastically changes for logarithmic binning and is applied for normalization in the pyirf 0.10 version. As of now I would not call each of them better or worse, just different. Although more thorough tests are ofc. needed.

@morcuended
Copy link
Member Author

Thanks a lot for the check @RuneDominik. In view of this, I think we can proceed with the merging.

@rlopezcoto
Copy link
Contributor

thanks a lot for the fix and the checks @morcuended and @RuneDominik! if you approve this, I think we can move on with this PR

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
ready for review Pull request is ready for review
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants