Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Chunk interpolation to select calibration data #2634

Open
wants to merge 16 commits into
base: main
Choose a base branch
from
Open

Conversation

ctoennis
Copy link

I will need a method to select calibration data for the strar tracker. I made some slides to decribe how it is supposed to work:
https://docs.google.com/presentation/d/1oxIcYSQvGnU7IQYy3fGdcv0qXiLpvaXR9YtmnDesj4Y/edit?usp=sharing

This comment has been minimized.

This comment has been minimized.

2 similar comments

This comment has been minimized.

Copy link

Passed

Analysis Details

1 Issue

  • Bug 0 Bugs
  • Vulnerability 0 Vulnerabilities
  • Code Smell 1 Code Smell

Coverage and Duplications

  • Coverage 98.00% Coverage (94.30% Estimated after merge)
  • Duplications 0.00% Duplicated Code (0.70% Estimated after merge)

Project ID: cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

View in SonarQube

@ctoennis
Copy link
Author

@maxnoe @kosack Can you have another look if there is something else to be changed?

@mexanick
Copy link
Contributor

@kosack @maxnoe this PR is needed to complete the pointing calibration (for the variance calibration application), can we advance it?


return np.nan

self._interpolators[tel_id][column] = interpolate_chunk
Copy link
Member

Choose a reason for hiding this comment

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

I think it's a bit weird here to store closures in a dict instead of just calling a method with tel_id and column that then checks if the necessary table exists.

Copy link
Author

Choose a reason for hiding this comment

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

ok, i'll change it to a method

Copy link
Contributor

Choose a reason for hiding this comment

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

The current implementation follows the PointingInterpolator approach, but of course can be changed.

Copy link
Author

Choose a reason for hiding this comment

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

so, _check_interpolators in the parent class assumes the current implementation. If i change it to what you describe i will need to make changes to the parent class or forgo this check for the ChunkInterpolator

Copy link
Author

Choose a reason for hiding this comment

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

Maybe keeping it like this is better.

Copy link
Member

@maxnoe maxnoe Nov 21, 2024

Choose a reason for hiding this comment

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

👍 Looks good, the __call__ method is missing on the interface though

I am not sure if putting the private methods you listed here to the base class will work out in the end, but you'll see how much code can be shared between the chunk interpolator and the linear interpolator.

Copy link
Author

Choose a reason for hiding this comment

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

ok, i'll work on the implementation now.

Copy link
Contributor

Choose a reason for hiding this comment

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

Shall we take an opportunity and stop using interp1d? It is legacy, andscipy suggests this: https://docs.scipy.org/doc/scipy/tutorial/interpolate/1D.html#tutorial-interpolate-1dsection

Also, if we call a base class LinearInterpolator, I guess we should enforce linearity (in the example above, e.g. cubic spline can be applied)

Copy link
Member

Choose a reason for hiding this comment

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

There is no replacement for linear 1d interpolation for interp1d. It is legacy, but not deprecated.

See scipy/scipy#21258

Copy link
Contributor

Choose a reason for hiding this comment

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

I also got tricked by "may be removed in the future"

@ctoennis
Copy link
Author

I am a bit stuck here with one of the tests. test_hdf5 is failing in pytest begause the data from the file is not loaded correctly. When i look in the test what x and y values the interpolators have i get some wrong values. However if i try to do the same outside of pytest it works. I used this code to test by myself:

import astropy.units as u
import numpy as np
import tables
from astropy.table import Table
from astropy.time import Time

from functools import partial
from ctapipe.core import Component, traits

from ctapipe.monitoring.interpolation import PointingInterpolator

from ctapipe.io import write_table

t0 = Time("2022-01-01T00:00:00")

table = Table(
    {"time": t0 + np.arange(0.0, 10.1, 2.0) * u.s, "azimuth": np.radians(np.linspace(0.0, 10.0, 6)) * u.rad, "altitude": np.radians(>)

path = "pointing.h5"

write_table(table, path, "/dl0/monitoring/telescope/pointing/tel_001")

with tables.open_file(path) as h5file:
    interpolator = PointingInterpolator(h5file)
    t = t0 + 1 * u.s
    alt, az = interpolator(tel_id=1, time=t)
    print(interpolator._interpolators[1]["alt"].y,interpolator._interpolators[1]["alt"].x)
    print(alt,az)

Has anyone an idea what is wrong here?

@@ -57,7 +173,7 @@ def test_invalid_input():
wrong_unit = Table(
{
"time": Time(1.7e9 + np.arange(3), format="unix"),
"azimuth": [1, 2, 3] * u.deg,
"azimuth": np.radians([1, 2, 3]) * u.rad,
Copy link
Member

Choose a reason for hiding this comment

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

Why did you change these tests here?

Copy link
Member

Choose a reason for hiding this comment

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

The test is explicitly to test the unit conversion, so please don't change to rad

Copy link
Author

@ctoennis ctoennis Nov 22, 2024

Choose a reason for hiding this comment

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

The PointingInterpolator is supposed to use radians for input:

157         expected_units = {"azimuth": u.rad, "altitude": u.rad}

Copy link
Contributor

Choose a reason for hiding this comment

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

right, and this test checks that correct exception is raised in case the units are not matching. Please revert back

@@ -72,8 +188,8 @@ def test_hdf5(tmp_path):
table = Table(
{
"time": t0 + np.arange(0.0, 10.1, 2.0) * u.s,
"azimuth": np.linspace(0.0, 10.0, 6) * u.deg,
Copy link
Member

Choose a reason for hiding this comment

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

same here

@@ -82,8 +198,8 @@ def test_hdf5(tmp_path):
with tables.open_file(path) as h5file:
interpolator = PointingInterpolator(h5file)
alt, az = interpolator(tel_id=1, time=t0 + 1 * u.s)
assert u.isclose(alt, 69 * u.deg)
assert u.isclose(az, 1 * u.deg)
assert u.isclose(alt, np.radians(69) * u.rad)
Copy link
Member

Choose a reason for hiding this comment

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

same here

@mexanick
Copy link
Contributor

@ctoennis the problem you're experiencing with the test is because you have class attributes in the base class that are shared between the class instances, including the _interpolators.

@ctoennis
Copy link
Author

@mexanick Thanks, i got the issue fixed.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants