Incorrect projection_weights results due to floating point rounding error

[goretkin]

I am getting incorrect results from the package with the latest version (Manifest.toml and Project.toml are attached). I haven't been able to figure out why, but projection_weights gives different answers from what it used to return.

To make things really fun, when I instrument the projection_weights_impl to print or track values, it seems to give the correct result.

import GeometryTypes: Point, FlexibleConvexHull
import EnhancedGJK

robot_piece = FlexibleConvexHull(Point{2,Float64}[[1.125, 2.0], [0.125, 2.0], [0.125, 2.8], [0.325, 3.0], [0.925, 3.0], [1.125, 2.8]])

world_peace = FlexibleConvexHull(Point{2,Float64}[Point(-0.025, 1.2321067811865476), Point(-0.025, 1.2821067811865474), Point(0.3267766952966369, 1.2821067811865474), Point(0.3267766952966369, 1.2321067811865476)])

result = EnhancedGJK.gjk(robot_piece, world_peace)
@show result.signed_distance

geomA = robot_piece
geomB = world_peace

cache = EnhancedGJK.CollisionCache(geomA, geomB)
poseA = poseB = EnhancedGJK.IdentityTransformation()
simplex = EnhancedGJK.transform_simplex(cache, poseA, poseB)
weights1 = EnhancedGJK.projection_weights(simplex)
weights2 = EnhancedGJK.projection_weights_reference(simplex)

@show weights1 weights2
# Geometries do not collide.
#=
import Plots

plt = Plots.plot(;aspect_ratio=1.0, reuse=false)

Plots.scatter!(plt, robot_piece.vertices, color=:red)
Plots.scatter!(plt, world_peace.vertices, color=:blue)
display(plt)

plt_gjk = Plots.plot(;aspect_ratio=1.0, reuse=false)
Plots.plot!(plt_gjk, Plots.Shape(collect(map(x->tuple(x...), result.simplex))))
=#

Running the above file gives the incorrect result

result.signed_distance = -0.7248609217812854
weights1 = [0.5, 0.5, 0.0]
weights2 = [1.0, 0.0, 0.0]
3-element StaticArrays.MArray{Tuple{3},Float64,1,3} with indices SOneTo(3):
 1.0
 0.0
 0.0

This is happening on both Julia 1.1.0 and 1.2rc.

Here are the Manifest.toml and Project.toml
BadGJK.zip

[tkoolen]

Starting to look at this. In 845e5ca I added a failing test case based on your example. One interesting thing is that it matters whether that test set is the first in the runtests.jl file or the last. There's definitely something funky going on.

[tkoolen]

Test also fails on Julia master (Version 1.3.0-DEV.435, JuliaLang/julia@75c10e435b). Also tried a debug build of Julia (make JULIA_PRECOMPILE=0 FORCE_ASSERTIONS=1 LLVM_DEBUG=Release LLVM_ASSERTIONS=1 -j 10), which didn't reveal anything.

[tkoolen]

Everything is OK with -O0 or -O1; bad with -O2. I suspect this is a SIMD optimization issue (I've seen something similar when I annotated a for loop with @simd while not satisfying the requirements for @simd.

[tkoolen]

So if you change the dot calls in projection_weights to mydot, defined by

@inline function mydot(x::SVector{2}, y::SVector{2})
    @inbounds return x[1] * y[1] + x[2] * y[2]
end

then everything is fine. But this function produces the exact same machine instructions as the StaticArrays dot implementation (which does have an @simd loop).

[tkoolen]

Maybe it's not a corruption/miscompilation issue. It could just be the fact that there's a floating point comparison to 0 here:

EnhancedGJK.jl/src/johnson_distance.jl

Line 96 in 59a9200

if d > 0

and with the SIMD loop you get a slightly different result for the dot products, which end up being close to zero for this example (e.g. 2.220446049250313e-16 instead of 0). Changing the 0 to 1e-15 fixes the issue also with the SIMD StaticArrays dot.

[tkoolen]

Also the mydot change from #17 (comment) does change the native code for johnson_distance, so it's not that mysterious.

[goretkin]

Thanks for figuring that out! I understand now that my debugging statements were preventing the SIMD instructions from being emitted. I did see the 2.220446049250313e-16 vs 0.0 at some point, but when I added more tracing, it disappeared. Maybe next time I'll think SIMD!

[tkoolen]

Thank you for the reproducer and the initial digging.

[rdeits]

Ahh, good catch. Sorry for my sloppy floating-point code!