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

Bug fix in full numeric pixel likelihood calculation #2388

Merged
merged 14 commits into from
Aug 31, 2023

Conversation

gschwefer
Copy link
Contributor

This fixes a bug in the implementation of the full numeric pixel likelihood from de Naurois 2009 , eq 21, in neg_log_likelihood_numeric. The approximate likelihood in neg_log_likelihood_approx is adjusted for consistency, which is important given they are combined in the neg_log_likelihood function. This was not caught in the tests because these also had a bug, which is now fixed too.

I will also implement the fix in the ImPACT branch where the overall pixel likelihood interface changes.

@gschwefer gschwefer added the bug label Aug 24, 2023
@gschwefer gschwefer self-assigned this Aug 24, 2023
@gschwefer
Copy link
Contributor Author

What exactly is isort complaining about here?

@Tobychev
Copy link
Contributor

What exactly is isort complaining about here?

Probably you haven't sorted your imports according to the standard scriptures of isort, the easiest fix is just to run isort locally and commit the changes it produces.

Do note that it's not actually isort itself that complains, it just modifies the files as part of the CI, but the files being modified by isort counts as a failure, so you never get to actually hear what it did. Thus the advice to just run it yourself on your own copy.

@@ -121,8 +120,9 @@ def neg_log_likelihood_numeric(
Width of single p.e. peak (:math:`σ_γ`).
pedestal: ndarray
Width of pedestal (:math:`σ_p`).
confidence: tuple(float, float), 0 < x < 1
Confidence interval of poisson integration.
confidence: loat, 0 < x < 1
Copy link
Contributor

Choose a reason for hiding this comment

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

probably mean float

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Fixed it

@maxnoe
Copy link
Member

maxnoe commented Aug 25, 2023

@gschwefer If you enable the pre-commit hooks, this will automatically run when you commit, see:

https://ctapipe.readthedocs.io/en/latest/getting_started/index.html#developing-a-new-feature-or-code-change

@maxnoe
Copy link
Member

maxnoe commented Aug 29, 2023

This fixes a bug in the implementation of the full numeric pixel likelihood

What is the actual bug here?

I see that you re-introduce terms that are constant under minimization of the likelihood, which is probably the reason why they were left out.

I guess using factorial is pretty slow, here is the optimized implementation of the poisson related functions in numba-stats:
https://github.com/HDembinski/numba-stats/blob/main/src/numba_stats/poisson.py

@gschwefer
Copy link
Contributor Author

Ok, so these are a couple of different things: The terms in the neg_log_likelihood_numeric that I reintroduce are actually not simple constants to the log likelihood because they are summed over in the n_signal summation and depend on n. Therefore, when the log is taken, they aren't just separate additive terms and therefore make an actual difference. This was not caught in the unit tests because an absolute value was missing in the tests of the relative differences that allowed large negative values to pass.

@gschwefer
Copy link
Contributor Author

The factorial I'm happy to switch out.

@maxnoe
Copy link
Member

maxnoe commented Aug 29, 2023

The factorial I'm happy to switch out.

Measure before :)

@gschwefer
Copy link
Contributor Author

The factorial I'm happy to switch out.

Measure before :)

Good hint...the explicit implementation is indeed faster
image

@maxnoe
Copy link
Member

maxnoe commented Aug 29, 2023

For arrays of k or lambda, numba_stats.poisson.logpmf was a lot faster (factor of 2 or more), maybe worth introducing this additional dependency.

@gschwefer
Copy link
Contributor Author

For arrays of k or lambda, numba_stats.poisson.logpmf was a lot faster (factor of 2 or more), maybe worth introducing this additional dependency.

Ok, maybe I just don't see it, but my version of numba_stats (1.3.0) doesn't seem to support arrays as input for lambda. And that's what really needed here...

@maxnoe
Copy link
Member

maxnoe commented Aug 29, 2023

Indeed, I only tested arrays for k, I'll open an issue

@kosack kosack merged commit 8f6d5aa into cta-observatory:main Aug 31, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants