Skip to content
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

test_velocity_in_ITRF_to_GCRS2 too tight for 32-bit Intel #411

Closed
bnavigator opened this issue Jul 22, 2020 · 10 comments
Closed

test_velocity_in_ITRF_to_GCRS2 too tight for 32-bit Intel #411

bnavigator opened this issue Jul 22, 2020 · 10 comments

Comments

@bnavigator
Copy link

In relation to #404: 32-bit intel is not precise enough for this test

[   93s] =================================== FAILURES ===================================
[   93s] ________________________ test_velocity_in_ITRF_to_GCRS2 ________________________
[   93s] 
[   93s]     def test_velocity_in_ITRF_to_GCRS2():
[   93s]         # TODO: Get test working with these vectors too, showing it works
[   93s]         # with a non-zero velocity vector, but in that case the test will
[   93s]         # have to be fancier in how it corrects.
[   93s]         # r = np.array([(1, 0, 0), (1, 1 / DAY_S, 0)]).T
[   93s]         # v = np.array([(0, 1, 0), (0, 1, 0)]).T
[   93s]     
[   93s]         ts = api.load.timescale(builtin=True)
[   93s]         t = ts.utc(2020, 7, 17, 8, 51, [0, 1])
[   93s]         r = np.array([(1, 0, 0), (1, 0, 0)]).T
[   93s]         v = np.array([(0, 0, 0), (0, 0, 0)]).T
[   93s]     
[   93s]         r, v = ITRF_to_GCRS2(t, r, v, True)
[   93s]     
[   93s]         # Rotate back to equinox-of-date before applying correction.
[   93s]         r = mxv(t.M, r)
[   93s]         v = mxv(t.M, v)
[   93s]     
[   93s]         r0, r1 = r.T
[   93s]         v0 = v[:,0]
[   93s]     
[   93s]         # Apply a correction: the instantaneous velocity does not in fact
[   93s]         # carry the position in a straight line, but in an arc around the
[   93s]         # origin; so use trigonometry to move the destination point to where
[   93s]         # linear motion would have carried it.
[   93s]         angvel = (t.gast[1] - t.gast[0]) / 24.0 * tau
[   93s]         r1 = mxv(rot_z(np.arctan(angvel) - angvel), r1)
[   93s]         r1 *= np.sqrt(1 + angvel*angvel)
[   93s]     
[   93s]         actual_motion = r1 - r0
[   93s]         predicted_motion = v0 / DAY_S
[   93s]     
[   93s]         relative_error = (length_of(actual_motion - predicted_motion)
[   93s]                           / length_of(actual_motion))
[   93s]     
[   93s] >       assert relative_error < 2e-12
[   93s] E       assert 3.866305785983307e-12 < 2e-12
[   93s] 
[   93s] skyfield/tests/test_positions.py:129: AssertionError
[   93s] =========================== short test summary info ============================
[   93s] FAILED skyfield/tests/test_positions.py::test_velocity_in_ITRF_to_GCRS2 - ass...
[   93s] ======================== 1 failed, 452 passed in 14.92s ========================
@brandon-rhodes
Copy link
Member

Interesting! I had expected IEEE 64-bit floats to behave the same whether the integers on a system were 32 bits or 64 bits. But maybe 32-bit processors did not have as many floating point registers, so some operations which a modern processor carries through at 80 bits precision drop to only 64 bits on an older processor?

I’ll see whether Travis CI can run Skyfield in a 32-bit environment, to give me a way to reproduce this problem. Alternately I guess I can look for a way to boot a 32-bit machine locally?

@bnavigator
Copy link
Author

These logs come from the openSUSE Build Service. I don't think they run 32-bit metal. Those are just virtual machines with 32-bit libraries installed.

@brandon-rhodes
Copy link
Member

The odd thing is: this bound isn’t even very tight! Because the velocity uses Earth's average angular motion rather than the true angular motion as of that date and time, the result has an error about 10⁴ above the noise floor of the floating point math itself. But improving the result by using the true angular velocity make Skyfield twice as slow, alas!

For my own reference: I do see the error you quote, by clicking on the i586 link at:
https://build.opensuse.org/package/show/devel:languages:python:numeric/python-skyfield

@brandon-rhodes
Copy link
Member

@bnavigator — It would be interesting to see by how much the assert inequality is missing its target. Do we have any hint as to why pytest’s assert introspection isn't working in the SUSE builds?

https://docs.pytest.org/en/3.0.1/assert.html

It should be on by default, unless SUSE ships a modified version of the library that turns it off?

@bnavigator
Copy link
Author

I don't understand?

[   93s] E       assert 3.866305785983307e-12 < 2e-12

what else could be rewritten? The error is off by a factor of ~2. So one bit.

FYI, the content of the link will change, because I am easing the tolerance for the package.

@brandon-rhodes
Copy link
Member

Thanks! My eye was down on the AssertionError — my eye ignored everything else as "the code where the error is", which I didn't expect would include the failed numbers, so thanks for pointing them out!

@brandon-rhodes
Copy link
Member

I have just committed 84951b3 which gives me my own Ubuntu 32-bit test container. Exactly as you suggested, my modern 64-bit processor was able to reproduce the error so long as it was running code compiled with a 32-bit compiler! It must be avoiding features that lead to higher precision when the modern features of a 64-bit processor can be assumed.

Now that I can reproduce the problem locally, I'll take a look — probably tomorrow, given my schedule — and see if there are any actual problems under a 32-bit processor that I can fix to restore the precision, since, again, this bound is not actually very tight compared to the machine precision available. If I can't make any headway, then I'll punt and simply reduce the precision here if a 32-bit Python is detected.

Either way, we'll hopefully get stock Skyfield where tests can pass on a 32-bit platform!

@bnavigator
Copy link
Author

bnavigator commented Jul 22, 2020

Some debugging:

def test_velocity_in_ITRF_to_GCRS2():
    # TODO: Get test working with these vectors too, showing it works
    # with a non-zero velocity vector, but in that case the test will
    # have to be fancier in how it corrects.
    # r = np.array([(1, 0, 0), (1, 1 / DAY_S, 0)]).T
    # v = np.array([(0, 1, 0), (0, 1, 0)]).T

    ts = api.load.timescale(builtin=True)
    t = ts.utc(2020, 7, 17, 8, 51, [0, 1])
    r = np.array([(1, 0, 0), (1, 0, 0)]).T
    v = np.array([(0, 0, 0), (0, 0, 0)]).T

    r, v = ITRF_to_GCRS2(t, r, v, True)

    # Rotate back to equinox-of-date before applying correction.
    r = mxv(t.M, r)
    v = mxv(t.M, v)

    r0, r1 = r.T
    v0 = v[:,0]

    # Apply a correction: the instantaneous velocity does not in fact
    # carry the position in a straight line, but in an arc around the
    # origin; so use trigonometry to move the destination point to where
    # linear motion would have carried it.
    angvel = (t.gast[1] - t.gast[0]) / 24.0 * tau
    r1 = mxv(rot_z(np.arctan(angvel) - angvel), r1)
    r1 *= np.sqrt(1 + angvel*angvel)

    actual_motion = r1 - r0
    predicted_motion = v0 / DAY_S
    motion_difference = actual_motion - predicted_motion

    print(f"r1 dtype: {r1.dtype}")
    info = np.finfo(r1.dtype)
    for attr in ["eps", "epsneg", "bits", "precision"]:
        print("{:>10s}: {}".format(attr, getattr(info, attr)))

    print("")
    print("                               "
          + ", ".join([" 1 23456789ABCDEF     "]*3))
    print("           r0                = "
          + ", ".join(["{: .15e}"]*3).format(*r0))
    print("           r1                = "
          + ", ".join(["{: .15e}"]*3).format(*r1))
    print("actual_motion = r1- r0       = "
          + ", ".join(["{: .15e}"]*3).format(*actual_motion))
    print("predicted_motion             = "
          + ", ".join(["{: .15e}"]*3).format(*predicted_motion))
    print("motion_difference            = "
          + ", ".join(["{: .15e}"]*3).format(*motion_difference))

    print("lenght_of(motion_difference) = "
          + "{: .15e}".format(length_of(motion_difference)))
    print("lenght_of(actual_motion)     = "
          + "{: .15e}".format(length_of(actual_motion)))

    relative_error = (length_of(motion_difference)
                      / length_of(actual_motion))

    print("relative_error               = "
          + "{: .15e}".format(relative_error))

    # fail for inspection
    assert False

x86_64:

[   12s] ----------------------------- Captured stdout call -----------------------------
[   12s] r1 dtype: float64
[   12s]        eps: 2.220446049250313e-16
[   12s]     epsneg: 1.1102230246251565e-16
[   12s]       bits: 64
[   12s]  precision: 15
[   12s] 
[   12s]                                 1 23456789ABCDEF     ,  1 23456789ABCDEF     ,  1 23456789ABCDEF     
[   12s]            r0                =  3.683658935252279e-01,  9.296808960537827e-01, -4.356432749264005e-15
[   12s]            r1                =  3.682981001133219e-01,  9.297077577230478e-01, -4.410426029179994e-15
[   12s] actual_motion = r1- r0       = -6.779341190599197e-05,  2.686166926513245e-05, -5.399327991598880e-17
[   12s] predicted_motion             = -6.779341190595400e-05,  2.686166926508038e-05, -3.113266460474117e-19
[   12s] motion_difference            = -3.797418109130479e-17,  5.206880933361635e-17, -5.368195326994138e-17
[   12s] lenght_of(motion_difference) =  8.387461738721598e-17
[   12s] lenght_of(actual_motion)     =  7.292116272773533e-05
[   12s] relative_error               =  1.150209544798091e-12

i586:

[   15s] ----------------------------- Captured stdout call -----------------------------
[   15s] r1 dtype: float64
[   15s]        eps: 2.220446049250313e-16
[   15s]     epsneg: 1.1102230246251565e-16
[   15s]       bits: 64
[   15s]  precision: 15
[   15s] 
[   15s]                                 1 23456789ABCDEF     ,  1 23456789ABCDEF     ,  1 23456789ABCDEF     
[   15s]            r0                =  3.683658935252279e-01,  9.296808960537826e-01, -4.356193568335524e-15
[   15s]            r1                =  3.682981001133219e-01,  9.297077577230479e-01, -4.410459963437534e-15
[   15s] actual_motion = r1- r0       = -6.779341190599197e-05,  2.686166926535449e-05, -5.426639510200982e-17
[   15s] predicted_motion             = -6.779341190595400e-05,  2.686166926508039e-05, -3.113556745723749e-19
[   15s] motion_difference            = -3.797418109130479e-17,  2.741066379950696e-16, -5.395503942743745e-17
[   15s] lenght_of(motion_difference) =  2.819355133751896e-16
[   15s] lenght_of(actual_motion)     =  7.292116272781712e-05
[   15s] relative_error               =  3.866305785983307e-12

more: https://build.opensuse.org/package/show/home:bnavigator:skyfield-debug/python-skyfield

Looks like error propagation starting from loss of significance by r1-r0 and getting worse with every difference you take.
You really should not expect motion_difference to be more precise than epsneg, even on 64-bit systems.

@brandon-rhodes
Copy link
Member

I have adjusted the test tolerances so that they should now all pass under 32-bit Linux! I do not have anything running under 32-bit in CI, but I now have a local command I can run before each Skyfield release (hopefully I remember!) that should hopefully keep the tests of future versions also running under 32-bit operating systems.

I plan to release a new Skyfield 1.25 within the next few hours that includes the improvement. Once you have the chance to try it out, let me know if you have any further problems that need to be tracked down for your openSUSE build!

@bnavigator
Copy link
Author

Can confirm, tests now pass using assay on i586

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants