Variance calibration #2636

ctoennis · 2024-11-04T10:18:19Z

This PR adds a trait to the CameraCalibrator called "image_calibration_type". This trait allows to switch the method of calibration between "charge" and "variance", with the former being the calibration method for the existing image extractors and the latter being the method to calibrate variance images.

I also made some small adjustment to the variance extractor so that it works in the CameraCalibrator.

…se with Calibrator

TjarkMiener · 2024-11-04T13:16:00Z

thanks for putting this together @ctoennis! In case I'm not missing something, I think we can (maybe should) apply this sqrt operation before writing the coefficients. Therefore we do not need the different cases here in the code. Would that work?

ctoennis · 2024-11-04T13:20:07Z

@TjarkMiener Yes, that should work. The calibration events would only need the sqrt coeficcients. Only thing is that this needs to be documented somewhere.

ctoennis · 2024-11-04T13:21:04Z

Though, the changes to the variance extractor should then still be kept so that it works with the calibrator.

mexanick · 2024-11-04T15:20:12Z

src/ctapipe/calib/camera/calibrator.py

@@ -61,6 +62,11 @@ class CameraCalibrator(TelescopeComponent):

    image_extractor_type: str
        The name of the ImageExtractor subclass to be used for image extraction
+
+    image_calibration_type: str


I don't think we need this. Instead, you can check the type of the image extractor.

if isinstance(VarianceExtractor, extractor): #apply squared gain calibration

OK, I'll do that

mexanick · 2024-11-04T15:20:26Z

src/ctapipe/calib/camera/calibrator.py

@@ -73,6 +79,15 @@ class CameraCalibrator(TelescopeComponent):
        help="Name of the ImageExtractor subclass to be used.",
    ).tag(config=True)

+    image_calibration_type = CaselessStrEnum(


mexanick · 2024-11-04T15:21:27Z

src/ctapipe/calib/camera/calibrator.py

+            if (
+                self.apply_peak_time_shift.tel[tel_id]
+                and remaining_shift is not None
+                and self.image_calibration_type == "charge"


Do we really need to alter this condition?

The variance extractor has no peak time, so applying the peak time shift otherwise throws an error

Just check if dl1.peak_time is not None instead of a specific extractor.

@maxnoe that makes sense

mexanick · 2024-11-04T15:21:39Z

src/ctapipe/calib/camera/calibrator.py

@@ -292,13 +311,21 @@ def _calibrate_dl1(self, event, tel_id):
            and dl1_calib.absolute_factor is not None
        ):
            if selected_gain_channel is None:
-                dl1.image *= dl1_calib.relative_factor / dl1_calib.absolute_factor
+                if self.image_calibration_type == "charge":


mexanick · 2024-11-04T15:22:32Z

src/ctapipe/calib/camera/tests/test_calibrator.py

+    calibrator(example_event)
+    image = example_event.dl1.tel[tel_id].image
+    assert image is not None
+    assert image.shape == (1764,)


from where this shape comes? What telescope is this?

That is from the example event. Not sure what telescope this is

FlashCam from the number of pixels

mexanick · 2024-11-04T15:24:08Z

src/ctapipe/image/extractor.py

@@ -1308,7 +1308,7 @@ def __call__(
        self, waveforms, tel_id, selected_gain_channel, broken_pixels
    ) -> DL1CameraContainer:
        container = DL1CameraContainer(
-            image=np.nanvar(waveforms, dtype="float32", axis=2),
+            image=np.nanvar(waveforms, dtype="float32", axis=2)[0],


We shall be able to work with multiple gains, why do you introduce this change?

np.nanvar was giving me an extra dimension on top of the gain. I think I can fix this with the keepdims parameter

TjarkMiener · 2024-11-04T15:24:54Z

src/ctapipe/calib/camera/calibrator.py

+                if self.image_calibration_type == "charge":
+                    dl1.image *= dl1_calib.relative_factor / dl1_calib.absolute_factor
+                elif self.image_calibration_type == "variance":
+                    dl1.image *= np.sqrt(
+                        dl1_calib.relative_factor / dl1_calib.absolute_factor
+                    )


You could check here if the extractor.__class__.__name__ is VarianceExtractor rather than passing the specific config parameter.

mexanick · 2024-11-05T10:28:54Z

src/ctapipe/image/extractor.py

@@ -1308,7 +1308,7 @@ def __call__(
        self, waveforms, tel_id, selected_gain_channel, broken_pixels
    ) -> DL1CameraContainer:
        container = DL1CameraContainer(
-            image=np.nanvar(waveforms, dtype="float32", axis=2),
+            image=np.nanvar(waveforms, dtype="float32", keepdims=False, axis=2),


A bit weird that you need to set keepdims=False explicitly, but ok.

i will check again if i make it work without that

Should work without it:

x = np.random.uniform(size=(2,1800, 1000)) In [1]: np.nanvar(x, axis=2).shape Out[2]: (2, 1800) In [2]: np.nanvar(x, axis=2, keepdims=True).shape Out[2]: (2, 1800, 1) In [3]: np.nanvar(x, axis=2, keepdims=False).shape Out[3]: (2, 1800)

mexanick · 2024-11-05T10:29:46Z

src/ctapipe/calib/camera/tests/test_calibrator.py

+    calibrator(example_event)
+    image = example_event.dl1.tel[tel_id].image
+    assert image is not None
+    assert image.shape == (


Add a test with the LST (two gains)

would it make sense to add a test for other extractor with the LST?

You could use a parametrized test here to test all telescope types. (@pytest.mark.parametrized). See for example the code in ctapipe/instrument/tests/test_telescope.py

mexanick · 2024-11-05T11:44:36Z

src/ctapipe/calib/camera/calibrator.py

@@ -292,13 +296,21 @@ def _calibrate_dl1(self, event, tel_id):
            and dl1_calib.absolute_factor is not None
        ):
            if selected_gain_channel is None:
-                dl1.image *= dl1_calib.relative_factor / dl1_calib.absolute_factor
+                if extractor.__class__.__name__ == "VarianceExtractor":
+                    dl1.image *= np.sqrt(


should be square, not square root

mexanick · 2024-11-05T11:45:00Z

src/ctapipe/calib/camera/calibrator.py

            else:
                corr = (
                    dl1_calib.relative_factor[selected_gain_channel, pixel_index]
                    / dl1_calib.absolute_factor[selected_gain_channel, pixel_index]
                )
-                dl1.image *= corr
+                if extractor.__class__.__name__ == "VarianceExtractor":
+                    dl1.image *= np.sqrt(corr)


same as l.300

maxnoe · 2024-11-06T08:30:57Z

docs/changes/2636.features.rst

@@ -0,0 +1 @@
+Makes changes to the CameraCalibrator in ctapipe.calib.camera.calibrator that allows it to correctly variance images generated with the VarianceExtractor


Please describe the actual change, not just say "there are changes". Also there is a verb missing here (calibrate?)

i added some explanation

mexanick · 2024-11-06T09:03:38Z

docs/changes/2636.features.rst

@@ -0,0 +1,3 @@
+Makes changes to the CameraCalibrator in ctapipe.calib.camera.calibrator that allows it to correctly calibrate variance images generated with the VarianceExtractor
+if the VarianceExtractor is used for the CameraCalibrator the element-wise square of the relative and absolute gain calibration factors are applied to the image


Please format according to ReST syntax. This will render as a single line.

I changed it

ctao-dpps-sonarqube · 2024-11-06T09:39:07Z

Analysis Details

2 Issues

0 Bugs
0 Vulnerabilities
2 Code Smells

Coverage and Duplications

97.30% Coverage (93.70% Estimated after merge)
0.00% Duplicated Code (0.70% Estimated after merge)

Project ID: cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

View in SonarQube

maxnoe · 2024-11-06T11:28:28Z

I think we can (maybe should) apply this sqrt operation before writing the coefficients. Therefore we do not need the different cases here in the code. Would that work?

I think this is the way to go, compute the appropriate calibration coefficients directly instead of adding special casing in the application.

Also: shouldn't it be square and not square root for the variance?

ctoennis · 2024-11-06T11:33:02Z

@maxnoe square is correct. If we apply the square operation when writing the calibration corrections then we will not need this PR

maxnoe · 2024-11-06T11:35:31Z

You still need the peak_time check, right?

ctoennis · 2024-11-06T11:37:10Z

Oh yes, that will be still needed

TjarkMiener · 2024-11-06T11:41:22Z

I discussed with @mexanick offline yesterday. We might not want to run the tool for computing camera calibration coefficients twice and store two separated files (one for comics; one for pointing calibration), just to square one column here. What do you think?

mexanick · 2024-11-06T13:09:04Z

at least for now, I would avoid recomputation of the calibrations just for variance. With the current code, we can always have this possibility to test, if needed.

kosack · 2024-11-07T10:32:52Z

src/ctapipe/calib/camera/calibrator.py

            else:
                corr = (
                    dl1_calib.relative_factor[selected_gain_channel, pixel_index]
                    / dl1_calib.absolute_factor[selected_gain_channel, pixel_index]
                )
-                dl1.image *= corr
+                if extractor.__class__.__name__ == "VarianceExtractor":


Please use this to check type, not a string comparison (which is slower and makes refactoring more difficult later)

Suggested change

if extractor.__class__.__name__ == "VarianceExtractor":

if isinstance(extractor, VarianceExtractor):

Makes sense, im changing it now

kosack · 2024-11-07T10:34:31Z

src/ctapipe/calib/camera/calibrator.py

@@ -292,13 +296,21 @@ def _calibrate_dl1(self, event, tel_id):
            and dl1_calib.absolute_factor is not None
        ):
            if selected_gain_channel is None:
-                dl1.image *= dl1_calib.relative_factor / dl1_calib.absolute_factor
+                if extractor.__class__.__name__ == "VarianceExtractor":


Suggested change

if extractor.__class__.__name__ == "VarianceExtractor":

if isinstance(extractor, VarianceExtractor):

kosack · 2024-11-07T10:37:34Z

src/ctapipe/calib/camera/tests/test_calibrator.py

+    calibrator(example_event)
+    image = example_event.dl1.tel[tel_id].image
+    assert image is not None
+    assert image.shape == (


You could use a parametrized test here to test all telescope types. (@pytest.mark.parametrized). See for example the code in ctapipe/instrument/tests/test_telescope.py

ctoennis · 2024-11-07T11:20:24Z

i will change the test to check all cameras

ctoennis · 2024-11-07T12:43:33Z

@kosack I used the camera_geometry fixture instead to check for all cameras. In the end I just need the geometry, so that should work well.

ctoennis · 2024-11-07T12:47:58Z

I think there is an issue with the recent change to the changelog check (#2637).

ctoennis · 2024-11-08T14:37:29Z

i switched it to a parametrized test for the different cameras now and caught the other test for name and switched it to isinstance

…data for test

maxnoe · 2024-11-08T16:35:28Z

src/ctapipe/calib/camera/calibrator.py

@@ -292,13 +296,21 @@ def _calibrate_dl1(self, event, tel_id):
            and dl1_calib.absolute_factor is not None
        ):
            if selected_gain_channel is None:


I think this section could use some refactoring to make it easier to follow.

I think the logic is the same if we simplify to:

if selected_gain_channel is None: calibration = dl1_calib.relative_factor / dl1_calib.absolute_factor else: calibration = ( dl1_calib.relative_factor[selected_gain_channel, pixel_index] / dl1_calib.absolute_factor[selected_gain_channel, pixel_index] ) if isinstance(extractor, VarianceExtractor): calibration = calibration**2 dl1.image *= calibration

Sounds good, I'll get to it when I arrive home

i put it in

Christoph Toennis added 4 commits November 1, 2024 16:54

Adding different method of applying calibration

7c328fa

Adding tests for variance calibration, fixing VarianceExtractor for u…

ceaf91b

…se with Calibrator

improving docustring

dd3a59f

Fixing error in VarianceExtractor test

4a5e93b

ctoennis added the calibration label Nov 4, 2024

ctoennis requested review from maxnoe and kosack November 4, 2024 10:18

ctoennis self-assigned this Nov 4, 2024

adding changelog

cb1929d

This comment has been minimized.

Sign in to view

mexanick requested changes Nov 4, 2024

View reviewed changes

TjarkMiener reviewed Nov 4, 2024

View reviewed changes

Christoph Toennis added 2 commits November 5, 2024 10:41

Addressing reviewer comments

e5d6fa6

Implementing change suggested by max

f24dfed

This comment has been minimized.

Sign in to view

mexanick reviewed Nov 5, 2024

View reviewed changes

Fixing variance calibration and adding test with LST camera

9861948

This comment has been minimized.

Sign in to view

mexanick requested review from TjarkMiener and mexanick November 6, 2024 08:25

mexanick previously approved these changes Nov 6, 2024

View reviewed changes

maxnoe reviewed Nov 6, 2024

View reviewed changes

Updating changelog

460857a

ctoennis dismissed mexanick’s stale review via 460857a November 6, 2024 08:37

Adding verb

13f83aa

mexanick reviewed Nov 6, 2024

View reviewed changes

This comment has been minimized.

Sign in to view

Formatting changelog

e1f4e13

This comment has been minimized.

Sign in to view

more formatting

29b9009

This comment has been minimized.

Sign in to view

Rename changelog entry and update formatting

9c7d4dc

kosack reviewed Nov 7, 2024

View reviewed changes

kosack requested changes Nov 7, 2024

View reviewed changes

Making test_dl1_variance_calib perform the test for all cameras

86b5db5

Christoph Toennis added 2 commits November 8, 2024 15:09

Switching to isinstance to test for variance extractor

554586f

parametrizing the test

6b3d336

ctoennis requested a review from kosack November 8, 2024 14:51

Testing all cameras through test_subarray fixture now, creating mock …

c7a8126

…data for test

maxnoe reviewed Nov 8, 2024

View reviewed changes

refactoring VarianceExtractor check

b9d7efb

		@@ -0,0 +1 @@
		Makes changes to the CameraCalibrator in ctapipe.calib.camera.calibrator that allows it to correctly variance images generated with the VarianceExtractor

		@@ -0,0 +1,3 @@
		Makes changes to the CameraCalibrator in ctapipe.calib.camera.calibrator that allows it to correctly calibrate variance images generated with the VarianceExtractor
		if the VarianceExtractor is used for the CameraCalibrator the element-wise square of the relative and absolute gain calibration factors are applied to the image

	if extractor.__class__.__name__ == "VarianceExtractor":
	if isinstance(extractor, VarianceExtractor):

Variance calibration #2636

Are you sure you want to change the base?

Variance calibration #2636

Conversation

ctoennis commented Nov 4, 2024

This comment has been minimized.

This comment has been minimized.

TjarkMiener commented Nov 4, 2024

ctoennis commented Nov 4, 2024

ctoennis commented Nov 4, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

This comment has been minimized.

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

This comment has been minimized.

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

This comment has been minimized.

This comment has been minimized.

This comment has been minimized.

ctao-dpps-sonarqube bot commented Nov 6, 2024

Analysis Details

2 Issues

Coverage and Duplications

maxnoe commented Nov 6, 2024

ctoennis commented Nov 6, 2024

maxnoe commented Nov 6, 2024

ctoennis commented Nov 6, 2024

TjarkMiener commented Nov 6, 2024

mexanick commented Nov 6, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ctoennis commented Nov 7, 2024

ctoennis commented Nov 7, 2024

ctoennis commented Nov 7, 2024

ctoennis commented Nov 8, 2024

maxnoe Nov 8, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

maxnoe Nov 8, 2024 •

edited

Loading