USNavalResearchLaboratory / eispac

Read the Docs
https://eispac.readthedocs.io/en/latest/
MIT License
23 stars 6 forks source link

NDcube "rebin" method applied to eiscube replaces uncertainty array with None #81

Open adolliou opened 1 year ago

adolliou commented 1 year ago

Dear all,

I tried to apply a spatial binning over the vertical axis before fitting the data. I tried the "rebin" method from the NDcube class, but it replaces the uncertainty array by none, making it impossible to perform a fitting afterwards. I suspect it might have something to do with the fact that uncertainty_type is not defined.

I wanted to know if there is a specific way to perform a binning with eispac ?

Thank you very much for your help.

As an example :

import eispac

if __name__ == '__main__':
    data_filepath = 'eis_20230404_065652.data.h5'
    template_filepath = 'fe_12_195_119.2c.template.h5'

    data_cube = eispac.read_cube(data_filepath, window=195.119)
    print(f'{data_cube.data.shape=}')
    print(f'{data_cube.uncertainty=}')
    shape = (2, 1, 1)
    print(f'{data_cube_binned.data.shape=}')
    data_cube_binned = data_cube.rebin(shape)
    print(f'{data_cube_binned.uncertainty=}')

returns :

data_cube.data.shape=(152, 25, 32) data_cube.uncertainty=UnknownUncertainty([[[ 34.06817081, 53.28593618, 49.9920177 , ..., 43.90691806, 43.88902026, -2112.19711304], ... [ 42.82511634, 40.77855076, 51.64254604, ..., 55.3329251 , 43.88857749, 53.80620185]]]) data_cube_binned.data.shape=(76, 25, 32) data_cube_binned.uncertainty=None

MJWeberg commented 1 year ago

Hello,

I am terribly sorry for the long delay! I somehow missed the notification for this issue.

It looks like the NDCube.rebin() method does NOT propagate uncertainties by default. Setting the keyword propagate_uncertainties=True should hopefully fix this issue. If not, you could alternatively use the EISCube.smooth_cube() method to apply a boxcar average to the EIS data using either pixel or real-world widths. At some point, we will probably depreciate smooth_cube(), since we originally added to a few years ago before NDCube had any native rebin function. However, we will need to test our code a bit more before then, to ensure the EISCube.wavelength array is copied correctly with NDCube.rebin

As as side note, the uncertainty type in EISCube was initially left undefined since, at the time, there was no appropriate type for instrumental uncertainties (being, as they are, neither the standard deviation nor variance of the actual data values). However, I see that Astropy now allows for arbitrary uncertainty arrays. We will take a look into them and see if we can use them for the EISCubes.

Thanks for using eispac! Please do not hesitate to let us know if you encounter any issues or have suggestions for how we can improve the code or documentation.

adolliou commented 1 year ago

Hi,

Thank you very much for your help!

Setting propagate_uncertainties=True to the rebin method returns the following warning : "UserWarning: self.uncertainty is of type UnknownUncertainty which does not support uncertainty propagation."

I will look into the smooth_cube function.

Bests,

DanRyanIrish commented 1 year ago

Hi @adolliou. This issue was mentioned in the ndcube matrix chat last week but I was on vacation. @MJWeberg is right in answering your first question. NDCube.rebin does not propagate uncertainties by default because it is slow, especially for large cubes. Setting propagate_uncertainties=True is the simple solution.

HOWEVER, this will only work so long as NONE of the following conditions true:

These conditions are checked in the NDCube.rebin code here.

A Quick Note on Uncertainty Types

Astropy currently supports two known uncertainty types: StdDevUncertainty and VarianceUncertainty. If a simple array is supplied to NDCube for the uncertainty attribute, then it is automatically converted to an UnknownUncertainty. I am not familiar with the eispac code-base, but based on @MJWeberg description above, it sounds like this is what's done. I believe this is why you are getting the warning you cited above.

Suggested Solutions

  1. (Simplest if allowed by eispac): If eispac does not force the uncertainty to be of type UnknownUncertainty, you could explicitly supply the uncertainties in a astropy.nddata.StdDevUncertainty or astropy.nddata.VarianceUncertainty object, and set propagate_uncertainties=True. However, according to @MJWeberg this might not be the strictly correct way if the uncertainties are not actually simple standard deviations or variances.
  2. (Advanced): The propagate_uncertainties= kwarg can be set to a function as well as a bool. This allows you to define how to propagate the uncertainties however you want. The default propagation function is defined here and can be used as a template. Alternatively you can implement the function however you want so long as it takes the same inputs as the default function and returns an uncertainty object.

I would recommend option 2, although it does require a little detailed knowledge of how NDCube.rebin works. Given that, I think it would be most helpful for users if the function was supplied by eispac. Perhaps it's something that @adolliou and @MJWeberg can work on together. I would be happy to advise since it will need a little ndcube expertise. Additionally, if you want to propagate unknown uncertainties, it would require a minor change to NDCube.rebin. I would be happy to work with either of you to implement that. I've raised an issue ndcube#639 describing what needs to be done.

MJWeberg commented 1 year ago

Thanks for checking in on this issue, @DanRyanIrish!

I have not had the time to work on this recently, but I would be happy to look into adding a proper propagation function in a week or two after I get back from a conference.

On this topic, would it make sense to add some kind of MeasuredUncertainty or EstimatedUncertainty class to either NDCube or nddata? As I am sure you are aware, instrumental uncertainties can have a variety of sources beyond just the statistical (e.g. CCD read noise, internal dispersion, calibration uncertainties, ect...). I would be surprised if EIS was the only instrument to have this issue. However, it might just be that we can treat most instrument uncertainties as StdDev for the purposes of propaganda and we are being unnecessarily particular here.