#641 2D channel and energy edges, with filtering options

[pstammler]

In this draft, a new structure is proposed for the extracted energies and channels. New features include the use of some datasets statistical errors and optional filtering regarding the quality of the data or the grouping of channels.

The approach for the grouping, however, is not complete yet. For now, it is working with channels grouped to get a minimum number of counts. Grouping with binning factors shall be added in future releases...

[sguillot]

Thank you Pierre. We will look at this important PR.
Can you tell us if this completely supersedes #503 ?

[thjsal]

I did not look yet into all the details, but it seems clear that this a currently backward-incompatible update. If keeping it like this, we should do a major version bump to X-PSI v4. Or modify all the functions so that they are backward-compatible.

[thjsal]

I think if you modify the code, it is not a "review". Thus, I suggest saying either "Changed" or "Modified" instead.

[pstammler]

The changelog file has been updated, with the ability to choose the structure of the edges.

[thjsal]

I see you changed "Review" to "Propose" but I think at least before merging this should be changed to something else, since it will not be a proposal anymore at that point. I think the changelog should describe what has changed. Not what has been proposed to be changed.

[pstammler]

Thank you Pierre. We will look at this important PR. Can you tell us if this completely supersedes #503 ?

This PR gives extended features from #503 and #624. So yes, it supersedes them.

[thjsal]

Thanks for the update @pstammler! The code appears now more backward-compatible to me. However, if executing the current Synthetic data notebook (https://xpsi-group.github.io/xpsi/Synthetic_data.html), I get the following error:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[19], line 1
----> 1 signal = CustomSignal(data = _data,
      2                         instrument = Instrument,
      3                         background = background,
      4                         interstellar = None,
      5                         cache = True,
      6                         prefix='Instrument')
      8 likelihood = xpsi.Likelihood(star = star, signals = signal,
      9                              num_energies=384,
     10                              threads=1,
     11                              externally_updated=False,
     12                              prior = prior)                             
     14 for h in hot_spot.objects:

Cell In[6], line 16, in CustomSignal.__init__(self, support, *args, **kwargs)
     15 def __init__(self, support = None, *args, **kwargs):
---> 16     super(CustomSignal, self).__init__(*args, **kwargs)

File ~/xpsi/xpsi_group/xpsi/xpsi/Signal.py:155, in Signal.__init__(self, data, instrument, background, interstellar, support, photosphere_prefix, cache, bounds, values, stokes, min_channel, max_channel, tolerance, *args, **kwargs)
    152     else:
    153         self._instrument.channels=self._instrument.channels[_np.where(self._data._quality[self._instrument.channels]==0)[0]]
--> 155 if self._data._grpdata == True:
    156     try:
    157         assert self._instrument.energy_edges.ndim == 2

AttributeError: 'SynthesiseData' object has no attribute '_grpdata'

[thjsal]

The general idea looks good to me. But maybe a good idea would be to mention an example of how the 2D energy edges should look like (e.g. in the docstring). If I understand right, they should be something like

energy_edges = [[E0_low, E1_low, E2_low], [E0_high,E1_high,E2_high]]
I also did not understand why this is marked as draft, so I clicked "Ready for a review".

[sguillot]

TO BE TESTED

Example scripts

• main.py
• TestRun_BB.py
• TestRun_Num.py
• TestRun_NumBeam.py
• TestRun_Pol.py
• TestRun_PolNum_split_1spot.py
• TestRun_PolNum_split_inference.py
• TestRun_AMXP.py

Notebooks

• XPSI 101
• Modelling Notebook --> CRASHES IN THE SYNTHETIC PART (See comment by @thjsal)
• Modelling Without Statistics Notebook
• Hot region complexity Notebook
• Global Surface emission
• Surface Radiation Field Tools
• Polarisation Notebook
• Synthesis Notebook
• Instrument Notebook

[pstammler]

Thanks for the update @pstammler! The code appears now more backward-compatible to me. However, if executing the current Synthetic data notebook (https://xpsi-group.github.io/xpsi/Synthetic_data.html), I get the following error:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[19], line 1
----> 1 signal = CustomSignal(data = _data,
      2                         instrument = Instrument,
      3                         background = background,
      4                         interstellar = None,
      5                         cache = True,
      6                         prefix='Instrument')
      8 likelihood = xpsi.Likelihood(star = star, signals = signal,
      9                              num_energies=384,
     10                              threads=1,
     11                              externally_updated=False,
     12                              prior = prior)                             
     14 for h in hot_spot.objects:

Cell In[6], line 16, in CustomSignal.__init__(self, support, *args, **kwargs)
     15 def __init__(self, support = None, *args, **kwargs):
---> 16     super(CustomSignal, self).__init__(*args, **kwargs)

File ~/xpsi/xpsi_group/xpsi/xpsi/Signal.py:155, in Signal.__init__(self, data, instrument, background, interstellar, support, photosphere_prefix, cache, bounds, values, stokes, min_channel, max_channel, tolerance, *args, **kwargs)
    152     else:
    153         self._instrument.channels=self._instrument.channels[_np.where(self._data._quality[self._instrument.channels]==0)[0]]
--> 155 if self._data._grpdata == True:
    156     try:
    157         assert self._instrument.energy_edges.ndim == 2

AttributeError: 'SynthesiseData' object has no attribute '_grpdata'

Maybe you should add the grpdata and the origdata parameters from the Data module in your SyntheticData class.

[pstammler]

Like this:

class SynthesiseData(xpsi.Data):
""" Custom data container to enable synthesis. """

def __init__(self, channels, phases, first, last, grpdata=False, origdata=None):

    self.channels = channels
    self._phases = phases
    self._grpdata=grpdata
    if self._grpdata==True:
        self._origdata=origdata
    else:
        self._origdata=None

    try:
        self._first = int(first)
        self._last = int(last)
    except TypeError:
        raise TypeError('The first and last channels must be integers.')
    if self._first >= self._last:
        raise ValueError('The first channel number must be lower than the '
                         'the last channel number.')

[pstammler]

The general idea looks good to me. But maybe a good idea would be to mention an example of how the 2D energy edges should look like (e.g. in the docstring). If I understand right, they should be something like

energy_edges = [[E0_low, E1_low, E2_low], [E0_high,E1_high,E2_high]] I also did not understand why this is marked as draft, so I clicked "Ready for a review".

Yes that's the general idea.
I marked this PR as a draft because of the potential impact of the changes in X-PSI, as you noted for the first version. Besides, there is only one kind of grouping available here.

[thjsal]

Like this:

class SynthesiseData(xpsi.Data): """ Custom data container to enable synthesis. """

def __init__(self, channels, phases, first, last, grpdata=False, origdata=None):

    self.channels = channels
    self._phases = phases
    self._grpdata=grpdata
    if self._grpdata==True:
        self._origdata=origdata
    else:
        self._origdata=None

    try:
        self._first = int(first)
        self._last = int(last)
    except TypeError:
        raise TypeError('The first and last channels must be integers.')
    if self._first >= self._last:
        raise ValueError('The first channel number must be lower than the '
                         'the last channel number.')

Or perhaps we can just add super().__init__(...) to the beginning of the init function of SynthesiseData class, so that it will inherit directly the properties from the __init__ function of the Data class. That would prevent some duplication of code. But this is still to be tested. Anyway, this means that the old Syntesise method syntax will not work anymore, and maybe it will be good to warn about that in the Changelog.

[sguillot]

Or perhaps we can just add super().__init__(...) to the beginning of the init function of SynthesiseData class, so that it will inherit directly the properties from the __init__ function of the Data class. That would prevent some duplication of code. But this is still to be tested. Anyway, this means that the old Syntesise method syntax will not work anymore, and maybe it will be good to warn about that in the Changelog.

class SynthesiseData(xpsi.Data):
    """ Custom data container to enable synthesis. """
    def __init__(self, channels, phases, first, last):
        counts = np.zeros((len(channels), len(phases)-1))
        super().__init__(counts, channels, phases, first, last, exposure_time=1.0)

This seems to work (but needs to be tested thoroughly). The xpsi.Data class requires a counts array with the right size, and an exposure_time. They are not used for the SyntheticData so, they are just initialise to zeros and 1.0.