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

Lhfit highsignal gaussian approximation #9

Merged
merged 6 commits into from Jan 6, 2021
Merged

Lhfit highsignal gaussian approximation #9

merged 6 commits into from Jan 6, 2021

Conversation

gabemery
Copy link

@gabemery gabemery commented Jan 6, 2021

Description

A Gaussian approximation of the pixel likelihood is added for high luminosity pixels. High luminosity pixels are identified as pixel requiring to consider a number of possible photo-electron above the parameter n_peak in the likelihood computation.

The goal is to reduce the computation time required for such pixels as the range of Poisson contributing terms for the exact likelihood computation increase with the expected signal. It also limits the number of pre-computed factorials required.

WIP :

  • Current implementation (7ff716e) only used for comparison. The Gaussian approximation is not used for the fitting procedure. Instead, the limit on the maximum value of number of photo-electron for the likelihood sum is increased by a factor 5. The comparison can thus be performed for pixels with n_peak < kmax < 5* n_peak

Related Issue

Related to #7

How Has This Been Tested?

Currently not tested

Types of changes

  • Bug fix (non-breaking change which fixes an issue)
  • New feature (non-breaking change which adds functionality)
  • Breaking change (fix or feature that would cause existing functionality to change)

Checklist:

  • My code follows the code style of this project.
  • My change requires a change to the documentation.
  • I have updated the documentation accordingly.
  • I have read the CONTRIBUTING document.
  • I have added tests to cover my changes.
  • All new and existing tests passed.

@gabemery gabemery added enhancement New feature or request in progress labels Jan 6, 2021
@gabemery gabemery requested a review from mexanick January 6, 2021 10:58
@gabemery gabemery self-assigned this Jan 6, 2021
@gabemery gabemery added this to In progress in LH-Fit via automation Jan 6, 2021
@pep8speaks
Copy link

pep8speaks commented Jan 6, 2021

Hello @gabemery! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

Line 639:37: E261 at least two spaces before inline comment
Line 639:38: E262 inline comment should start with '# '
Line 639:80: E501 line too long (84 > 79 characters)
Line 656:80: E501 line too long (106 > 79 characters)

Line 65:80: E501 line too long (82 > 79 characters)
Line 66:80: E501 line too long (89 > 79 characters)
Line 512:80: E501 line too long (103 > 79 characters)
Line 586:13: E265 block comment should start with '# '
Line 586:80: E501 line too long (82 > 79 characters)
Line 620:80: E501 line too long (87 > 79 characters)
Line 623:80: E501 line too long (87 > 79 characters)
Line 632:9: E303 too many blank lines (2)
Line 634:80: E501 line too long (94 > 79 characters)
Line 636:80: E501 line too long (94 > 79 characters)
Line 637:80: E501 line too long (105 > 79 characters)
Line 640:80: E501 line too long (85 > 79 characters)
Line 643:80: E501 line too long (85 > 79 characters)
Line 646:80: E501 line too long (114 > 79 characters)

Comment last updated at 2021-01-06 15:46:44 UTC

Copy link

@mexanick mexanick left a comment

Choose a reason for hiding this comment

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

I thought you would leave the non-saturated processing intact and change to gaussian approximation only for saturated events. In fact, in the code I see that you're now computing gaussian, poissonian and full sum for every event, if I'm not mistaken... I would actually start with only doing this for saturated events (nb, the approximation is expected to work better for brighter events/pixels). Then, check the convergence on less bright events in bins of n_photopeaks (for faint events I would expect the approximation to perform poorly).

Also, please make the PEP8 bot happy again :)

@@ -65,6 +65,7 @@
from ..io.io import dl1_params_lstcam_key

logger = logging.getLogger(__name__)
logger.setLevel(DEBUG)
Copy link

Choose a reason for hiding this comment

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

This is a bit dangerous... I would much prefer the logging level to be passed as a parameter, otherwise one needs to modify a source code to change it.

Copy link
Author

Choose a reason for hiding this comment

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

Ok I will try

#May need rework after accelaration addition
photoelectron_peak = np.arange(n_peaks, dtype=np.int)
# May need rework after accelaration addition
photoelectron_peak = np.arange(5*n_peaks, dtype=np.int) # the 5* is more testing only
Copy link

Choose a reason for hiding this comment

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

Please expand comment on what's this for

@@ -574,32 +576,33 @@ def log_pdf(self, charge, t_cm, x_cm, y_cm, width,
it = it + 1
kmin[kmin<0] = 0
kmax = np.ceil(kmax)
if len(kmin)==0 or len(kmax)==0:
kmin,kmax=0,len(self.log_k)
mask = (kmax <= len(self.log_k) / 5) # /5 for testing only
Copy link

Choose a reason for hiding this comment

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

Is this hardcoded constant related to l511? Please avoid such things. If you need, use a local variable

Copy link
Author

Choose a reason for hiding this comment

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

Yes it is the same 5, I can replace it. But it should disappear before merging as it is only used to have a common range for the previous likelihood computation and the new Gaussian approximation.

LH-Fit automation moved this from In progress to Review in progress Jan 6, 2021
@gabemery
Copy link
Author

gabemery commented Jan 6, 2021

The approximation cannot be employed event wise. This is due to low signal pixels on the side of the shower. It thus need to have a criteria based on each pixels. Additionally, the likelihood uses a modeled value for the estimation of the Poisson probability of k photo-electron and not the measured value.

Here you are a bit confused by the naming of the variables. Both log_poisson and log_gaussian are part of the sum likelihood. the Gaussian approximated log likelihood is log_pdf_HL. So the sum is computed for all pixels (limiting kmax to 5* n_peak, and thus should be accurate as long as this limit is not reached) while the Gaussian approximation is computed only for pixels requiring kmax > n_peak.

The next step would be to ignore pixels requiring kmax > n_peak in the estimation of the sum range, and to compute only the Gaussian approximation for such pixels. At the same time, the sum would be only computed for pixels requiring kmax < n_peak.

…tor used to test the compatibility of the new gaussian approximation with the previous computation
@mexanick
Copy link

mexanick commented Jan 6, 2021

Apparently a typo was introduced at some point. As reported by Travis:

____________ ERROR collecting lstchain/reco/tests/test_r0_to_dl1.py ____________
997../../../miniconda/envs/lst-dev/lib/python3.6/site-packages/_pytest/python.py:578: in _importtestmodule
998    mod = import_path(self.fspath, mode=importmode)
999../../../miniconda/envs/lst-dev/lib/python3.6/site-packages/_pytest/pathlib.py:531: in import_path
1000    importlib.import_module(module_name)
1001../../../miniconda/envs/lst-dev/lib/python3.6/importlib/__init__.py:126: in import_module
1002    return _bootstrap._gcd_import(name[level:], package, level)
1003<frozen importlib._bootstrap>:994: in _gcd_import
1004    ???
1005<frozen importlib._bootstrap>:971: in _find_and_load
1006    ???
1007<frozen importlib._bootstrap>:955: in _find_and_load_unlocked
1008    ???
1009<frozen importlib._bootstrap>:665: in _load_unlocked
1010    ???
1011../../../miniconda/envs/lst-dev/lib/python3.6/site-packages/_pytest/assertion/rewrite.py:170: in exec_module
1012    exec(co, module.__dict__)
1013lstchain/reco/tests/test_r0_to_dl1.py:3: in <module>
1014    from lstchain.reco.r0_to_dl1 import rescale_dl1_charge
1015E     File "/home/travis/build/cta-sst-1m/cta-lstchain/lstchain/reco/r0_to_dl1.py", line 668
1016E       )
1017E       ^
1018E   SyntaxError: invalid syntax

Please fix

@mexanick
Copy link

mexanick commented Jan 6, 2021

Ok, I will approve and merge it for the moment, but there's some refactoring required...

LH-Fit automation moved this from Review in progress to Reviewer approved Jan 6, 2021
LH-Fit automation moved this from Reviewer approved to Review in progress Jan 6, 2021
Copy link

@mexanick mexanick left a comment

Choose a reason for hiding this comment

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

____________ ERROR collecting lstchain/reco/tests/test_r0_to_dl1.py ____________
996lstchain/reco/tests/test_r0_to_dl1.py:3: in
from travis again:

997    from lstchain.reco.r0_to_dl1 import rescale_dl1_charge
998lstchain/reco/r0_to_dl1.py:68: in <module>
999    logger.setLevel(DEBUG)
1000E   NameError: name 'DEBUG' is not defined

Copy link

@mexanick mexanick left a comment

Choose a reason for hiding this comment

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

Travis build failure is now unrelated to this PR

LH-Fit automation moved this from Review in progress to Reviewer approved Jan 6, 2021
@mexanick mexanick merged commit 929645e into cta-sst-1m:develop Jan 6, 2021
LH-Fit automation moved this from Reviewer approved to Done Jan 6, 2021
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request in progress
Projects
Development

Successfully merging this pull request may close these issues.

None yet

3 participants