Identity and Mapping transforms always return units (wrongly) in v0.24

[chris-simpson]

There seems to have been an unintended consequence of PR #457 with regard to the logic for units. This code

>>> in_frame = Frame2D(name="in_frame")
>>> out_frame = Frame2D(name="out_frame")
>>> wcs = gWCS([(in_frame, models.Identity(2)),
                (out_frame, None)])
>>> wcs(10, 10)

returns (10.0, 10.0) under gwcs 0.22.1 (as it should, since with_units is False by default), but (<Quantity 10. pix>, <Quantity 10. pix>) under 0.24.

The cause is the logic for adding units, which is if not input_is_quantity and transform.uses_quantity while the documentation for astropy.modeling.models.Model states that uses_quantity is True if "this model has been created with Quantity objects or if there are no parameters", and Identity models fall into the second category. Mapping models also have the same issue. (The Projection models have no parameters but are unlikely to be used in isolation.)

This looks fixable with the logic

if not input_is_quantity and transform.uses_quantity and transform.parameters.size

(pick your own idiom for checking for a non-empty ndarray) since that won't add units if the inputs don't have them and uses_quantity is only True because there are no parameters. But this may also have a knock-on effect (or I may have missed some additional logic) so I'm not going to submit a PR; at least not yet.

Edit: Tabular1D and Tabular2D also have no parameters.

[chris-simpson]

OK, there's quite a bit more needed and I don't understand some of the code so I'm struggling to write a PR to fix this bug. All these things should return output that's the same as the input:

wcs(10, 10)
wcs(10 * u.pix, 10 * u.pix)
wcs.invert(10, 10)
wcs.invert(10 * u.pix, 10 * u.pix)

but the problem is the last one. This is prevented from knowing about its units because they're stripped by these lines at the start of _call_backward():

        if is_high_level(*args, low_level_wcs=self):
            args = high_level_objects_to_values(*args, low_level_wcs=self)

There are no similar lines at the start of _call_forward() so I don't understand what they're for since I believe the forward and backward transforms should be treated equally. If I remove them, the behaviour of the above calls is correct.

[WilliamJamieson]

There are more issues. Tablular models with units on the "points" argument do actually need units and are indistinguishable using the method you suggest from ones without units on the "points".

[WilliamJamieson]

The core issue at play here is that uses_quantity doesn't actually work that well. Especially, when we consider people writing custom models. E.G. what about a custom model built off the function $f(x) = p* e^x$, for some parameter $p$. uses_quantity will indicate that such a model uses units, but passing units in this case will cause unit related errors (unless the unit is "dimensionless").

The larger issue is that there is not a reliable way to determine if a model requires units or not.

[WilliamJamieson]

        if is_high_level(*args, low_level_wcs=self):
            args = high_level_objects_to_values(*args, low_level_wcs=self)

There are no similar lines at the start of _call_forward() so I don't understand what they're for since I believe the forward and backward transforms should be treated equally. If I remove them, the behaviour of the above calls is correct.

There is an API difference __call__ does not accept things like astropy.coordinates objects while invert does. That call turns the astropy.coordinates objects into arrays.

Moreover, Quantities are considered "high-level" objects by APE 14 (which defines our API conventions) and the with_units on invert controls whether a "high-level" or "low-level" object is returned. Thus in order to get units out, we should expect to have to pass with_units=True to invert.

Now lets consider your example

from gwcs import coordinate_frames as cf
from gwcs import wcs
from astropy.modeling import models
from astropy import units as u

in_frame = cf.Frame2D(name="in_frame")
out_frame = cf.Frame2D(name="out_frame")

gwcs = wcs.WCS(
    [
        (in_frame, models.Identity(2)),
        (out_frame, None),
    ]
)

print(f"1. {gwcs(10, 10)=}")
print(f"2. {gwcs(10*u.pix, 10*u.pix)=}")
print(f"3. {gwcs.invert(10, 10)=}")
print(f"4. {gwcs.invert(10*u.pix, 10*u.pix)=}")

Outputs:

1. gwcs(10, 10)=(<Quantity 10. pix>, <Quantity 10. pix>)
2. gwcs(10*u.pix, 10*u.pix)=(<Quantity 10. pix>, <Quantity 10. pix>)
3. gwcs.invert(10, 10)=(<Quantity 10. pix>, <Quantity 10. pix>)
4. gwcs.invert(10*u.pix, 10*u.pix)=(<Quantity 10. pix>, <Quantity 10. pix>)

It is my understanding that you believe the first and third cases are wrong. However, the 4th case is also wrong per the API for invert.

That is the expected behavior is

1. gwcs(10, 10)=(10.0, 10.0)
2. gwcs(10*u.pix, 10*u.pix)=(<Quantity 10. pix>, <Quantity 10. pix>)
3. gwcs.invert(10, 10)=(10.0, 10.0)
4. gwcs.invert(10*u.pix, 10*u.pix)=(10.0, 10.0)

There are two code blocks that where the issue originates

gwcs/gwcs/wcs/_wcs.py

Lines 207 to 211 in bfedbdf

	input_is_quantity = any(isinstance(a, u.Quantity) for a in args)
	if not input_is_quantity and transform.uses_quantity:
	args = self._add_units_input(args, from_frame)
	if not transform.uses_quantity and input_is_quantity:
	args = self._remove_units_input(args, from_frame)

gwcs/gwcs/wcs/_wcs.py

Lines 337 to 342 in bfedbdf

	# Validate that the input type matches what the transform expects
	input_is_quantity = any(isinstance(a, u.Quantity) for a in args)
	if not input_is_quantity and transform.uses_quantity:
	args = self._add_units_input(args, self.output_frame)
	if not transform.uses_quantity and input_is_quantity:
	args = self._remove_units_input(args, self.output_frame)

Which both are (literally) identical logic. In all cases the root cause of the issue is self._add_units_input being called erroneously. Making your fix at both points does indeed correct the behavior to match the API for this case. However, it breaks down for the general model. For example it will break the behavior of:

gwcs/gwcs/examples.py

Lines 258 to 279 in bfedbdf

	def gwcs_stokes_lookup():
	transform = models.Tabular1D(
	[0, 1, 2, 3] * u.pix,
	[1, 2, 3, 4] * u.one,
	method="nearest",
	fill_value=np.nan,
	bounds_error=False,
	)
	frame = cf.StokesFrame()
	detector_frame = cf.CoordinateFrame(
	name="detector",
	naxes=1,
	axes_order=(0,),
	axes_type=("pixel",),
	axes_names=("x",),
	unit=(u.pix,),
	)
	return wcs.WCS(
	forward_transform=transform, output_frame=frame, input_frame=detector_frame
	)

much more completely because an actual error is raised as I mentioned in my earlier comments.

Ultimately, the issue is that model.uses_quantity is not quite telling us the correct thing.

[chris-simpson]

There is an API difference call does not accept things like astropy.coordinates objects while invert does. That call turns the astropy.coordinates objects into arrays.

Why do __call__() and invert() behave differently? I note that __call__() handles to_frame and from_frame in _call_forward() but invert() has no such code in _call_backward(). I'm struggling to see a reason for handling forward and backward transforms differently; indeed, I feel I should be able to do gw(*coords, from_frame=gw.output_frame, to_frame=gw.input_frame) and get the same result as gw.invert(*coords). I appreciate that, where the transform has no analytical inverse and so the numerical_inverse must be used, there may be limitations, but this seems fine to me.

Here's what I think the logic should be when evaluating a gWCS object from a from_frame to a to_frame via a transform.

1. Convert high-level objects in the input to Quantities
2. If the transform has units and the input doesn't, attach the from_frame's units to the input
3. If the transform doesn't have units (including the case of a numerical inverse) and the input does, convert the input to the from_frame's units and strip the units
4. Evaluate the transform with the (modified if necessary) inputs
5. If the output doesn't have units, attach the to_frame's units to the output; otherwise convert the output to the to_frame's units
6. Strip the units in appropriate circumstances, e.g., no input_units and with_units=False, or convert the output into high-level objects

That looks to me to be totally agnostic to the direction of travel. I'm going to work on this and see where it leads.

The issue of whether the transform has units or not is still unclear since, as you note, uses_quantity doesn't do what we want. But it may be possible to get this information another way by determining whether its parameters have units defined. Or just by attempting to evaluate it and catching the appropriate exception.

[nden]

Clearly we need to reassess how the code works with and w/o units. @chris-simpson Is there a way for us to run your tests downstream on PRs, similar to dkist, jwst and roman? If not, can you run nightly tests with gwcs/main so we catch these as early as possible?

[WilliamJamieson]

The issue of whether the transform has units or not is still unclear since, as you note, uses_quantity doesn't do what we want. But it may be possible to get this information another way by determining whether its parameters have units defined.

This is precisely what uses_quantity does. The issue that one has with this approach is that in general it is fairly easy to define a model that violates any attempt at analysis via this method. The Tabular1D model for example has no parameters but depending on an instance's construction could require or not require units. If we take custom defined models into account one can get even more perverse.

Or just by attempting to evaluate it and catching the appropriate exception.

This is the method that I came up with last week to handle this issue. I will open a PR to this effect once I am happy with my changes.

[WilliamJamieson]

Why do __call__() and invert() behave differently? I note that __call__() handles to_frame and from_frame in _call_forward() but invert() has no such code in _call_backward(). I'm struggling to see a reason for handling forward and backward transforms differently; indeed, I feel I should be able to do gw(*coords, from_frame=gw.output_frame, to_frame=gw.input_frame) and get the same result as gw.invert(*coords). I appreciate that, where the transform has no analytical inverse and so the numerical_inverse must be used, there may be limitations, but this seems fine to me.

Here's what I think the logic should be when evaluating a gWCS object from a from_frame to a to_frame via a transform.

1. Convert high-level objects in the input to Quantities
2. If the transform has units and the input doesn't, attach the from_frame's units to the input
3. If the transform doesn't have units (including the case of a numerical inverse) and the input does, convert the input to the from_frame's units and strip the units
4. Evaluate the transform with the (modified if necessary) inputs
5. If the output doesn't have units, attach the to_frame's units to the output; otherwise convert the output to the to_frame's units
6. Strip the units in appropriate circumstances, e.g., no input_units and with_units=False, or convert the output into high-level objects

That looks to me to be totally agnostic to the direction of travel. I'm going to work on this and see where it leads.

I believe the reason for why they act differently is that they arose organically and their API has been adjusted as needs arose. The APE 14 API is much more well defined and self consistent. For "low-level" inputs the those are

• gwcs/gwcs/api.py

  Line 83 in bfedbdf

  def pixel_to_world_values(self, *pixel_arrays):

• gwcs/gwcs/api.py

  Line 100 in bfedbdf

  def array_index_to_world_values(self, *index_arrays):

• gwcs/gwcs/api.py

  Line 111 in bfedbdf

  def world_to_pixel_values(self, *world_arrays):

• gwcs/gwcs/api.py

  Line 127 in bfedbdf

  def world_to_array_index_values(self, *world_arrays):

For "high-level" inputs those are (inherited from astropy):

• pixel_to_world
• array_index_to_world
• world_to_pixel
• world_to_array_index

My suggestion is that if possible you should call one of these functions rather than the __call__ or invert.

[jehturner]

Clearly we need to reassess how the code works with and w/o units. @chris-simpson Is there a way for us to run your tests downstream on PRs, similar to dkist, jwst and roman? If not, can you run nightly tests with gwcs/main so we catch these as early as possible?

Hi Nadia!

My apologies that we have not been doing a good job of testing upstream changes early (especially when you've pinged us on some of these issues before and we haven't had the bandwidth to comment much). I'll try to get a job set up for this (probably just on the main branch for now).

Thanks,

James.

[chris-simpson]

My suggestion is that if possible you should call one of these functions rather than the call or invert.

That's a significant refactor of our code. I don't see any reason why the two directions of flow can't be homogenized as I suggest. I was going to work on a fix for uses_quantity first but if you're going to do that, then I'll get straight onto having a go at this refactor.

Edit: No I won't. We have too many pressing issues. We'll need to fix to an older version of gwcs. Once this issue and #559 are fixed I'll try to deal with things. There seems to also be vestigial code, as in the transform() method which calls _call_forward() even when it's doing a backward transform, and so breaks if the inverse isn't defined. I'm not sure whether this method is needed (apart from in a couple of tests).

[jehturner]

@nden, I see what you're doing with the downstream tests now. Our complete DRAGONS tests run on 100s of GB of data on an internal Jenkins server, but there are also some more basic unit tests that run on GitHub actions. The basic tests related to WCS live in the astrodata sub-package, which we are just finishing splitting out from DRAGONS and has just been accepted as an affiliated package. It would therefore make sense if you're able to incorporate the new astrodata into your downstream testing, which @teald hopes to work on in the next few days. That would have shown up at least some of the problems reported by Chris. We can also try to run Jenkins against your dev branch periodically, but it probably doesn't make sense to tie those tests into PRs, since other people can't see the logs (just a pass/fail status). Chris thinks the AstroData tests will already be a good sanity check though and it's obviously important for affiliated packages to work together...

[WilliamJamieson]

Edit: No I won't. We have too many pressing issues. We'll need to fix to an older version of gwcs. Once this issue and #559 are fixed I'll try to deal with things. There seems to also be vestigial code, as in the transform() method which calls _call_forward() even when it's doing a backward transform, and so breaks if the inverse isn't defined. I'm not sure whether this method is needed (apart from in a couple of tests).

I'm moving the API discussion to a separate issue as it is really as separate issue from the bug originally reported here. I also feel that it is important in its own right.