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

Variance calibration #2636

Open
wants to merge 18 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 7 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions docs/changes/2636.features.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
Makes changes to the CameraCalibrator in ctapipe.calib.camera.calibrator that allows it to correctly variance images generated with the VarianceExtractor
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i added some explanation

18 changes: 15 additions & 3 deletions src/ctapipe/calib/camera/calibrator.py
Original file line number Diff line number Diff line change
Expand Up @@ -218,7 +218,7 @@
calibration_monitoring_id=r1.calibration_monitoring_id,
)

def _calibrate_dl1(self, event, tel_id):

Check failure on line 221 in src/ctapipe/calib/camera/calibrator.py

View check run for this annotation

CTAO-DPPS-SonarQube / ctapipe Sonarqube Results

src/ctapipe/calib/camera/calibrator.py#L221

Refactor this function to reduce its Cognitive Complexity from 30 to the 15 allowed.
waveforms = event.dl0.tel[tel_id].waveform
if self._check_dl0_empty(waveforms):
return
Expand Down Expand Up @@ -283,7 +283,11 @@
)

# correct non-integer remainder of the shift if given
if self.apply_peak_time_shift.tel[tel_id] and remaining_shift is not None:
if (
self.apply_peak_time_shift.tel[tel_id]
and remaining_shift is not None
and dl1.peak_time is not None
):
dl1.peak_time -= remaining_shift

# Calibrate extracted charge
Expand All @@ -292,13 +296,21 @@
and dl1_calib.absolute_factor is not None
):
if selected_gain_channel is None:
Copy link
Member

@maxnoe maxnoe Nov 8, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i put it in

dl1.image *= dl1_calib.relative_factor / dl1_calib.absolute_factor
if extractor.__class__.__name__ == "VarianceExtractor":
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

dl1.image *= np.sqrt(
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should be square, not square root

dl1_calib.relative_factor / dl1_calib.absolute_factor
)
else:
dl1.image *= dl1_calib.relative_factor / dl1_calib.absolute_factor
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":
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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):

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Makes sense, im changing it now

dl1.image *= np.sqrt(corr)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same as l.300

else:
dl1.image *= corr

# handle invalid pixels
if self.invalid_pixel_handler is not None:
Expand Down
18 changes: 18 additions & 0 deletions src/ctapipe/calib/camera/tests/test_calibrator.py
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
GlobalPeakWindowSum,
LocalPeakWindowSum,
NeighborPeakWindowSum,
VarianceExtractor,
)
from ctapipe.image.reducer import NullDataVolumeReducer, TailCutsDataVolumeReducer

Expand Down Expand Up @@ -130,6 +131,23 @@ def test_check_dl0_empty(example_event, example_subarray):
assert (event.dl1.tel[tel_id].image == 2).all()


def test_dl1_variance_calib(example_event, example_subarray):
# test the calibration of variance images
tel_id = list(example_event.r0.tel)[0]
calibrator = CameraCalibrator(
subarray=example_subarray,
image_extractor=VarianceExtractor(subarray=example_subarray),
apply_waveform_time_shift=False,
)
calibrator(example_event)
image = example_event.dl1.tel[tel_id].image
assert image is not None
assert image.shape == (
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add a test with the LST (two gains)

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ok

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

1,
1764,
)


def test_dl1_charge_calib(example_subarray):
# copy because we mutate the camera, should not affect other tests
subarray = deepcopy(example_subarray)
Expand Down
2 changes: 1 addition & 1 deletion src/ctapipe/image/extractor.py
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i will check again if i make it work without that

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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)

)
container.meta["ExtractionMethod"] = str(VarianceType.WAVEFORM)
return container
Expand Down