Repository: cta-observatory/ctapipe
Pull request #2543 (2024-04-24 15:14:41)

Variance extractor
This pull request is to add a new ImageExtractor subclass to generate images from the variance of waveforms that will be needed for the pointing calibration using stars. This request is related to issue #2542 .

kosack: If I understand correctly, what you are basically doing here is computing an event-wise pedestal variance. In the data model (see the original top-level data model document), there are a few places you can put pedestals:

- _dl1.event.telescope.pedestal_: pedestals measured for each event (usually in a time window as @maxnoe mentioned, since you can for example just use time slices outside the trace integration window, even for Cherenkov events). This is also computed by some cameras directly in hardware as dl0 information.
- _dl1.monitoring.telescope.pedestal_: is the pedestal computed from an ensemble of pedestal events, or from Cherenkov events with the shower cut out  (this is also in the DL0.monitoring.telescope.pedestal, computed by the cameras) computed every few seconds.
- dl1.service.telescope.pedestal: pedestals measured from some standard run that are applied for e.g. RTA

So I think what you are computing is the first one, and you should store it as pedestal information:  It is in fact a second measurement of the pedestal and its variance, and could be compared to the other methods (even if the use case you are describing here is for doing star tracking, star tracking is computed from the pedestal variances, which can be over many time scales)
fwerner: In addition to what has been said I'd suggest to clarify the purpose of this extractor: it is calculating the waveform _sample_ variance, which should be reflected in the name/docstring to communicate that it has a complicated relationship to charge.
ctoennis: The thing we want to calculate is the variance of the waveforms specifically of interleaved pedestal events and not the variance of the pedestal in units of charge. 

@kosack : What we get is a sort of image based on the variance of the waveform in each pixel. What we get is not a pedestal value, so putting it in a pedestal container is not the best option i think. 

I think it is best to write a container appropriate for these variance images similar to the DL1CameraContainer
ctoennis: I just added a custom continer for the variance images. I also compared the variance calculation using just a plain variance and a variance of integrated chunks of the waveforms using the FixedWindowSum extractor in the plot below.
![comparemethod](https://github.com/cta-observatory/ctapipe/assets/63591332/76a85d09-0db9-4f04-b476-3238126129e1)

The x axis shows the number of images used in calculating an average variance value that is shown on the y axis. 

Though the methods are generally statistically equivalent the average variance becomes stable with fewer images when using a simple variance over the variance from integrated chunks.

maxnoe: > Though the methods are generally statistically equivalent the average variance becomes stable with fewer images when using a simple variance over the variance from integrated chunks.

That's not surprising, since you have just more samples. But you change *of what you compute the variance*
ctoennis: Hi all, can we restart the discussion here? I made a small update to the extractor to contain timing information that will later be needed in the corresponding subclass of StatisticsExtractor that i am working on with Tjark. 

There was a decision needed as to where the variance images will be attached to. Should this find a place somewhere in dl1.event.telescope.pedestal or something similar?
maxnoe: Hi @ctoennis the lint step in the CI is failing due to some typos and wrong formatting.

To check locallly, make sure to install the pre-commit hooks:
https://ctapipe.readthedocs.io/en/latest/developer-guide/getting-started.html#installing-ctapipe-in-development-mode
ctao-dpps-sonarqube[bot]: ![Passed](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/QualityGateBadge/passed.svg?sanitize=true)

# Analysis Details
## 1 Issue
- ![Bug](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/bug.svg?sanitize=true) 0 Bugs
- ![Vulnerability](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/vulnerability.svg?sanitize=true) 0 Vulnerabilities
- ![Code Smell](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/code_smell.svg?sanitize=true) 1 Code Smell

## Coverage and Duplications
- ![Coverage](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/CoverageChart/100.svg?sanitize=true) 100.00% Coverage (93.60% Estimated after merge)
- ![Duplications](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/Duplications/3.svg?sanitize=true) 0.00% Duplicated Code (0.60% Estimated after merge)

**Project ID:** cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

[View in SonarQube](https://sonar-cta-dpps.zeuthen.desy.de/dashboard?id=cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs&pullRequest=2543)


ctao-dpps-sonarqube[bot]: ![Passed](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/QualityGateBadge/passed.svg?sanitize=true)

# Analysis Details
## 1 Issue
- ![Bug](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/bug.svg?sanitize=true) 0 Bugs
- ![Vulnerability](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/vulnerability.svg?sanitize=true) 0 Vulnerabilities
- ![Code Smell](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/code_smell.svg?sanitize=true) 1 Code Smell

## Coverage and Duplications
- ![Coverage](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/CoverageChart/100.svg?sanitize=true) 100.00% Coverage (94.20% Estimated after merge)
- ![Duplications](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/Duplications/3.svg?sanitize=true) 0.00% Duplicated Code (0.60% Estimated after merge)

**Project ID:** cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

[View in SonarQube](https://sonar-cta-dpps.zeuthen.desy.de/dashboard?id=cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs&pullRequest=2543)


ctoennis: i fixed all the isses with the spelling, code checks etc.
ctoennis: i can remove trigger time. But i am not sure how to proceed on the difference to the ImageExtractor API. I could make it that the variance image is stored in a DL1ImageContainer and that it complies with the API or make it a separate class that doesn't inherit from the ImageExtractor. 

I think what i generate here is a type of image as i try to look for stars in the camera FOV with them. But if it makes more sense to you it can be separate from the ImageExtractors.
maxnoe: The question is if it is possible to adhere to the API that is expected from the image extractors: if yes, no problem, let it be an image extractor. But that means for example the call above should work and correctly write out the data to the HDF5 file.

If you need the different API, it should be its own thing and we need to think about how to integrate it into the process tool for example.
kosack: > The question is if it is possible to adhere to the API that is expected from the image extractors: if yes, no problem, let it be an image extractor. But that means for example the call above should work and correctly write out the data to the HDF5 file.

I think what this class does is not the same as an `ImageExtractor`, or am I missing something? `ImageExtractors` take a "movie" (a set of waveforms per pixel for one event) for a single event and convert them to  set of time-independent dependent images that must include the  _peak time_ and _integral_ over a window, but can also optionally inlcude any other parameters of the waveform like the pulse width or variabce.    Generally, there is only one ImageExtractor that is run at a time.    That doesn't seem to be what you have here.
maxnoe: > I think what this class does is not the same as an ImageExtractor, or am I missing something?

You could think of this as an `ImageExtractor` that uses the variance of each waveform to produce the image and doesn't produce a `peak_time`.
mexanick: I believe, it is an ImageExtractor, the `__call__()` signature is of course wrong here. It should copy the one of the base extractor:
```
@abstractmethod
    def __call__(
        self, waveforms, tel_id, selected_gain_channel, broken_pixels
    ) -> DL1CameraContainer:
```

Only `waveforms` are then used, the other parameters are simply ignored.

The only issue, that has to be resolved, is the return type. `DL1CameraContainer` in its current shape contains a number of fields that are not needed in this case (actually, everything except `image`). 

I'd propose to add an intermediate container:

```python
class CameraImageContainer(Container):
    """
    Storage of output of generic camera images (any reduction of waveforms).
    """

    image = Field(
        None,
        "Numpy array of camera image, after waveform extraction."
        "Shape: (n_pixel) if n_channels is 1 or data is gain selected"
        "else: (n_channels, n_pixel)",
    )

class DL1CameraContainer(CameraImageContainer):
    """
    Storage of output of camera calibration e.g the final calibrated
    image in intensity units and the pulse time.
    """
    peak_time = Field(
        None,
        "Numpy array containing position of the peak of the pulse as determined by "
        "the extractor."
        "Shape: (n_pixel) if n_channels is 1 or data is gain selected"
        "else: (n_channels, n_pixel)",
    )
    # rest of the DL1CameraContainer fields
    
class VarianceContainer(CameraImageContainer):
    """
    Storage of output of camera variance image.

    Stores the variance of waveform in each pixel, composed as an image.
    """
    variance_method = Field(
        VarianceType.SIMPLE,
        "Method by which the variance was calculated"
        "This can either be a plain variance"
        "or a variance of integrated samples",
        type=VarianceType,
    )
    # I don't think you need anything else here. The trigger time you extract from elsewhere in the event tree.
```
and change the return type of base `ImageExtractor` to `CameraImageContainer`
maxnoe: If it's just that the `VarianceExtractor` doesn't provide a `peak_time`, I think that's fine, just leave it `None` (the default value).

Adding an additional field to the `DL1CameraContainer` that's an "image_type" sounds like a good idea.

I would avoid having extractors that have different return types.
mexanick: Following a private discussion with @maxnoe , the new Container class is not needed. 
@ctoennis just use the `DL1CameraContainer` as a return type, setting the extractor type as a container meta keyword. You only need to fill `Image` field and leave all others as defaults (`None`). In this case, the serialization will be done correctly and only the image that you need will be stored. I.e. something like this (assuming you change VarianceType enum to WAVEFORM):

```python
def __call__(
        self, waveforms, tel_id, selected_gain_channel, broken_pixels
    ) -> DL1CameraContainer:
    container = DL1CameraContainer(image=np.nanvar(waveforms, dtype="float32", axis=2))
    container.meta["ExtractionMethod"] = str(VarianceType.WAVEFORM)
    return container
```
    
maxnoe: Since this is also somewhat connected to #1243 

I.e. we could introduce a new class e.g. `WaveformProperties` that runs after the `CameraCalibrator` and computes stuff like variance, FWHM, etc.
mexanick: > Since this is also somewhat connected to #1243
> 
> I.e. we could introduce a new class e.g. `WaveformProperties` that runs after the `CameraCalibrator` and computes stuff like variance, FWHM, etc.

While this might be a good idea, I'd save it for a separate refactoring PR and merge the current (small) addition as is.  

Patch:
diff --git a/docs/changes/2543.feature.rst b/docs/changes/2543.feature.rst
new file mode 100644
index 00000000..81544985
--- /dev/null
+++ b/docs/changes/2543.feature.rst
@@ -0,0 +1,2 @@
+A new ImageExtractor called ``VarianceExtractor`` was added
+An Enum class was added to containers.py that is used in the metadata of the VarianceExtractor output
diff --git a/src/ctapipe/containers.py b/src/ctapipe/containers.py
index 9a5b39f4..6dbf7404 100644
--- a/src/ctapipe/containers.py
+++ b/src/ctapipe/containers.py
@@ -166,6 +166,15 @@ class EventType(enum.Enum):
     UNKNOWN = 255
 
 
+class VarianceType(enum.Enum):
+    """Enum of variance types used for the VarianceContainer"""
+
+    # Simple variance of waveform
+    WAVEFORM = 0
+    # Variance of integrated samples of a waveform
+    INTEGRATED = 1
+
+
 class PixelStatus(enum.IntFlag):
     """
     Pixel status information
@@ -510,7 +519,6 @@ class DL1CameraContainer(Container):
             "pass only was returned."
         ),
     )
-
     parameters = Field(
         None, description="Image parameters", type=ImageParametersContainer
     )
diff --git a/src/ctapipe/image/extractor.py b/src/ctapipe/image/extractor.py
index ce772744..88e7e683 100644
--- a/src/ctapipe/image/extractor.py
+++ b/src/ctapipe/image/extractor.py
@@ -13,6 +13,7 @@ __all__ = [
     "NeighborPeakWindowSum",
     "BaselineSubtractedNeighborPeakWindowSum",
     "TwoPassWindowSum",
+    "VarianceExtractor",
     "extract_around_peak",
     "extract_sliding_window",
     "neighbor_average_maximum",
@@ -33,7 +34,10 @@ from numba import float32, float64, guvectorize, int64, njit, prange
 from scipy.ndimage import convolve1d
 from traitlets import Bool, Int
 
-from ctapipe.containers import DL1CameraContainer
+from ctapipe.containers import (
+    DL1CameraContainer,
+    VarianceType,
+)
 from ctapipe.core import TelescopeComponent
 from ctapipe.core.traits import (
     BoolTelescopeParameter,
@@ -1297,6 +1301,19 @@ class TwoPassWindowSum(ImageExtractor):
         )
 
 
+class VarianceExtractor(ImageExtractor):
+    """Calculate the variance over samples in each waveform."""
+
+    def __call__(
+        self, waveforms, tel_id, selected_gain_channel, broken_pixels
+    ) -> DL1CameraContainer:
+        container = DL1CameraContainer(
+            image=np.nanvar(waveforms, dtype="float32", axis=2),
+        )
+        container.meta["ExtractionMethod"] = str(VarianceType.WAVEFORM)
+        return container
+
+
 def deconvolution_parameters(
     camera: CameraDescription,
     upsampling: int,


Test patch:
diff --git a/src/ctapipe/image/tests/test_extractor.py b/src/ctapipe/image/tests/test_extractor.py
index 067a0780..2b7fd1a5 100644
--- a/src/ctapipe/image/tests/test_extractor.py
+++ b/src/ctapipe/image/tests/test_extractor.py
@@ -19,6 +19,7 @@ from ctapipe.image.extractor import (
     NeighborPeakWindowSum,
     SlidingWindowMaxSum,
     TwoPassWindowSum,
+    VarianceExtractor,
     __filtfilt_fast,
     adaptive_centroid,
     deconvolve,
@@ -281,6 +282,17 @@ def test_extract_around_peak_charge_expected():
     assert_equal(charge, n_samples)
 
 
+def test_variance_extractor(toymodel):
+    _, subarray, _, _, _, _ = toymodel
+    # make dummy data with known variance
+    rng = np.random.default_rng(0)
+    var_data = rng.normal(2.0, 5.0, size=(2, 1855, 5000))
+    extractor = ImageExtractor.from_name("VarianceExtractor", subarray=subarray)
+
+    variance = extractor(var_data, 0, None, None).image
+    np.testing.assert_allclose(variance, np.var(var_data, axis=2), rtol=1e-3)
+
+
 @pytest.mark.parametrize("toymodels", camera_toymodels)
 def test_neighbor_average_peakpos(toymodels, request):
     waveforms, subarray, tel_id, _, _, _ = request.getfixturevalue(toymodels)
@@ -408,6 +420,9 @@ def test_extractors(Extractor, toymodels, request):
             extractor(waveforms, tel_id, selected_gain_channel, broken_pixels)
         return
 
+    if Extractor is VarianceExtractor:
+        return
+
     dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels)
     assert dl1.is_valid
     if dl1.image.ndim == 1:
@@ -423,7 +438,7 @@ def test_extractors(Extractor, toymodels, request):
 @pytest.mark.parametrize("Extractor", extractors)
 def test_integration_correction_off(Extractor, toymodels, request):
     # full waveform extractor does not have an integration correction
-    if Extractor is FullWaveformSum:
+    if Extractor in (FullWaveformSum, VarianceExtractor):
         return
 
     (
@@ -714,8 +729,11 @@ def test_dtype(Extractor, subarray):
     n_channels, n_pixels, _ = waveforms.shape
     broken_pixels = np.zeros((n_channels, n_pixels), dtype=bool)
     dl1 = extractor(waveforms, tel_id, selected_gain_channel, broken_pixels)
+
+    if Extractor is not VarianceExtractor:
+        assert dl1.peak_time.dtype == np.float32
+
     assert dl1.image.dtype == np.float32
-    assert dl1.peak_time.dtype == np.float32
 
 
 def test_global_peak_window_sum_with_pixel_fraction(subarray):