I am trying to clean up our floating point warning handling:
And an upcoming PR to remove most floating point error clearing. There
are some things I am unsure about, though.
Part of why it got so confusing, is that GCC seemed to have
fixed/changed their behaviour for comparisons with NaN. In GCC 7 it
did not give the warning, but GCC 8 does (correctly).
Comparison with NaN
IEEE says that the default comparisons should give warnings for
comparison with NaN (except == and !=). And notes that an alternative
should be provided (C99 does this with `isless`, etc.).
We currently break this by suppressing invalid value warnings for all
comparisons (e.g. also `NaN > 0.`).
We can easily do either version (aside possibly compiler issues).
Making it give warnings had one test case fail for `quantile`, which
uses the pattern:
if not (np.all(q >= 0) and np.all(q <= 1)):
raise ValueError("bad q")
This would additionally (and first) give an "invalid value" warning and
require `np.errstate(invalid="ignore") to suppress it.
I dislike diverging from IEEE, but Python also does not warn for :
float("nan") >= 0
and presumably the user either explicitly created the NaN or has seen a
warning earlier during computation when the NaN was first created.
(IEEE does not distinguish creating a new NaN with `0./0.` from a
comparison with `NaN > 0.` . So we can't easily make this settable
via `np.errstate` or so.)
So, should we ignore the warning here?
Some compilers may get flags wrong. How much effort do we want to
spend on details few users will notice?
My current problem is `1 % 0` and `divmod(1, 0)`. The MacOS/clang CI
does not set the correct "invalid value" warning flag. (The remainder
is NaN, so a new NaN is created and that should be indicated but the
C99 `fmod` does not set it.)
I propose dropping any special concern for signalling NaNs. Which
means they raise almost always. Although, rarely we might suppress the
warning if we do it manually for normal NaNs .
We have two tests which check for behaviour on signalling NaNs. I could
not find having any logic to them besides someone being surprised at
signalling NaN behaviour at the time – not based on use-cases.
Even functions like `isnan` give a warning for signalling NaNs!
The "fix" for anyone having sNaN's is to convert them to qNaNs as early
as possible. Which e.g. `np.positive(arr, out=arr)` should probably
do. If this becomes an issue, maybe we could have an explicit ufunc.
 Mainly it seems SSE2 does not provide some non-error comparisons.
So trying to avoid manually clearing errors might make some SSE code
considerable slower (e.g. `isfinite`, `np.min`).
 Probably Python just does not check the CPU warning flags
There will be a NumPy Community meeting Wednesday July 7th at
20:00 UTC. Everyone is invited and encouraged to
join in and edit the work-in-progress meeting topics and notes at:
DNV Energy USA is looking for an experienced software engineer to help
accelerate the renewable energy transition. Do you know any software
engineers interested in clean energy?
Would you mind sharing the following link with your network?
Mark A. Mikofski
I'm writing some helper Cythpm functions for scipy.linalg which is kinda
performant and usable. And there is still quite some wiggle room for more.
In many linalg routines there is a lot of performance benefit if the
structure can be discovered in a cheap and reliable way at the outset. For
example if symmetric then eig can delegate to eigh or if triangular then
triangular solvers can be used in linalg.solve and lstsq so forth
Here is the Cythonized version for Jupyter notebook to paste to discover
the lower/upper bandwidth of square array A that competes well with A != 0
just to use some low level function (note the latter returns an array hence
more cost is involved) There is a higher level supervisor function that
checks C-contiguousness otherwise specializes to different versions of it
Then another cell
# cython: language_level=3
# cython: linetrace=True
# cython: binding = True
# distutils: define_macros=CYTHON_TRACE=1
# distutils: define_macros=CYTHON_TRACE_NOGIL=1
cimport numpy as cnp
import numpy as np
ctypedef fused np_numeric_t:
cpdef inline (int, int) band_check_internal(np_numeric_t[:, ::1]A):
cdef Py_ssize_t n = A.shape, lower_band = 0, upper_band = 0, r, c
cdef np_numeric_t zero = 0
for r in xrange(n):
# Only bother if outside the existing band:
for c in xrange(r-lower_band):
if A[r, c] != zero:
lower_band = r - c
for c in xrange(n - 1, r + upper_band, -1):
if A[r, c] != zero:
upper_band = c - r
return lower_band, upper_band
Final cell for use-case ---------------
# Make arbitrary lower-banded array
n = 50 # array size
k = 3 # k'th subdiagonal
R = np.zeros([n, n], dtype=np.float32)
R[[x for x in range(n)], [x for x in range(n)]] = 1
R[[x for x in range(n-1)], [x for x in range(1,n)]] = 1
R[[x for x in range(1,n)], [x for x in range(n-1)]] = 1
R[[x for x in range(k,n)], [x for x in range(n-k)]] = 2
Some very haphazardly put together metrics
2.59 µs ± 84.7 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit np.linalg.solve(R, zzz)
824 µs ± 6.24 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit R != 0.
1.65 µs ± 43.1 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
So the worst case cost is negligible in general (note that the given code
is slower as it uses the fused type however if I go with tempita standalone
version is faster)
1) This is missing np.half/float16 functionality since any arithmetic with
float16 is might not be reliable including nonzero check. IS it safe to
view it as np.uint16 and use that specialization? I'm not sure about the
sign bit hence the question. I can leave this out since almost all linalg
suite rejects this datatype due to well-known lack of supprt.
2) Should this be in NumPy or SciPy linalg? It is quite relevant to be on
SciPy but then again this stuff is purely about array structures. But if
the opinion is for NumPy then I would need a volunteer because NumPy
codebase flies way above my head.
All feedback welcome
The PR https://github.com/numpy/numpy/pull/19211 proposes to extend
argmin and argmax with a `keepdims=False` keyword-only argument.
This is a standard argument in NumPy, so it is a small API addition.
The PR also proposes to add:
in the C-API. We have barely extended the C-API in a very long time,
so if anyone has concerns, we could pull that out again .
Otherwise, this should go in soon, and we will have `keepdims` for both
of those functions in the next release :).
 I do not see this is much of a maintenance concern, since the
original function is just a one-line wrapper of the new one.
The API is fairly large and it probably is not used much. So it doesn't
feel important to add to me. Overally, I just don't have a preference.