How do i sort mathematical vectors by strict weak ordering for a map?

Question

I try to write a std::map< Vector3D, double > where colinear (parallel or anti-parallel) vectors should share the same key.

As a compare function I use the following function (with a 1e-9 tolerance in isEqualEnough()) that I created with the help of Using a (mathematical) vector in a std::map

struct Vector3DComparator 
{ 
    bool operator() (const Vector3D& lhsIn, const Vector3D& rhsIn) const
    {
        Vector3D lhs = lhsIn.absolute(); // make all members positive
        Vector3D rhs = rhsIn.absolute(); 

        if ((lhs.z < rhs.z)) 
            return true;

        if ((isEqualEnough(lhs.z, rhs.z)) 
            && (lhs.y < rhs.y)) 
            return true;

        if ((isEqualEnough(lhs.z, rhs.z)) 
            && (isEqualEnough(lhs.y, rhs.y))
            && (lhs.x < rhs.x))
            return true;

        return false;
    }
};

When I insert the normals of a cube into my map I should get 3 different values (because I don't care about the direction) but I get 4:

• x=1 y=0 z=0
• x=0 y=1 z=0
• x=0 y=0 z=1
• x=0 y=2.2e-16 z=1

The compare function is somehow wrong but whenever I try to fix it I get an assert telling me "Expression: invalid comparator".

Anyone spots the mistake?

Answer 1

It is mathematically impossible to use a tolerance with relational operators and yield a strict weak ordering. Any kind of convergence criterion will fail to satisfy ordering algorithms and data structures requirements. The reason is very simple: the incompatibility of two values, using a tolerance, doesn't yield an equivalence relation since it is not transitive. You may have almostEqual(a, b) and almostEqual(b, c) and yet ~almostEqual(a, c). Try this using a=1.0; b=2.0; c=3.0; tolerance=1.5;. You may look at this answer: Is floating-point == ever OK?.

You may still define an equivalence relation on floats using truncation, floor, roof, or round kind of functions. Let's define for example less3(a, b) if and only if floor(a * 8) < floor(b * 8) assuming a and b are binary floats and are not NAN and multiplications doesn't yield both the same signed infinite; this compares a and b using 3 bits of precision (0.125 in decimal). Now define equiv3(a, b) if and only if !less3(a, b) && ~less3(b, a). It can be shown that eqiv3(a, b) yields an appropriate equivalence relation. Since less3 is an order relation and equiv3 is an equivalence relation, then less3 is a strict weak order on floats (excluding NANs). Furthermore, in the case a * 8 == +INF && b * 8 == +INF || a * 8 == -INF && b * 8 == -INF you may fallback with ordinary < operator on floats.

Answer 2

Combining the answer from Julien Villemure-Fréchette with the link posted by @alterigel I made my comparison function work:

struct Vector3DComparator 
{ 
    bool operator() (const Vector3D& lhsIn, const Vector3D& rhsIn) const
    {
        int p = 100000; // precision factor

        Vector3D lhs = lhsIn.absolute(); // make all members positive
        Vector3D rhs = rhsIn.absolute();

        auto lhsTied = std::tie((int)(lhs.x * p), (int)(lhs.y * p), (int)(lhs.z * p));
        auto rhsTied = std::tie((int)(rhs.x * p), (int)(rhs.y * p), (int)(rhs.z * p));
        return lhsTied < rhsTied;
    }
};

Please note: This code contains bad style like c-style casts, bad naming and more. My functions and classes are different to the ones posted here. I just stripped everything away to make it easy to understand.

EDIT:

I noticed two more mistakes:

First: It didn't always work for nearly identical vectors. Based on @tetorea last comment on my question I changed the function to always get very similar values compared. I use the dot product for that as it is ±1 (or at least close to) for parallel vectors.

Second: .absolute() was not working because with this function two vectors (-1,1,0) and (1,1,0) were considered parallel which they are clearly not.

In the code below you can find the corrected version:

struct Vector3DComparator 
{ 
    bool operator() (const Vector3D& lhs, const Vector3D& rhs) const
    {
        if (isEqualEnough(fabs(lhs * rhs), 1.0)) // dot product
            return false;

        int p = 100000; // precision factor

        auto lhsTied = std::tie((int)(lhs.x * p), (int)(lhs.y * p), (int)(lhs.z * p));
        auto rhsTied = std::tie((int)(rhs.x * p), (int)(rhs.y * p), (int)(rhs.z * p));

        return lhsTied < rhsTied;
    }
};