Hi, On Sat, Oct 18, 2014 at 6:17 PM, Nathaniel Smith <n...@pobox.com> wrote: > Okay! I think I now actually understand what was going on with > np.gradient! The discussion has been pretty confused, and I'm worried > that what's in master right now is not the right solution, which is a > problem because we're supposed to cut 1.9.1 tomorrow. > > Background: > > np.gradient computes gradients using a finite difference method. > Wikipedia has a good explanation of how this works: > https://en.wikipedia.org/wiki/Numerical_differentiation > The key point is that there are multiple formulas one can use that all > converge to the derivative, e.g. > > (f(x + h) - f(x)) / h > > or > > (f(x + h) - f(x - h)) / 2h > > The first is the textbook definition of the derivative. As h -> 0, the > error in that approximation shrinks like O(h). The second formula also > converges to the derivative as h -> 0, but it converges like O(h^2), > i.e., much faster. And there's are many many formulas like this with > different trade-offs: > https://en.wikipedia.org/wiki/Finite_difference_coefficient > In practice, given a grid of values of f(x) and a grid stepsize of h, > all of these formulas come down to doing a simple convolution of the > data with certain coefficients (e.g. [-1/h, 1/h] or [-1/2h, 0, 1/2h] > for the two formulas above). > > Traditionally np.gradient has used the second formula above, with its > quadratically diminishing error term, for interior points in the grid. > For points on the boundary of the grid, though, this formula has a > problem, b/c it requires you look at points to both the left and the > right of the current point -- but for boundary points one of these > doesn't exist. In such situations np.gradient has traditionally used > the first formula instead (the "forward (or backward) finite > difference approximation with linear error", where > "forward"/"backward" means that it works on the boundary). As the web > page linked above shows, though, there's an easy alternative formula > that works on
Did you lose some text here? > The change in numpy 1.9.0 that has caused all the fuss, is that we > switched to using the quadratically accurate approximation instead of > the linearly accurate approximation for boundary points. > > This had two effects: > > 1) it produced (hopefully small) changes in the calculation for the > boundaries. In general these should have made the calculations more > accurate, but apparently this did cause problems for some people who > wanted perfectly reproducible output. > > 2) it tickled a nasty bug in how gradient and masked arrays work > together. It turns out gradient applied to masked arrays has always > been seriously buggy; there's some horribleness in the masked array > implementation that means the output array that np.gradient is writing > into ends up sharing a mask array with the input. So as we go along > calculating the output, we also mutate the input. This can definitely > corrupt the input, for example here: > https://gist.github.com/njsmith/551738469b74d175e039 > ...and I strongly suspect it's possible to find cases where it means > you get the wrong answer. > > For mysterious reasons involving the order in which values were > calculated or something, the version of np.gradient in 1.9.0 tickled > this bug in a much worse manner, with early mutations to the input > array ended up affecting later calculations using the same input > array, causing cascading errors. This is the cause of the massive > problems that matplotlib encountered in their test suite: > https://github.com/matplotlib/matplotlib/issues/3598#issuecomment-58859663 > > There's a proposed fix for the masked array issues at: > https://github.com/numpy/numpy/pull/5203 > For the reasons described above, I would be very wary about using > *any* currently released version of np.gradient on masked arrays. A > workaround would be to make a local copy of the definition of > np.gradient, and immediately after the 'out' array is allocated do > 'out.mask = out.mask.copy()'. > > Once this bug is fixed, the "new" gradient code produces results which > are identical to the "old" gradient code on the interior; only the > boundary is different. > > ------ > > So here are my concerns: > > - We decided to revert the changes to np.gradient in 1.9.1 (at least > by default). I'm not sure how much of that decision was based on the > issues with masked arrays, though, which turns out to be a different > issue entirely. The previous discussion conflated the two issues > above. > > - We decided to gate the more-accurate boundary calculation behind a > kwarg called "edge_order=2". AFAICT now that I actually understand > what this code is doing, this is a terrible name -- we're using a > 3-coefficient kernel to compute a quadratically accurate approximation > to a first-order derivative. There probably exist other kernels that > are also quadratically accurate. "Order 2" simply doesn't refer to > this in any unique or meaningful way. And it will be even more > confusing if we ever add the option to select which kernel to use on > the interior, where "quadratically accurate" is definitely not enough > information to uniquely define a kernel. Maybe we want something like > edge_accuracy=2 or edge_accuracy="quadratic" or something? I'm not > sure. > > - If edge_order=2 escapes into a release tomorrow then we'll be stuck with it. > > Some possible options for 1.9.1: > - Just revert the entire np.gradient code to match 1.8, and put off a > real fix until 1.9. This at least avoids getting stuck with a poor > API. > > - Put the masked array bug fix into 1.9.1, and leave the gradient code > the same as 1.9.0; put off designing a proper API for selecting the > old behavior for 1.9.2. (Or 1.10 I guess.) This doesn't solve the pure > reproduceability problems, but it might be enough to make matplotlib > work again at least? > > - Delay 1.9.1 to rehash the issues and make a decision on what we want > to support long-term. Excellent summary, thanks for doing all the work to get to the bottom of the problem. Is there any significant disadvantage to delaying the 1.9.1 release for a few days? I would appreciate a short delay to see if it's possible to agree a workaround for the OSX Accelerate bug : https://github.com/numpy/numpy/pull/5205 - but in any case, unless there is a compelling reason to release tomorrow, it seems reasonable to take the time to agree on a good solution to the np.gradient API. Cheers, Matthew _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
