Implementation of errors and masks - #202

@eteq I think for errors and masks (and I would focus on the errors first as it is more obvious there) are not specific to spectrum1d. They are a base level NDData problem which is the same ranging from 1d spectra, images, cubes, to n-dimensonal simulations. There are a couple of ways to incorporate them in nddata (that I thought about) and this ones seems to be the most straight forward. I think implementing them in Spectrum1d would inadvertently lead to a lot of code duplication.

We should try and make progress with this. First I think that first we need to separate the issue of masks and errors, as they are quite different, so I propose splitting the pull request.

In the case of errors, I have several issues:

• I have major reservations about the current implementation. In particular, the division and multiplication are only approximations, which break down once errors are significant (e.g. 50%), and I think the results should on the contrary be as accurate as possible.
• I think we should have already several different classes in the initial pull request, to see how things are going to work. For instance, we could have SDError, VarianceError, PDFError, and some kind of log-normal error.
• Related to the previous point, what happens when you add two datasets that have errors with different classes?
• Do we really want to inherit from ndarray? What are the benefits, as opposed to storing the array in a hidden attribute? For example, do all the methods on arrays, e.g. sum, mean, etc. even make sense here.

I guess what I'm trying to say is that we're going to have to put a lot of thought into how to deal with the errors properly. I suggest that we create a wiki page with a list of requirements, code examples, how we want things to behave, then work on the implementation.

For masks, I think that we also want to think about this more, in particular, I kind of agree with @eteq that the arithmetic cannot be done at the level of the mask objects, because they don't know what the user wants (for errors, things are well formulated mathematically, but not for masks).

@astrofrog I agree that this is one of the most important things and we should put some effort into this.
I do not believe masks and errors are that different issues as both of these need to handle arithmetic operations in conjunction with nddata or nddata-based (spectrum1d, ccdimage, ...) objects. The idea is that if the arithemtic it is then passed through to functions of error and mask that know what they are doing.
Anyways I'll address the points in order of appearance.

1. do you mean numerically unstable? How would you make them stable?
2. More is always better ;-) But I think the current SDError pretty well shows what happens. In case of addition, error_add is called and does whatever you want.
3. That should not work. Same with mismatching units and mask types. We came to that conclusion in Spectrum1D for units and I think that should hold here as well. I think it is prudent to have some convenience functions that transform errors from one to the other. so if you have datasets with mismatching errors for example:
   nddata1; nddata2 # are our two nddata sets with one VarianceError and one SDError nddata1 + nddata2 # raises an exception nddata1.error = VarianceError2SDError(nddata1.error) #could also be a method; I don't care nddata1 + nddata2 # all is peachy
4. I have a slight preference for ndarray, as many functions that we use (like broadcast, methods like shape, ndim) automatically work without us having to implement it. But if people have good arguments for not doing that then, we can certainly inherit from object. I think it is important that we have a common superclass that we can do isinstance checking.

I agree we need to put some thought into this, but I think one good thing is to have a pull request that shows one implementation (here's one I prepared earlier ;-) ). I have thought about this problem in a bit of detail as I was trying to implement it for Spectrum1D. I thought of another approach using functions pointers hidden in the meta data, which seemed inferior.

Last but not least: masks. Of course they know what the user wants: If the user wants a BoolMask then it will behave like masked arrays in numpy. Ignoring values in sums and so on. If the user specifies other masks then it will do what is written in these masks. Masks are the same for Spectrum1D, Images, IFUs, simulation data ,.... the same as errors.

A couple of quick notes (need to think about the rest more)

1. Ignore me - I was thinking of error propagation in general, but after double checking, it does seem that even the equation for division is exact. It's when you get to things like power, etc. that things are no longer exact.
2. I think we should have at least two to show how interactions happen (or don't). I think it's also important to show two different types to figure out what is going to be common and can go in a base Error class.

First, regarding splitting: I agree with @astrofrog - they aren't the same thing, if nothing else to be evidenced by the fact that very different topics are coming up here. However, I also agree with @wkerzendorf that we should treat them the same in the sense of how they interact with NDData (e.g. if we say the error objects always get added if the image objects are being added, then we should treat the mask objects the same). But we can split the pull request while still agreeing these should be consistent (so that the specific details of each don't drag down the other).

1. (no longer relevant, but github "cleverly" makes a list starting with 2 start at 1...)
2. This is tied into 3 - just the SD one to start is fine if they don't interact (as @wkerzendorf is suggesting), but if they do interact then we should have at least one other as an example of how they should do it.
3. Sorry to keep causing trouble :) But I think this is not the same as the situation in Spectrum1D - that was about what to do if the wavelength values don't match, because there's lots of ways to combine/interpolate mis-matched spectra. In this case, for errors (but not for masks, as I am pointing out above), there are completely deterministic ways of combining them. E.g. if you add two values where one has S.D. error and one is var errors, it's pretty obvious how they are supposed to be combined. I agree that in many more complicated cases we want to throw an exception (if you try to combine a laplace distributed error with a lognormal error, lets say), but when the way to combine is unambiguous, we want it to just work.
4. I agree we do not want to use ndarray, because it implies there errors are always arrays, which may not be true. What If I want an NDData object where the error is just the poisson distribution centered around the mean value (the value in data)? It should be possible to do this without having to keep a second array in-sync with data, and you can't do that if you subclass from ndarray (at least not at all easily). I agree with @wkerzendorf that a base class with a few basic functions/properties is a very good idea, though (e.g., var, which gives the second-moment of the distribution).

(see my example above for a case where we don't know how the user wants to combine masks without making them decide. People use mask flags for all kinds of wacky things!)

@eteq regarding the arithmetic of nddata objects with different errors and masks: I believe for now it is easiest to throw an exception if they do not match. It is very similar to Spectrum1D's behaviour: e.g. Spectrum1D throws a tantrum if the units are not the same. Both units are frequency or wavelength units, so one can convert the into each other (like sd and var). The question remains what units should the end product have, the units from the first and second one.
I'm not saying that's impossible to do, but it is ambiguous.

Anyways, the whole point is that this implementation allows for such magic to happen in the background. Var-errors can be converted to sd errors on the fly if that is what we want to do. For now, I would throw an exception and worry about making it automagical later.

(4) ...screw you github, I can make my own numbering scheme and there's nothing you can do about it ;-). That is fine. However this reconfirms me in my assumption that error_add and mask_add should have access to the data, masks and errors of both operand and self nddata.

I think in general for the first release we want bool masks and sd errors to lay a foundation that others can build upon.

Just a couple of comments:

• I do think that we should not sub-class from ndarray, because we actually don't expect the dimensions of the errors to match the dimensions of the data. For example, as @eteq pointed out, errors could be stored as scalar values (constant value over the image), we might want to be able to have them be stored as functions, or we might want to define a full PDF for each pixel, or a single PDF over the image, etc. You mentioned that we would have to re-implement some array methods, but most methods would need to be reimplemented to make sense. For example, what does sum mean on an SDError? Is it the sum of the values? Or should they be added in quadrature like standard deviations? Therefore I believe that at least errors should not inherit from ndarray. Note that masks are different, as the point of a mask is to have a 1-to-1 correspondance with the data, so in that case one could argue that subclassing from ndarray makes sense.
• I do think we want to show at least two types of errors, maybe just SDError and VarianceError, to show the mechanism that would check for compatibility with arithmetic operations. I agree we can implement most of the magic later, but we need to ensure that the API is magic-friendly to start with, and the best way is to show a single example (and SDError and VarianceError should be particularly simple to implement).

@wkerzendorf - as I said above, I've come around to your point regarding masks, but I'm not quite following what you are saying about errors... I agree that the point is to have the magic happen in the background, when it's unambiguous -

(4) You raise a good point that there is a need for the masks to have access to the data, though... this is why I was saying before that it might be better to implement this at the e.g. Spectrum1D or Image level, as then there's no need to cross-connect the data and the erors/mask. But I can also see the virtue of this design... maybe the thing to do is use weakrefs to access the data so they don't prevent the data from being garbage-collected?

But I still agree with @astrofrog that if we are going to allow error-conversion, there should be at least two to make it clear who the conversion is supposed to work.

I also still think it would be better if you split this into two pull requests, especially if you are going to split it into a base class for each (mask and errors), which I think you should. Then it will be much easier to review the PRs in detail without having tons of confusing threads in the diffs.

@eteq's comments:
First, It's never unambiguous. Imagine a + b = c(where a,b nddata objects). a has errortype x, b has errortype y. x can be converted to y and vice versa. The ambiguity (which will never go away, unless x=y) , what error type does c have? x or y?

Second: sorry you misunderstood me. I meant the error_add (and friends....) function needs to have access to all the data and not the error object. This way we go around the problematic garbage collection (but good point, which we should think about in this interconnected hell). I believe that the masks and errors do exactley the same no matter what kind of object it is (a spectrum, an ifu, a bird, a plane, ...) so we need to solve this here at ground zero.

Third: Agreed, we need to check if in this framework we can make the conversion-magic happen. This is the last thing on my list (next ones are two pull requests and subclass from object (probably)).

Irregardless of the pull requests, when thinking about masks and errors there are many similarities that we should keep in mind (they are both tag-alongs of nddata, they need to be taken care of in the right way when doing arithmetic, ...).

Finally, I think that when talking about masks, it often seems to be more like flags (@eteq mentioned, saturated, defect, ...) maybe we should call them flags. What do people think? In the end we'll probably inherit the new NA value from numpy 2.0 anyways (which might make boolmasks obsolete; only boolmasks).

@wkerzendorf
First: If you're doing c = a + b, it's the __add__ method of a that gets called, so the error object type to take is the one from a if it's compatible with b. The point I'm making is that it's safe to automatically convert as long as no information is lost. If e.g. a has SDError and b has VarError, c.error.var should still give the same number whether or not c.error is SDError or VarError. If the mapping ever requires a user decision, then it should give up and throw an exception.

second: I don't understand what you're saying here, still. Error.__add__ always has the error information, because it's a method of the error object. Isn't that your whole reason for implementing this at the error/mask level instead of in nddata subclass operators? So that e.g. NDData.__add__ can just do newdata = self.data + other.data, newerror = self.error + other.error and so on?

And you are absolutely right that there are similarities we should keep in mind - I promise we will be looking at both PRs if we split them :) But I think when we have the API ironed out and are getting into the detailed code-review stage, it will be far less confusing if there are separate, smaller PRs.

As for masks/flags, this is a hopeless ambiguity - I've heard both terms used to mean the other one many times. I think mask makes slightly more sense because "mask" is generally more commonly used for pixel-like data (usually there's some agreed definition of a "bad" pixel), and "flags" is more common in catalog settings. But it's really just semantics at this point. But you can't count on using the numpy 2.0 NA value for a long time yet - as you say, we may want to internally use that eventually if numpy 2.x is installed, but we have to continue supporting older versions of numpy (and probably will for many years, based on how slowly astronomers can be to upgrade).

It seems we agree at least on the basics. Which is a huge leap forward in astropy progress (doesn't happen to often ;-) ).

@eteq

1. what error type does c have SDError or VarError. You say the one from a, but why not the one from b. It's ambiguous and just a choice that we make at the implementation time. anyways, we'll worry about that later
2. I'll implement this and you'll see what I mean.
3. I'm looking forward to your critical eyes perusing both my pull requests (that sounds wrong ;-) ).
4. well you're thinking too much along the line of images and spectra, but it can be any array based stuff. anyways, there's a s**t storm coming our way anyways when talking about arithmetic with boolmasks (what to do?: ignore the masked values, add them up but flag them as bad in the final product, etc.). But for now I'm just pretending that everything is fine and will try to live in blissful ignorance.

Good to see progress is being made!

Should we close this PR in favor of #223 and #224 ?

I have left this open to see if there's a few things missing, I guess no one complained that there's fundamental differences. So I close this one.