ctools

Implemented here:
https://cta-gitlab.irap.omp.eu/sadeh/ctools/tree/cssens_source_min_event_cuts

And yet another new feature in this pull request:
- in case edisp is set, only start using it in practice once the test flux is close to convergence.
This saves considerable CPU time.

Specifically, the modification follows (as part of the source flux calculation loop):

# if the flux has converged, check if the original value of
# edisp was set, and is not what has fo far been used
if edisp_orig and not self['edisp'].boolean():
    # set edisp on and continue the calculation
    self['edisp'] = True

    value = 'set edisp on after initial convergence withouth it'
    self._log_value(gammalib.TERSE, 'Converged result', value)
else:
    # finish the calculation after convergence
    break

Sadeh Iftach wrote:

Hi Jürgen,

I’m not sure I follow. Which duplicated issue are you referring to?
I’d like to keep track.

This was feature #2130 that was there since a while asking to incorporate minimum counts into the sensitivity calculation. Since no one was working on that I closed this feature now and keep this one alive.

In any case, once you go through your other updates, please check out
https://cta-gitlab.irap.omp.eu/sadeh/ctools/tree/cssens_source_min_event_cuts

To summarize all changes:
1. Add optional absolute minimal event count on test source
2. Add optional minimal threshold of test source counts with respect to non-source counts
3. Numerical stability updates for source prefactor and for the _predict_flux regression function
4. Avoid extremely slow unnecessary energy dispersion calculations
5. Add optional random seed initialization for event generation

I know this is a lot for a single pull request. I hope you can integrate the content.

No worries, I will go over the changes. I probably will propose shorter names for the new parameters as they are rather long compared to typical ctools names and a bit cryptic.

Knödlseder Jürgen wrote:

A duplicated feature was created with issue #3668. Since some code changes were done on this feature I close this issue now.

Closing was in fact not my intention here, sorry for the confusion. I closed the duplicated feature.

I implemented the following code to replace the median function:

    def _median(self, array):
        """ 
        Compute median value for an array

        Parameters
        ----------
        array : list of floats
            Array

        Returns
        -------
        median : float
            Median value of array
        """ 
        # Initialise median
        median = 0.0

        # Get number of elements in array
        n = len(array)

        # Continue only if there are elements in the array
        if n > 0:

            # Sort the array
            wrk_array = sorted(array)

            # Get median
            if n % 2 != 0:
                median = wrk_array[int(n/2)]
            else:
                median = (wrk_array[int((n-1)/2)] + wrk_array[int(n/2)]) / 2.0

        # Return median
        return median

0.0
Traceback (most recent call last):
  File "./test.py", line 19, in <module>
    print(statistics.median(array))
  File "/Users/jurgen/anaconda3/envs/ctools-1.7.4/lib/python3.8/statistics.py", line 430, in median
    raise StatisticsError("no median for empty data")
statistics.StatisticsError: no median for empty data

Here a few tests with different arrays:

The code seems to work as expected.

Median is a very low-level function. Shouldn’t it be implemented somewhere else so that it can be used by other ctools/cscripts in the future? For example GMath or GNdarray?

Tibaldo Luigi wrote:

Median is a very low-level function. Shouldn’t it be implemented somewhere else so that it can be used by other ctools/cscripts in the future? For example GMath or GNdarray?

Fully agree. The best place would probably be the GNdarray class, which should then also get mean() and std() methods. I created an issue for that: #3693.

Thanks, Jürgen.

Of course, please feel free to make any changes to variable names etc. to fit the ctools conventions.
At the moment I understand that I have no action items for myself, so please let me know if you’d like me to modify any code.

An outstanding open question is the choice of default settings.
In my opinion, the minimal 10 evt source count can be the new default (rather than the current default of 0). This complies with what everyone expects with respect to the “official” cta results (produced with EventDisplay), and shown on the public website.

The second condition on the background counts is not relevant in most cases (except for very long exposures). In my opinion it is also less rigorously defined, especially considering the full-enclosure IRFs.
I’m agnostic about wether or not this one should be turned on by default.

I agree that setting the default value to 10 is maybe what we want to do, although Gernot mentioned recently that these criteria may be revisited.

I’m basically done with the code integration, I did a bit of refactoring (introducing a method cssens._simulate_events that avoids having a loop over source and background in cssens._sim_evt_excess), modified slightly the logging (so that values all fit in the expected length and that detailed information is only written when chatter is in EXPLICIT) and changed the parameters to shorter names:

mincounts, i, h, 0,,, "Minimum number of source counts" 
bkgexcess, r, h, 0.0,,, "Background uncertainty fraction" 
bkgrad,    r, h, 0.33,,, "Background radius (deg)" 
seed,      i, h, 1,,, "Initial random number generator seed" 

I also updated the reference manual to add the new parameters. The documentation of all tools is in doc/source/users/reference_manual/ which is also used to build the help text that you get with cssens --help.

Great!

One final note, not really critical...
I tried to explicitly use “non-source” rather than “background” terminology for

bkgexcess, r, h, 0.0,,, "Background uncertainty fraction" 
bkgrad,    r, h, 0.33,,, "Background radius (deg)" 

I understand, but I don’t think that a user will understand what non-source means. We may clarify in the documentation that “background” means everything else than the source (a bit like for the DM community).

Great.
Thanks again for the quick integration of this pull request.