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

Add light curve upper limits #1456

Merged
merged 8 commits into from Jul 10, 2018
Merged

Add light curve upper limits #1456

merged 8 commits into from Jul 10, 2018

Conversation

@bkhelifi
Copy link
Contributor

@bkhelifi bkhelifi commented Jul 6, 2018

Hi,
After the work made for the CTA consortium meeting, I made this pull request for only the UL management.
As it is, the LC still produces NaN as the rest of the modifications within spectrum, proposed in the PR #1424, is not pushed.
B.

@bkhelifi bkhelifi added the feature label Jul 6, 2018
@cdeil cdeil assigned cdeil and unassigned registerrier Jul 6, 2018
@cdeil cdeil added feature and removed feature labels Jul 6, 2018
@cdeil cdeil added this to the 0.8 milestone Jul 6, 2018
@cdeil cdeil added this to To do in gammapy.time via automation Jul 6, 2018
@cdeil cdeil changed the title LC: add upper limites Add light curve upper limits Jul 6, 2018
Copy link
Member

@cdeil cdeil left a comment

@bkhelifi - Thank you for splitting this out into a separate PR!

Can you resolve the test fails (see inline comments)?
Or would you prefer I do it and merge?

@@ -37,6 +38,9 @@ def test_excess_error():
assert_allclose(excess_error(n_on=10, n_off=20, alpha=0.1), 3.1937439)
assert_allclose(excess_error(n_on=4, n_off=9, alpha=0.5), 2.5)

def test_Helene_ULs():
assert_allclose(Helene_ULs(excess=50, error=40, conf_level=99.73), 119.708038)

This comment has been minimized.

@cdeil

cdeil Jul 6, 2018
Member

I get a different value in the test?!
@bkhelifi - does the test pass for you locally?

python setup.py test -V -t gammapy/stats/tests/test_poisson.py
_____________________________________________________________________________ test_Helene_ULs _____________________________________________________________________________

    def test_Helene_ULs():
>       assert_allclose(Helene_ULs(excess=50, error=40, conf_level=99.73), 119.708038)
E       AssertionError: 
E       Not equal to tolerance rtol=1e-07, atol=0
E       
E       (mismatch 100.0%)
E        x: array(162.727686)
E        y: array(119.708038)

gammapy/stats/tests/test_poisson.py:42: AssertionError
============================================================= 1 failed, 18 passed, 1 xfailed in 1.14 seconds ==============================================================

This comment has been minimized.

@cdeil

cdeil Jul 6, 2018
Member

It looks like the paper has plenty of "reference values" on what the result should be?
Could you pick one of those and leave a one line comment in the test saying where to find the test case in the paper (e.g. Table 1, 2 or 3?)

This comment has been minimized.

@cdeil

cdeil Jul 6, 2018
Member

Using this command, I see that the elif excess < 0.: case is never executed by the tests:

python setup.py test -V -t gammapy/stats/tests/test_poisson.py --coverage
open htmlcov/index.html 

Please add at least one test case, preferably using a test case given in the paper.

This comment has been minimized.

@bkhelifi

bkhelifi Jul 9, 2018
Author Contributor

I did a stupid typo. Thanks a lot to have re-check that. Corrected.
I have added a test to make a wider coverage of the code


# Store measurements in a dict and return that
return useinterval, OrderedDict([
('time_min', Time(tmin, format='mjd')),
('time_max', Time(tmax, format='mjd')),
('flux', flux * u.Unit('1 / (s cm2)')),
('flux_err', flux_err * u.Unit('1 / (s cm2)')),

('flux_ul', flux_ul * u.Unit('1 / (s cm2)')),
('is_ul', sigma <= sigma_ul_thres)

This comment has been minimized.

@cdeil

cdeil Jul 6, 2018
Member

@bkhelifi - You're missing a comma here:

$ python setup.py test -V -t gammapy/time/tests/test_lightcurve.py
``

________________________________________________________________________ test_lightcurve_estimator ________________________________________________________________________

@requires_data('gammapy-extra')
@requires_dependency('scipy')
def test_lightcurve_estimator():
    spec_extract = spec_extraction()
    lc_estimator = LightCurveEstimator(spec_extract)

    # param
    intervals = []
    for obs in spec_extract.obs_list:
        intervals.append([obs.events.time[0], obs.events.time[-1]])

    model = PowerLaw(
        index=2.3 * u.Unit(''),
        amplitude=3.4e-11 * u.Unit('1 / (cm2 s TeV)'),
        reference=1 * u.TeV,
    )

    lc = lc_estimator.light_curve(
        time_intervals=intervals,
        spectral_model=model,
      energy_range=[0.5, 100] * u.TeV,
    )

gammapy/time/tests/test_lightcurve.py:206:


gammapy/time/lightcurve.py:672: in light_curve
useinterval, row = self.compute_flux_point(time_interval, spectral_model, energy_range, sigma_ul_thres)


self = <gammapy.time.lightcurve.LightCurveEstimator object at 0xb147e2da0>
time_interval = [, ]
spectral_model = PowerLaw(), energy_range = <Quantity [ 0.5, 100. ] TeV>, sigma_ul_thres = 3.0

def compute_flux_point(self, time_interval, spectral_model, energy_range, sigma_ul_thres=3.):
    """Compute one flux point for one time interval.

        Parameters
        ----------
        time_interval : `~astropy.time.Time`
            Time interval (2-element array, or a tuple of Time objects)
        spectral_model : `~gammapy.spectrum.models.SpectralModel`
            Spectral model
        energy_range : `~astropy.units.Quantity`
            True energy range to evaluate integrated flux (true energy)
        sigma_ul_thres : float
            Threshold on detection significance to compute Flux Upper Limits

        Returns
        -------
        useinterval : bool
            Is True if the time_interval produce a valid flux point
        measurements : dict
            Dictionary with flux point measurement in the time interval
        """
    tmin, tmax = time_interval[0], time_interval[1]
    livetime = 0
    alpha_mean = 0.
    alpha_mean_backup = 0.
    measured_excess = 0
    predicted_excess = 0
    n_on = 0
    n_off = 0
    useinterval = False

    # Loop on observations
    for t_index, obs in enumerate(self.obs_list):

        # discard observations not matching the time interval
        obs_start = obs.events.observation_time_start
        obs_stop = obs.events.observation_time_end
        if (tmin < obs_start and tmax < obs_start) or (tmin > obs_stop):
            continue
        useinterval = True
        # get ON and OFF evt list
        on, off, spec, e_reco = self._create_and_filter_onofflists(t_index=t_index, energy_range=energy_range,
                                                                   interval=[tmin, tmax], extra=True)
        n_on_obs = len(on.table)
        n_off_obs = len(off.table)

        # compute effective livetime (for the interval)
        if tmin >= obs_start and tmax <= obs_stop:
            # interval included in obs
            livetime_to_add = (tmax - tmin).to('s')
        elif tmin >= obs_start and tmax >= obs_stop:
            # interval min above tstart from obs
            livetime_to_add = (obs_stop - tmin).to('s')
        elif tmin <= obs_start and tmax <= obs_stop:
            # interval min below tstart from obs
            livetime_to_add = (tmax - obs_start).to('s')
        elif tmin <= obs_start and tmax >= obs_stop:
            # obs included in interval
            livetime_to_add = (obs_stop - obs_start).to('s')
        else:
            livetime_to_add = 0 * u.s

        # Take into account dead time
        livetime_to_add *= (1. - obs.observation_dead_time_fraction)

        # Compute excess
        obs_measured_excess = n_on_obs - spec.alpha * n_off_obs

        # Compute the expected excess in the range given by the user
        # but must respect the energy threshold of the observation
        # (to match the energy range of the measured excess)
        # We use the effective livetime and the right energy threshold
        e_idx = np.where(np.logical_and.reduce(
            (e_reco >= spec.lo_threshold,  # threshold
             e_reco <= spec.hi_threshold,  # threshold
             e_reco >= energy_range[0],  # user
             e_reco <= energy_range[-1])  # user
        ))[0]
        counts_predictor = CountsPredictor(
            livetime=livetime_to_add,
            aeff=spec.aeff,
            edisp=spec.edisp,
            model=spectral_model
        )
        counts_predictor.run()
        counts_predicted_excess = counts_predictor.npred.data.data[e_idx[:-1]]

        obs_predicted_excess = np.sum(counts_predicted_excess)

        # compute effective normalisation between ON/OFF (for the interval)
        livetime += livetime_to_add
        alpha_mean += spec.alpha * n_off_obs
        alpha_mean_backup += spec.alpha * livetime_to_add
        measured_excess += obs_measured_excess
        predicted_excess += obs_predicted_excess
        n_on += n_on_obs
        n_off += n_off_obs

    # Fill time interval information
    if useinterval:
        int_flux = spectral_model.integral(energy_range[0], energy_range[1])

        if n_off > 0.:
            alpha_mean /= n_off
        if livetime > 0.:
            alpha_mean_backup /= livetime
        if alpha_mean == 0.:  # use backup if necessary
            alpha_mean = alpha_mean_backup

        flux = measured_excess / predicted_excess.value
        flux *= int_flux

        # Gaussian errors, TODO: should be improved
        flux_err = int_flux / predicted_excess.value
        delta_excess = excess_error(n_on=n_on, n_off=n_off, alpha=alpha_mean)
        flux_err *= delta_excess

        flux_ul = -1
        sigma = significance_on_off(n_on=n_on, n_off=n_off, alpha=alpha_mean, method='lima')
        if sigma <= sigma_ul_thres:
            flux_ul = Helene_ULs(n_on - alpha_mean * n_off, delta_excess, 99.73) #3 sigma ULs
            flux_ul *= int_flux / predicted_excess.value
    else:
        flux = 0
        flux_err = 0
        flux_ul = -1
        sigma = 0.

    # Store measurements in a dict and return that
    return useinterval, OrderedDict([
        ('time_min', Time(tmin, format='mjd')),
        ('time_max', Time(tmax, format='mjd')),
        ('flux', flux * u.Unit('1 / (s cm2)')),
        ('flux_err', flux_err * u.Unit('1 / (s cm2)')),
        ('flux_ul', flux_ul * u.Unit('1 / (s cm2)')),
        ('is_ul', sigma <= sigma_ul_thres)
      ('livetime', livetime * u.s),
        ('alpha', alpha_mean),
        ('n_on', n_on),
        ('n_off', n_off),
        ('measured_excess', measured_excess),
        ('expected_excess', predicted_excess),
    ])

E TypeError: 'tuple' object is not callable

gammapy/time/lightcurve.py:832: TypeError
=================================================================== 1 failed, 16 passed in 6.83 seconds ===================================================================

This comment has been minimized.

@bkhelifi

bkhelifi Jul 9, 2018
Author Contributor

Of course! My configuration has not run all the tests, in particular the one of the class LightCurve. I missed indeed bugs, showing that the test functions are useful!
Corrected

Parameters
----------
excess : double

This comment has been minimized.

@cdeil

cdeil Jul 6, 2018
Member

Write always "float" instead of "double", that's what the type is called in Python.

This comment has been minimized.

@bkhelifi

bkhelifi Jul 9, 2018
Author Contributor

Sorry, a too fast copy/paste. Corrected

integral = (1. + erf(value_new))/2.

while integral < integral2:
value_old = value_new

This comment has been minimized.

@cdeil

cdeil Jul 6, 2018
Member

PyCharm shows me that this value_old is unused; you can delete this line and the result will be the same.

This code is really hard to understand.
It looks like the excess > 0 and < 0 are the same or a similar algorithm?
If yes, can you put the algorithm in a helper functino and call it twice, instead of duplicating the ~ 10 lines of the algorithm?

This comment has been minimized.

@bkhelifi

bkhelifi Jul 9, 2018
Author Contributor

The comment on value_old is accepted and its removal is done.
For the different cases on excess, the algo is slightly different and a difference should be made. I have however simplify a bit the code as there are common ligns. Concerning its complexity, the code is coming from a C lib and optimised in this respect. I have not changed it and I do not plan to optimise it in a short term.

@@ -413,6 +414,84 @@ def _significance_direct_on_off(n_on, n_off, alpha):

return significance

def Helene_ULs(excess, error, conf_level=95.45):
"""Compute Upper limit based on Helene et al., NIM 212 (1983) 319

This comment has been minimized.

@cdeil

cdeil Jul 6, 2018
Member

We usually give a URL to the paper in the docstring.
I think it's this one?
http://adsabs.harvard.edu/abs/1984NIMPA.228..120H
Can you please add it?

This comment has been minimized.

@bkhelifi

bkhelifi Jul 9, 2018
Author Contributor

Done

@cdeil
cdeil approved these changes Jul 10, 2018
Copy link
Member

@cdeil cdeil left a comment

Some more pair coding with @bkhelifi just now. Should be OK to merge if CI passes.

@cdeil cdeil merged commit fd19724 into gammapy:master Jul 10, 2018
1 of 2 checks passed
1 of 2 checks passed
continuous-integration/travis-ci/pr The Travis CI build is in progress
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
gammapy.time automation moved this from To do to Done Jul 10, 2018
@bkhelifi bkhelifi deleted the bkhelifi:LC_Err branch Oct 27, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
gammapy.time
  
Done
Linked issues

Successfully merging this pull request may close these issues.

None yet

3 participants