Add algorithm to compute image parameters via fitting the light distribution #2275

clara-escanuela · 2023-03-01T17:39:09Z

No description provided.

clara-escanuela · 2023-03-01T17:43:07Z

This is my first step towards this new image-fitting method. Although I would like to get suggestions about:

How to make ellipsoid.py more general and include an option for PDFs
I want to now integrate this into HillasReconstructor. How would you proceed?

clara-escanuela · 2023-03-02T07:20:49Z

I also think that the parameter boundaries should be calculated in ellipsoid.py and not being optional. So maybe I could include another function in ellipsoid.py that finds the limits of each parameter

maxnoe · 2023-03-02T09:23:51Z

Let's do one thing after another. Adding new image parameters is already a large enough pull request, let's not mix it with changing the HillasReconstructor

maxnoe · 2023-03-08T13:07:18Z

ctapipe/containers.py

@@ -243,6 +243,86 @@ class HillasParametersContainer(BaseHillasParametersContainer):
    width_uncertainty = Field(nan * u.deg, "uncertainty of width", unit=u.deg)
    psi = Field(nan * u.deg, "rotation angle of ellipse", unit=u.deg)

+class BaseImageFitParametersContainer(Container):


Since we are planning to drop the camera frame stuff altogether, for new stuff I wouldn't at it. So just having the TelescopeFrame parameter container is fine.

maxnoe · 2023-03-08T13:09:10Z

ctapipe/image/toymodel.py

@@ -369,3 +371,72 @@ def pdf(self, x, y):
        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)
+
+
+class SkewedCauchy(ImageModel):


It's not a SkewedCauchy, is it? It's a SkewedGaussian for the longitudinal and a normal cauchy for trans, right?

Yes. I will clarify

kosack

not a complete review, but a few initial comments

ctapipe/image/ellipsoid.py

kosack · 2023-04-17T15:10:50Z

ctapipe/image/ellipsoid.py

+        Camera geometry
+    cleaned_image: ndarray
+        Charge for each pixel, cleaning mask should be applied
+    pdf: str


Might be useful to use an Enum here and elsewhere for pdf (for speed and correctness and to make it easier to include in a config file later).

kosack · 2023-04-17T15:11:21Z

ctapipe/image/ellipsoid.py

+    x0 : dict
+       seeds of the fit
+    pdf: str
+       name of the PDF, options = "Gaussian", "Laplace", "Skewed"


same here, use an Enum

kosack · 2023-04-17T15:11:52Z

ctapipe/image/ellipsoid.py

+        0.99, x0["skewness"] + 0.3
+    )
+
+    if pdf == "Gaussian":


this then becomes, e.g. PDFType.Gaussian

Is that what you were asking for?

clara-escanuela · 2023-05-06T14:43:14Z

Here are some examples of the performance of the fit when applying a 2D Gaussian:

clara-escanuela · 2023-05-09T06:39:25Z

Maybe this is also interesting. For these tests, I am using an array of 14 MSTs and the HillasReconstructor, telescope multiplicity cut >= 5 triggered telescopes. Although these plots were made with a slightly older version of this code.

maxnoe · 2023-05-09T12:48:29Z

@clara-escanuela Could you also include please the standard hillas parameters based approach without applying a centroid cut?

maxnoe · 2023-05-09T12:50:40Z

As discussed via slack, I don't think that the centroid is a good criterion, the leakage_* are better parameters.

Thinking about it: since you get a predicted image from the model, and the model pdf is normalized, you should very easily be able to compute the expected containment from the image fit as (predicted_image.sum() / normalization) , which might be an interesting parameter.

gaia-verna · 2023-05-09T12:52:19Z

Maybe this is also interesting. For these tests, I am using an array of 14 MSTs and the HillasReconstructor, telescope multiplicity cut >= 5 triggered telescopes. Although these plots were made with a slightly older version of this code.

Just for clarity, you considered:

the MST sub-array of CTA-South
MC gammas from Prod5
an on-axis observation @ 20° ZA

and the 2D fit is applied to all the shower images, right?

maxnoe · 2023-05-09T12:54:37Z

And when you do the cut on the centroid: is it the "new" centroid from the image fit for the fit or the same hillas centroid for both?

clara-escanuela · 2023-05-10T08:25:36Z

And when you do the cut on the centroid: is it the "new" centroid from the image fit for the fit or the same hillas centroid for both?

The centroid from Hillas in both cases. I wanted to compare the exact same events.

clara-escanuela · 2023-05-10T08:32:55Z

@clara-escanuela Could you also include please the standard hillas parameters based approach without applying a centroid cut?

Here I am plotting Hillas when I do not apply any cut on the centroid distance (light blue), Hillas when I apply a cut on centroid (need to change definition) at 0.8m (dark blue), and 2D Fit with no cut on centroid. I fit all images and not only truncated. It would be interesting to see how this plot looks like when fitting only truncated images to check what it is the actual impact on only truncated images.

maxnoe · 2023-09-08T14:03:12Z

ctapipe/containers.py

+    goodness_of_fit = Field(
+        nan, "measure of goodness of fit, mean likelihood subtracted to the likelihood"
+    )
+    n_pix_fit = Field(nan, "number of pixels used in the fit")


Please avoid such unneeded abbreviations.

maxnoe · 2023-09-08T14:04:12Z

ctapipe/containers.py

@@ -216,6 +217,37 @@ class CameraHillasParametersContainer(BaseHillasParametersContainer):
    psi = Field(nan * u.deg, "rotation angle of ellipse", unit=u.deg)


+class ImageFitParametersContainer(CameraHillasParametersContainer):


We added support for the toy models in telescope frame now, so it should be possible (and should be the default) that this happens also in telescope frame (i.e. using fov_lon / fov_lat in degrees instead of x, y in m)

ctapipe/image/ellipsoid.py

maxnoe · 2023-09-08T14:05:26Z

ctapipe/image/ellipsoid.py

+    Seed parameters of the fit.
+    """
+    unit = geometry.pix_x.unit
+    hillas = hillas_parameters(geometry, image)  # compute Hillas parameters


Why compute hillas parameters again instead of just taking them as input?

maxnoe · 2023-09-08T14:05:51Z

ctapipe/image/ellipsoid.py

+    return initial_guess
+
+
+def dilation(n, cleaned_mask, geometry):


Pretty sure this already exists in ctapipe.image.cleaning, no need to redefine

actually, we only have dilate which does a single dilation. THe loop for multiple is not there, but this code should in any case go into cleaning.py (maybe call it dilate_multiple)

or use functools.reduce to do the dilate(dilate(dilate()) call chaining

maxnoe · 2023-09-08T14:06:10Z

ctapipe/image/ellipsoid.py

+    -------
+    list of boundaries
+    """
+    hillas = hillas_parameters(geometry, cleaned_image)


Again a recomputation of parameters...

maxnoe · 2023-09-08T14:06:36Z

ctapipe/image/ellipsoid.py

+    hillas = hillas_parameters(geometry, cleaned_image)
+
+    unit = geometry.pix_x.unit
+    camera_radius = geometry.guess_radius()


Is this memoized? If not it might be expensive to call this in the event loop

it is current not memoized, so it will be quite slow to call this frequently. I didn't add a @lazyproperty to this call, since this is only one way to compute the edge of the camera (not sure it works well in all cases). It was just to give a quick way of computing a reasonable value. Ideally we would have a geom.radius property that is supplied externally like the pixel positions or telescope focal length, etc, and only guessed if none is provided.

In the mean-time, we could just add:

@lazyproperty def radius(self): return self.guess_radius()

maxnoe · 2023-09-08T14:06:58Z

ctapipe/image/ellipsoid.py

+    camera_radius = geometry.guess_radius()
+
+    cogx_min, cogx_max = np.sign(hillas.x) * min(
+        np.abs(hillas.x - u.Quantity(0.2, unit)), camera_radius


magic number 0.2?

maxnoe · 2023-09-08T14:07:58Z

ctapipe/image/ellipsoid.py

+
+from ..containers import ImageFitParametersContainer
+
+PED_TABLE = {


Where do these come from? Why these keys?

These seem to be used in the likelihood to provide the pedestal, so this is not the right way to load them: it will break if you use different simulation productions, or even sims with different cameras (there is no "SST-Camera" for example, it's call "SSTCam"), and clearly won't work for real data. Note that even between prod5 and prod6, the pedestal values are not the same

However, I would have two questions:

since the waveforms and images derived from them are already pedestal-subtracted, do you really need the pedestal? Maybe the likelihood function needs modification for this case?

if you need the pedestal, you should get it from the calibration container (but have to check if it is actually filled for simulations)

I do not think the pedestal width is available at the moment. I used those values that I found in Dan's code as he also had to load the pedestal width like that.

maxnoe · 2023-09-08T14:10:14Z

ctapipe/image/toymodel.py

@@ -369,3 +399,88 @@ def pdf(self, x, y):
        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)
+
+
+class SkewedGaussianLaplace(ImageModel):


Is there a physical motivation for this? Does it fit better?

We thought this could better estimate the width. But, after some tests, I do not think it works better than the 2D Skewed Gaussian so I think I could remove it

kosack · 2023-09-08T15:02:47Z

ctapipe/image/ellipsoid.py

+    "DUMMY": 0,
+    "testcam": 0,
+}
+SPE_WIDTH = 0.5


I guess this is also not the correct value, and you should get it from the appropriate Container.

I also think this is not available. At least I could not find the SPE width

maxnoe

I left some more inline comments.

One main one, that also should simplify the implementation a lot is that we shouldn't add new features using the CameraFrame, since we plan to remove image parameters in CameraFrame anyways very soon.

Second comment is that there are seemingly multiple PDFType choices allowed, but I see several parts of the code that should only work with the skewed version.

maxnoe · 2023-11-23T16:10:57Z

ctapipe/containers.py

+    goodness_of_fit = Field(
+        nan, "measure of goodness of fit, mean likelihood subtracted to the likelihood"
+    )
+    n_pixels = Field(nan, "number of pixels used in the fit")


Positive integer values should have defaults of the correct type, e.g. -1

maxnoe · 2023-11-23T16:11:29Z

ctapipe/containers.py

+    )
+    n_pixels = Field(nan, "number of pixels used in the fit")
+    free_parameters = Field(nan, "number of free parameters")
+    is_valid = Field(False, "True if the fit is valid")


What's the difference between is_valid and is_accurate?

maxnoe · 2023-11-23T16:12:43Z

ctapipe/image/ellipsoid.py

+
+def create_initial_guess(geometry, hillas, size):
+    """
+


The comment below was only for the "Parameters" and "Returns" sections, not for the initial line. So this empty line should be removed.

maxnoe · 2023-11-23T16:13:13Z

ctapipe/image/ellipsoid.py

+
+    initial_guess["skewness"] = hillas.skewness
+
+    if (hillas.width.to_value(unit) == 0) or (hillas.length.to_value(unit) == 0):


For checking against 0, the unit doesn't matter, so just doing .value != 0 is fine and faster.

maxnoe · 2023-11-23T16:13:46Z

ctapipe/image/ellipsoid.py

+
+def boundaries(geometry, image, dilated_mask, hillas, pdf):
+    """
+


No empty line here

ctapipe/image/ellipsoid.py

maxnoe · 2023-11-23T16:16:55Z

ctapipe/image/ellipsoid.py

+        like_array[like_array > 0]
+        - mean_poisson_likelihood_gaussian(
+            pdf_dict[pdf](
+                pars[0] * unit,


This won't work for different kinds of PDFs that have different numbers of parameters, does it?

ctapipe/image/ellipsoid.py

codecov · 2023-12-15T15:06:14Z

Codecov Report

All modified and coverable lines are covered by tests ✅

❗ No coverage uploaded for pull request base (main@3c5b4f3). Click here to learn what that means.

Additional details and impacted files

@@           Coverage Diff           @@
##             main    #2275   +/-   ##
=======================================
  Coverage        ?   92.54%           
=======================================
  Files           ?      236           
  Lines           ?    20288           
  Branches        ?        0           
=======================================
  Hits            ?    18776           
  Misses          ?     1512           
  Partials        ?        0

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

clara-escanuela linked an issue Mar 1, 2023 that may be closed by this pull request

Alternative image parametrization based on gaussian fitting #2270

Open

maxnoe reviewed Mar 8, 2023

View reviewed changes

clara-escanuela force-pushed the image-fitting branch from 2e09fb1 to 30abb64 Compare March 16, 2023 10:32

clara-escanuela force-pushed the image-fitting branch from 2b36b13 to 90b5045 Compare April 3, 2023 15:05

clara-escanuela force-pushed the image-fitting branch from 232377e to c4f0c76 Compare April 15, 2023 17:57

kosack requested changes Apr 17, 2023

View reviewed changes

clara-escanuela marked this pull request as ready for review May 5, 2023 17:46

maxnoe reviewed Sep 8, 2023

View reviewed changes

ctapipe/image/ellipsoid.py Show resolved Hide resolved

maxnoe reviewed Sep 8, 2023

View reviewed changes

kosack reviewed Sep 8, 2023

View reviewed changes

clara-escanuela added 17 commits September 13, 2023 10:38

add docs/changes/2275.*.rst file

02774fe

added more tests and made some changes in the method

20b0b0d

solved issues with tests

36eae9a

solved issues with tests

2f6e677

added simple test on boundaries

242d90a

added one test

35648fa

added a test on the skewness of the Skewed and Cauchy

9f6897f

replaced Cauchy function by Laplace, seems to be a better match

03d9877

added PDFType instances and removed not needed line in dilation

6b5d8ed

solved problem with skewness boundaries

22e1247

removed no needed containers and repeated lines of code

e9a9d4a

reduce code duplications

198832d

added one test and adapted code

6b1b70d

improved documentation

0a239ca

added goodness of fit

0043aad

removed unnecesary function and repetitions

c9bef50

removed unnecesary function and repetitions

329a772

clara-escanuela force-pushed the image-fitting branch from 269cd82 to 329a772 Compare September 13, 2023 08:48

clara-escanuela added 8 commits September 13, 2023 10:51

fixed single pixel unit test

09beba4

removed dilation function

e07d947

fixed unit tests for toymodel

790ab1f

added telescope frame

f81e7d5

changed variable names

fd22cf0

removed Laplace model

dec2e6d

test changes

6a08f04

more unit tests

879c18f

maxnoe changed the title ~~Image-fitting algorithm~~ Add algorithm to compute image parameters via fitting the light distribution Oct 22, 2023

maxnoe requested changes Nov 23, 2023

View reviewed changes

removed camera frame

12d3738

		@@ -216,6 +217,37 @@ class CameraHillasParametersContainer(BaseHillasParametersContainer):
		psi = Field(nan * u.deg, "rotation angle of ellipse", unit=u.deg)


		class ImageFitParametersContainer(CameraHillasParametersContainer):

		return initial_guess


		def dilation(n, cleaned_mask, geometry):


		from ..containers import ImageFitParametersContainer

		PED_TABLE = {


		initial_guess["skewness"] = hillas.skewness

		if (hillas.width.to_value(unit) == 0) or (hillas.length.to_value(unit) == 0):


		def boundaries(geometry, image, dilated_mask, hillas, pdf):
		"""

Add algorithm to compute image parameters via fitting the light distribution #2275

Are you sure you want to change the base?

Add algorithm to compute image parameters via fitting the light distribution #2275

Conversation

clara-escanuela commented Mar 1, 2023

clara-escanuela commented Mar 1, 2023

clara-escanuela commented Mar 2, 2023

maxnoe commented Mar 2, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kosack left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

clara-escanuela commented May 6, 2023

clara-escanuela commented May 9, 2023 • edited Loading

maxnoe commented May 9, 2023

maxnoe commented May 9, 2023

gaia-verna commented May 9, 2023

maxnoe commented May 9, 2023

clara-escanuela commented May 10, 2023 • edited Loading

clara-escanuela commented May 10, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kosack Sep 8, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kosack Sep 8, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

maxnoe left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

codecov bot commented Dec 15, 2023

Codecov Report

clara-escanuela commented May 9, 2023 •

edited

Loading

clara-escanuela commented May 10, 2023 •

edited

Loading

kosack Sep 8, 2023 •

edited

Loading

kosack Sep 8, 2023 •

edited

Loading