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.

Sign up for GitHub

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

VaryingCelestialTransform fix and unit tests #285

Merged
merged 6 commits into from
Aug 7, 2023
Merged

VaryingCelestialTransform fix and unit tests #285

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
7a9e33b
Fixed coordinate reversal in VaryingCelestialTransformSlit2D._map_tra…
svensk-pojke Aug 3, 2023
1b9c515
Added changelog fragment.
svensk-pojke Aug 3, 2023
0eefab6
Merge branch 'main' into time_varying_transform_fix
svensk-pojke Aug 3, 2023
3c03c2e
Apply suggestions from code review
svensk-pojke Aug 4, 2023
b53f41d
Refactored vct unit tests.
svensk-pojke Aug 4, 2023
46e6977
Fixed some inaccurate comments.
svensk-pojke Aug 5, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions changelog/285.bugfix.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
Fixed coordinate reversal in VaryingCelestialTransformSlit2D._map_transform().

Rewrote cvt unit tests and ensured both imager and slit classes have 1D - 3D LUT table support and full round-trip
tests for the entire pixel table. There is lots of code duplication in the vct classes and unit tests.
This will be remediated with the refactor of the vct classes in issue 279.

Minor cleanup and reformatting.
63 changes: 35 additions & 28 deletions dkist/wcs/models.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -205,7 +205,7 @@ def _map_transform(self, x, y, z, crpix, cdelt, lon_pole, inverse=False):
# We need to broadcast the arrays together so they are all the same shape
bx, by, bz = np.broadcast_arrays(x, y, z, subok=True)
# Convert the z coordinate into an index to the lookup tables
ind = self.sanitize_index(bz)
zind = self.sanitize_index(bz)

# Generate output arrays (ignore units for simplicity)
if isinstance(bx, u.Quantity):
Expand All @@ -217,12 +217,13 @@ def _map_transform(self, x, y, z, crpix, cdelt, lon_pole, inverse=False):

# We now loop over every unique value of z and compute the transform.
# This means we make the minimum number of calls possible to the transform.
for zind in np.unique(ind):
z_range = np.unique(zind)
for zzind in z_range:
# Scalar parameters are reshaped to be length one arrays by modeling
sct = self.transform_at_index(zind, crpix[0], cdelt[0], lon_pole[0])
sct = self.transform_at_index(zzind, crpix[0], cdelt[0], lon_pole[0])

# Call this transform for all values of x, y where z == zind
mask = ind == zind
mask = zind == zzind
if inverse:
xx, yy = sct.inverse(bx[mask], by[mask])
else:
Expand Down Expand Up @@ -323,8 +324,10 @@ def _map_transform(self, x, y, z, q, crpix, cdelt, lon_pole, inverse=False):

# We now loop over every unique value of z and compute the transform.
# This means we make the minimum number of calls possible to the transform.
for qq in np.unique(qind):
for zz in np.unique(zind):
z_range = np.unique(zind) # raster
q_range = np.unique(qind) # other: maps or meas
for zz in z_range:
for qq in q_range:
# Scalar parameters are reshaped to be length one arrays by modeling
sct = self.transform_at_index((zz, qq), crpix[0], cdelt[0], lon_pole[0])

Expand Down Expand Up @@ -416,15 +419,17 @@ def _map_transform(self, x, y, m, z, q, crpix, cdelt, lon_pole, inverse=False):

# We now loop over every unique value of z and compute the transform.
# This means we make the minimum number of calls possible to the transform.
for mm in np.unique(mind):
for qq in np.unique(qind):
for zz in np.unique(zind):
m_range = np.unique(mind) # raster
q_range = np.unique(qind) # maps
z_range = np.unique(zind) # meas
for mm in m_range:
for zz in z_range:
for qq in q_range:
# Scalar parameters are reshaped to be length one arrays by modeling
sct = self.transform_at_index((mm, zz, qq), crpix[0], cdelt[0], lon_pole[0])

# Call this transform for all values of x, y where z == zind q == qind and m == mind
mask = np.logical_and(mind == mm, zind == zz, qind == qq)
# TODO: Figure out what to do here...
mask = np.logical_and(np.logical_and(mind == mm, zind == zz), qind == qq)
if inverse:
xx, yy = sct.inverse(bx[mask], by[mask])
else:
Expand Down Expand Up @@ -476,7 +481,6 @@ def inverse(self):
return ivct


# TODO: Test
class InverseVaryingCelestialTransform3D(BaseVaryingCelestialTransform3D):
n_inputs = 5
n_outputs = 2
Expand All @@ -497,6 +501,7 @@ def evaluate(self, lon, lat, m, z, q, crpix, cdelt, lon_pole, **kwargs):

class BaseVaryingCelestialTransformSlit(BaseVaryingCelestialTransform, ABC):
def _map_transform(self, x, y, crpix, cdelt, lon_pole, inverse=False):
# x: along slit, y: raster position, both in pixels
# We need to broadcast the arrays together so they are all the same shape
bx, by = np.broadcast_arrays(x, y, subok=True)
# Convert the z coordinate into an index to the lookup tables
Expand All @@ -512,14 +517,15 @@ def _map_transform(self, x, y, crpix, cdelt, lon_pole, inverse=False):

# We now loop over every unique value of y and compute the transform.
# This means we make the minimum number of calls possible to the transform.
for yyind in np.unique(yind):
y_range = np.unique(yind)
for yyind in y_range:
# Scalar parameters are reshaped to be length one arrays by modeling
sct = self.transform_at_index(yyind, crpix[0], cdelt[0], lon_pole[0])

# Call this transform for all values of x, y where y == yind
mask = yind == yyind
if inverse:
# Note this isn't use as the inverse of this model uses the
# Note this isn't used as the inverse of this model uses the
# standard 2D inverse.
xx, yy = sct.inverse(bx[mask], by[mask]) # pragma: no cover
else:
Expand Down Expand Up @@ -636,10 +642,12 @@ def _map_transform(self, x, y, z, crpix, cdelt, lon_pole, inverse=False):

# We now loop over every unique value of y and z and compute the transform.
# This means we make the minimum number of calls possible to the transform.
for yyind in np.unique(yind):
for zzind in np.unique(zind):
y_range = np.unique(yind)
z_range = np.unique(zind)
for yyind in y_range:
for zzind in z_range:
# Scalar parameters are reshaped to be length one arrays by modeling
sct = self.transform_at_index((zzind, yyind), crpix[0], cdelt[0], lon_pole[0])
sct = self.transform_at_index((yyind, zzind), crpix[0], cdelt[0], lon_pole[0])

# Call this transform for all values of x, y where z == zind
mask = np.logical_and(zind == zzind, yind == yyind)
Expand Down Expand Up @@ -689,7 +697,6 @@ def __init__(self, *args, **kwargs):
self.inputs = ("along_slit", "meas_num", "raster", "repeat")
self.outputs = ("lon", "lat")

# TODO: Will this work for 2D if num_meas == 1?
if len(self.table_shape) != 3:
raise ValueError("This model can only be constructed with a three dimensional lookup table.")

Expand Down Expand Up @@ -730,19 +737,20 @@ def _map_transform(self, x, m, y, z, crpix, cdelt, lon_pole, inverse=False):

# We now loop over every unique value of y and z and compute the transform.
# This means we make the minimum number of calls possible to the transform.
# TODO: why this particular order?
for mmind in np.unique(mind):
for yyind in np.unique(yind):
for zzind in np.unique(zind):
m_range = np.unique(mind)
y_range = np.unique(yind)
z_range = np.unique(zind)
for mmind in m_range:
for yyind in y_range:
for zzind in z_range:
# Scalar parameters are reshaped to be length one arrays by modeling
sct = self.transform_at_index((mmind, zzind, yyind), crpix[0], cdelt[0], lon_pole[0])
sct = self.transform_at_index((mmind, yyind, zzind), crpix[0], cdelt[0], lon_pole[0])

# Call this transform for all values of x, y where m == mind, y == yind and z == zind
mask = np.logical_and(mmind == mind, zind == zzind, yind == yyind)
mask = np.logical_and(np.logical_and(mmind == mind, zind == zzind), yind == yyind)
if inverse:
# TODO: figure out what to do here...
# Note this isn't use as the inverse of this model uses the
# standard 2D inverse.
# Note this isn't used as the inverse of this model uses the
# standard 3D inverse.
xx, yy = sct.inverse(bx[mask], by[mask]) # pragma: no cover
else:
xx, yy = sct(bx[mask], by[mask])
Expand All @@ -759,7 +767,6 @@ def _map_transform(self, x, m, y, z, crpix, cdelt, lon_pole, inverse=False):

return x_out, y_out

# TODO: Test
class InverseVaryingCelestialTransformSlit3D(BaseVaryingCelestialTransform3D):
n_inputs = 5
n_outputs = 1
Expand Down
Loading