sunpy
Added new features to the ndcube.__add__ method
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.
A changelog file needs to be added.
And your branch needs to be updated with the latest version of main.
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].
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.
@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)
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:
- 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.
- 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
- 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.
- 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.
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.
@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.
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?
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.
I am not sure how I feel about adding methods to the API for every arithmetic operation. How about functions instead? import ndcube; ndcube.add(cube1, cube2)?
I am not sure how I feel about adding methods to the API for every arithmetic operation. How about functions instead?
import ndcube; ndcube.add(cube1, cube2)?
I'm open to considering that. It might go well in that analysis subpackage I've proposed. I can also see an argument for making it a method as the first argument must always be the NDCube and it's/they are operation(s) valid for any cube. It needs more thought.
Either way, a method can be converted to a function and vice versa quite easily. So for now, I think we should proceed with a NDCube.add method until a longer-term design decision is made. We should discuss this on the sunpy call.
I'm not sure at first glance why this test is failing. I have made some inline suggestions on chages you should make to the code. (You can accept the changes on GitHub, and then pull them you your local branch.)
Once you've made those changes, and assuming the test is still failing, you figure out each step the code takes when the test is being run, and find where the problem is. You can try printing out things at various points in the code and running pytest. See this sunpy explainer: https://docs.sunpy.org/en/latest/dev_guide/contents/tests.html#id2
Alternatively, you can use the pytest --pdb as explained in that sunpy explainer. This will pause the code in an interactive session where the error is raised and you can explore the values of current variables and try running specific lines of code.
This is looking good. The test_cube_add_uncertainty should also check the other attributes are as expected, i.e., the data and the unit. For this, it might be better to not initiate your fixture with np.random.rand, but a known set of values, e.g. np.ones.
Notes for the current stage of development: 1, The 'TypeError: None is not a valid Unit' occurs whenever an NDCube/NDData object with a unit enters the propagate method. Their values become None after they enter the propagate method. Their values remain correct before entering the propagate method. It might happen due to how the unit is handled or converted inside the propagate method.
The stack trace output of pdb w tells me:
the test function performs the addition, the dunder add method is called, inside the dunder add method, the propagate method is called, it then calls the _propagate_add method within the astropy nduncertainty file, which then calls the _propagate_add_sub method, which calls the method .to() for the result_unit_sq, then, inside the .to() method, it tries to convert the unit to the unit of the argument, but finds that unit of result_unit_sq is None. the __call__ method inside the core file is called and raises a TypeError.
result_unit_sq is not set correctly, in _propagate_add_sub() .
result_unit_sq is None, but other objects do have values:
> c:\users\junya\anaconda3\envs\ndcube-dev\lib\site-packages\astropy\nddata\nduncertainty.py(717)_propagate_add_sub()
-> .to(result_unit_sq)
(Pdb) locals()
{'self': StdDevUncertainty([[0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,
0.2],
[0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,
0.2],
[0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,
0.2],
[0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,
0.2],
[0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,
0.2],
[0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,
0.2],
[0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,
0.2],
[0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,
0.2],
[0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,
0.2],
[0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2,
0.2]]), 'other_uncert': StdDevUncertainty([[0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
0.1],
[0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
0.1],
[0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
0.1],
[0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
0.1],
[0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
0.1],
[0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
0.1],
[0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
0.1],
[0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
0.1],
[0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
0.1],
[0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
0.1]]), 'result_data': array([[ 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11.,
12.],
[ 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23.,
24.],
[ 25., 26., 27., 28., 29., 30., 31., 32., 33., 34., 35.,
36.],
[ 37., 38., 39., 40., 41., 42., 43., 44., 45., 46., 47.,
48.],
[ 49., 50., 51., 52., 53., 54., 55., 56., 57., 58., 59.,
60.],
[ 61., 62., 63., 64., 65., 66., 67., 68., 69., 70., 71.,
72.],
[ 73., 74., 75., 76., 77., 78., 79., 80., 81., 82., 83.,
84.],
[ 85., 86., 87., 88., 89., 90., 91., 92., 93., 94., 95.,
[ 85., 86., 87., 88., 89., 90., 91., 92., 93., 94., 95.,
96.],
[ 97., 98., 99., 100., 101., 102., 103., 104., 105., 106., 107.,
[ 97., 98., 99., 100., 101., 102., 103., 104., 105., 106., 107.,
108.],
108.],
[109., 110., 111., 112., 113., 114., 115., 116., 117., 118., 119.,
120.]]), 'correlation': 0, 'subtract': False, 'to_variance': <ufunc 'square'>, 'from_variance': <ufunc 'sqrt'>, 'correlation_sign': 1, 'result_unit_sq': No [109., 110., 111., 112., 113., 114., 115., 116., 117., 118., 119.,
120.]]), 'correlation': 0, 'subtract': False, 'to_variance': <ufunc 'square'>, 'from_variance': <ufunc 'sqrt'>, 'correlation_sign': 1, 'result_unit_sq': None}
Then, when checking how result_unit_sq is used, it shows this:
(Pdb) l
712 ):
713 # If the other uncertainty has a unit and this unit differs
714 # from the unit of the result convert it to the results unit
715 other = (
716 to_variance(other_uncert.array << other_uncert.unit)
717 -> .to(result_unit_sq)
718 .value
719 )
720 else:
721 other = to_variance(other_uncert.array)
722 else:
This is not the same as expected, since it entered the conditional scenario when the two objects' units are different, but in the test case written, the two units are set the same.
TODO: 1) look more into the stack trace that pytest --pdb, ll, w, returns. 2) using pytest --pdb, u, check how result_unit_sq is assigned.
2, The test cases testing the case when neither of them has a unit, and that when only the NDData has a unit, have been implemented and passed. TODO: add the test case for when only NDCube has a unit.
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
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.
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 @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, or3? 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.)
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.
@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
Hi @DanRyanIrish , I have pushed four new commits. Any suggestions would be appreciated. 1, I added three different conditional blocks for the three with-mask scenarios, and to set the uncertainty and mask value for the addition result.
2, I added a new method named combine_uncertainty(), since the code for uncertainty setting is used repeatedly. But I am not sure whether I should add a new method inside the NDCube class.
3, I wrote the tests for the three different with-mask scenarios. But I think the uncertainty testing is not written in the right way. I might need some suggestions on this, thank you!
4, To work on the __mul__ method, I wonder whether I should create a new PR or add a new method within the NDCube class in this branch, please let me know what you think, thank you!
Test Scenarios with mask
-
NDCube has mask but NDData does not
-
NDCube does not have mask but NDData does
-
Both NDCube and NDData have mask handle_mask is True handle_mask is False