ENH: Set symmetric threshold for identifying unit quaternions in qform calculation #1182
Conversation
Codecov ReportBase: 92.17% // Head: 92.16% // Decreases project coverage by
Additional details and impacted files@@ Coverage Diff @@
## master #1182 +/- ##
==========================================
- Coverage 92.17% 92.16% -0.01%
==========================================
Files 97 97
Lines 12341 12337 -4
Branches 2535 2534 -1
==========================================
- Hits 11375 11371 -4
Misses 645 645
Partials 321 321
Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here. ☔ View full report at Codecov. |
nibabel/tests/test_quaternions.py
Outdated
| @@ -69,6 +69,20 @@ def test_fillpos(): | |||
| assert wxyz[0] == 0.0 | |||
|
|
|||
|
|
|||
| @pytest.mark.parametrize('dtype', ('f4', 'f8')) | |||
| def test_fillpos_precision(dtype): | |||
There was a problem hiding this comment.
Is the idea that one of these 50 random vectors will generate a small negative difference from zero? Worth testing that explicitly?
There was a problem hiding this comment.
Would this be easier to read with some example of a positive and negative error? Here we are leaving it to chance whether the test fails, are we not?
There was a problem hiding this comment.
Just to be clear, you're suggesting two vectors that are exemplars of positive and negative error? Would random generation work, as long as we verify that we get them?
def make_unit_vector(shape, dtype, constraint=lambda x: True, tries=50):
for _ in range(tries):
vec = rand.uniform(low=-1, high=1, size=shape).astype(dtype)
vec /= np.sqrt(vec @ vec)
if constraint(vec):
return vec
short = make_unit_vector((3,), dtype, lambda x: x @ x < 1)
long = make_unit_vector((3,), dtype, lambda x: x @ x > 1)
assert nq.fillpositive(short, w2_thresh)[0] == 0
assert nq.fillpositive(long, w2_thresh)[0] == 0There was a problem hiding this comment.
It seems to me cleaner to find some examplar vectors that we know give positive and negative < threshold differences, and a > threshold difference, and then test some random vectors - rather than hoping (even with reasonable expectation) that we'll get suitable vectors.
This comment was marked as resolved.
This comment was marked as resolved.
effigies
force-pushed
the
enh/qform_calculation
branch
2 times, most recently
from
5009ad5 to
9afb5dd
Compare
| minus = baseline * nptype(1 - np.finfo(dtype).eps) | ||
|
|
||
| assert nq.fillpositive(plus)[0] == 0.0 | ||
| assert nq.fillpositive(minus)[0] == 0.0 |
There was a problem hiding this comment.
@matthew-brett How's this for a deterministic test?
There was a problem hiding this comment.
Nice - but - consider pusing to threshold to show where it fails?
There was a problem hiding this comment.
Sure. It actually fails at 2*eps, since the error compounds. Updated.
| # Create random vectors, normalize to unit length, and count on floating point | ||
| # error to result in magnitudes larger/smaller than one | ||
| # This is to simulate cases where a unit quaternion with w == 0 would be encoded | ||
| # as xyz with small error, and we want to recover the w of 0 |
There was a problem hiding this comment.
I still like this as a bit of a stress test that will hit more realistic rounding errors, and added a comment to explain more. Can drop it if it seems excessive.
|
@matthew-brett Gentle bump on this, if you get a few minutes. |
| minus = baseline * nptype(1 - np.finfo(dtype).eps) | ||
|
|
||
| assert nq.fillpositive(plus)[0] == 0.0 | ||
| assert nq.fillpositive(minus)[0] == 0.0 |
There was a problem hiding this comment.
Nice - but - consider pusing to threshold to show where it fails?
nibabel/tests/test_quaternions.py
Outdated
| elif magnitude > 1: | ||
| neg_error = True | ||
|
|
||
| assert pos_error, 'Did not encounter a case where 1 - |xyz| > 0' |
There was a problem hiding this comment.
But do you still want to fail the tests here, on a random criterion?
There was a problem hiding this comment.
Good call. Just went ahead and checked the frequencies:
In [26]: vec = np.array([np.sqrt(np.sum(norm(gen_vec('f4'))**2)) for _ in range(10000)])
In [27]: np.unique(np.sign(1 - vec), return_counts=True)
Out[27]: (array([-1., 0., 1.], dtype=float32), array([ 81, 6946, 2973]))Seems negative errors are much less likely and we might have failed out for no good reason. I had assumed that the majority would be non-zero error and roughly split in either direction.
Also remove explicit check that we randomly generated positive and negative errors. Failing this check is unlikely, but not a bug.
effigies
force-pushed
the
enh/qform_calculation
branch
from
9afb5dd to
aa4b017
Compare
The
quaternions.fillpositivefunction is intended to augment a partial quaternion to make a unit 4-vector by findingw, given(x, y, z), such that|(w, x, y, z)| == 1. Floating point precision error means that a(w, x, y, z)selected to be a unit vector can still sometimes result in|(w, x, y, z)| != 1by a small amount. Currently if1 - |(x, y, z)| > -3*eps, we call that close enough to zero thatw = 0will still produce a unit quaternion.This is asymmetric, as if
1 - 3*eps < |(x, y, z)| < 1, we will getw = sqrt(1 - |(x, y, z)|) > 0, which can result in a surprising failure to preserve an affine that encodes a rotation in only one plane (see #1181 for an example). This is not exactly wrong, but the asymmetry is unsatisfying and inconsistent withfslhdat the least.I'm starting this PR with a (failing) test to assert that a passing
fillpositive()unit 3-vector, up to the precision allowed by its dtype, will result in a 0w.Following failing tests, I will use the change proposed in #1181 (comment) to satisfy the constraint imposed by the test.
Closes #1181.