Skip to content

Commit

Permalink
removed Laplace model
Browse files Browse the repository at this point in the history
  • Loading branch information
@clara-escanuela
clara-escanuela committed Sep 14, 2023
1 parent fd22cf0 commit dec2e6d
Show file tree
Hide file tree
Showing 3 changed files with 1 addition and 178 deletions.
3 changes: 1 addition & 2 deletions ctapipe/image/ellipsoid.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@
mean_poisson_likelihood_gaussian,
neg_log_likelihood_approx,
)
from ctapipe.image.toymodel import SkewedGaussian, SkewedGaussianLaplace
from ctapipe.image.toymodel import SkewedGaussian

from ..containers import CameraImageFitParametersContainer, ImageFitParametersContainer

Expand Down Expand Up @@ -227,7 +227,6 @@ def image_fit_parameters(
pdf_dict = {
PDFType.gaussian: SkewedGaussian,
PDFType.skewed: SkewedGaussian,
PDFType.laplace: SkewedGaussianLaplace,
}

unit = geom.pix_x.unit
Expand Down
90 changes: 0 additions & 90 deletions ctapipe/image/tests/test_ellipsoid.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -209,18 +209,6 @@ def test_truncated(prod5_lst):
assert conc_fit > conc_hillas
assert conc_fit > 0.4

# Laplace
result = image_fit_parameters(
geom, image, n_row=2, cleaned_mask=clean_mask, pdf=PDFType("laplace")
)

assert result.length.value > hillas.length.value

conc_fit = concentration_parameters(geom, cleaned_image, result).core

assert conc_fit > conc_hillas
assert conc_fit > 0.4


def test_percentage(prod5_lst):
geom = prod5_lst.camera.geometry
Expand Down Expand Up @@ -274,9 +262,6 @@ def test_with_toy_mst_tel(prod5_mst_flashcam):
model_skewed = toymodel.SkewedGaussian(
x=x, y=y, width=width, length=length, psi=psi, skewness=0.5
)
model_laplace = toymodel.SkewedGaussianLaplace(
x=x, y=y, width=width, length=length, psi=psi, skewness=0.5
)

image, signal, noise = model_gaussian.generate_image(
geom, intensity=intensity, nsb_level_pe=0, rng=rng
Expand Down Expand Up @@ -320,24 +305,6 @@ def test_with_toy_mst_tel(prod5_mst_flashcam):
result.psi.to_value(u.deg) - psi.deg
) == approx(180.0, abs=2)

image, signal, noise = model_laplace.generate_image(
geom, intensity=intensity, nsb_level_pe=0, rng=rng
)
clean_mask = np.array(signal) > 0
result = image_fit_parameters(
geom, signal, n_row=0, cleaned_mask=clean_mask, pdf=PDFType("laplace")
)

if result.is_valid or result.is_accurate:
assert u.isclose(result.x, x, rtol=0.1)
assert u.isclose(result.y, y, rtol=0.1)

assert u.isclose(result.width, width, rtol=0.15)
assert u.isclose(result.length, length, rtol=0.1)
assert (result.psi.to_value(u.deg) == approx(psi.deg, abs=2)) or abs(
result.psi.to_value(u.deg) - psi.deg
) == approx(180.0, abs=2)


def test_with_toy_lst(prod5_lst):
rng = np.random.default_rng(42)
Expand All @@ -362,9 +329,6 @@ def test_with_toy_lst(prod5_lst):
model_skewed = toymodel.SkewedGaussian(
x=x, y=y, width=width, length=length, psi=psi, skewness=0.5
)
model_laplace = toymodel.SkewedGaussianLaplace(
x=x, y=y, width=width, length=length, psi=psi, skewness=0.5
)

image, signal, noise = model_gaussian.generate_image(
geom, intensity=intensity, nsb_level_pe=0, rng=rng
Expand Down Expand Up @@ -408,24 +372,6 @@ def test_with_toy_lst(prod5_lst):
result.psi.to_value(u.deg) - psi.deg
) == approx(180.0, abs=2)

image, signal, noise = model_laplace.generate_image(
geom, intensity=intensity, nsb_level_pe=0, rng=rng
)
clean_mask = np.array(signal) > 0
result = image_fit_parameters(
geom, signal, n_row=0, cleaned_mask=clean_mask, pdf=PDFType("laplace")
)

if result.is_valid or result.is_accurate:
assert u.isclose(result.x, x, rtol=0.1)
assert u.isclose(result.y, y, rtol=0.1)

assert u.isclose(result.width, width, rtol=0.15)
assert u.isclose(result.length, length, rtol=0.1)
assert (result.psi.to_value(u.deg) == approx(psi.deg, abs=2)) or abs(
result.psi.to_value(u.deg) - psi.deg
) == approx(180.0, abs=2)


def test_skewness(prod5_lst):
rng = np.random.default_rng(42)
Expand All @@ -449,9 +395,6 @@ def test_skewness(prod5_lst):
model_skewed = toymodel.SkewedGaussian(
x=x, y=y, width=width, length=length, psi=psi, skewness=skew
)
model_laplace = toymodel.SkewedGaussianLaplace(
x=x, y=y, width=width, length=length, psi=psi, skewness=skew
)

image, signal, noise = model_skewed.generate_image(
geom, intensity=intensity, nsb_level_pe=0, rng=rng
Expand Down Expand Up @@ -482,34 +425,6 @@ def test_skewness(prod5_lst):

assert signal.sum() == result.intensity

image, signal, noise = model_laplace.generate_image(
geom, intensity=intensity, nsb_level_pe=0, rng=rng
)
clean_mask = np.array(signal) > 0
result = image_fit_parameters(
geom, signal, n_row=0, cleaned_mask=clean_mask, pdf=PDFType("laplace")
)

if result.is_valid or result.is_accurate:
assert u.isclose(result.x, x, rtol=0.1)
assert u.isclose(result.y, y, rtol=0.1)

assert u.isclose(result.width, width, rtol=0.1)
assert u.isclose(result.length, length, rtol=0.15)

psi_same = result.psi.to_value(u.deg) == approx(psi.deg, abs=1)
psi_opposite = abs(result.psi.to_value(u.deg) - psi.deg) == approx(
180.0, abs=1
)
assert psi_same or psi_opposite

if psi_same:
assert result.skewness == approx(skew, abs=0.3)
else:
assert result.skewness == approx(-skew, abs=0.3)

assert signal.sum() == result.intensity


def test_gaussian_skewness(prod5_lst):
rng = np.random.default_rng(42)
Expand Down Expand Up @@ -554,11 +469,6 @@ def test_single_pixel():

clean_mask = np.array(image) > 0

with pytest.raises(ImageFitParameterizationError):
image_fit_parameters(
geom, image, n_row=2, cleaned_mask=clean_mask, pdf=PDFType("laplace")
)

with pytest.raises(ImageFitParameterizationError):
image_fit_parameters(
geom, image, n_row=2, cleaned_mask=clean_mask, pdf=PDFType("skewed")
Expand Down
86 changes: 0 additions & 86 deletions ctapipe/image/toymodel.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,6 @@
"WaveformModel",
"Gaussian",
"SkewedGaussian",
"SkewedGaussianLaplace",
"ImageModel",
"obtain_time_image",
]
Expand Down Expand Up @@ -404,88 +403,3 @@ 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 self.dist.pdf(r)


class SkewedGaussianLaplace(ImageModel):
"""A shower image that has a skewness along the major axis and follows the Laplace distribution along the
transverse axis"""

@u.quantity_input
def __init__(self, x, y, length, width, psi, skewness, amplitude=None):
"""Create 2D function with a Skewed Gaussian in the longitudinal direction
and a Laplace function modelling the transverse direction of the shower.
See https://en.wikipedia.org/wiki/Skew_normal_distribution ,
https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.skewnorm.html , and
https://en.wikipedia.org/wiki/Laplace_distribution for details.
Parameters
----------
x : u.Quantity[cx]
position of the centroid in x-axis of the shower in camera coordinates
y : u.Quantity[cy]
position of the centroid in y-axis of the shower in camera coordinates
width: u.Quantity[length]
width of shower (minor axis)
length: u.Quantity[length]
length of shower (major axis)
psi : convertable to `astropy.coordinates.Angle`
rotation angle about the centroid (0=x-axis)
skewness: float
skewness of the shower in longitudinal direction
amplitude : float
normalization amplitude
Returns
-------
model : ndarray
Skewed Gaussian * Laplace distribution
"""
self.x = x
self.y = y
self.width = width
self.length = length
self.psi = psi
self.skewness = skewness
self.amplitude = amplitude
self.unit = self.x.unit

def _moments_to_parameters(self):
"""Returns loc and scale from mean, std and skewness."""
# see https://en.wikipedia.org/wiki/Skew_normal_distribution#Estimation
skew23 = np.abs(self.skewness) ** (2 / 3)
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(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
def pdf(self, x, y):
"""2d probability for photon electrons in the camera plane."""
mu = u.Quantity([self.x, self.y]).to_value(self.unit)

rotation = linalg.rotation_matrix_2d(-Angle(self.psi))
pos = np.column_stack([x.to_value(self.unit), y.to_value(self.unit)])
long, trans = rotation @ (pos - mu).T

a, loc, scale = self._moments_to_parameters()

if self.amplitude is None:
self.amplitude = 1 / (
np.sqrt(2 * np.pi)
* np.sqrt(2)
* scale
* (self.width.to_value(self.unit))
)

trans_pdf = np.exp(
-np.sqrt(trans**2) * np.sqrt(2) / self.width.to_value(self.unit)
)
skew_pdf = np.exp(-1 / 2 * ((long - loc) / scale) ** 2) * (
1 + scipy.special.erf(a / np.sqrt(2) * (long - loc) / scale)
)
return self.amplitude * trans_pdf * skew_pdf

0 comments on commit dec2e6d

Please sign in to comment.