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?

Your comparator should not be modifying the two objects. What exactly does `lhs.absolute();` do? Also, take a look at this: https://stackoverflow.com/questions/6218812/implementing-comparison-operators-via-tuple-and-tie-a-good-idea — alter_igel, Nov 22 '18 at 20:02
@alterigel absolute() makes x,y and z positive. Also he is not really changing them, I just messed up my minimal example. I will edit it. — jaba, Nov 22 '18 at 20:07
maybe I'm missing something obvious... but your comparator only test if lhs is inferior to rhs (=sorting the key), not if they are equal. Where is the code to compute your map key? — tetorea, Nov 22 '18 at 20:25
@tetorea The map key is the lhs vector. And if they are equal enough they should share the same key. I insert into the map with myMap[vector.absolute()] += faceArea; — jaba, Nov 22 '18 at 20:28
@alterigel I like the solution from the link but then I am not able to make values identical in the map that are really close because tie will compare the double values without precision. This will cause what I have now that 1e-16 and 0 are different. — jaba, Nov 22 '18 at 20:30
as I see it, your comparator is the third argument of your map initialization, right? If so, it is used to sort your keys, not testing if they are equal enough to ignore the newly inserted key. — tetorea, Nov 22 '18 at 20:39
just checked the documentation and realized the comparator is also used to check if keys are equal... sorry! — tetorea, Nov 22 '18 at 20:48
your comparator should return false at the beginning if the x,y and z values are equal enough. then you can test for the ordering. — tetorea, Nov 22 '18 at 20:55

Julien Villemure-Fréchette · Accepted Answer · 2018-11-22T21:01:31.640

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.

jaba · Answer 2 · 2018-11-26T19:01:06.103

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;
    }
};

watch out for operator precedence, you probably meant to write `(int)(lhs.x * prec)` instead of `(int)lhs.x * prec` which is equivalent to `floor(lhs.x) * prec` — Julien Villemure-Fréchette, Nov 22 '18 at 21:29
@JulienVillemure-Fréchette Thats correct. I will edit my answer. — jaba, Nov 22 '18 at 21:40
Another note: While testing int was to small for my values. Either use long, unsigned long long or whatever or even better normalize the vectors if you only care for the direction like me. — jaba, Nov 22 '18 at 21:49
A detail: using unsigned something in such situations is not necessarily a good idea as overflows become much more difficult to detect — Damien, Nov 23 '18 at 06:16

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

2 Answers2