diff --git a/ctapipe/image/tests/test_toy.py b/ctapipe/image/tests/test_toy.py index 85f762421af..3060f0f9fdf 100644 --- a/ctapipe/image/tests/test_toy.py +++ b/ctapipe/image/tests/test_toy.py @@ -25,7 +25,7 @@ def test_intensity(seed, monkeypatch, prod5_lst): width = 0.05 * u.m length = 0.15 * u.m intensity = 200 - psi = "30d" + psi = 30 * u.deg # make sure we set a fixed seed for this test even when testing the # API without giving the rng @@ -69,7 +69,7 @@ def test_skewed(prod5_lst): width = 0.05 * u.m length = 0.15 * u.m intensity = 50 - psi = "30d" + psi = 30 * u.deg skewness = 0.3 model = SkewedGaussian( @@ -94,7 +94,7 @@ def test_compare(prod5_lst): width = 0.05 * u.m length = 0.15 * u.m intensity = 50 - psi = "30d" + psi = 30 * u.deg skewed = SkewedGaussian(x=x, y=y, width=width, length=length, psi=psi, skewness=0) normal = Gaussian(x=x, y=y, width=width, length=length, psi=psi) diff --git a/ctapipe/image/toymodel.py b/ctapipe/image/toymodel.py index bb9a9f5abf6..03598790748 100644 --- a/ctapipe/image/toymodel.py +++ b/ctapipe/image/toymodel.py @@ -18,15 +18,16 @@ >>> print(image.shape) (400,) """ -import numpy as np -from ctapipe.utils import linalg -from ctapipe.image.hillas import camera_to_shower_coordinates -import astropy.units as u -from astropy.coordinates import Angle -from scipy.stats import multivariate_normal, skewnorm, norm -from scipy.ndimage import convolve1d from abc import ABCMeta, abstractmethod + +import astropy.units as u +import numpy as np from numpy.random import default_rng +from scipy.ndimage import convolve1d +from scipy.stats import multivariate_normal, norm, skewnorm + +from ctapipe.image.hillas import camera_to_shower_coordinates +from ctapipe.utils import linalg __all__ = [ "WaveformModel", @@ -40,15 +41,13 @@ TOYMODEL_RNG = default_rng(0) -@u.quantity_input( - x=u.m, - y=u.m, - centroid_x=u.m, - centroid_y=u.m, - psi=u.deg, - time_gradient=u.ns / u.m, - time_intercept=u.ns, -) +def _obtain_image_no_units( + x, y, centroid_x, centroid_y, psi, time_gradient, time_intercept +): + longitudinal, _ = camera_to_shower_coordinates(x, y, centroid_x, centroid_y, psi) + return longitudinal * time_gradient + time_intercept + + def obtain_time_image(x, y, centroid_x, centroid_y, psi, time_gradient, time_intercept): """Create a pulse time image for a toymodel shower. Assumes the time development occurs only along the longitudinal (major) axis of the shower, and scales @@ -56,19 +55,19 @@ def obtain_time_image(x, y, centroid_x, centroid_y, psi, time_gradient, time_int Parameters ---------- - x : u.Quantity[length] + x : u.Quantity[angle] X camera coordinate to evaluate the time at. Usually the array of pixel X positions - y : u.Quantity[length] + y : u.Quantity[angle] Y camera coordinate to evaluate the time at. Usually the array of pixel Y positions - centroid_x : u.Quantity[length] + centroid_x : u.Quantity[angle] X camera coordinate for the centroid of the shower - centroid_y : u.Quantity[length] + centroid_y : u.Quantity[angle] Y camera coordinate for the centroid of the shower psi : convertible to `astropy.coordinates.Angle` rotation angle about the centroid (0=x-axis) - time_gradient : u.Quantity[time/length] + time_gradient : u.Quantity[time/angle] Rate at which the time changes with distance along the shower axis time_intercept : u.Quantity[time] Pulse time at the shower centroid @@ -79,11 +78,27 @@ def obtain_time_image(x, y, centroid_x, centroid_y, psi, time_gradient, time_int Pulse time in nanoseconds at (x, y) """ - longitudinal, _ = camera_to_shower_coordinates(x, y, centroid_x, centroid_y, psi) - longitudinal_m = longitudinal.to_value(u.m) - time_gradient_ns_m = time_gradient.to_value(u.ns / u.m) - time_intercept_ns = time_intercept.to_value(u.ns) - return longitudinal_m * time_gradient_ns_m + time_intercept_ns + # camera frame + if x.unit == u.m: + return _obtain_image_no_units( + x=x.to_value(u.m), + y=y.to_value(u.m), + centroid_x=centroid_x.to_value(u.m), + centroid_y=centroid_y.to_value(u.m), + psi=psi.to_value(u.rad), + time_gradient=time_gradient.to_value(u.ns / u.m), + time_intercept=time_intercept.to_value(u.ns), + ) + + return _obtain_image_no_units( + x=x.to_value(u.deg), + y=y.to_value(u.deg), + centroid_x=centroid_x.to_value(u.dg), + centroid_y=centroid_y.to_value(u.deg), + psi=psi.to_value(u.rad), + time_gradient=time_gradient.to_value(u.ns / u.deg), + time_intercept=time_intercept.to_value(u.ns), + ) class WaveformModel: @@ -183,11 +198,9 @@ def from_camera_readout(cls, readout, gain_channel=0): class ImageModel(metaclass=ABCMeta): - @u.quantity_input(x=u.m, y=u.m) @abstractmethod def pdf(self, x, y): - """Probability density function. - """ + """Probability density function.""" def generate_image(self, camera, intensity=50, nsb_level_pe=20, rng=None): """Generate a randomized DL1 shower image. @@ -244,7 +257,6 @@ def expected_signal(self, camera, intensity): class Gaussian(ImageModel): - @u.quantity_input(x=u.m, y=u.m, length=u.m, width=u.m) def __init__(self, x, y, length, width, psi): """Create 2D Gaussian model for a shower image in a camera. @@ -264,32 +276,36 @@ def __init__(self, x, y, length, width, psi): a `scipy.stats` object """ + self.unit = x.unit self.x = x self.y = y self.width = width self.length = length self.psi = psi - @u.quantity_input(x=u.m, y=u.m) - def pdf(self, x, y): - """2d probability for photon electrons in the camera plane""" aligned_covariance = np.array( - [[self.length.to_value(u.m) ** 2, 0], [0, self.width.to_value(u.m) ** 2]] + [ + [self.length.to_value(self.unit) ** 2, 0], + [0, self.width.to_value(self.unit) ** 2], + ] ) # rotate by psi angle: C' = R C R+ rotation = linalg.rotation_matrix_2d(self.psi) rotated_covariance = rotation @ aligned_covariance @ rotation.T + self.dist = multivariate_normal( + mean=u.Quantity([self.x, self.y]).to_value(self.unit), + cov=rotated_covariance, + ) - return multivariate_normal( - mean=[self.x.to_value(u.m), self.y.to_value(u.m)], cov=rotated_covariance - ).pdf(np.column_stack([x.to_value(u.m), y.to_value(u.m)])) + def pdf(self, x, y): + """2d probability for photon electrons in the camera plane""" + X = np.column_stack([x.to_value(self.unit), y.to_value(self.unit)]) + return self.dist.pdf(X) class SkewedGaussian(ImageModel): - """A shower image that has a skewness along the major axis. - """ + """A shower image that has a skewness along the major axis.""" - @u.quantity_input(x=u.m, y=u.m, length=u.m, width=u.m) def __init__(self, x, y, length, width, psi, skewness): """Create 2D skewed Gaussian model for a shower image in a camera. Skewness is only applied along the main shower axis. @@ -313,6 +329,7 @@ def __init__(self, x, y, length, width, psi, skewness): a `scipy.stats` object """ + self.unit = x.unit self.x = x self.y = y self.width = width @@ -320,6 +337,12 @@ def __init__(self, x, y, length, width, psi, skewness): self.psi = psi self.skewness = skewness + a, loc, scale = self._moments_to_parameters() + self.long_dist = skewnorm(a=a, loc=loc, scale=scale) + self.trans_dist = norm(loc=0, scale=self.width.to_value(self.unit)) + self.rotation = linalg.rotation_matrix_2d(-self.psi) + self.mu = u.Quantity([self.x, self.y]).to_value(u.m) + def _moments_to_parameters(self): """Returns loc and scale from mean, std and skewnewss.""" # see https://en.wikipedia.org/wiki/Skew_normal_distribution#Estimation @@ -327,26 +350,17 @@ def _moments_to_parameters(self): delta = np.sign(self.skewness) * np.sqrt( (np.pi / 2 * skew23) / (skew23 + (0.5 * (4 - np.pi)) ** (2 / 3)) ) - a = delta / np.sqrt(1 - delta ** 2) - scale = self.length.to_value(u.m) / np.sqrt(1 - 2 * delta ** 2 / np.pi) + a = delta / np.sqrt(1 - delta**2) + scale = self.length.to_value(self.unit) / np.sqrt(1 - 2 * delta**2 / np.pi) loc = -scale * delta * np.sqrt(2 / np.pi) return a, loc, scale - @u.quantity_input(x=u.m, y=u.m) def pdf(self, x, y): """2d probability for photon electrons in the camera plane.""" - mu = u.Quantity([self.x, self.y]).to_value(u.m) - - rotation = linalg.rotation_matrix_2d(-Angle(self.psi)) - pos = np.column_stack([x.to_value(u.m), y.to_value(u.m)]) - long, trans = rotation @ (pos - mu).T - - trans_pdf = norm(loc=0, scale=self.width.to_value(u.m)).pdf(trans) - - a, loc, scale = self._moments_to_parameters() - - return trans_pdf * skewnorm(a=a, loc=loc, scale=scale).pdf(long) + pos = np.column_stack([x.to_value(self.unit), y.to_value(self.unit)]) + long, trans = self.rotation @ (pos - self.mu).T + return self.trans_dist.pdf(trans) * self.long_dist.pdf(long) class RingGaussian(ImageModel): @@ -355,17 +369,17 @@ class RingGaussian(ImageModel): """ - @u.quantity_input(x=u.m, y=u.m, radius=u.m, sigma=u.m) def __init__(self, x, y, radius, sigma): + self.unit = x.unit self.x = x self.y = y self.sigma = sigma self.radius = radius + self.dist = norm( + self.radius.to_value(self.unit), self.sigma.to_value(self.unit) + ) - @u.quantity_input(x=u.m, y=u.m) def pdf(self, x, y): """2d probability for photon electrons in the camera plane.""" - r = np.sqrt((x - self.x) ** 2 + (y - self.y) ** 2) - - return norm(self.radius.to_value(u.m), self.sigma.to_value(u.m)).pdf(r) + return self.dist.pdf(r)