Object3d.unique ordering

[viljarjf]

Hello!

I have been playing around with diffsims' ReciprocalLatticeVector.unique, which uses Miller.unique, which again uses Object3d.unique.
Specifically, I needed to be able to re-construct the initial array from the unique ones, so I specified "return_inverse=True".
This led to some issues with the ordering of the elements, since the returned data from Object3d.unique sorts the indices of the unique elements, without the necessary changes to the array of indices to reconstruct the original.

In the process of addressing this, I also fixed the symmetry-aware Miller.unique, so initial order is preserved here as well. However, that broke some tests, in particular the final assertion in orix/tests/test_vector3d/test_miller.py::TestMiller::test_unique. This line checks the indices of the first 10 unique elements in a large Miller object, and all these indices have fairly large and unexpectedly disorganized values (65, 283, 278 ect.)
By changing the "idx" and "inv" to preserve the order in the input when using Object3d.unique and Miller.unique, the indices in the mentioned test becomes [0, 1, 2, ....] instead.

Another test that broke was orix/tests/test_quaternion/test_orientation_region.py::test_get_large_cell_normals. I did not investigate this one further, but from what I can see it is a lot of raw data that gets compared to computed data. The order of the pre-computed data might very well be intentional, and not a consequence of Object3d.unique.

I personally think it better that the order of the input is preserved, but as this seems to be an old implementation I worry it will break things other places too.

Edit: To illustrate:

from numpy import array_equal
from orix.vector import Vector3d

# Make some random vector with duplicates
rnd = Vector3d.random(30)
rnd.data = rnd.data.round(5)
v = Vector3d.stack([rnd.get_random_sample(10, replace=True) for _ in range(10)]).flatten()

# Get unique elements and reconstruction arrays
u, idx, inv = v.unique(return_index=True, return_inverse=True)

# Check that there are duplicates
assert v.size > u.size

# Check that all elements are present in both, ignoring order
def unordered_equal(vec1: Vector3d, vec2: Vector3d) -> bool:
    assert vec1.shape == vec2.shape
    # Use tuples so all three elements are checked
    vec1_tuples = [tuple(v.data.squeeze().tolist()) for v in vec1]
    vec2_tuples = [tuple(v.data.squeeze().tolist()) for v in vec2]
    # Check both ways
    return all(v1 in vec2_tuples for v1 in vec1_tuples) and all(v2 in vec1_tuples for v2 in vec2_tuples) 

# Passes
assert unordered_equal(v[idx], u)
assert unordered_equal(v, u[inv])
assert unordered_equal(v[idx][inv], v)

# idx and inv are self-consistent
assert array_equal(v[idx][inv].data, v.data)

# Fails 
assert array_equal(v[idx].data, u.data)
assert array_equal(v.data, u[inv].data)

[argerlt]

@viljarjf , can you check if perhaps this is related to #595? I suspect it is, but cannot tell for sure.

[viljarjf]

Hi @argerlt , I don't think it is related. The issue arises from

orix/orix/_base.py

Lines 253 to 255 in ca1db33

	_, idx, inv = np.unique(data, axis=0, return_index=True, return_inverse=True)
	obj = self.__class__(data[np.sort(idx), : self.dim])
	obj._data = data[np.sort(idx)]

, where a sort is performed on the indices before accessing the data, without retaining said sort. The inverse of that sort needs to be performed on the inverse indices, to retain the correct properties. The following passes all tests in my initial example, but breaks a few tests as mentioned:

        _, idx, inv = np.unique(data, axis=0, return_index=True, return_inverse=True, sorted=False)
        order = np.argsort(idx)
        inv_order = np.argsort(order)
        idx = idx[order]
        inv = inv_order[inv]
        obj = self.__class__(data[idx, : self.dim])
        obj._data = data[idx]

Using your fix for the flatten order did not pass the example tests, although I have not tried unique on non-flat arrays yet so there might be even more issues when including that...

[hakonanes]

The order of the pre-computed data might very well be intentional, and not a consequence of Object3d.unique.

@viljarjf, as long as the fundamental zones look the same, I cannot see how the order of zone normals (quaternions) should matter.

I personally think it better that the order of the input is preserved, but as this seems to be an old implementation I worry it will break things other places too.

I think a clean solution would be to add your implementation alongside the existing one, with a unique(, ...legacy: bool = True) flag. We can then start using your implementation internally, and with time make it the default. How does that sound?

[argerlt]

@viljarjf I played with this and I think your solution is the right one. If you submit a PR, I can help with the failing tests, but as far as i can tell, the tests are just compensating for the scrambling.

I also second @hakonanes suggestion of a legacy option, it gives internal and external code something to to fall back on if your update breaks something, which I believe is the main concern you are expressing here?

I personally think it better that the order of the input is preserved, but as this seems to be an old implementation I worry it will break things other places too.

It might be worth checking with @bm424 since that was the last person who touched the actual logic in this function back in 2018 (commit for reference ), but what you are saying makes a lot of sense, and I'm all for imitating numpy function behavior when we can.