From e9296fe329029f5c491dd2c963bdda8ffd55409e Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Thu, 31 Oct 2024 13:02:09 +0100 Subject: [PATCH 01/33] adding PSF model --- src/ctapipe/image/psf_model.py | 147 +++++++++++++++++++++++++++++++++ 1 file changed, 147 insertions(+) create mode 100644 src/ctapipe/image/psf_model.py diff --git a/src/ctapipe/image/psf_model.py b/src/ctapipe/image/psf_model.py new file mode 100644 index 00000000000..9f6f97ffeb2 --- /dev/null +++ b/src/ctapipe/image/psf_model.py @@ -0,0 +1,147 @@ +""" +Models for the Point Spread Functions of the different telescopes +""" + +__all__ = ["PSFModel", "ComaModel"] + +from abc import abstractmethod + +import numpy as np +from scipy.stats import laplace, laplace_asymmetric + +from ctapipe.core.component import non_abstract_children +from ctapipe.core.plugins import detect_and_import_plugins + + +class PSFModel: + def __init__(self, **kwargs): + """ + Base component to describe image distortion due to the optics of the different cameras. + """ + + @classmethod + def from_name(cls, name, **kwargs): + """ + Obtain an instance of a subclass via its name + + Parameters + ---------- + name : str + Name of the subclass to obtain + + Returns + ------- + Instance + Instance of subclass to this class + """ + requested_subclass = cls.non_abstract_subclasses()[name] + return requested_subclass(**kwargs) + + @classmethod + def non_abstract_subclasses(cls): + """ + Get a dict of all non-abstract subclasses of this class. + + This method is using the entry-point plugin system + to also check for registered plugin implementations. + + Returns + ------- + subclasses : dict[str, type] + A dict mapping the name to the class of all found, + non-abstract subclasses of this class. + """ + if hasattr(cls, "plugin_entry_point"): + detect_and_import_plugins(cls.plugin_entry_point) + + subclasses = {base.__name__: base for base in non_abstract_children(cls)} + return subclasses + + @abstractmethod + def pdf(self, *args): + pass + + @abstractmethod + def update_location(self, *args): + pass + + @abstractmethod + def update_model_parameters(self, *args): + pass + + +class ComaModel(PSFModel): + """ + PSF model, describing pure coma aberrations PSF effect + """ + + def __init__( + self, + asymmetry_params=[0.49244797, 9.23573115, 0.15216096], + radial_scale_params=[0.01409259, 0.02947208, 0.06000271, -0.02969355], + az_scale_params=[0.24271557, 7.5511501, 0.02037972], + ): + """ + PSF model, describing pure coma aberrations PSF effect + + param list asymmetry_params Parameters describing the dependency of the asymmetry of the psf on the distance to the center of the camera + param list radial_scale_params Parameters describing the dependency of the radial scale on the distance to the center of the camera + param list radial_scale_params Parameters describing the dependency of the azimuthal scale scale on the distance to the center of the camera + """ + + self.asymmetry_params = asymmetry_params + self.radial_scale_params = radial_scale_params + self.az_scale_params = az_scale_params + + self.update_location(0.0, 0.0) + + def k_func(self, x): + return ( + 1 + - self.asymmetry_params[0] * np.tanh(self.asymmetry_params[1] * x) + - self.asymmetry_params[2] * x + ) + + def sr_func(self, x): + return ( + self.radial_scale_params[0] + - self.radial_scale_params[1] * x + + self.radial_scale_params[2] * x**2 + - self.radial_scale_params[3] * x**3 + ) + + def sf_func(self, x): + return self.az_scale_params[0] * np.exp( + -self.az_scale_params[1] * x + ) + self.az_scale_params[2] / (self.az_scale_params[2] + x) + + def pdf(self, r, f): + return laplace_asymmetric.pdf(r, *self.radial_pdf_params) * laplace.pdf( + f, *self.azimuthal_pdf_params + ) + + def update_model_parameters(self, model_params): + if not ( + len(model_params["asymmetry_params"]) == 3 + and len(model_params["radial_scale_params"]) == 4 + and len(model_params["az_scale_params"]) == 3 + ): + raise ValueError( + "asymmetry_params and az_scale_params needs to have length 3 and radial_scale_params length 4" + ) + + self.asymmetry_params = model_params["asymmetry_params"] + self.radial_scale_params = model_params["radial_scale_params"] + self.az_scale_params = model_params["az_scale_params"] + k = self.k_func(self.radial_pdf_params[1]) + sr = self.sr_func(self.radial_pdf_params[1]) + sf = self.sf_func(self.radial_pdf_params[1]) + self.radial_pdf_params = (k, self.radial_pdf_params[1], sr) + self.azimuthal_pdf_params = (self.azimuthal_pdf_params[0], sf) + + def update_location(self, r, f): + k = self.k_func(r) + sr = self.sr_func(r) + sf = self.sf_func(r) + self.radial_pdf_params = (k, r, sr) + self.azimuthal_pdf_params = (f, sf) From dc46bfb90be9225178d116a5419f7dd16d72475b Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Fri, 1 Nov 2024 14:41:02 +0100 Subject: [PATCH 02/33] adding tests and documentation --- src/ctapipe/image/psf_model.py | 47 +++++++++++++++++++++-- src/ctapipe/image/tests/test_psf_model.py | 27 +++++++++++++ 2 files changed, 70 insertions(+), 4 deletions(-) create mode 100644 src/ctapipe/image/tests/test_psf_model.py diff --git a/src/ctapipe/image/psf_model.py b/src/ctapipe/image/psf_model.py index 9f6f97ffeb2..cba6fc2098a 100644 --- a/src/ctapipe/image/psf_model.py +++ b/src/ctapipe/image/psf_model.py @@ -82,11 +82,16 @@ def __init__( az_scale_params=[0.24271557, 7.5511501, 0.02037972], ): """ - PSF model, describing pure coma aberrations PSF effect + PSF model, describing purely the effect of coma aberration on the PSF - param list asymmetry_params Parameters describing the dependency of the asymmetry of the psf on the distance to the center of the camera - param list radial_scale_params Parameters describing the dependency of the radial scale on the distance to the center of the camera - param list radial_scale_params Parameters describing the dependency of the azimuthal scale scale on the distance to the center of the camera + Parameters + ---------- + asymmetry_params: list + Parameters describing the dependency of the asymmetry of the psf on the distance to the center of the camera + radial_scale_params : list + Parameters describing the dependency of the radial scale on the distance to the center of the camera + radial_scale_params : list + Parameters describing the dependency of the azimuthal scale scale on the distance to the center of the camera """ self.asymmetry_params = asymmetry_params @@ -116,11 +121,36 @@ def sf_func(self, x): ) + self.az_scale_params[2] / (self.az_scale_params[2] + x) def pdf(self, r, f): + """ + Calculates the value of the psf at a given location + + Parameters + ---------- + r : float + distance to the center of the camera in meters + f : float + polar angle in radians + + Returns + ---------- + psf : float + value of the PSF at the specified location + """ return laplace_asymmetric.pdf(r, *self.radial_pdf_params) * laplace.pdf( f, *self.azimuthal_pdf_params ) def update_model_parameters(self, model_params): + """ + Updates the model parameters for the psf + + Parameters + ---------- + model_params : dict + dictionary with the model parameters + needs to have the keys `asymmetry_params`, `radial_scale_params` and `az_scale_params` + The values need to be lists of length 3, 4 and 3 respectively + """ if not ( len(model_params["asymmetry_params"]) == 3 and len(model_params["radial_scale_params"]) == 4 @@ -140,6 +170,15 @@ def update_model_parameters(self, model_params): self.azimuthal_pdf_params = (self.azimuthal_pdf_params[0], sf) def update_location(self, r, f): + """ + Updates the location on the camera focal plane from where the psf is calculated + Parameters + ---------- + r : float + distance to the center of the camera in meters + f : float + polar angle in radians + """ k = self.k_func(r) sr = self.sr_func(r) sf = self.sf_func(r) diff --git a/src/ctapipe/image/tests/test_psf_model.py b/src/ctapipe/image/tests/test_psf_model.py new file mode 100644 index 00000000000..21ecd22fe34 --- /dev/null +++ b/src/ctapipe/image/tests/test_psf_model.py @@ -0,0 +1,27 @@ +""" +This module contains the ctapipe.image.psf_model unit tests +""" +import numpy as np +import pytest + +from ctapipe.image.psf_model import PSFModel + + +def test_psf(): + psf = PSFModel.from_name("ComaModel") + wrong_parameters = { + "asymmetry_params": [0.49244797, 9.23573115, 0.15216096], + "radial_scale_params": [0.01409259, 0.02947208, 0.06000271, -0.02969355], + "az_scale_params": [0.24271557, 7.5511501], + } + with pytest.raises( + ValueError, + match="asymmetry_params and az_scale_params needs to have length 3 and radial_scale_params length 4", + ): + psf.update_model_parameters(wrong_parameters) + + +def test_asymptotic_behavior(): + psf_coma = PSFModel.from_name("ComaModel") + psf_coma.update_location(1.0, 0.0) + assert np.isclose(psf_coma.psf(10.0, 0.0), 0.0) From 3a005018994abb041e3f2f9f6d358d2d80751a60 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Sun, 10 Nov 2024 13:29:02 +0100 Subject: [PATCH 03/33] adding changelog, fixing test --- docs/changes/2641.feature.rst | 4 ++++ src/ctapipe/image/psf_model.py | 4 ++-- src/ctapipe/image/tests/test_psf_model.py | 2 +- 3 files changed, 7 insertions(+), 3 deletions(-) create mode 100644 docs/changes/2641.feature.rst diff --git a/docs/changes/2641.feature.rst b/docs/changes/2641.feature.rst new file mode 100644 index 00000000000..18e979ae82f --- /dev/null +++ b/docs/changes/2641.feature.rst @@ -0,0 +1,4 @@ +Adds ``ctapipe.image.psf_model`` with the parent class for psf models called ``PSFModel`` and a psf model based on pure coma abbaration called ``ComaModel`` +- All psf models have a function ``PSFModel.update_location`` that updates the location in the focal plane where the PSF is evaluated from +- The function ``PSFModel.update_model_parameters`` allows to update the parameters of the PSF model +- The function ``PSFModel.pdf`` gives the value of the PSF in a given location diff --git a/src/ctapipe/image/psf_model.py b/src/ctapipe/image/psf_model.py index cba6fc2098a..92b81744e83 100644 --- a/src/ctapipe/image/psf_model.py +++ b/src/ctapipe/image/psf_model.py @@ -90,7 +90,7 @@ def __init__( Parameters describing the dependency of the asymmetry of the psf on the distance to the center of the camera radial_scale_params : list Parameters describing the dependency of the radial scale on the distance to the center of the camera - radial_scale_params : list + az_scale_params : list Parameters describing the dependency of the azimuthal scale scale on the distance to the center of the camera """ @@ -126,7 +126,7 @@ def pdf(self, r, f): Parameters ---------- - r : float + r : float distance to the center of the camera in meters f : float polar angle in radians diff --git a/src/ctapipe/image/tests/test_psf_model.py b/src/ctapipe/image/tests/test_psf_model.py index 21ecd22fe34..d6c953ea483 100644 --- a/src/ctapipe/image/tests/test_psf_model.py +++ b/src/ctapipe/image/tests/test_psf_model.py @@ -24,4 +24,4 @@ def test_psf(): def test_asymptotic_behavior(): psf_coma = PSFModel.from_name("ComaModel") psf_coma.update_location(1.0, 0.0) - assert np.isclose(psf_coma.psf(10.0, 0.0), 0.0) + assert np.isclose(psf_coma.pdf(10.0, 0.0), 0.0) From c260bfb54700162604788bce7076e2279bb8c089 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Sun, 10 Nov 2024 13:36:28 +0100 Subject: [PATCH 04/33] renaming changelog --- docs/changes/{2641.feature.rst => 2643.feature.rst} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename docs/changes/{2641.feature.rst => 2643.feature.rst} (100%) diff --git a/docs/changes/2641.feature.rst b/docs/changes/2643.feature.rst similarity index 100% rename from docs/changes/2641.feature.rst rename to docs/changes/2643.feature.rst From f4c95a37372333a720178cece23307ad9a4e5148 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Mon, 11 Nov 2024 15:29:42 +0100 Subject: [PATCH 05/33] Moving the psf models to ctapipe.instrument.optics --- docs/changes/2643.feature.rst | 2 +- src/ctapipe/image/psf_model.py | 186 ------------------ src/ctapipe/instrument/optics.py | 180 +++++++++++++++++ .../tests/test_psf_model.py | 2 +- 4 files changed, 182 insertions(+), 188 deletions(-) delete mode 100644 src/ctapipe/image/psf_model.py rename src/ctapipe/{image => instrument}/tests/test_psf_model.py (94%) diff --git a/docs/changes/2643.feature.rst b/docs/changes/2643.feature.rst index 18e979ae82f..336b0b65f30 100644 --- a/docs/changes/2643.feature.rst +++ b/docs/changes/2643.feature.rst @@ -1,4 +1,4 @@ -Adds ``ctapipe.image.psf_model`` with the parent class for psf models called ``PSFModel`` and a psf model based on pure coma abbaration called ``ComaModel`` +Adds psf models to ``ctapipe.instrument.optics`` with the parent class ``PSFModel`` and a psf model based on pure coma abbaration called ``ComaModel`` - All psf models have a function ``PSFModel.update_location`` that updates the location in the focal plane where the PSF is evaluated from - The function ``PSFModel.update_model_parameters`` allows to update the parameters of the PSF model - The function ``PSFModel.pdf`` gives the value of the PSF in a given location diff --git a/src/ctapipe/image/psf_model.py b/src/ctapipe/image/psf_model.py deleted file mode 100644 index 92b81744e83..00000000000 --- a/src/ctapipe/image/psf_model.py +++ /dev/null @@ -1,186 +0,0 @@ -""" -Models for the Point Spread Functions of the different telescopes -""" - -__all__ = ["PSFModel", "ComaModel"] - -from abc import abstractmethod - -import numpy as np -from scipy.stats import laplace, laplace_asymmetric - -from ctapipe.core.component import non_abstract_children -from ctapipe.core.plugins import detect_and_import_plugins - - -class PSFModel: - def __init__(self, **kwargs): - """ - Base component to describe image distortion due to the optics of the different cameras. - """ - - @classmethod - def from_name(cls, name, **kwargs): - """ - Obtain an instance of a subclass via its name - - Parameters - ---------- - name : str - Name of the subclass to obtain - - Returns - ------- - Instance - Instance of subclass to this class - """ - requested_subclass = cls.non_abstract_subclasses()[name] - return requested_subclass(**kwargs) - - @classmethod - def non_abstract_subclasses(cls): - """ - Get a dict of all non-abstract subclasses of this class. - - This method is using the entry-point plugin system - to also check for registered plugin implementations. - - Returns - ------- - subclasses : dict[str, type] - A dict mapping the name to the class of all found, - non-abstract subclasses of this class. - """ - if hasattr(cls, "plugin_entry_point"): - detect_and_import_plugins(cls.plugin_entry_point) - - subclasses = {base.__name__: base for base in non_abstract_children(cls)} - return subclasses - - @abstractmethod - def pdf(self, *args): - pass - - @abstractmethod - def update_location(self, *args): - pass - - @abstractmethod - def update_model_parameters(self, *args): - pass - - -class ComaModel(PSFModel): - """ - PSF model, describing pure coma aberrations PSF effect - """ - - def __init__( - self, - asymmetry_params=[0.49244797, 9.23573115, 0.15216096], - radial_scale_params=[0.01409259, 0.02947208, 0.06000271, -0.02969355], - az_scale_params=[0.24271557, 7.5511501, 0.02037972], - ): - """ - PSF model, describing purely the effect of coma aberration on the PSF - - Parameters - ---------- - asymmetry_params: list - Parameters describing the dependency of the asymmetry of the psf on the distance to the center of the camera - radial_scale_params : list - Parameters describing the dependency of the radial scale on the distance to the center of the camera - az_scale_params : list - Parameters describing the dependency of the azimuthal scale scale on the distance to the center of the camera - """ - - self.asymmetry_params = asymmetry_params - self.radial_scale_params = radial_scale_params - self.az_scale_params = az_scale_params - - self.update_location(0.0, 0.0) - - def k_func(self, x): - return ( - 1 - - self.asymmetry_params[0] * np.tanh(self.asymmetry_params[1] * x) - - self.asymmetry_params[2] * x - ) - - def sr_func(self, x): - return ( - self.radial_scale_params[0] - - self.radial_scale_params[1] * x - + self.radial_scale_params[2] * x**2 - - self.radial_scale_params[3] * x**3 - ) - - def sf_func(self, x): - return self.az_scale_params[0] * np.exp( - -self.az_scale_params[1] * x - ) + self.az_scale_params[2] / (self.az_scale_params[2] + x) - - def pdf(self, r, f): - """ - Calculates the value of the psf at a given location - - Parameters - ---------- - r : float - distance to the center of the camera in meters - f : float - polar angle in radians - - Returns - ---------- - psf : float - value of the PSF at the specified location - """ - return laplace_asymmetric.pdf(r, *self.radial_pdf_params) * laplace.pdf( - f, *self.azimuthal_pdf_params - ) - - def update_model_parameters(self, model_params): - """ - Updates the model parameters for the psf - - Parameters - ---------- - model_params : dict - dictionary with the model parameters - needs to have the keys `asymmetry_params`, `radial_scale_params` and `az_scale_params` - The values need to be lists of length 3, 4 and 3 respectively - """ - if not ( - len(model_params["asymmetry_params"]) == 3 - and len(model_params["radial_scale_params"]) == 4 - and len(model_params["az_scale_params"]) == 3 - ): - raise ValueError( - "asymmetry_params and az_scale_params needs to have length 3 and radial_scale_params length 4" - ) - - self.asymmetry_params = model_params["asymmetry_params"] - self.radial_scale_params = model_params["radial_scale_params"] - self.az_scale_params = model_params["az_scale_params"] - k = self.k_func(self.radial_pdf_params[1]) - sr = self.sr_func(self.radial_pdf_params[1]) - sf = self.sf_func(self.radial_pdf_params[1]) - self.radial_pdf_params = (k, self.radial_pdf_params[1], sr) - self.azimuthal_pdf_params = (self.azimuthal_pdf_params[0], sf) - - def update_location(self, r, f): - """ - Updates the location on the camera focal plane from where the psf is calculated - Parameters - ---------- - r : float - distance to the center of the camera in meters - f : float - polar angle in radians - """ - k = self.k_func(r) - sr = self.sr_func(r) - sf = self.sf_func(r) - self.radial_pdf_params = (k, r, sr) - self.azimuthal_pdf_params = (f, sf) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index ddec040fc5f..2bb799050ed 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -3,11 +3,16 @@ """ import logging +from abc import abstractmethod from enum import Enum, auto, unique import astropy.units as u import numpy as np from astropy.table import QTable +from scipy.stats import laplace, laplace_asymmetric + +from ctapipe.core.component import non_abstract_children +from ctapipe.core.plugins import detect_and_import_plugins from ..compat import StrEnum from ..utils import get_table_dataset @@ -18,6 +23,8 @@ __all__ = [ "OpticsDescription", "FocalLengthKind", + "PSFModel", + "ComaModel", ] @@ -263,3 +270,176 @@ def __repr__(self): def __str__(self): return self.name + + +class PSFModel: + def __init__(self, **kwargs): + """ + Base component to describe image distortion due to the optics of the different cameras. + """ + + @classmethod + def from_name(cls, name, **kwargs): + """ + Obtain an instance of a subclass via its name + + Parameters + ---------- + name : str + Name of the subclass to obtain + + Returns + ------- + Instance + Instance of subclass to this class + """ + requested_subclass = cls.non_abstract_subclasses()[name] + return requested_subclass(**kwargs) + + @classmethod + def non_abstract_subclasses(cls): + """ + Get a dict of all non-abstract subclasses of this class. + + This method is using the entry-point plugin system + to also check for registered plugin implementations. + + Returns + ------- + subclasses : dict[str, type] + A dict mapping the name to the class of all found, + non-abstract subclasses of this class. + """ + if hasattr(cls, "plugin_entry_point"): + detect_and_import_plugins(cls.plugin_entry_point) + + subclasses = {base.__name__: base for base in non_abstract_children(cls)} + return subclasses + + @abstractmethod + def pdf(self, *args): + pass + + @abstractmethod + def update_location(self, *args): + pass + + @abstractmethod + def update_model_parameters(self, *args): + pass + + +class ComaModel(PSFModel): + """ + PSF model, describing pure coma aberrations PSF effect + """ + + def __init__( + self, + asymmetry_params=[0.49244797, 9.23573115, 0.15216096], + radial_scale_params=[0.01409259, 0.02947208, 0.06000271, -0.02969355], + az_scale_params=[0.24271557, 7.5511501, 0.02037972], + ): + """ + PSF model, describing purely the effect of coma aberration on the PSF + + Parameters + ---------- + asymmetry_params: list + Parameters describing the dependency of the asymmetry of the psf on the distance to the center of the camera + radial_scale_params : list + Parameters describing the dependency of the radial scale on the distance to the center of the camera + az_scale_params : list + Parameters describing the dependency of the azimuthal scale scale on the distance to the center of the camera + """ + + self.asymmetry_params = asymmetry_params + self.radial_scale_params = radial_scale_params + self.az_scale_params = az_scale_params + + self.update_location(0.0, 0.0) + + def k_func(self, x): + return ( + 1 + - self.asymmetry_params[0] * np.tanh(self.asymmetry_params[1] * x) + - self.asymmetry_params[2] * x + ) + + def sr_func(self, x): + return ( + self.radial_scale_params[0] + - self.radial_scale_params[1] * x + + self.radial_scale_params[2] * x**2 + - self.radial_scale_params[3] * x**3 + ) + + def sf_func(self, x): + return self.az_scale_params[0] * np.exp( + -self.az_scale_params[1] * x + ) + self.az_scale_params[2] / (self.az_scale_params[2] + x) + + def pdf(self, r, f): + """ + Calculates the value of the psf at a given location + + Parameters + ---------- + r : float + distance to the center of the camera in meters + f : float + polar angle in radians + + Returns + ---------- + psf : float + value of the PSF at the specified location + """ + return laplace_asymmetric.pdf(r, *self.radial_pdf_params) * laplace.pdf( + f, *self.azimuthal_pdf_params + ) + + def update_model_parameters(self, model_params): + """ + Updates the model parameters for the psf + + Parameters + ---------- + model_params : dict + dictionary with the model parameters + needs to have the keys `asymmetry_params`, `radial_scale_params` and `az_scale_params` + The values need to be lists of length 3, 4 and 3 respectively + """ + if not ( + len(model_params["asymmetry_params"]) == 3 + and len(model_params["radial_scale_params"]) == 4 + and len(model_params["az_scale_params"]) == 3 + ): + raise ValueError( + "asymmetry_params and az_scale_params needs to have length 3 and radial_scale_params length 4" + ) + + self.asymmetry_params = model_params["asymmetry_params"] + self.radial_scale_params = model_params["radial_scale_params"] + self.az_scale_params = model_params["az_scale_params"] + k = self.k_func(self.radial_pdf_params[1]) + sr = self.sr_func(self.radial_pdf_params[1]) + sf = self.sf_func(self.radial_pdf_params[1]) + self.radial_pdf_params = (k, self.radial_pdf_params[1], sr) + self.azimuthal_pdf_params = (self.azimuthal_pdf_params[0], sf) + + def update_location(self, r, f): + """ + Updates the location on the camera focal plane from where the psf is calculated + Parameters + ---------- + r : float + distance to the center of the camera in meters + f : float + polar angle in radians + """ + k = self.k_func(r) + sr = self.sr_func(r) + sf = self.sf_func(r) + self.radial_pdf_params = (k, r, sr) + self.azimuthal_pdf_params = (f, sf) diff --git a/src/ctapipe/image/tests/test_psf_model.py b/src/ctapipe/instrument/tests/test_psf_model.py similarity index 94% rename from src/ctapipe/image/tests/test_psf_model.py rename to src/ctapipe/instrument/tests/test_psf_model.py index d6c953ea483..5baae1182d3 100644 --- a/src/ctapipe/image/tests/test_psf_model.py +++ b/src/ctapipe/instrument/tests/test_psf_model.py @@ -4,7 +4,7 @@ import numpy as np import pytest -from ctapipe.image.psf_model import PSFModel +from ctapipe.instrument.optics import PSFModel def test_psf(): From 31a56bb022b8e743b7427190c834ef97c74496f8 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Mon, 11 Nov 2024 16:04:36 +0100 Subject: [PATCH 06/33] Updating docustrings --- src/ctapipe/instrument/optics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index 2bb799050ed..f9843fd4cec 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -407,7 +407,7 @@ def update_model_parameters(self, model_params): ---------- model_params : dict dictionary with the model parameters - needs to have the keys `asymmetry_params`, `radial_scale_params` and `az_scale_params` + needs to have the keys asymmetry_params, radial_scale_params and az_scale_params The values need to be lists of length 3, 4 and 3 respectively """ if not ( From 317609b167cf272b596cd5fd197b87829f5040e8 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Mon, 11 Nov 2024 16:33:25 +0100 Subject: [PATCH 07/33] upate docustring --- src/ctapipe/instrument/optics.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index f9843fd4cec..ccaafcfc264 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -431,6 +431,7 @@ def update_model_parameters(self, model_params): def update_location(self, r, f): """ Updates the location on the camera focal plane from where the psf is calculated + Parameters ---------- r : float @@ -438,6 +439,7 @@ def update_location(self, r, f): f : float polar angle in radians """ + k = self.k_func(r) sr = self.sr_func(r) sf = self.sf_func(r) From 163b0a9f5200a831f6fc3e16426d61f13e793658 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Mon, 11 Nov 2024 16:51:02 +0100 Subject: [PATCH 08/33] Updating and renaming from_name method --- src/ctapipe/instrument/optics.py | 26 +++---------------- .../instrument/tests/test_psf_model.py | 4 +-- 2 files changed, 6 insertions(+), 24 deletions(-) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index ccaafcfc264..f2a8a7f0df4 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -12,7 +12,6 @@ from scipy.stats import laplace, laplace_asymmetric from ctapipe.core.component import non_abstract_children -from ctapipe.core.plugins import detect_and_import_plugins from ..compat import StrEnum from ..utils import get_table_dataset @@ -279,7 +278,7 @@ def __init__(self, **kwargs): """ @classmethod - def from_name(cls, name, **kwargs): + def subclass_from_name(cls, name, **kwargs): """ Obtain an instance of a subclass via its name @@ -293,28 +292,11 @@ def from_name(cls, name, **kwargs): Instance Instance of subclass to this class """ - requested_subclass = cls.non_abstract_subclasses()[name] - return requested_subclass(**kwargs) - - @classmethod - def non_abstract_subclasses(cls): - """ - Get a dict of all non-abstract subclasses of this class. - - This method is using the entry-point plugin system - to also check for registered plugin implementations. - - Returns - ------- - subclasses : dict[str, type] - A dict mapping the name to the class of all found, - non-abstract subclasses of this class. - """ - if hasattr(cls, "plugin_entry_point"): - detect_and_import_plugins(cls.plugin_entry_point) subclasses = {base.__name__: base for base in non_abstract_children(cls)} - return subclasses + + requested_subclass = subclasses[name] + return requested_subclass(**kwargs) @abstractmethod def pdf(self, *args): diff --git a/src/ctapipe/instrument/tests/test_psf_model.py b/src/ctapipe/instrument/tests/test_psf_model.py index 5baae1182d3..5a0f1a21bdc 100644 --- a/src/ctapipe/instrument/tests/test_psf_model.py +++ b/src/ctapipe/instrument/tests/test_psf_model.py @@ -8,7 +8,7 @@ def test_psf(): - psf = PSFModel.from_name("ComaModel") + psf = PSFModel.subclass_from_name("ComaModel") wrong_parameters = { "asymmetry_params": [0.49244797, 9.23573115, 0.15216096], "radial_scale_params": [0.01409259, 0.02947208, 0.06000271, -0.02969355], @@ -22,6 +22,6 @@ def test_psf(): def test_asymptotic_behavior(): - psf_coma = PSFModel.from_name("ComaModel") + psf_coma = PSFModel.subclass_from_name("ComaModel") psf_coma.update_location(1.0, 0.0) assert np.isclose(psf_coma.pdf(10.0, 0.0), 0.0) From 5b807fec268e9f981b78f88f7fdb4bb6a6cd147d Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Wed, 13 Nov 2024 20:49:35 +0100 Subject: [PATCH 09/33] updating doumentation with comamodel math --- docs/references.bib | 11 +++++++++ src/ctapipe/instrument/optics.py | 39 +++++++++++++++++++++++--------- 2 files changed, 39 insertions(+), 11 deletions(-) diff --git a/docs/references.bib b/docs/references.bib index aa40c827f29..65146670f51 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -136,3 +136,14 @@ @article{bright-star-catalog year={1991}, url={https://heasarc.gsfc.nasa.gov/W3Browse/star-catalog/bsc5p.html}, } + +@article{startracker, + author = {Abe, K. and others} + title = {Star tracking for pointing determination of Imaging Atmospheric Cherenkov Telescopes - Application to the Large-Sized Telescope of the Cherenkov Telescope Array}, + doi= "10.1051/0004-6361/202347128", + url= "https://doi.org/10.1051/0004-6361/202347128", + journal = {A&A}, + year = 2023, + volume = 679, + pages = "A90", +} diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index f2a8a7f0df4..b6cad715979 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -324,20 +324,39 @@ def __init__( ): """ PSF model, describing purely the effect of coma aberration on the PSF + for reference see :cite:p:`startracker` + uses an asymmetric laplacian for the radial part: + .. math:: f_{R}(r, K) = \begin{cases}\frac{1}{S_{R}(K+K^{-1})}e^{-K\frac{r-L}{S_{R}}}, r\ge L\\ \frac{1}{S_{R}(K+K^{-1})}e^{\frac{r-L}{KS_{R}}}, r < L\end{cases} + and a symmetric laplacian in azimuthal direction: + .. math:: f_{\Phi}(\phi) = \frac{1}{2S_\phi}e^{-|\frac{\phi-\phi_0}{S_\phi}|} Parameters ---------- asymmetry_params: list Parameters describing the dependency of the asymmetry of the psf on the distance to the center of the camera + Used to calculate a pdf asymmetry parameter K of the asymmetric radial laplacian of the psf as a function of the distance r to the optical axis: + .. math:: K(r) = 1 - asym_0 \tanh(asym_1 r) - asym_2 r radial_scale_params : list Parameters describing the dependency of the radial scale on the distance to the center of the camera + Used to calculate width Sr of the asymmetric radial laplacian in the PSF as a of function the distance r to the optical axis: + .. math:: S_{R}(r) & = b_1 - b_2\,r + b_3\,r^2 + b_4\,r^3 az_scale_params : list - Parameters describing the dependency of the azimuthal scale scale on the distance to the center of the camera + Parameters describing the dependency of the azimuthal scale on the distance to the center of the camera + Used to calculate width Sf of the azimuthal laplacian in the PSF as a function of the angle phi: + .. math:: S_{\phi}(r) & = a_1\,\exp{(-a_2\,r)}+\frac{a_3}{a_3+r} """ - - self.asymmetry_params = asymmetry_params - self.radial_scale_params = radial_scale_params - self.az_scale_params = az_scale_params + if ( + len(asymmetry_params) == 3 + and len(radial_scale_params) == 4 + and len(az_scale_params) == 3 + ): + self.asymmetry_params = asymmetry_params + self.radial_scale_params = radial_scale_params + self.az_scale_params = az_scale_params + else: + raise ValueError( + "asymmetry_params and az_scale_params needs to have length 3 and radial_scale_params length 4" + ) self.update_location(0.0, 0.0) @@ -349,12 +368,7 @@ def k_func(self, x): ) def sr_func(self, x): - return ( - self.radial_scale_params[0] - - self.radial_scale_params[1] * x - + self.radial_scale_params[2] * x**2 - - self.radial_scale_params[3] * x**3 - ) + return np.polyval(self.radial_scale_params, x) def sf_func(self, x): return self.az_scale_params[0] * np.exp( @@ -413,6 +427,9 @@ def update_model_parameters(self, model_params): def update_location(self, r, f): """ Updates the location on the camera focal plane from where the psf is calculated + This means it updates the parameters L and phi0 in the laplacian functions of the psf: + .. math:: f_{\Phi}(\phi) = \frac{1}{2S_\phi}e^{-|\frac{\phi-\phi_0}{S_\phi}|} + .. math:: f_{R}(r, K) = \begin{cases}\frac{1}{S_{R}(K+K^{-1})}e^{-K\frac{r-L}{S_{R}}}, r\ge L\\ \frac{1}{S_{R}(K+K^{-1})}e^{\frac{r-L}{KS_{R}}}, r < L\end{cases} Parameters ---------- From 3697b35964df5ad803fa91e59726f01cc61eeda8 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Wed, 13 Nov 2024 21:04:42 +0100 Subject: [PATCH 10/33] Fixing issues in documentation --- docs/references.bib | 15 +++++++-------- src/ctapipe/instrument/optics.py | 8 ++++---- 2 files changed, 11 insertions(+), 12 deletions(-) diff --git a/docs/references.bib b/docs/references.bib index 65146670f51..453dc9d85a6 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -138,12 +138,11 @@ @article{bright-star-catalog } @article{startracker, - author = {Abe, K. and others} - title = {Star tracking for pointing determination of Imaging Atmospheric Cherenkov Telescopes - Application to the Large-Sized Telescope of the Cherenkov Telescope Array}, - doi= "10.1051/0004-6361/202347128", - url= "https://doi.org/10.1051/0004-6361/202347128", - journal = {A&A}, - year = 2023, - volume = 679, - pages = "A90", + author={Abe, K. and others}, + title={Star tracking for pointing determination of Imaging Atmospheric Cherenkov Telescopes - Application to the Large-Sized Telescope of the Cherenkov Telescope Array}, + doi={10.1051/0004-6361/202347128}, + journal={A\&A}, + year={2023}, + volume={679}, + pages={A90}, } diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index b6cad715979..d6338944eb8 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -334,7 +334,7 @@ def __init__( ---------- asymmetry_params: list Parameters describing the dependency of the asymmetry of the psf on the distance to the center of the camera - Used to calculate a pdf asymmetry parameter K of the asymmetric radial laplacian of the psf as a function of the distance r to the optical axis: + Used to calculate a pdf asymmetry parameter K of the asymmetric radial laplacian of the psf as a function of the distance r to the optical axis .. math:: K(r) = 1 - asym_0 \tanh(asym_1 r) - asym_2 r radial_scale_params : list Parameters describing the dependency of the radial scale on the distance to the center of the camera @@ -342,8 +342,8 @@ def __init__( .. math:: S_{R}(r) & = b_1 - b_2\,r + b_3\,r^2 + b_4\,r^3 az_scale_params : list Parameters describing the dependency of the azimuthal scale on the distance to the center of the camera - Used to calculate width Sf of the azimuthal laplacian in the PSF as a function of the angle phi: - .. math:: S_{\phi}(r) & = a_1\,\exp{(-a_2\,r)}+\frac{a_3}{a_3+r} + Used to calculate width Sf of the azimuthal laplacian in the PSF as a function of the angle phi: + .. math:: S_{\phi}(r) & = a_1\,\exp{(-a_2\,r)}+\frac{a_3}{a_3+r} """ if ( len(asymmetry_params) == 3 @@ -427,7 +427,7 @@ def update_model_parameters(self, model_params): def update_location(self, r, f): """ Updates the location on the camera focal plane from where the psf is calculated - This means it updates the parameters L and phi0 in the laplacian functions of the psf: + This means it updates the parameters L and phi0 in the laplacian functions of the psf .. math:: f_{\Phi}(\phi) = \frac{1}{2S_\phi}e^{-|\frac{\phi-\phi_0}{S_\phi}|} .. math:: f_{R}(r, K) = \begin{cases}\frac{1}{S_{R}(K+K^{-1})}e^{-K\frac{r-L}{S_{R}}}, r\ge L\\ \frac{1}{S_{R}(K+K^{-1})}e^{\frac{r-L}{KS_{R}}}, r < L\end{cases} From e8d29b834cadd85a5267004eb58d6dc453d0a524 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Thu, 14 Nov 2024 08:58:01 +0100 Subject: [PATCH 11/33] update to docustrings --- src/ctapipe/instrument/optics.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index d6338944eb8..28915968dd8 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -325,9 +325,9 @@ def __init__( """ PSF model, describing purely the effect of coma aberration on the PSF for reference see :cite:p:`startracker` - uses an asymmetric laplacian for the radial part: + uses an asymmetric laplacian for the radial part .. math:: f_{R}(r, K) = \begin{cases}\frac{1}{S_{R}(K+K^{-1})}e^{-K\frac{r-L}{S_{R}}}, r\ge L\\ \frac{1}{S_{R}(K+K^{-1})}e^{\frac{r-L}{KS_{R}}}, r < L\end{cases} - and a symmetric laplacian in azimuthal direction: + and a symmetric laplacian in azimuthal direction .. math:: f_{\Phi}(\phi) = \frac{1}{2S_\phi}e^{-|\frac{\phi-\phi_0}{S_\phi}|} Parameters @@ -338,11 +338,11 @@ def __init__( .. math:: K(r) = 1 - asym_0 \tanh(asym_1 r) - asym_2 r radial_scale_params : list Parameters describing the dependency of the radial scale on the distance to the center of the camera - Used to calculate width Sr of the asymmetric radial laplacian in the PSF as a of function the distance r to the optical axis: + Used to calculate width Sr of the asymmetric radial laplacian in the PSF as a of function the distance r to the optical axis .. math:: S_{R}(r) & = b_1 - b_2\,r + b_3\,r^2 + b_4\,r^3 az_scale_params : list Parameters describing the dependency of the azimuthal scale on the distance to the center of the camera - Used to calculate width Sf of the azimuthal laplacian in the PSF as a function of the angle phi: + Used to calculate width Sf of the azimuthal laplacian in the PSF as a function of the angle phi .. math:: S_{\phi}(r) & = a_1\,\exp{(-a_2\,r)}+\frac{a_3}{a_3+r} """ if ( From 9ddc35c505c7c28d824c71db4ef50a54a47f0bc5 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Thu, 14 Nov 2024 09:28:19 +0100 Subject: [PATCH 12/33] fixing issue in docustring --- src/ctapipe/instrument/optics.py | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index 28915968dd8..e21e0793e23 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -326,8 +326,11 @@ def __init__( PSF model, describing purely the effect of coma aberration on the PSF for reference see :cite:p:`startracker` uses an asymmetric laplacian for the radial part + .. math:: f_{R}(r, K) = \begin{cases}\frac{1}{S_{R}(K+K^{-1})}e^{-K\frac{r-L}{S_{R}}}, r\ge L\\ \frac{1}{S_{R}(K+K^{-1})}e^{\frac{r-L}{KS_{R}}}, r < L\end{cases} + and a symmetric laplacian in azimuthal direction + .. math:: f_{\Phi}(\phi) = \frac{1}{2S_\phi}e^{-|\frac{\phi-\phi_0}{S_\phi}|} Parameters @@ -335,15 +338,21 @@ def __init__( asymmetry_params: list Parameters describing the dependency of the asymmetry of the psf on the distance to the center of the camera Used to calculate a pdf asymmetry parameter K of the asymmetric radial laplacian of the psf as a function of the distance r to the optical axis + .. math:: K(r) = 1 - asym_0 \tanh(asym_1 r) - asym_2 r + radial_scale_params : list Parameters describing the dependency of the radial scale on the distance to the center of the camera Used to calculate width Sr of the asymmetric radial laplacian in the PSF as a of function the distance r to the optical axis + .. math:: S_{R}(r) & = b_1 - b_2\,r + b_3\,r^2 + b_4\,r^3 + az_scale_params : list Parameters describing the dependency of the azimuthal scale on the distance to the center of the camera Used to calculate width Sf of the azimuthal laplacian in the PSF as a function of the angle phi + .. math:: S_{\phi}(r) & = a_1\,\exp{(-a_2\,r)}+\frac{a_3}{a_3+r} + """ if ( len(asymmetry_params) == 3 @@ -427,8 +436,12 @@ def update_model_parameters(self, model_params): def update_location(self, r, f): """ Updates the location on the camera focal plane from where the psf is calculated - This means it updates the parameters L and phi0 in the laplacian functions of the psf - .. math:: f_{\Phi}(\phi) = \frac{1}{2S_\phi}e^{-|\frac{\phi-\phi_0}{S_\phi}|} + This means it updates the parameters :math:`\phi_0` in the azimuthal laplacian function of the psf + + .. math:: f_{\Phi}(\phi) = \frac{1}{2 S_{\phi}}e^{-|\frac{\phi-\phi_0}{S_{\phi}}|} + + and L in the radial laplacian function of the psf + .. math:: f_{R}(r, K) = \begin{cases}\frac{1}{S_{R}(K+K^{-1})}e^{-K\frac{r-L}{S_{R}}}, r\ge L\\ \frac{1}{S_{R}(K+K^{-1})}e^{\frac{r-L}{KS_{R}}}, r < L\end{cases} Parameters From 2a3ce78d41ad2edaf1cd4210ad9631b6df37a877 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Thu, 14 Nov 2024 09:41:58 +0100 Subject: [PATCH 13/33] removing the update_location method --- src/ctapipe/instrument/optics.py | 42 ++++++++++---------------------- 1 file changed, 13 insertions(+), 29 deletions(-) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index e21e0793e23..022cd5d9b15 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -384,22 +384,31 @@ def sf_func(self, x): -self.az_scale_params[1] * x ) + self.az_scale_params[2] / (self.az_scale_params[2] + x) - def pdf(self, r, f): + def pdf(self, r, f, r0, f0): """ Calculates the value of the psf at a given location Parameters ---------- r : float - distance to the center of the camera in meters + distance to the center of the camera in meters, location at where the PSF is evaluated f : float - polar angle in radians - + polar angle in radians, location at where the PSF is evaluated + r0 : float + distance to the center of the camera in meters, location from where the PSF is evaluated + f0 : float + polar angle in radians, location from where the PSF is evaluated Returns ---------- psf : float value of the PSF at the specified location """ + k = self.k_func(r0) + sr = self.sr_func(r0) + sf = self.sf_func(r0) + self.radial_pdf_params = (k, r0, sr) + self.azimuthal_pdf_params = (f0, sf) + return laplace_asymmetric.pdf(r, *self.radial_pdf_params) * laplace.pdf( f, *self.azimuthal_pdf_params ) @@ -432,28 +441,3 @@ def update_model_parameters(self, model_params): sf = self.sf_func(self.radial_pdf_params[1]) self.radial_pdf_params = (k, self.radial_pdf_params[1], sr) self.azimuthal_pdf_params = (self.azimuthal_pdf_params[0], sf) - - def update_location(self, r, f): - """ - Updates the location on the camera focal plane from where the psf is calculated - This means it updates the parameters :math:`\phi_0` in the azimuthal laplacian function of the psf - - .. math:: f_{\Phi}(\phi) = \frac{1}{2 S_{\phi}}e^{-|\frac{\phi-\phi_0}{S_{\phi}}|} - - and L in the radial laplacian function of the psf - - .. math:: f_{R}(r, K) = \begin{cases}\frac{1}{S_{R}(K+K^{-1})}e^{-K\frac{r-L}{S_{R}}}, r\ge L\\ \frac{1}{S_{R}(K+K^{-1})}e^{\frac{r-L}{KS_{R}}}, r < L\end{cases} - - Parameters - ---------- - r : float - distance to the center of the camera in meters - f : float - polar angle in radians - """ - - k = self.k_func(r) - sr = self.sr_func(r) - sf = self.sf_func(r) - self.radial_pdf_params = (k, r, sr) - self.azimuthal_pdf_params = (f, sf) From a1f0b1e99aaff2399a413741ebeec2843549631b Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Thu, 14 Nov 2024 09:54:23 +0100 Subject: [PATCH 14/33] fixing more issues with docustrings --- src/ctapipe/instrument/optics.py | 4 +--- src/ctapipe/instrument/tests/test_psf_model.py | 3 +-- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index 022cd5d9b15..d6ad71ff18e 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -322,7 +322,7 @@ def __init__( radial_scale_params=[0.01409259, 0.02947208, 0.06000271, -0.02969355], az_scale_params=[0.24271557, 7.5511501, 0.02037972], ): - """ + r""" PSF model, describing purely the effect of coma aberration on the PSF for reference see :cite:p:`startracker` uses an asymmetric laplacian for the radial part @@ -367,8 +367,6 @@ def __init__( "asymmetry_params and az_scale_params needs to have length 3 and radial_scale_params length 4" ) - self.update_location(0.0, 0.0) - def k_func(self, x): return ( 1 diff --git a/src/ctapipe/instrument/tests/test_psf_model.py b/src/ctapipe/instrument/tests/test_psf_model.py index 5a0f1a21bdc..2bc52628046 100644 --- a/src/ctapipe/instrument/tests/test_psf_model.py +++ b/src/ctapipe/instrument/tests/test_psf_model.py @@ -23,5 +23,4 @@ def test_psf(): def test_asymptotic_behavior(): psf_coma = PSFModel.subclass_from_name("ComaModel") - psf_coma.update_location(1.0, 0.0) - assert np.isclose(psf_coma.pdf(10.0, 0.0), 0.0) + assert np.isclose(psf_coma.pdf(10.0, 0.0, 1.0, 0.0), 0.0) From 5041e52490727ef444322434e9744619c511916e Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Thu, 14 Nov 2024 17:10:33 +0100 Subject: [PATCH 15/33] Refactoring as `TelescopeComponent` --- src/ctapipe/instrument/optics.py | 147 +++++++----------- .../instrument/tests/test_psf_model.py | 26 ++-- 2 files changed, 74 insertions(+), 99 deletions(-) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index d6ad71ff18e..cc464e733d5 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -11,7 +11,8 @@ from astropy.table import QTable from scipy.stats import laplace, laplace_asymmetric -from ctapipe.core.component import non_abstract_children +from ctapipe.core import TelescopeComponent +from ctapipe.core.traits import List from ..compat import StrEnum from ..utils import get_table_dataset @@ -271,56 +272,69 @@ def __str__(self): return self.name -class PSFModel: - def __init__(self, **kwargs): +class PSFModel(TelescopeComponent): + def __init__(self, subarray, **kwargs): """ Base component to describe image distortion due to the optics of the different cameras. """ - - @classmethod - def subclass_from_name(cls, name, **kwargs): - """ - Obtain an instance of a subclass via its name - - Parameters - ---------- - name : str - Name of the subclass to obtain - - Returns - ------- - Instance - Instance of subclass to this class - """ - - subclasses = {base.__name__: base for base in non_abstract_children(cls)} - - requested_subclass = subclasses[name] - return requested_subclass(**kwargs) + super().__init__(subarray=subarray, **kwargs) @abstractmethod def pdf(self, *args): pass - @abstractmethod - def update_location(self, *args): - pass - - @abstractmethod - def update_model_parameters(self, *args): - pass - class ComaModel(PSFModel): - """ + r""" PSF model, describing pure coma aberrations PSF effect + asymmetry_params describe the dependency of the psf on the distance to the center of the camera + Used to calculate a pdf asymmetry parameter K of the asymmetric radial laplacian of the psf as + a function of the distance r to the optical axis + + .. math:: K(r) = 1 - asym_0 \tanh(asym_1 r) - asym_2 r + + radial_scale_params describes the dependency of the radial scale on the distance to the center of the camera + Used to calculate width Sr of the asymmetric radial laplacian in the PSF as a of function the distance r to the optical axis + + .. math:: S_{R}(r) & = b_1 - b_2\,r + b_3\,r^2 + b_4\,r^3 + + az_scale_params Describes the dependency of the azimuthal scale on the distance to the center of the camera + Used to calculate the width Sf of the azimuthal laplacian in the PSF as a function of the angle :math:`phi` + + .. math:: S_{\phi}(r) & = a_1\,\exp{(-a_2\,r)}+\frac{a_3}{a_3+r} """ + asymmetry_params = List( + help=( + "Describes the dependency of the psf on the distance" + "to the center of the camera. Used to calculate a pdf" + "asymmetry parameter K of the asymmetric radial laplacian" + "of the psf as a function of the distance r to the optical axis" + ) + ).tag(config=True) + + radial_scale_params = List( + help=( + "Describes the dependency of the radial scale on the" + "distance to the center of the camera Used to calculate" + "width Sr of the asymmetric radial laplacian in the PSF" + "as a of function the distance r to the optical axis" + ) + ).tag(config=True) + + az_scale_params = List( + help=( + "Describes the dependency of the azimuthal scale on the" + "distance to the center of the camera. Used to calculate" + "the the width Sf of the azimuthal laplacian in the PSF" + "as a function of the angle phi" + ) + ).tag(config=True) + def __init__( self, - asymmetry_params=[0.49244797, 9.23573115, 0.15216096], - radial_scale_params=[0.01409259, 0.02947208, 0.06000271, -0.02969355], - az_scale_params=[0.24271557, 7.5511501, 0.02037972], + subarray, + **kwargs, ): r""" PSF model, describing purely the effect of coma aberration on the PSF @@ -335,37 +349,11 @@ def __init__( Parameters ---------- - asymmetry_params: list - Parameters describing the dependency of the asymmetry of the psf on the distance to the center of the camera - Used to calculate a pdf asymmetry parameter K of the asymmetric radial laplacian of the psf as a function of the distance r to the optical axis - - .. math:: K(r) = 1 - asym_0 \tanh(asym_1 r) - asym_2 r - - radial_scale_params : list - Parameters describing the dependency of the radial scale on the distance to the center of the camera - Used to calculate width Sr of the asymmetric radial laplacian in the PSF as a of function the distance r to the optical axis - - .. math:: S_{R}(r) & = b_1 - b_2\,r + b_3\,r^2 + b_4\,r^3 - - az_scale_params : list - Parameters describing the dependency of the azimuthal scale on the distance to the center of the camera - Used to calculate width Sf of the azimuthal laplacian in the PSF as a function of the angle phi - - .. math:: S_{\phi}(r) & = a_1\,\exp{(-a_2\,r)}+\frac{a_3}{a_3+r} - + subarray: ctapipe.instrument.SubarrayDescription + Description of the subarray. """ - if ( - len(asymmetry_params) == 3 - and len(radial_scale_params) == 4 - and len(az_scale_params) == 3 - ): - self.asymmetry_params = asymmetry_params - self.radial_scale_params = radial_scale_params - self.az_scale_params = az_scale_params - else: - raise ValueError( - "asymmetry_params and az_scale_params needs to have length 3 and radial_scale_params length 4" - ) + super().__init__(subarray=subarray, **kwargs) + self.check_model_parameters() def k_func(self, x): return ( @@ -411,31 +399,12 @@ def pdf(self, r, f, r0, f0): f, *self.azimuthal_pdf_params ) - def update_model_parameters(self, model_params): - """ - Updates the model parameters for the psf - - Parameters - ---------- - model_params : dict - dictionary with the model parameters - needs to have the keys asymmetry_params, radial_scale_params and az_scale_params - The values need to be lists of length 3, 4 and 3 respectively - """ + def check_model_parameters(self): if not ( - len(model_params["asymmetry_params"]) == 3 - and len(model_params["radial_scale_params"]) == 4 - and len(model_params["az_scale_params"]) == 3 + len(self.asymmetry_params) == 3 + and len(self.radial_scale_params) == 4 + and len(self.az_scale_params) == 3 ): raise ValueError( "asymmetry_params and az_scale_params needs to have length 3 and radial_scale_params length 4" ) - - self.asymmetry_params = model_params["asymmetry_params"] - self.radial_scale_params = model_params["radial_scale_params"] - self.az_scale_params = model_params["az_scale_params"] - k = self.k_func(self.radial_pdf_params[1]) - sr = self.sr_func(self.radial_pdf_params[1]) - sf = self.sf_func(self.radial_pdf_params[1]) - self.radial_pdf_params = (k, self.radial_pdf_params[1], sr) - self.azimuthal_pdf_params = (self.azimuthal_pdf_params[0], sf) diff --git a/src/ctapipe/instrument/tests/test_psf_model.py b/src/ctapipe/instrument/tests/test_psf_model.py index 2bc52628046..c0da8f5ce4e 100644 --- a/src/ctapipe/instrument/tests/test_psf_model.py +++ b/src/ctapipe/instrument/tests/test_psf_model.py @@ -7,20 +7,26 @@ from ctapipe.instrument.optics import PSFModel -def test_psf(): - psf = PSFModel.subclass_from_name("ComaModel") - wrong_parameters = { - "asymmetry_params": [0.49244797, 9.23573115, 0.15216096], - "radial_scale_params": [0.01409259, 0.02947208, 0.06000271, -0.02969355], - "az_scale_params": [0.24271557, 7.5511501], - } +def test_psf(example_subarray): with pytest.raises( ValueError, match="asymmetry_params and az_scale_params needs to have length 3 and radial_scale_params length 4", ): - psf.update_model_parameters(wrong_parameters) + PSFModel.from_name( + "ComaModel", + subarray=example_subarray, + asymmetry_params=[0.49244797, 9.23573115, 0.15216096], + radial_scale_params=[0.01409259, 0.02947208, 0.06000271, -0.02969355], + az_scale_params=[0.24271557, 7.5511501], + ) -def test_asymptotic_behavior(): - psf_coma = PSFModel.subclass_from_name("ComaModel") +def test_asymptotic_behavior(example_subarray): + psf_coma = PSFModel.from_name( + "ComaModel", + subarray=example_subarray, + asymmetry_params=[0.49244797, 9.23573115, 0.15216096], + radial_scale_params=[0.01409259, 0.02947208, 0.06000271, -0.02969355], + az_scale_params=[0.24271557, 7.5511501, 0.02037972], + ) assert np.isclose(psf_coma.pdf(10.0, 0.0, 1.0, 0.0), 0.0) From 30233b048456e74fea48445da9cbc024760dd364 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Fri, 15 Nov 2024 10:27:22 +0100 Subject: [PATCH 16/33] Moving some documentation --- src/ctapipe/instrument/optics.py | 47 ++++++++++++++++++++++---------- 1 file changed, 32 insertions(+), 15 deletions(-) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index cc464e733d5..a971f40f7b2 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -273,14 +273,30 @@ def __str__(self): class PSFModel(TelescopeComponent): - def __init__(self, subarray, **kwargs): - """ - Base component to describe image distortion due to the optics of the different cameras. - """ - super().__init__(subarray=subarray, **kwargs) + """ + Base component to describe image distortion due to the optics of the different cameras. + """ @abstractmethod - def pdf(self, *args): + def pdf(self, r, f, r0, f0, *args): + """ + Calculates the value of the psf at a given location + + Parameters + ---------- + r : float + distance to the center of the camera in meters, location at where the PSF is evaluated + f : float + polar angle in radians, location at where the PSF is evaluated + r0 : float + distance to the center of the camera in meters, location from where the PSF is evaluated + f0 : float + polar angle in radians, location from where the PSF is evaluated + Returns + ---------- + psf : float + value of the PSF at the specified location + """ pass @@ -302,6 +318,8 @@ class ComaModel(PSFModel): Used to calculate the width Sf of the azimuthal laplacian in the PSF as a function of the angle :math:`phi` .. math:: S_{\phi}(r) & = a_1\,\exp{(-a_2\,r)}+\frac{a_3}{a_3+r} + + for reference see :cite:p:`startracker` """ asymmetry_params = List( @@ -337,20 +355,19 @@ def __init__( **kwargs, ): r""" - PSF model, describing purely the effect of coma aberration on the PSF - for reference see :cite:p:`startracker` + PSF model, describing purely the effect of coma aberration on the PSF uses an asymmetric laplacian for the radial part - .. math:: f_{R}(r, K) = \begin{cases}\frac{1}{S_{R}(K+K^{-1})}e^{-K\frac{r-L}{S_{R}}}, r\ge L\\ \frac{1}{S_{R}(K+K^{-1})}e^{\frac{r-L}{KS_{R}}}, r < L\end{cases} + .. math:: f_{R}(r, K) = \begin{cases}\frac{1}{S_{R}(K+K^{-1})}e^{-K\frac{r-L}{S_{R}}}, r\ge L\\ \frac{1}{S_{R}(K+K^{-1})}e^{\frac{r-L}{KS_{R}}}, r < L\end{cases} - and a symmetric laplacian in azimuthal direction + and a symmetric laplacian in azimuthal direction - .. math:: f_{\Phi}(\phi) = \frac{1}{2S_\phi}e^{-|\frac{\phi-\phi_0}{S_\phi}|} + .. math:: f_{\Phi}(\phi) = \frac{1}{2S_\phi}e^{-|\frac{\phi-\phi_0}{S_\phi}|} - Parameters - ---------- - subarray: ctapipe.instrument.SubarrayDescription - Description of the subarray. + Parameters + ---------- + subarray: ctapipe.instrument.SubarrayDescription + Description of the subarray. """ super().__init__(subarray=subarray, **kwargs) self.check_model_parameters() From 6dfc537fe373619fb269ffba73f23c045cabe83c Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Fri, 15 Nov 2024 13:46:05 +0100 Subject: [PATCH 17/33] Fixing some docustring --- src/ctapipe/instrument/optics.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index a971f40f7b2..48a60fa56dc 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -358,16 +358,16 @@ def __init__( PSF model, describing purely the effect of coma aberration on the PSF uses an asymmetric laplacian for the radial part - .. math:: f_{R}(r, K) = \begin{cases}\frac{1}{S_{R}(K+K^{-1})}e^{-K\frac{r-L}{S_{R}}}, r\ge L\\ \frac{1}{S_{R}(K+K^{-1})}e^{\frac{r-L}{KS_{R}}}, r < L\end{cases} + .. math:: f_{R}(r, K) = \begin{cases}\frac{1}{S_{R}(K+K^{-1})}e^{-K\frac{r-L}{S_{R}}}, r\ge L\\ \frac{1}{S_{R}(K+K^{-1})}e^{\frac{r-L}{KS_{R}}}, r < L\end{cases} - and a symmetric laplacian in azimuthal direction + and a symmetric laplacian in azimuthal direction - .. math:: f_{\Phi}(\phi) = \frac{1}{2S_\phi}e^{-|\frac{\phi-\phi_0}{S_\phi}|} + .. math:: f_{\Phi}(\phi) = \frac{1}{2S_\phi}e^{-|\frac{\phi-\phi_0}{S_\phi}|} - Parameters - ---------- - subarray: ctapipe.instrument.SubarrayDescription - Description of the subarray. + Parameters + ---------- + subarray: ctapipe.instrument.SubarrayDescription + Description of the subarray. """ super().__init__(subarray=subarray, **kwargs) self.check_model_parameters() From 36c8748277e13b710d0aa1c8e4d8f5425a1e185b Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Wed, 27 Nov 2024 09:46:12 +0100 Subject: [PATCH 18/33] fixing issue in reference bibtex file --- docs/references.bib | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/references.bib b/docs/references.bib index 1dab43771dc..dde248295fb 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -146,7 +146,7 @@ @article{startracker year={2023}, volume={679}, pages={A90}, - +} @manual{ctao-software-standards, title={Software Programming Standards}, From 3518a8b2637677eadc83edd7e1a63f3219368f40 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Mon, 2 Dec 2024 17:28:08 +0100 Subject: [PATCH 19/33] implementing reviewer comments --- docs/changes/2643.feature.rst | 4 +-- src/ctapipe/instrument/optics.py | 2 +- .../instrument/tests/test_psf_model.py | 33 ++++++++++++++----- 3 files changed, 27 insertions(+), 12 deletions(-) diff --git a/docs/changes/2643.feature.rst b/docs/changes/2643.feature.rst index 336b0b65f30..0b799d5623e 100644 --- a/docs/changes/2643.feature.rst +++ b/docs/changes/2643.feature.rst @@ -1,4 +1,2 @@ -Adds psf models to ``ctapipe.instrument.optics`` with the parent class ``PSFModel`` and a psf model based on pure coma abbaration called ``ComaModel`` -- All psf models have a function ``PSFModel.update_location`` that updates the location in the focal plane where the PSF is evaluated from -- The function ``PSFModel.update_model_parameters`` allows to update the parameters of the PSF model +Adds psf models to ``ctapipe.instrument.optics`` with the parent class ``PSFModel`` and a psf model based on pure coma aberration called ``ComaModel`` - The function ``PSFModel.pdf`` gives the value of the PSF in a given location diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index 48a60fa56dc..7251bfe9c9d 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -312,7 +312,7 @@ class ComaModel(PSFModel): radial_scale_params describes the dependency of the radial scale on the distance to the center of the camera Used to calculate width Sr of the asymmetric radial laplacian in the PSF as a of function the distance r to the optical axis - .. math:: S_{R}(r) & = b_1 - b_2\,r + b_3\,r^2 + b_4\,r^3 + .. math:: S_{R}(r) & = b_1 + b_2\,r + b_3\,r^2 + b_4\,r^3 az_scale_params Describes the dependency of the azimuthal scale on the distance to the center of the camera Used to calculate the width Sf of the azimuthal laplacian in the PSF as a function of the angle :math:`phi` diff --git a/src/ctapipe/instrument/tests/test_psf_model.py b/src/ctapipe/instrument/tests/test_psf_model.py index c0da8f5ce4e..6883c3b7171 100644 --- a/src/ctapipe/instrument/tests/test_psf_model.py +++ b/src/ctapipe/instrument/tests/test_psf_model.py @@ -7,7 +7,22 @@ from ctapipe.instrument.optics import PSFModel -def test_psf(example_subarray): +@pytest.fixture(scope="session") +def asymmetry_params(): + return [0.49244797, 9.23573115, 0.15216096] + + +@pytest.fixture(scope="session") +def radial_scale_params(): + return [0.01409259, -0.02947208, 0.06000271, -0.02969355] + + +@pytest.fixture(scope="session") +def az_scale_params(): + return [0.24271557, 7.5511501, 0.02037972] + + +def test_psf(example_subarray, asymmetry_params, radial_scale_params): with pytest.raises( ValueError, match="asymmetry_params and az_scale_params needs to have length 3 and radial_scale_params length 4", @@ -15,18 +30,20 @@ def test_psf(example_subarray): PSFModel.from_name( "ComaModel", subarray=example_subarray, - asymmetry_params=[0.49244797, 9.23573115, 0.15216096], - radial_scale_params=[0.01409259, 0.02947208, 0.06000271, -0.02969355], - az_scale_params=[0.24271557, 7.5511501], + asymmetry_params=asymmetry_params, + radial_scale_params=radial_scale_params, + az_scale_params=[0.0], ) -def test_asymptotic_behavior(example_subarray): +def test_asymptotic_behavior( + example_subarray, asymmetry_params, radial_scale_params, az_scale_params +): psf_coma = PSFModel.from_name( "ComaModel", subarray=example_subarray, - asymmetry_params=[0.49244797, 9.23573115, 0.15216096], - radial_scale_params=[0.01409259, 0.02947208, 0.06000271, -0.02969355], - az_scale_params=[0.24271557, 7.5511501, 0.02037972], + asymmetry_params=asymmetry_params, + radial_scale_params=radial_scale_params, + az_scale_params=az_scale_params, ) assert np.isclose(psf_coma.pdf(10.0, 0.0, 1.0, 0.0), 0.0) From b7d6498225a4d1ea93024e764711059091c4f1e1 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Tue, 3 Dec 2024 13:27:09 +0100 Subject: [PATCH 20/33] changing tests --- src/ctapipe/instrument/optics.py | 24 ++++++----- .../instrument/tests/test_psf_model.py | 40 +++++++------------ 2 files changed, 28 insertions(+), 36 deletions(-) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index 7251bfe9c9d..efce24db019 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -9,6 +9,7 @@ import astropy.units as u import numpy as np from astropy.table import QTable +from erfa.ufunc import p2s as cartesian_to_spherical from scipy.stats import laplace, laplace_asymmetric from ctapipe.core import TelescopeComponent @@ -387,25 +388,28 @@ def sf_func(self, x): -self.az_scale_params[1] * x ) + self.az_scale_params[2] / (self.az_scale_params[2] + x) - def pdf(self, r, f, r0, f0): + def pdf(self, x, y, x0, y0): """ Calculates the value of the psf at a given location Parameters ---------- - r : float - distance to the center of the camera in meters, location at where the PSF is evaluated - f : float - polar angle in radians, location at where the PSF is evaluated - r0 : float - distance to the center of the camera in meters, location from where the PSF is evaluated - f0 : float - polar angle in radians, location from where the PSF is evaluated + x : float + x-coordinate of the point on the focal plane where the psf is evaluated in meters + y : float + y-coordinate of the point on the focal plane where the psf is evaluated in meters + x0 : float + x-coordinate of the point source on the focal plane in meters + y0 : float + y-coordinate of the point source on the focal plane in meters Returns ---------- psf : float - value of the PSF at the specified location + value of the PSF at the specified location with the specified position of the point source """ + f, _, r = cartesian_to_spherical((x, y, 0.0)) + f0, _, r0 = cartesian_to_spherical((x0, y0, 0.0)) + print(f, r, x, y) k = self.k_func(r0) sr = self.sr_func(r0) sf = self.sf_func(r0) diff --git a/src/ctapipe/instrument/tests/test_psf_model.py b/src/ctapipe/instrument/tests/test_psf_model.py index 6883c3b7171..372c650b50e 100644 --- a/src/ctapipe/instrument/tests/test_psf_model.py +++ b/src/ctapipe/instrument/tests/test_psf_model.py @@ -8,21 +8,18 @@ @pytest.fixture(scope="session") -def asymmetry_params(): - return [0.49244797, 9.23573115, 0.15216096] - - -@pytest.fixture(scope="session") -def radial_scale_params(): - return [0.01409259, -0.02947208, 0.06000271, -0.02969355] - - -@pytest.fixture(scope="session") -def az_scale_params(): - return [0.24271557, 7.5511501, 0.02037972] +def coma_psf(example_subarray): + psf = PSFModel.from_name( + "ComaModel", + subarray=example_subarray, + asymmetry_params=[0.5, 10, 0.15], + radial_scale_params=[0.015, -0.1, 0.06, -0.03], + az_scale_params=[0.25, 7.5, 0.02], + ) + return psf -def test_psf(example_subarray, asymmetry_params, radial_scale_params): +def test_psf(example_subarray): with pytest.raises( ValueError, match="asymmetry_params and az_scale_params needs to have length 3 and radial_scale_params length 4", @@ -30,20 +27,11 @@ def test_psf(example_subarray, asymmetry_params, radial_scale_params): PSFModel.from_name( "ComaModel", subarray=example_subarray, - asymmetry_params=asymmetry_params, - radial_scale_params=radial_scale_params, + asymmetry_params=[0.0, 0.0, 0.0], + radial_scale_params=[0.0, 0.0, 0.0], az_scale_params=[0.0], ) -def test_asymptotic_behavior( - example_subarray, asymmetry_params, radial_scale_params, az_scale_params -): - psf_coma = PSFModel.from_name( - "ComaModel", - subarray=example_subarray, - asymmetry_params=asymmetry_params, - radial_scale_params=radial_scale_params, - az_scale_params=az_scale_params, - ) - assert np.isclose(psf_coma.pdf(10.0, 0.0, 1.0, 0.0), 0.0) +def test_asymptotic_behavior(coma_psf): + assert np.isclose(coma_psf.pdf(10.0, 0.0, 1.0, 0.0), 0.0) From 4af9bdb53e576ef0613d8b2d3855b39c55a0deed Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Tue, 3 Dec 2024 13:46:49 +0100 Subject: [PATCH 21/33] fixing sign issue in psf test parameters --- src/ctapipe/instrument/optics.py | 1 - src/ctapipe/instrument/tests/test_psf_model.py | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index efce24db019..77b0e3e2f2d 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -409,7 +409,6 @@ def pdf(self, x, y, x0, y0): """ f, _, r = cartesian_to_spherical((x, y, 0.0)) f0, _, r0 = cartesian_to_spherical((x0, y0, 0.0)) - print(f, r, x, y) k = self.k_func(r0) sr = self.sr_func(r0) sf = self.sf_func(r0) diff --git a/src/ctapipe/instrument/tests/test_psf_model.py b/src/ctapipe/instrument/tests/test_psf_model.py index 372c650b50e..4150724ccfc 100644 --- a/src/ctapipe/instrument/tests/test_psf_model.py +++ b/src/ctapipe/instrument/tests/test_psf_model.py @@ -13,7 +13,7 @@ def coma_psf(example_subarray): "ComaModel", subarray=example_subarray, asymmetry_params=[0.5, 10, 0.15], - radial_scale_params=[0.015, -0.1, 0.06, -0.03], + radial_scale_params=[0.015, -0.1, 0.06, 0.03], az_scale_params=[0.25, 7.5, 0.02], ) return psf From 2f9900f3ed6e5a689c049bc0ea35a6d79ce06db7 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Tue, 3 Dec 2024 20:04:56 +0100 Subject: [PATCH 22/33] using trait validation --- src/ctapipe/core/traits.py | 2 + src/ctapipe/instrument/optics.py | 58 +++++++++++-------- .../instrument/tests/test_psf_model.py | 29 +++++++++- 3 files changed, 61 insertions(+), 28 deletions(-) diff --git a/src/ctapipe/core/traits.py b/src/ctapipe/core/traits.py index b787c3fad22..f4791d7bb69 100644 --- a/src/ctapipe/core/traits.py +++ b/src/ctapipe/core/traits.py @@ -45,6 +45,7 @@ "Unicode", "flag", "observe", + "validate", ] import logging @@ -83,6 +84,7 @@ def __repr__(self): Tuple = traitlets.Tuple observe = traitlets.observe flag = traitlets.config.boolean_flag +validate = traitlets.validate class AstroQuantity(TraitType): diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index 77b0e3e2f2d..d395816ec30 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -13,7 +13,7 @@ from scipy.stats import laplace, laplace_asymmetric from ctapipe.core import TelescopeComponent -from ctapipe.core.traits import List +from ctapipe.core.traits import List, TraitError, validate from ..compat import StrEnum from ..utils import get_table_dataset @@ -279,25 +279,26 @@ class PSFModel(TelescopeComponent): """ @abstractmethod - def pdf(self, r, f, r0, f0, *args): + def pdf(self, x, y, x0, y0, *args): """ Calculates the value of the psf at a given location Parameters ---------- - r : float - distance to the center of the camera in meters, location at where the PSF is evaluated - f : float - polar angle in radians, location at where the PSF is evaluated - r0 : float - distance to the center of the camera in meters, location from where the PSF is evaluated - f0 : float - polar angle in radians, location from where the PSF is evaluated + x : float + x-coordinate of the point on the focal plane where the psf is evaluated + y : float + y-coordinate of the point on the focal plane where the psf is evaluated + x0 : float + x-coordinate of the point source on the focal plane + y0 : float + y-coordinate of the point source on the focal plane Returns ---------- psf : float - value of the PSF at the specified location + value of the PSF at the specified location with the specified position of the point source """ + pass @@ -371,7 +372,6 @@ def __init__( Description of the subarray. """ super().__init__(subarray=subarray, **kwargs) - self.check_model_parameters() def k_func(self, x): return ( @@ -395,13 +395,13 @@ def pdf(self, x, y, x0, y0): Parameters ---------- x : float - x-coordinate of the point on the focal plane where the psf is evaluated in meters + x-coordinate of the point on the focal plane where the psf is evaluated y : float - y-coordinate of the point on the focal plane where the psf is evaluated in meters + y-coordinate of the point on the focal plane where the psf is evaluated x0 : float - x-coordinate of the point source on the focal plane in meters + x-coordinate of the point source on the focal plane y0 : float - y-coordinate of the point source on the focal plane in meters + y-coordinate of the point source on the focal plane Returns ---------- psf : float @@ -419,12 +419,20 @@ def pdf(self, x, y, x0, y0): f, *self.azimuthal_pdf_params ) - def check_model_parameters(self): - if not ( - len(self.asymmetry_params) == 3 - and len(self.radial_scale_params) == 4 - and len(self.az_scale_params) == 3 - ): - raise ValueError( - "asymmetry_params and az_scale_params needs to have length 3 and radial_scale_params length 4" - ) + @validate("asymmetry_params") + def _check_asymmetry_params(self, proposal): + if not len(proposal["value"]) == 3: + raise TraitError("asymmetry_params needs to have length 3") + return proposal["value"] + + @validate("radial_scale_params") + def _check_radial_scale_params(self, proposal): + if not len(proposal["value"]) == 4: + raise TraitError("radial_scale_params needs to have length 4") + return proposal["value"] + + @validate("az_scale_params") + def _check_az_scale_params(self, proposal): + if not len(proposal["value"]) == 3: + raise TraitError("az_scale_params needs to have length 3") + return proposal["value"] diff --git a/src/ctapipe/instrument/tests/test_psf_model.py b/src/ctapipe/instrument/tests/test_psf_model.py index 4150724ccfc..c7f88772537 100644 --- a/src/ctapipe/instrument/tests/test_psf_model.py +++ b/src/ctapipe/instrument/tests/test_psf_model.py @@ -4,6 +4,7 @@ import numpy as np import pytest +from ctapipe.core.traits import TraitError from ctapipe.instrument.optics import PSFModel @@ -21,16 +22,38 @@ def coma_psf(example_subarray): def test_psf(example_subarray): with pytest.raises( - ValueError, - match="asymmetry_params and az_scale_params needs to have length 3 and radial_scale_params length 4", + TraitError, + match="az_scale_params needs to have length 3", ): PSFModel.from_name( "ComaModel", subarray=example_subarray, asymmetry_params=[0.0, 0.0, 0.0], - radial_scale_params=[0.0, 0.0, 0.0], + radial_scale_params=[0.0, 0.0, 0.0, 0.0], az_scale_params=[0.0], ) + with pytest.raises( + TraitError, + match="radial_scale_params needs to have length 4", + ): + PSFModel.from_name( + "ComaModel", + subarray=example_subarray, + asymmetry_params=[0.0, 0.0, 0.0], + radial_scale_params=[0.0, 0.0, 0.0], + az_scale_params=[0.0, 0.0, 0.0], + ) + with pytest.raises( + TraitError, + match="asymmetry_params needs to have length 3", + ): + PSFModel.from_name( + "ComaModel", + subarray=example_subarray, + asymmetry_params=[0.0, 0.0, 0.0, 0.0], + radial_scale_params=[0.0, 0.0, 0.0, 0.0], + az_scale_params=[0.0, 0.0, 0.0], + ) def test_asymptotic_behavior(coma_psf): From c0f27cf03e71697433560a7be45b736ce1ed2461 Mon Sep 17 00:00:00 2001 From: Christoph Toennis Date: Wed, 4 Dec 2024 10:06:45 +0100 Subject: [PATCH 23/33] moving traitlets validate import statement --- src/ctapipe/core/traits.py | 2 -- src/ctapipe/instrument/optics.py | 3 ++- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/ctapipe/core/traits.py b/src/ctapipe/core/traits.py index f4791d7bb69..b787c3fad22 100644 --- a/src/ctapipe/core/traits.py +++ b/src/ctapipe/core/traits.py @@ -45,7 +45,6 @@ "Unicode", "flag", "observe", - "validate", ] import logging @@ -84,7 +83,6 @@ def __repr__(self): Tuple = traitlets.Tuple observe = traitlets.observe flag = traitlets.config.boolean_flag -validate = traitlets.validate class AstroQuantity(TraitType): diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index d395816ec30..eba72e16cc6 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -11,9 +11,10 @@ from astropy.table import QTable from erfa.ufunc import p2s as cartesian_to_spherical from scipy.stats import laplace, laplace_asymmetric +from traitlets import validate from ctapipe.core import TelescopeComponent -from ctapipe.core.traits import List, TraitError, validate +from ctapipe.core.traits import List, TraitError from ..compat import StrEnum from ..utils import get_table_dataset From c3c0af83d10d3843c9823770ae82db0d23652940 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Thu, 23 Jan 2025 19:02:08 +0100 Subject: [PATCH 24/33] Make psf work on nd inputs with units --- src/ctapipe/instrument/optics.py | 56 +++++++++++++------------------- 1 file changed, 23 insertions(+), 33 deletions(-) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index eba72e16cc6..e5addac0256 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -9,15 +9,14 @@ import astropy.units as u import numpy as np from astropy.table import QTable -from erfa.ufunc import p2s as cartesian_to_spherical from scipy.stats import laplace, laplace_asymmetric from traitlets import validate -from ctapipe.core import TelescopeComponent -from ctapipe.core.traits import List, TraitError - from ..compat import StrEnum +from ..core import TelescopeComponent +from ..core.traits import List, TraitError from ..utils import get_table_dataset +from ..utils.quantities import all_to_value from .warnings import warn_from_name logger = logging.getLogger(__name__) @@ -279,6 +278,7 @@ class PSFModel(TelescopeComponent): Base component to describe image distortion due to the optics of the different cameras. """ + @u.quantity_input(x=u.m, y=u.m, x0=u.m, y0=u.m) @abstractmethod def pdf(self, x, y, x0, y0, *args): """ @@ -286,23 +286,29 @@ def pdf(self, x, y, x0, y0, *args): Parameters ---------- - x : float + x : u.Quantity[length] x-coordinate of the point on the focal plane where the psf is evaluated - y : float + y : u.Quantity[length] y-coordinate of the point on the focal plane where the psf is evaluated - x0 : float + x0 : u.Quantity[length] x-coordinate of the point source on the focal plane - y0 : float + y0 : u.Quantity[length] y-coordinate of the point source on the focal plane Returns ---------- - psf : float + psf : np.ndarray value of the PSF at the specified location with the specified position of the point source """ pass +def _cartesian_to_polar(x, y): + r = np.sqrt(x**2 + y**2) + phi = np.arctan2(y, x) + return r, phi + + class ComaModel(PSFModel): r""" PSF model, describing pure coma aberrations PSF effect @@ -389,36 +395,20 @@ def sf_func(self, x): -self.az_scale_params[1] * x ) + self.az_scale_params[2] / (self.az_scale_params[2] + x) + @u.quantity_input(x=u.m, y=u.m, x0=u.m, y0=u.m) def pdf(self, x, y, x0, y0): - """ - Calculates the value of the psf at a given location + x, y, x0, y0 = all_to_value(x, y, x0, y0, unit=u.m) + r, f = _cartesian_to_polar(x, y) + r0, f0 = _cartesian_to_polar(x0, y0) - Parameters - ---------- - x : float - x-coordinate of the point on the focal plane where the psf is evaluated - y : float - y-coordinate of the point on the focal plane where the psf is evaluated - x0 : float - x-coordinate of the point source on the focal plane - y0 : float - y-coordinate of the point source on the focal plane - Returns - ---------- - psf : float - value of the PSF at the specified location with the specified position of the point source - """ - f, _, r = cartesian_to_spherical((x, y, 0.0)) - f0, _, r0 = cartesian_to_spherical((x0, y0, 0.0)) k = self.k_func(r0) sr = self.sr_func(r0) sf = self.sf_func(r0) - self.radial_pdf_params = (k, r0, sr) - self.azimuthal_pdf_params = (f0, sf) - return laplace_asymmetric.pdf(r, *self.radial_pdf_params) * laplace.pdf( - f, *self.azimuthal_pdf_params - ) + radial_pdf = laplace_asymmetric.pdf(r, k, r0, sr) + polar_pdf = laplace.pdf(f, f0, sf) + + return radial_pdf * polar_pdf @validate("asymmetry_params") def _check_asymmetry_params(self, proposal): From f8aef5d31ffaadbb5bb9b263a367cc5913c46725 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Thu, 23 Jan 2025 19:37:30 +0100 Subject: [PATCH 25/33] Improve notation with respect to paper, fix polyval evaluation (wrong order of parameters) --- src/ctapipe/instrument/optics.py | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index e5addac0256..a0d6b88d190 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -280,7 +280,7 @@ class PSFModel(TelescopeComponent): @u.quantity_input(x=u.m, y=u.m, x0=u.m, y0=u.m) @abstractmethod - def pdf(self, x, y, x0, y0, *args): + def pdf(self, x, y, x0, y0) -> np.ndarray: """ Calculates the value of the psf at a given location @@ -381,19 +381,15 @@ def __init__( super().__init__(subarray=subarray, **kwargs) def k_func(self, x): - return ( - 1 - - self.asymmetry_params[0] * np.tanh(self.asymmetry_params[1] * x) - - self.asymmetry_params[2] * x - ) + c1, c2, c3 = self.asymmetry_params + return 1 - c1 * np.tanh(c2 * x) + c3 * x def sr_func(self, x): - return np.polyval(self.radial_scale_params, x) + return np.polyval(self.radial_scale_params[::-1], x) def sf_func(self, x): - return self.az_scale_params[0] * np.exp( - -self.az_scale_params[1] * x - ) + self.az_scale_params[2] / (self.az_scale_params[2] + x) + a1, a2, a3 = self.az_scale_params + return a1 * np.exp(-a2 * x) + a3 / (a3 + x) @u.quantity_input(x=u.m, y=u.m, x0=u.m, y0=u.m) def pdf(self, x, y, x0, y0): From cc2457b87846d8be6ececa91a493d1fc467983a1 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Thu, 23 Jan 2025 19:44:05 +0100 Subject: [PATCH 26/33] More changes to bring notation closer to the paper --- src/ctapipe/instrument/optics.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index a0d6b88d190..f84e8fce041 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -380,29 +380,29 @@ def __init__( """ super().__init__(subarray=subarray, **kwargs) - def k_func(self, x): + def _k(self, r): c1, c2, c3 = self.asymmetry_params - return 1 - c1 * np.tanh(c2 * x) + c3 * x + return 1 - c1 * np.tanh(c2 * r) + c3 * r - def sr_func(self, x): - return np.polyval(self.radial_scale_params[::-1], x) + def _s_r(self, r): + return np.polyval(self.radial_scale_params[::-1], r) - def sf_func(self, x): + def _s_phi(self, r): a1, a2, a3 = self.az_scale_params - return a1 * np.exp(-a2 * x) + a3 / (a3 + x) + return a1 * np.exp(-a2 * r) + a3 / (a3 + r) @u.quantity_input(x=u.m, y=u.m, x0=u.m, y0=u.m) def pdf(self, x, y, x0, y0): x, y, x0, y0 = all_to_value(x, y, x0, y0, unit=u.m) - r, f = _cartesian_to_polar(x, y) - r0, f0 = _cartesian_to_polar(x0, y0) + r, phi = _cartesian_to_polar(x, y) + r0, phi0 = _cartesian_to_polar(x0, y0) - k = self.k_func(r0) - sr = self.sr_func(r0) - sf = self.sf_func(r0) + k = self._k(r0) + s_r = self._s_r(r0) + s_phi = self._s_phi(r0) - radial_pdf = laplace_asymmetric.pdf(r, k, r0, sr) - polar_pdf = laplace.pdf(f, f0, sf) + radial_pdf = laplace_asymmetric.pdf(r, k, r0, s_r) + polar_pdf = laplace.pdf(phi, phi0, s_phi) return radial_pdf * polar_pdf From 72459342b868b1a9e34f32fbb83b905c9756f3c4 Mon Sep 17 00:00:00 2001 From: "mykhailo.dalchenko" Date: Tue, 28 Jan 2025 10:36:58 +0100 Subject: [PATCH 27/33] Enforcing the edge case of the PSF model - The PSF Model is only applicable under assumption that polar axis is perpendicular to the radial axis - Now this is enforced, the PSF model further away will be truncated - Another limitation is that polar angle can't be defined at the center of the camera. - This is fixed by setting the polar component to a constant very close to the center of the camera --- docs/changes/2643.feature.rst | 2 +- src/ctapipe/instrument/optics.py | 103 +++++++++++------- .../instrument/tests/test_psf_model.py | 22 ++-- 3 files changed, 78 insertions(+), 49 deletions(-) diff --git a/docs/changes/2643.feature.rst b/docs/changes/2643.feature.rst index 0b799d5623e..b378d09e15a 100644 --- a/docs/changes/2643.feature.rst +++ b/docs/changes/2643.feature.rst @@ -1,2 +1,2 @@ -Adds psf models to ``ctapipe.instrument.optics`` with the parent class ``PSFModel`` and a psf model based on pure coma aberration called ``ComaModel`` +Adds psf models to ``ctapipe.instrument.optics`` with the parent class ``PSFModel`` and a psf model based on pure coma aberration called ``ComaPSFModel`` - The function ``PSFModel.pdf`` gives the value of the PSF in a given location diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index f84e8fce041..26647202715 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -14,7 +14,7 @@ from ..compat import StrEnum from ..core import TelescopeComponent -from ..core.traits import List, TraitError +from ..core.traits import Float, List, TraitError from ..utils import get_table_dataset from ..utils.quantities import all_to_value from .warnings import warn_from_name @@ -25,7 +25,7 @@ "OpticsDescription", "FocalLengthKind", "PSFModel", - "ComaModel", + "ComaPSFModel", ] @@ -309,90 +309,107 @@ def _cartesian_to_polar(x, y): return r, phi -class ComaModel(PSFModel): +class ComaPSFModel(PSFModel): r""" - PSF model, describing pure coma aberrations PSF effect - asymmetry_params describe the dependency of the psf on the distance to the center of the camera - Used to calculate a pdf asymmetry parameter K of the asymmetric radial laplacian of the psf as - a function of the distance r to the optical axis + PSF model describing pure coma aberrations PSF effect. - .. math:: K(r) = 1 - asym_0 \tanh(asym_1 r) - asym_2 r + The PSF is described by a combination of an asymmetric Laplacian for the radial part and a symmetric Laplacian in the polar direction. - radial_scale_params describes the dependency of the radial scale on the distance to the center of the camera - Used to calculate width Sr of the asymmetric radial laplacian in the PSF as a of function the distance r to the optical axis - - .. math:: S_{R}(r) & = b_1 + b_2\,r + b_3\,r^2 + b_4\,r^3 - - az_scale_params Describes the dependency of the azimuthal scale on the distance to the center of the camera - Used to calculate the width Sf of the azimuthal laplacian in the PSF as a function of the angle :math:`phi` + Attributes + ---------- + asymmetry_params : list + Describes the dependency of the PSF on the distance to the center of the camera. + Used to calculate a PDF asymmetry parameter K of the asymmetric radial Laplacian + of the PSF as a function of the distance r to the optical axis. + .. math:: K(r) = 1 - \text{asym}_0 \tanh(\text{asym}_1 r) - \text{asym}_2 r + + radial_scale_params : list + Describes the dependency of the radial scale on the distance to the center of the camera. + Used to calculate width Sr of the asymmetric radial Laplacian in the PSF as a function of the distance :math:`r` to the optical axis. + .. math:: S_{R}(r) = b_1 + b_2\,r + b_3\,r^2 + b_4\,r^3 + + phi_scale_params : list + Describes the dependency of the polar angle (:math:`\phi`) scale on the distance to the center of the camera. + Used to calculate the width Sf of the polar Laplacian in the PSF as a function of the distance :math:`r` to the optical axis. + .. math:: S_{\phi}(r) = a_1\,\exp{(-a_2\,r)}+\frac{a_3}{a_3+r} - .. math:: S_{\phi}(r) & = a_1\,\exp{(-a_2\,r)}+\frac{a_3}{a_3+r} + Parameters + ---------- + subarray : ctapipe.instrument.SubarrayDescription + Description of the subarray. - for reference see :cite:p:`startracker` + References + ---------- + For reference, see :cite:p:`startracker` """ asymmetry_params = List( help=( - "Describes the dependency of the psf on the distance" - "to the center of the camera. Used to calculate a pdf" - "asymmetry parameter K of the asymmetric radial laplacian" - "of the psf as a function of the distance r to the optical axis" + "Describes the dependency of the PSF on the distance" + "to the center of the camera. Used to calculate a PDF" + "asymmetry parameter K of the asymmetric radial Laplacian" + "of the PSF as a function of the distance r to the optical axis" ) ).tag(config=True) radial_scale_params = List( help=( "Describes the dependency of the radial scale on the" - "distance to the center of the camera Used to calculate" - "width Sr of the asymmetric radial laplacian in the PSF" - "as a of function the distance r to the optical axis" + "distance to the center of the camera. Used to calculate" + "width Sr of the asymmetric radial Laplacian in the PSF" + "as a function of the distance r to the optical axis" ) ).tag(config=True) - az_scale_params = List( + phi_scale_params = List( help=( - "Describes the dependency of the azimuthal scale on the" + "Describes the dependency of the polar scale on the" "distance to the center of the camera. Used to calculate" - "the the width Sf of the azimuthal laplacian in the PSF" - "as a function of the angle phi" + "the width Sf of the polar Laplacian in the PSF" + "as a function of the distance r to the optical axis" ) ).tag(config=True) + pixel_width = Float( + default_value=0.05, + help="Width of a pixel in the camera in meters", + ).tag(config=True) + def __init__( self, subarray, **kwargs, ): r""" - PSF model, describing purely the effect of coma aberration on the PSF - uses an asymmetric laplacian for the radial part + PSF model, describing purely the effect of coma aberration on the PSF + uses an asymmetric Laplacian for the radial part .. math:: f_{R}(r, K) = \begin{cases}\frac{1}{S_{R}(K+K^{-1})}e^{-K\frac{r-L}{S_{R}}}, r\ge L\\ \frac{1}{S_{R}(K+K^{-1})}e^{\frac{r-L}{KS_{R}}}, r < L\end{cases} - and a symmetric laplacian in azimuthal direction + and a symmetric Laplacian in azimuthal direction .. math:: f_{\Phi}(\phi) = \frac{1}{2S_\phi}e^{-|\frac{\phi-\phi_0}{S_\phi}|} Parameters ---------- - subarray: ctapipe.instrument.SubarrayDescription + subarray : ctapipe.instrument.SubarrayDescription Description of the subarray. """ super().__init__(subarray=subarray, **kwargs) def _k(self, r): c1, c2, c3 = self.asymmetry_params - return 1 - c1 * np.tanh(c2 * r) + c3 * r + return 1 - c1 * np.tanh(c2 * r) - c3 * r def _s_r(self, r): return np.polyval(self.radial_scale_params[::-1], r) def _s_phi(self, r): - a1, a2, a3 = self.az_scale_params + a1, a2, a3 = self.phi_scale_params return a1 * np.exp(-a2 * r) + a3 / (a3 + r) @u.quantity_input(x=u.m, y=u.m, x0=u.m, y0=u.m) - def pdf(self, x, y, x0, y0): + def pdf(self, x, y, x0, y0) -> np.ndarray: x, y, x0, y0 = all_to_value(x, y, x0, y0, unit=u.m) r, phi = _cartesian_to_polar(x, y) r0, phi0 = _cartesian_to_polar(x0, y0) @@ -404,6 +421,16 @@ def pdf(self, x, y, x0, y0): radial_pdf = laplace_asymmetric.pdf(r, k, r0, s_r) polar_pdf = laplace.pdf(phi, phi0, s_phi) + # Phi is not defined at the center + at_center = np.isclose(r0, 0, atol=self.pixel_width) + polar_pdf = np.where(at_center, 1 / (2 * s_phi), polar_pdf) + # Polar PDF is valid under approximation that the polar axis is orthogonal to the radial axis + # Thus, we limit the PDF to a chord of 6 pixels or covering ~30deg around the radial axis, whichever is smaller + chord_length = min(6 * self.pixel_width, 0.5 * r0) + dphi = np.arcsin(chord_length / (2 * r0)) + polar_pdf[phi < phi0 - dphi] = 0 + polar_pdf[phi > phi0 + dphi] = 0 + return radial_pdf * polar_pdf @validate("asymmetry_params") @@ -418,8 +445,8 @@ def _check_radial_scale_params(self, proposal): raise TraitError("radial_scale_params needs to have length 4") return proposal["value"] - @validate("az_scale_params") - def _check_az_scale_params(self, proposal): + @validate("phi_scale_params") + def _check_phi_scale_params(self, proposal): if not len(proposal["value"]) == 3: - raise TraitError("az_scale_params needs to have length 3") + raise TraitError("phi_scale_params needs to have length 3") return proposal["value"] diff --git a/src/ctapipe/instrument/tests/test_psf_model.py b/src/ctapipe/instrument/tests/test_psf_model.py index c7f88772537..fa5f51819bd 100644 --- a/src/ctapipe/instrument/tests/test_psf_model.py +++ b/src/ctapipe/instrument/tests/test_psf_model.py @@ -1,6 +1,8 @@ """ This module contains the ctapipe.image.psf_model unit tests """ + +import astropy.units as u import numpy as np import pytest @@ -11,11 +13,11 @@ @pytest.fixture(scope="session") def coma_psf(example_subarray): psf = PSFModel.from_name( - "ComaModel", + "ComaPSFModel", subarray=example_subarray, asymmetry_params=[0.5, 10, 0.15], radial_scale_params=[0.015, -0.1, 0.06, 0.03], - az_scale_params=[0.25, 7.5, 0.02], + phi_scale_params=[0.25, 7.5, 0.02], ) return psf @@ -23,38 +25,38 @@ def coma_psf(example_subarray): def test_psf(example_subarray): with pytest.raises( TraitError, - match="az_scale_params needs to have length 3", + match="phi_scale_params needs to have length 3", ): PSFModel.from_name( - "ComaModel", + "ComaPSFModel", subarray=example_subarray, asymmetry_params=[0.0, 0.0, 0.0], radial_scale_params=[0.0, 0.0, 0.0, 0.0], - az_scale_params=[0.0], + phi_scale_params=[0.0], ) with pytest.raises( TraitError, match="radial_scale_params needs to have length 4", ): PSFModel.from_name( - "ComaModel", + "ComaPSFModel", subarray=example_subarray, asymmetry_params=[0.0, 0.0, 0.0], radial_scale_params=[0.0, 0.0, 0.0], - az_scale_params=[0.0, 0.0, 0.0], + phi_scale_params=[0.0, 0.0, 0.0], ) with pytest.raises( TraitError, match="asymmetry_params needs to have length 3", ): PSFModel.from_name( - "ComaModel", + "ComaPSFModel", subarray=example_subarray, asymmetry_params=[0.0, 0.0, 0.0, 0.0], radial_scale_params=[0.0, 0.0, 0.0, 0.0], - az_scale_params=[0.0, 0.0, 0.0], + phi_scale_params=[0.0, 0.0, 0.0], ) def test_asymptotic_behavior(coma_psf): - assert np.isclose(coma_psf.pdf(10.0, 0.0, 1.0, 0.0), 0.0) + assert np.isclose(coma_psf.pdf(*([10.0, 0.0, 1.0, 0.0] * u.m)), 0.0) From 70438689b6d0f58e691a842be6a9ad52af3ccd8b Mon Sep 17 00:00:00 2001 From: "mykhailo.dalchenko" Date: Tue, 4 Feb 2025 14:54:24 +0100 Subject: [PATCH 28/33] Fixing division by 0 warning at the camera center --- src/ctapipe/instrument/optics.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index 26647202715..a3654708aca 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -427,9 +427,10 @@ def pdf(self, x, y, x0, y0) -> np.ndarray: # Polar PDF is valid under approximation that the polar axis is orthogonal to the radial axis # Thus, we limit the PDF to a chord of 6 pixels or covering ~30deg around the radial axis, whichever is smaller chord_length = min(6 * self.pixel_width, 0.5 * r0) - dphi = np.arcsin(chord_length / (2 * r0)) - polar_pdf[phi < phi0 - dphi] = 0 - polar_pdf[phi > phi0 + dphi] = 0 + if r0 != 0: + dphi = np.arcsin(chord_length / (2 * r0)) + polar_pdf[phi < phi0 - dphi] = 0 + polar_pdf[phi > phi0 + dphi] = 0 return radial_pdf * polar_pdf From f2f1e33e605eb1fb514ce608c12b94c8c87d7471 Mon Sep 17 00:00:00 2001 From: "mykhailo.dalchenko" Date: Tue, 4 Feb 2025 15:10:04 +0100 Subject: [PATCH 29/33] Exposing PSFModel and ComaPSFModel at the ``instrument`` module level --- src/ctapipe/instrument/__init__.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/ctapipe/instrument/__init__.py b/src/ctapipe/instrument/__init__.py index b7c540d35ad..c6e28a4fe91 100644 --- a/src/ctapipe/instrument/__init__.py +++ b/src/ctapipe/instrument/__init__.py @@ -1,7 +1,14 @@ from .atmosphere import get_atmosphere_profile_functions from .camera import CameraDescription, CameraGeometry, CameraReadout, PixelShape from .guess import guess_telescope -from .optics import FocalLengthKind, OpticsDescription, ReflectorShape, SizeType +from .optics import ( + ComaPSFModel, + FocalLengthKind, + OpticsDescription, + PSFModel, + ReflectorShape, + SizeType, +) from .subarray import SubarrayDescription, UnknownTelescopeID from .telescope import TelescopeDescription from .trigger import SoftwareTrigger @@ -23,4 +30,6 @@ "SizeType", "SoftwareTrigger", "FromNameWarning", + "PSFModel", + "ComaPSFModel", ] From f0a686e15131264c14379c0664df5157ab2e9949 Mon Sep 17 00:00:00 2001 From: "mykhailo.dalchenko" Date: Tue, 4 Feb 2025 16:46:41 +0100 Subject: [PATCH 30/33] Fix formulas in doc --- src/ctapipe/instrument/optics.py | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index a3654708aca..10054138384 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -321,16 +321,19 @@ class ComaPSFModel(PSFModel): Describes the dependency of the PSF on the distance to the center of the camera. Used to calculate a PDF asymmetry parameter K of the asymmetric radial Laplacian of the PSF as a function of the distance r to the optical axis. + .. math:: K(r) = 1 - \text{asym}_0 \tanh(\text{asym}_1 r) - \text{asym}_2 r radial_scale_params : list Describes the dependency of the radial scale on the distance to the center of the camera. Used to calculate width Sr of the asymmetric radial Laplacian in the PSF as a function of the distance :math:`r` to the optical axis. + .. math:: S_{R}(r) = b_1 + b_2\,r + b_3\,r^2 + b_4\,r^3 phi_scale_params : list Describes the dependency of the polar angle (:math:`\phi`) scale on the distance to the center of the camera. Used to calculate the width Sf of the polar Laplacian in the PSF as a function of the distance :math:`r` to the optical axis. + .. math:: S_{\phi}(r) = a_1\,\exp{(-a_2\,r)}+\frac{a_3}{a_3+r} Parameters @@ -345,34 +348,34 @@ class ComaPSFModel(PSFModel): asymmetry_params = List( help=( - "Describes the dependency of the PSF on the distance" - "to the center of the camera. Used to calculate a PDF" - "asymmetry parameter K of the asymmetric radial Laplacian" + "Describes the dependency of the PSF on the distance " + "to the center of the camera. Used to calculate a PDF " + "asymmetry parameter :math:`K` of the asymmetric radial Laplacian " "of the PSF as a function of the distance r to the optical axis" ) ).tag(config=True) radial_scale_params = List( help=( - "Describes the dependency of the radial scale on the" - "distance to the center of the camera. Used to calculate" - "width Sr of the asymmetric radial Laplacian in the PSF" + "Describes the dependency of the radial scale on the " + "distance to the center of the camera. Used to calculate " + "width :math:`S_R` of the asymmetric radial Laplacian in the PSF " "as a function of the distance r to the optical axis" ) ).tag(config=True) phi_scale_params = List( help=( - "Describes the dependency of the polar scale on the" - "distance to the center of the camera. Used to calculate" - "the width Sf of the polar Laplacian in the PSF" + "Describes the dependency of the polar scale on the " + "distance to the center of the camera. Used to calculate " + "the width :math:`S_\phi` of the polar Laplacian in the PSF " "as a function of the distance r to the optical axis" ) ).tag(config=True) pixel_width = Float( default_value=0.05, - help="Width of a pixel in the camera in meters", + help="Width of a pixel of the camera in meters", ).tag(config=True) def __init__( @@ -386,7 +389,7 @@ def __init__( .. math:: f_{R}(r, K) = \begin{cases}\frac{1}{S_{R}(K+K^{-1})}e^{-K\frac{r-L}{S_{R}}}, r\ge L\\ \frac{1}{S_{R}(K+K^{-1})}e^{\frac{r-L}{KS_{R}}}, r < L\end{cases} - and a symmetric Laplacian in azimuthal direction + and a symmetric Laplacian in polar direction .. math:: f_{\Phi}(\phi) = \frac{1}{2S_\phi}e^{-|\frac{\phi-\phi_0}{S_\phi}|} From b01fa777b44efba66b43652f3432757d449d2182 Mon Sep 17 00:00:00 2001 From: "mykhailo.dalchenko" Date: Tue, 4 Feb 2025 17:50:06 +0100 Subject: [PATCH 31/33] Update docs, remove redundant constructor --- src/ctapipe/instrument/optics.py | 30 ++++++++---------------------- 1 file changed, 8 insertions(+), 22 deletions(-) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index 10054138384..bdbb6c28031 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -314,6 +314,14 @@ class ComaPSFModel(PSFModel): PSF model describing pure coma aberrations PSF effect. The PSF is described by a combination of an asymmetric Laplacian for the radial part and a symmetric Laplacian in the polar direction. + It uses an asymmetric Laplacian for the radial part + + .. math:: f_{R}(r, K) = \begin{cases}\frac{1}{S_{R}(K+K^{-1})}e^{-K\frac{r-L}{S_{R}}}, r\ge L\\ \frac{1}{S_{R}(K+K^{-1})}e^{\frac{r-L}{KS_{R}}}, r < L\end{cases} + + and a symmetric Laplacian in polar direction + + .. math:: f_{\Phi}(\phi) = \frac{1}{2S_\phi}e^{-|\frac{\phi-\phi_0}{S_\phi}|} + Attributes ---------- @@ -378,28 +386,6 @@ class ComaPSFModel(PSFModel): help="Width of a pixel of the camera in meters", ).tag(config=True) - def __init__( - self, - subarray, - **kwargs, - ): - r""" - PSF model, describing purely the effect of coma aberration on the PSF - uses an asymmetric Laplacian for the radial part - - .. math:: f_{R}(r, K) = \begin{cases}\frac{1}{S_{R}(K+K^{-1})}e^{-K\frac{r-L}{S_{R}}}, r\ge L\\ \frac{1}{S_{R}(K+K^{-1})}e^{\frac{r-L}{KS_{R}}}, r < L\end{cases} - - and a symmetric Laplacian in polar direction - - .. math:: f_{\Phi}(\phi) = \frac{1}{2S_\phi}e^{-|\frac{\phi-\phi_0}{S_\phi}|} - - Parameters - ---------- - subarray : ctapipe.instrument.SubarrayDescription - Description of the subarray. - """ - super().__init__(subarray=subarray, **kwargs) - def _k(self, r): c1, c2, c3 = self.asymmetry_params return 1 - c1 * np.tanh(c2 * r) - c3 * r From 16562ddf44ee807ba0685b538ab2f8496311b089 Mon Sep 17 00:00:00 2001 From: "mykhailo.dalchenko" Date: Thu, 6 Feb 2025 13:20:35 +0100 Subject: [PATCH 32/33] Update docstrings addressing @Tobychev comments --- src/ctapipe/instrument/optics.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index bdbb6c28031..d14783dc33e 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -313,15 +313,17 @@ class ComaPSFModel(PSFModel): r""" PSF model describing pure coma aberrations PSF effect. - The PSF is described by a combination of an asymmetric Laplacian for the radial part and a symmetric Laplacian in the polar direction. - It uses an asymmetric Laplacian for the radial part + The PSF is described by a product of an asymmetric Laplacian for the radial part and a symmetric Laplacian in the polar direction. + Explicitly, the radial part is given by - .. math:: f_{R}(r, K) = \begin{cases}\frac{1}{S_{R}(K+K^{-1})}e^{-K\frac{r-L}{S_{R}}}, r\ge L\\ \frac{1}{S_{R}(K+K^{-1})}e^{\frac{r-L}{KS_{R}}}, r < L\end{cases} + .. math:: f_{R}(r, K) = \begin{cases}\frac{1}{S_{R}(K+K^{-1})}e^{-K\frac{r-r_0}{S_{R}}}, r\ge r_0\\ \frac{1}{S_{R}(K+K^{-1})}e^{\frac{r-r_0}{KS_{R}}}, r < r_0\end{cases} - and a symmetric Laplacian in polar direction + and the polar part is given by .. math:: f_{\Phi}(\phi) = \frac{1}{2S_\phi}e^{-|\frac{\phi-\phi_0}{S_\phi}|} + The parameters :math:`K`, :math:`S_{R}`, and :math:`S_{\phi}` are functions of the distance :math:`r` to the optical axis. + Their detailed description is provided in the attributes section. Attributes ---------- @@ -330,7 +332,7 @@ class ComaPSFModel(PSFModel): Used to calculate a PDF asymmetry parameter K of the asymmetric radial Laplacian of the PSF as a function of the distance r to the optical axis. - .. math:: K(r) = 1 - \text{asym}_0 \tanh(\text{asym}_1 r) - \text{asym}_2 r + .. math:: K(r) = 1 - c_0 \tanh(c_1 r) - c_2 r radial_scale_params : list Describes the dependency of the radial scale on the distance to the center of the camera. From 3a9248644ef524a2d48ddb7b1c0082cbaaa36132 Mon Sep 17 00:00:00 2001 From: "mykhailo.dalchenko" Date: Thu, 6 Feb 2025 14:08:11 +0100 Subject: [PATCH 33/33] Fix SonarQube-detected code smells --- src/ctapipe/instrument/optics.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/ctapipe/instrument/optics.py b/src/ctapipe/instrument/optics.py index d14783dc33e..7ddbd744803 100644 --- a/src/ctapipe/instrument/optics.py +++ b/src/ctapipe/instrument/optics.py @@ -427,18 +427,18 @@ def pdf(self, x, y, x0, y0) -> np.ndarray: @validate("asymmetry_params") def _check_asymmetry_params(self, proposal): - if not len(proposal["value"]) == 3: + if len(proposal["value"]) != 3: raise TraitError("asymmetry_params needs to have length 3") return proposal["value"] @validate("radial_scale_params") def _check_radial_scale_params(self, proposal): - if not len(proposal["value"]) == 4: + if len(proposal["value"]) != 4: raise TraitError("radial_scale_params needs to have length 4") return proposal["value"] @validate("phi_scale_params") def _check_phi_scale_params(self, proposal): - if not len(proposal["value"]) == 3: + if len(proposal["value"]) != 3: raise TraitError("phi_scale_params needs to have length 3") return proposal["value"]