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

Added new features to the ndcube.__add__ method #794

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

Conversation

PCJY
Copy link
Contributor

@PCJY PCJY commented Dec 11, 2024

PR Description

This PR aims to fix issue #734 created by @DanRyanIrish .

write down the scenarios.

(list, to see what has been covered)

Testing scenarios without mask

redrafting the tests.

  • One and only one of them has a unit. [Error]

  • Both NDCube and NDData have unit
    Both have uncertainty.
    NDCube has uncertainty.
    NDData has uncertainty.
    Neither has uncertainty.

  • Neither has a unit.
    Both have uncertainty.
    NDCube has uncertainty.
    NDData has uncertainty.
    Neither has uncertainty.

name them again:
test_cube_add_cube_unit_mask_nddata_unc_unit_mask

Handling Masks

Determing data result of operation

The operation_ignores_mask kwarg determines the resulting data value when adding objects, and accounting for mask values. Below are the different scenarios for the addition of a single pixel NDCube, named cube with a data value of 1, and a single pixel NDData object, named nddata with a data value of 2:

my draft
what this means:
when OIM is T, does it ignore the mask after the addition or during the addition?
during the addition, because we want to see the arithmetic operation results here but not the mask. if it deals with the mask after the addition, there would be no need to check the arithmetic values here
no need to do boolean operations here for the masks?
this is why mask handling and arithmetic operation are separate with each other.
data and mask are two separate things: First data, then mask;
addition is done on both the data and the mask.

cube.data = 1, nddata.data = 2

cube.mask nddata.mask operation_ignores_mask resulting data value
F F T 3
F F F 3
T F T 3
T F F 2
F T T 3
F T F 1
T T T 3
T T F None

What the distinct cases are:
When OIM is T, the result is always 3, i.e. the actual result of the addition of the two data values.
When OIM is F, the result is always the addition result of any value with its corresponding mask being F (when there is not one, the result is None).

Determining the new mask value produced by the operation

The handle_mask kwarg takes a function which determines the new mask value resulting from the addition. While this function can be anything that takes two boolean arrays and outputs a resulting boolean array of the shape, the most commonly used are expected to by numpy.logical_and and numpy.logical_or. Since the user supplies the function, the resulting mask can be implemented something like:

new_mask = handle_mask(self.mask, value.mask) if handle_mask else None

My understanding:
handle_mask is a function,
as long as it has a value, it will be True, and does the operation on the two masks,
otherwise, the new_mask value would be None, meaning:
1), the user did not set anything for the handle_mask kwarg, should an error be raised here?
2), there is no need to set any value for the handle_mask kwarg.

ndcube/ndcube.py Outdated
Comment on lines 930 to 931
# addition
new_data = self.data + value_data
Copy link
Member

Choose a reason for hiding this comment

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

The addition should be done as part of the masked array addition. You've already done this below, you just need to extract the added data from the results as well as the mask.

ndcube/ndcube.py Outdated
Comment on lines 950 to 952
return self._new_instance(
data=new_data, uncertainty=new_uncertainty, mask=new_mask
)
Copy link
Member

Choose a reason for hiding this comment

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

Instead of having a separate return here for the NDData case, I think we should build a dictionary of kwargs that we can give self._new_instance, here. So, you can create an empty kwargs dictionary at the start of the method, and add the new data, uncertainty, etc. in the relevant place, e.g.

kwargs["uncertainty"] = new_uncertainty

Then the final line of the method would become

return self._new_instance(**kwargs)

Copy link
Member

Choose a reason for hiding this comment

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

Let me know if this doesn't make sense

@nabobalis nabobalis added this to the 2.4.0 milestone Dec 18, 2024
ndcube/ndcube.py Outdated
if self.uncertainty is not None and value.uncertainty is not None:
new_uncertainty = self.uncertainty.propagate(
np.add, value.uncertainty, correlation=0
np.add, value.uncertainty, result_data = value.data, correlation=0
Copy link
Member

Choose a reason for hiding this comment

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

The result_data needs to be the result of the operation. So, assuming you moved the addition of the datas using the masked array to before the uncertainty propagation, you could do:

Suggested change
np.add, value.uncertainty, result_data = value.data, correlation=0
np.add, value.uncertainty, result_data = kwargs["data"], correlation=0

ndcube/ndcube.py Outdated
Comment on lines 1061 to 1070
# combine mask
self_ma = np.ma.MaskedArray(self.data, mask=self.mask)
value_ma = np.ma.MaskedArray(value_data, mask=value.mask)

# addition
result_ma = self_ma + value_ma
new_mask = result_ma.mask

# extract new mask and new data
kwargs["mask"] = result_ma.mask
kwargs["data"] = result_ma.data
Copy link
Member

Choose a reason for hiding this comment

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

As mentioned in above comment, I think it makes sense to do this before the uncertainty propagation so you can use the kwargs["data"] value in that propagation.

ndcube/ndcube.py Outdated
kwargs["data"] = result_ma.data

# return the new NDCube instance
return self._new_instance(**kwargs)
Copy link
Member

Choose a reason for hiding this comment

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

Move this line to the end of the method and use the kwargs approach when handling the other cases, e.g. Quantity. So, for example, L1082 would become:

kwargs["data"] = self.data + value.to_value(cube_unit)

@DanRyanIrish
Copy link
Member

A changelog file needs to be added.

And your branch needs to be updated with the latest version of main.

@PCJY
Copy link
Contributor Author

PCJY commented Jan 18, 2025

Hi @DanRyanIrish, as we have discussed in our project meetings, below are the issues we encountered and may need further discussions with others in the community:

The issue is mainly around how NumPy handles masks when performing an addition for two NumPy.MaskedArray.
We think the expected outcome for an addition should be: the sum of any value that is not masked by its individual mask.
E.g. [1] ([T]) + [2] ([F]) = [2].

However, from experimentation, it can be seen that
NumPy returns in this way:
[1] ([T]) + [2] ([F]) = [1].
Screenshot 2025-01-18 220004

I find this confusing because even if it does combine the mask and then apply it on the result, it should be:
[1] ([T]) + [2] ([F]) = [-].

Please correct me if there is anything wrong in my understanding.

@PCJY
Copy link
Contributor Author

PCJY commented Jan 18, 2025

@DanRyanIrish, Secondly, we also encountered some issues around the propagate method:
it ignores the mask of the objects that are passed in, and still takes into account the uncertainties of the masked elements when it should not have done so.
Following your guidance and suggestions, this issue is currently being worked on by setting the corresponding entries that should be masked of the uncertainty array to be 0, before passed in to the propagate method.
A clearer example of the issue was implemented as shown in the code below with the screenshot of its output attached.

from ndcube import NDCube
import numpy as np
from astropy.nddata import StdDevUncertainty
from astropy.wcs import WCS

data = np.array([[1, 2], [3, 4]])  
uncertainty = StdDevUncertainty(np.array([[0.1, 0.2], [0.3, 0.4]])) 
mask = np.array([[False, True], [False, False]])  
wcs1 = WCS(naxis=2) 
wcs1.wcs.ctype = ["HPLT-TAN", "HPLN-TAN"]

cube = NDCube(data, wcs=wcs1, uncertainty=uncertainty, mask=mask)
print(cube)

def add_operation(cube1, cube2):
    """
    Example function to add two data arrays with uncertainty propagation.
    """
    result_data = cube1.data + cube2.data 
    # Propagate the uncertainties using the NDCube objects
    propagated_uncertainty = cube1.uncertainty.propagate(
        np.add, cube2, result_data=result_data, correlation = 0
    )
    return result_data, propagated_uncertainty

# adding the cube to itself
result_data, propagated_uncertainty = add_operation(cube, cube)

print("Original Data:\n", cube.data)
print("Original Uncertainty:\n", cube.uncertainty.array)
print("Result Data (after addition):\n", result_data)
print("Propagated Uncertainty:\n", propagated_uncertainty.array)

Screenshot 2025-01-18 222146

@DanRyanIrish
Copy link
Member

DanRyanIrish commented Jan 20, 2025

Hi @PCJY. I think the first thing we need to do is decide what behaviours we want to implement in the case where at least one of the NDCube and NDData have a mask. I think we need to get some feedback from other users on this decision. I propose the following scheme (@Cadair, thoughts on this?):

Firstly, if an object has no mask, that is equivalent to all pixels being unmasked.
Secondly, for a given pixel in both objects:

  1. If both are unmasked, the resultant
    i. data value is the sum of both pixels
    ii. mask value is False
    iii. uncertainty value is the propagation of the two uncertainties. If one or other object doesn't have uncertainty, the uncertainty of that component is assumed to be 0.
  2. If it is masked in one object, but not the other, the resultant:
    i. data value is equal to the unmasked value
    ii. mask value is False
    iii. uncertainty value is the same as the unmasked pixel
  3. If both pixels are masked, this is where is gets ambiguous. I propose, in order to remain consistent with the above:
    i. The operation is not performed and the data, mask and uncertainty values remain the same as the left-hand operand, i.e. the NDCube.

Alternatives for parts of the scheme could include:
2. If it is masked in one object, but not the other, the resultant:
ii. mask value is True.

  1. If both pixels are masked:
    i. The operation IS performed as normal but the mask value is True.

Once we agree on a scheme, the way forward on your uncertainty questions will become clear.

@Cadair what are your thoughts on this scheme. I also think we should bring this up at the sunpy weekly meeting to get other thoughts.

@DanRyanIrish
Copy link
Member

I find this confusing because even if it does combine the mask and then apply it on the result, it should be: [1] ([T]) + [2] ([F]) = [-].

This is where I also find numpy masked arrays counter-intuitive. However, the logic is as follows:

  • If one pixel is masked, retain the data of the left-hand operand and set the mask value to True.
    Notice that order of the operands matters. Because you did [1] ([T]) + [2] ([F]), the results is [1] ([T]), which is displayed as [--]. I would expect if you did the operation the other way around ([2] ([F]) + [1] ([T])), the result would be [2] ([T]).

Notice that this is not the same as the scheme I've proposed in my previous comment, in part because it's confusing, as you've found.

@DanRyanIrish
Copy link
Member

DanRyanIrish commented Jan 20, 2025

@PCJY, until we agree a way forward with the mask, you should proceed by implementing in the case for which neither object has a mask. So no need for masked arrays.

@Cadair
Copy link
Member

Cadair commented Jan 21, 2025

I propose the following scheme

I haven't thought too much each of these individual cases, but the fact there is a list is enough to make me think we probably need a way for the user to choose. This is obviously not possible with the __add__ operator, so we would need to have a NDCube.add method (and presumably subtract) which accepted kwargs.

Is this also not a problem for other operators as well?

@DanRyanIrish
Copy link
Member

I think this is a good idea. As well as add and subtract, I think we would also need multiply and divide methods.

As far as I can see, this ambiguity only arises when there are masks involved. So we could still implement the dunder methods as wrappers around the above methods, but have the raise and error/return NotImplemented if the non-NDCube operand has a mask, and require users use the NDCube.add method instead.

@DanRyanIrish DanRyanIrish changed the title Added new features to the ndcube._add_ method Added new features to the ndcube.__add__ method Feb 25, 2025
@DanRyanIrish
Copy link
Member

Here is an example of including fixtures in parameterisations: https://github.com/sunpy/ndcube/blob/main/ndcube/tests/test_ndcube.py#L48
The trick is to add this indirect kwarg at the end of the parameterisation given the name of the input variable(s) that has to be retrieved from conftext.py

Copy link
Member

@DanRyanIrish DanRyanIrish left a comment

Choose a reason for hiding this comment

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

The test function names are now much clearer, and it's easier to see the cases that are tested.

A couple more comments on your test writing:

  1. At least one of your tests is somewhat circular, i.e. you use the result of the operation you're testing to form the expected values. This means that an error might make it into the expected values and the test will pass when it should fail. See inline comments for a better way to structure tests. These are for information only. There is an even better way to structure these tests. See item below.

  2. ndcube provides a test helper function that checks if two cubes are equal. This means your test only needs to construct the expected cube, and then pass it and the result of the addition to that helper function. See below for an example of how to do this for one test

def test_cube_add_cube_unit_unc_nddata_unit_unc(ndc, value):
    output_cube = ndc + value # perform the addition
    # Construct expected cube
    expected_unit = u.ct
    expected_data = ((ndc.data * ndc.unit) + (value.data * value.unit)).to_value(expected_unit)
    expected_uncertainty = ndc.uncertainty.propagate(
                            operation=np.add,
                            other_nddata=value,
                            result_data=expected.data*expected.unit,
                            correlation=0,
    )
    expected_cube = NDCube(expected_data, ndc.wcs, uncertainty=expected_uncertainty, unit=expected_unit)
    # Assert output cube is same as expected cube
    assert_cubes_equal(output_cube, expected_cube)

All your tests in this PR should be written like this to ensure the output cubes are fully tested. It should also make them easier to read.

Comment on lines 1177 to 1180
expected_uncertainty = ndc.uncertainty.propagate(
operation=np.add,
other_nddata=value,
result_data=new_cube.data*new_cube.unit,
Copy link
Member

Choose a reason for hiding this comment

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

This test is a bit circular. You shouldn't use outputs of the test to form expected values. Otherwise the test may pass when it shouldn't.

Suggested change
expected_uncertainty = ndc.uncertainty.propagate(
operation=np.add,
other_nddata=value,
result_data=new_cube.data*new_cube.unit,
expected_unit = u.ct
expected_data = (ndc.data * ndc.unit).to_value(expected_unit) + (value.data * value.unit).to_value(expected_unit)
expected_uncertainty = ndc.uncertainty.propagate(
operation=np.add,
other_nddata=value,
result_data=expected_data*expected_unit,

Copy link
Contributor Author

@PCJY PCJY Feb 28, 2025

Choose a reason for hiding this comment

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

Hi @DanRyanIrish , thank you for the suggestions. I checked the code for the assert_cubes_equal method. I am unsure whether it checks the values, type and units of the uncertainty attributes as well.
It looks like it only checks whether the shapes of uncertainties are the same.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@DanRyanIrish Or maybe I can
either: use the assert_cube_equal method together with a few more lines that checks the values, type and units of the uncertainty attributes,
or: adding a few more lines into the assert_cube_equal method?

Copy link
Member

@DanRyanIrish DanRyanIrish Feb 28, 2025

Choose a reason for hiding this comment

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

Hi @PCJY. Well spotted. You're right. You should include a check_uncertainty_values=False kwarg to assert_cube_equal and make it check those aspects of the uncertainty if set to True. So the code here would be replaced by something like:

if check_uncertainty_values:
    # Check output and expected uncertainty are of same type. Remember they could be None.
    # If the uncertainties are not None,...
    # Check units, shape, and values of the uncertainty.
elif test_input.uncertainty:
        assert test_input.uncertainty.array.shape == expected_cube.uncertainty.array.shape

Then you can set check_uncertainty_values to True when you call assert_cubes_equal in your tests, and that should do what you need it to do.

Copy link
Contributor Author

@PCJY PCJY Feb 28, 2025

Choose a reason for hiding this comment

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

@DanRyanIrish I see what you mean, thank you for this suggestion! I will implement this.

Comment on lines 1183 to 1184
assert np.allclose(new_cube.data, ndc.data + value.data)
assert new_cube.unit == u.ct
Copy link
Member

Choose a reason for hiding this comment

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

Using above suggested changes:

Suggested change
assert np.allclose(new_cube.data, ndc.data + value.data)
assert new_cube.unit == u.ct
assert np.allclose(new_cube.data, expected_data)
assert new_cube.unit == expected_unit

@DanRyanIrish
Copy link
Member

Hi @PCJY. The tests are looking good now. However, could you please rename them test_arithmetic_add_..., rather than test_cube_add..., as the test actually uses the arithmetic operator +, and only tests the NDCube.add method indirectly.

else:
assert (test_input.uncertainty is None and expected_cube.uncertainty is None)
pass
Copy link
Member

Choose a reason for hiding this comment

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

Now the way you have it, you could do this:

Suggested change
pass
assert type(test_input.uncertainty) is type(expected_cube.uncertainty)
assert np.allclose(test_input.uncertainty.array, expected_cube.uncertainty.array), \
f"Expected uncertainty: {expected_cube.uncertainty}, but got: {test_input.uncertainty.array}"

and then remove lines 135-139. Up to you though.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@DanRyanIrish Thank you for your suggestion. I re-understood what the code aims to do, and changed it again.

@DanRyanIrish
Copy link
Member

Hi @PCJY. I've edited your mask discussion in the PR description. I've created a table for all the scenarios for determining the result of adding the masked data. Can you fill in the table with the result of each case? Either 1, 2, or 3? Once you've done that, we can look at the table and, like before, determine how to implement them, and if we have some effectively duplicated cases.

@PCJY
Copy link
Contributor Author

PCJY commented Mar 4, 2025

Hi @PCJY. I've edited your mask discussion in the PR description. I've created a table for all the scenarios for determining the result of adding the masked data. Can you fill in the table with the result of each case? Either 1, 2, or 3? Once you've done that, we can look at the table and, like before, determine how to implement them, and if we have some effectively duplicated cases.

Hi @DanRyanIrish, thank you for providing the logic structures. I have filled in the results, please feel free to tell me whether I made any mistakes. (I showed my understanding process in the comments as well.)

@DanRyanIrish
Copy link
Member

DanRyanIrish commented Mar 6, 2025

Hi @PCJY Regarding you question in the PR description:

My understanding:
handle_mask is a function, as long as it has a value, it will be True, and does the operation on the two masks, otherwise, the new_mask value would be None, meaning:
1), the user did not set anything for the handle_mask kwarg, should an error be raised here?
2), there is no need to set any value for the handle_mask kwarg.

Yes, that's right. Regarding your two options, handle_mask will have a default value. (We can decide what that should be later.) But that means the user doesn't have to set it and there's no need to raise an error. The default value will trigger the default mask handling behaviour.

@DanRyanIrish
Copy link
Member

DanRyanIrish commented Mar 6, 2025

@PCJY: Regarding your comments on the distinct cases in the PR description:

When OIM is T, the result is always 3, i.e. the actual result of the addition of the two data values.

Yes :)

When OIM is F, the result is always the addition result of any value with its corresponding mask being F

Yes!

(when there is not one, the result is None).

Yet to be decided

So I think a way to implement this would be:

self_data, value_data = self.data, value.data # May require a copy
self_mask, value_mask = self.mask, value.mask # May require handling/converting of cases when masks aren't boolean arrays but are None, True, or False.
if not operation_ignores_mask:
    no_op_value = 0 # Value to set masked values since we are doing addition. (Would need to be 1 if we were doing multiplication.)
    idx = np.logical_and(self_mask, np.logical_not(value_mask))
    self_data[idx] = no_op_value
    idx = np.logical_and(value_mask, np.logical_not(self_mask))
    value_data[idx] = no_op_value
    idx = np.logical_and(self_mask, value_mask)
    #self_data[idx], value_data[idx] = ?, ? # Handle case when both values are masked here. # We are yet to decide the best behaviour here.
else:
    pass # If operation ignores mask, nothing needs to be done. This line not needed in actual code.  Only here for clarity.
# Perform addition
new_data = self_data + value_data
# Calculate new mask.
new_mask = handle_mask(self_mask, value_mask) if handle_mask else None

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

Successfully merging this pull request may close these issues.

4 participants