-
Notifications
You must be signed in to change notification settings - Fork 13
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
Inaccuracy of some operations (e.g. stacked average) when underlying data is float32 #217
Comments
@chuaxr we do not do any explicit typecasting in |
@spencerkclark To follow up on our discussion yesterday (see /home/xrc/WRF/aospy_code/float_issue.ipynb for the actual calculations):
Either way, it seems that typecasting all relevant variables to float would circumvent these problems. However, typecasting the entire dataset is inefficient given that all variables in the dataset are typecasted when only one might be used. Perhaps variables could be cast to float in _sel_var in data_loader.py? Re your comment: my main aim is to avoid my results being contaminated by floating point issues, not that I would like to use float32 throughout. Separately, the output of reg.av is float64, although it is not clear when exactly the typecasting occurs. |
@chuaxr thanks, this seems much clearer now. The note in the numpy documentation is particularly helpful. It does seem that one way to address this would be in
In my reading of the xarray documentation for
However that does not seem to be the case:
So perhaps it seems that the only immediate fix would indeed be casting all variables to
It's remarkable it makes such a difference in your case (even in this numpy example the difference in means is really small), but your notebook makes it pretty apparent. @spencerahill do you have any thoughts here? |
@spencerkclark My issue can be reproduced as follows:
|
@chuaxr I just remembered an issue in xarray directly related to this: pydata/xarray#1346 The issue is that xarray uses another library (called bottleneck) which provides a set of optimized functions for computing things like means (and numpy does not use this library), and it has issues with the So long story short, I think for our purposes casting up to |
@chuaxr @spencerkclark thanks for taking what looks like a pretty deep dive to get this sorted out. Admittedly I have only skimmed the above discussion, but it seems that casting to float64 does indeed seem to be the best solution for this particular case, especially in light of the xarray/bottleneck problem that @spencerkclark nicely remembered. The question in my mind then is whether or not we want to do this automatically for all calculations, i.e. not rely on users including it themselves within their Either way, we should add a Warning in our docs about the potential for this problem with float32 data. |
Thanks for weighing in @spencerahill; I think this is the way to go. Considering that this impacts reduction operations, which can be used on any variable (including model-native ones), dealing with this within
Very much agreed. This is a tricky issue to diagnose! |
By the way @chuaxr, thanks for reporting this. I'm sorry that this puzzling issue impacted some of the figures you showed Yi, but I'm glad we got to the bottom of it (and it's something to keep in mind even outside of using |
Thanks @spencerkclark for your help in getting to the bottom of this, it makes sense that my problems are consistent with other users using xarray. I never meant to imply that this was a fault of aospy per se. Rather, I raised the issue here because I figured that a workaround would likely need to be placed within aospy. It's quite likely that this problem would have (or perhaps already have) manifested at some point, aospy or no. |
@chuaxr no worries even if you did imply that! We want to get to the bottom of any issue affecting any of our users, even if technically as in this case they don't come from our code. Because regardless, we need to develop a solution. We wouldn't have caught this one without your help. So keep the bug reports/feature requests/usage questions/pretty much anything else coming! 😄 |
If no data is missing in the underlying dataset, then the stacked average (#208) and reg.av should yield the same answer. (In my case, doing a reg.av automatically results in a time average, since all data is output within the same year).
However, this is not the case.
Using the full data file (1000 time steps at 20s, ~5h on a 96x96 horizontal grid), I obtain the following:
With only the first hour of data, the difference isn't as large:
Adding a typecast statement in stack_mean (see code below) appeared to fix that (see plot below), presumably because the model output is in float32 and reg.av does a typecast somewhere while the stacked average did not. Since any errors here contaminate all my analyses, I would very much appreciate knowing if there are other fixes I should make. For example, I may need to average over days at a very high frequency (e.g. 20s). In that case, would numpy float precision be enough, or would I have to average it within the model to make the size manageable?
Same as the first plot, except with modified stack_mean
The text was updated successfully, but these errors were encountered: