Nderror implementation - #223

(I haven't looked at this yet, but @wkerzendorf, FYI, it's useful to use github's automatic issue-linking so we can see the extra bonus-info github provides - e.g., if you just do #202 when leaving a github comment, it creates a link, and shows up in the other place. e.g.: #202 )

@wkerzendorf - are you going to have time to work on this, or shall I start a branch with your commits in my fork, implement my own comments, and open a new PR?

@astrofrog I have some time to work on it, but you can also start a new PR and close this one if you think it's easier.

So here is what I was thinking when implementing it: the errror functions (arithmetic, interpolation, convolution) should have access to results, flags masks and errors. With very simple errors that we're implementing right now, that is not necessarily needed (well if masked values have 0 then you get a nasty output on the console when dividing by them afaik). But more complex errors can have a different behaviour according to what the mask or flag is at that point. I would just pass it everything, thus later on we don't have to change it if someone wants a more complex error.
Does that make sense?

I don't see how we can implement methods like interpolation and convolution in NDError classes, as the implementation will be different depending on the dimensionality of the data.

I would vote to remove the convolution and interpolation methods until we actually need them. After all, when finding the errors when interpolating and convolving, we might be able to get by with just combining errors with the standard operations. Let's keep the error class as simple as possible.

I don't quite understand what you mean: what happens to the error if you convolve your nddata?

How do we see interacting with these Error classes when we need to access subsections of the error array? For example, in the photutils affiliated package, we have an aperture_photometry() function. It currently accepts separate data, mask, error arrays. But, we would like to let it accept Image objects (derived from NDData). Image.error would be some sort of NDError-derived object.

Currently the function simply accepts an error array, which we dicate is the 1-sigma error. In the course of the function, we do things like taking a slice of the error array, multiplying it by another array, and summing the array.

What would we do instead of these things if handed an NDError-derived object instead of a simple array?

well it will always get an NDError derived object. As with arithmetic I envisage that if slicing happens (which is a __getitem__ call) this will be handed down to the array.
pseudo code

class SDError(object):
   def __getitem__(self, index):
          return self.__class__(self.error.__getitem__(index))

This just gives back a new sliced error array (well I guess more of a new view if I remember my numpy correctly).

One thing to think about is how slicing works if the dimensions don't match up. I guess some fancy broadcasting magic.

In the simplest case where error.ndim==0 that should be easy.

... and then I guess there would be methods in SDError for multiplying by an ndarray, and for summing itself? It seems like there are a lot more methods that need to be implemented.

What is the overhead of having to create a new object every time I want to get a slice? I know we don't have to copy the internal array, but an object still needs to be created. In aperture_photometry() we're doing these slices a lot.

@kbarbary So what ever NDData does like a sum of itself in either axes SDError needs to be able to do on some level or raise an NotImplementedError. multiplying with ndarray is not one of the things that nddata does, afaik.

I understand that, this is a lot to do for NDError, but I don't have a better idea. That this is a lot also stems from fact that NDData is one of the pilars that astropy will stand on.

@kbarbary when you do slices you are automatically creating a new object anyways. What we need to make sure is that NDError doesn't perform lots of checks each time it is instantiated

I'm not thinking about operations on NDData that get passed to NDError. In aperture_photometry() we would be accessing the error attribute of an NDData object directly and Doing Things with it. These Things include slicing, and basically taking a weighted sum. Currently we do the weighted sum by multiplying by an array of floats and then summing. NDError will have to support operations that NDData doesn't necessarily support.

Regarding slicing : I guess that's true, but this is creating another new object in addition to the new ndarray. Still, you're right, it probably isn't much overhead as long as there aren't a lot of checks.

I see the draw of having these classes, but it seems like a lot of effort and complication and I don't know the use cases for non-Gaussian or Poisson distributed errors. Perhaps someone can give a good example??

Personally, my solution to this would be to say "we only support Gaussian-distributed errors", just have NDData.error be an ndarray. Poisson error could also be represented with something analagous to a gain parameter.

I mean poisson distributed error are when we only have a few counts in x-rays for example. I think you're thinking very much along the line of nddata being a ccd image, which is just one of the many use cases.

Anyway, I don't think it is that complicated as you think: You're right many operations that photutils does can only be done (I guess) using gaussian errors. For example, (and I think that might be the example that you are talking about) you want to fit a 2d-gaussian using a chi-squared method. In that case nothing, but gaussian errors, make sense and in that case you can easily access the array. I think many cases that you're thinking of can only be done with gaussian errors and thus at the beginning of the function it will just say isinstance(nddata.error, SDError) and if not true we either raise a warning and don't take errors into account or raise an exception and tell the user that for this function any other error but SDError doesn't make sense.

I think that is a very sensible approach actually: In the SDError class we check if the shape is the same or if it's a scalar. Then in either case you can access the data directly.

what do you think?

Sorry I meant "I don't see the use cases for non-[Gaussian or Poisson] distributed errors". (You're right I'm mainly thinking about CCD data, but of course I know Poisson errors are a common case.) Are there others besides Gaussian and Poisson you have in mind?

You can store a whole probability distribution function as an error.

Yes I understand that, but is it a common thing to do? What is a specific use case in astronomy? Perhaps people do this for spectra? I've never seen or heard of it for imaging data.

I'm pushing on this because supporting this will most likely make the code much more complex, so we'd better have a good reason to do it.

My aim in this is to give the users a framework that they can build things with that we never envisioned, and I think errors is a big part of it as no other package properly supports that. But as you rightly pointed out I think that other errors unless they have specific functions don't make sense.

So you can just require your image class to have SDErrors and things become incredibly easy.

Regarding storing whole PDFs, people don't usually do it because they can't, but I'm sure there are use cases for them (and it's the right thing to do in some cases, e.g. if datasets are combined non-linearly). One doesn't necessarily need to be able to combine all different kinds of errors, but normal and poisson errors are only a very specific case and I do see the value in having classes to store errors.

I don't disagree with @astrofrog - I'm not necessarily against an error class. What I was getting at a month ago is that for a given implementation, we should carefully think through what the use of such a class would look like as part of NDData-derived classes (e.g., an Image class, using photutils.aperture functions as a real-world example). If such use turns out to be more complicated, less transparent, or significantly slower than simply using a numpy.ndarray to represent errors, I feel we are doing something wrong.

see #348