Implement rounding consistently to 12 decimal places in orix, eg. in unique() functions

[harripj]

Description of the change

A discussed in #281 there appears to be a bug in Miller.unique which subsequently causes a failing test in that PR. As As @hakonanes has confirmed in that PR, this appears to be a floating point bug, and rounding, even to a large precision (14 dp), corrects the bug.

So whilst implementing #281 has shown the bug and corrects it, I am not sure whether we should rely completely on this or whether we should manually implement some rounding in unique to be sure we always get the expected behaviour.

Some thoughts on this would be great.

This PR also corrects the test and removes an unneeded line as discussed in #281.

Update:
Rounding, for example in unique functions, has been consistently fixed to 12 dp for the repository.

Progress of the PR

• [n/a] Docstrings for all functions
• Unit tests with pytest for all lines
• Clean code style by running black via pre-commit

Example code

From @hakonanes in #281
#281 (comment)

For reviewers

• The PR title is short, concise, and will make sense 1 year later.
• [n/a] New functions are imported in corresponding __init__.py.
• New features, API changes, and deprecations are mentioned in the
  unreleased section in CHANGELOG.rst.

[hakonanes]

manually implement some rounding in unique to be sure we always get the expected behaviour

np.finfo(float).eps returns 2.220446049250313e-16 on my machine. Rounding to 14 decimals fixed the bug for the example in #281. I'm OK with rounding to 14 decimals. But Miller indices in orix aren't restricted to integers, and can be cartesian coordinates xyz. So we shouldn't round too much. np.finfo(np.float32).eps is 1.1920929e-07, this might be a bit too much.

[harripj]

From my testing the reason between the two versions ('stable' and 'paddy') in #281 (comment) is that at this line in Miller.unique

orix/orix/vector/miller.py

Line 675 in ff6e541

data = v2.data

the min and max differences between the two data arrays are -1.7763568394002505e-15 and 1.7763568394002505e-15, which then flows through the rest of the function.

[harripj]

In Rotation.unique we round to 6 dp., which I suspect is for the same reason:

orix/orix/quaternion/rotation.py

Lines 149 to 159 in ff6e541

	else:
	abcd = np.stack(
	[
	rotation.a.data,
	rotation.b.data,
	rotation.c.data,
	rotation.d.data,
	rotation.improper,
	],
	axis=-1,
	).round(6)

[harripj]

The behaviour is the same with the next test

orix/orix/tests/test_miller.py

Line 218 in ff6e541

assert m3.size == 345

The output size is 205 and stable rounding the data from 0 dp all the way down to 14 dp.

[harripj]

The indices test also needed correcting

orix/orix/tests/test_miller.py

Line 220 in ff6e541

assert np.allclose(idx[:10], [74, 448, 434, 409, 374, 330, 278, 447, 412, 432])

However the actual vectors returned before the fix and after are the same, ie. m2.unit[idx[:10]] did not change.

[harripj]

From @hakonanes #281:

np.finfo(float).eps returns 2.220446049250313e-16 on my machine. Rounding to 14 decimals fixed the bug for the example in #281. I'm OK with rounding to 14 decimals. But Miller indices in orix aren't restricted to integers, and can be cartesian coordinates xyz. So we shouldn't round too much. np.finfo(np.float32).eps is 1.1920929e-07, this might be a bit too much.

I agree, if we decide to round in unique functions (we already do in Rotation.unique) then I suggest we keep the same value. We could offer this as precision as a function argument, but I suggest we keep it the same for all functions. In Rotation.unique it is currently 6, but as @hakonanes suggested maybe something in the range 12-14 would be more appropriate.

[harripj]

The tests for Rotation.unique also pass at 12 and 14 dp precision.

[hakonanes]

This issue reminds me of mtex-toolbox/mtex#1261, where rounding to 12 decimals of vectors fixed an issue in TSL IPF color keys for some Laue classes.

We shouldn't round too much, but not too little so that we have to fix this bug again. I suggest 12 decimals for arrays of 64 bit floats, in all rounding cases as you say. If users don't set array data types themselves, they should always be 64 bit floats.

We could offer this as precision as a function argument

We could, but I don't see a use case for it?

[harripj]

I suggest 12 decimals for arrays of 64 bit floats, in all rounding cases as you say.

I agree.

There was a nasty rounding case in the tests. In sample_S2_cube_mesh:

orix/orix/sampling/S2_sampling.py

Lines 88 to 92 in ff6e541

	elif grid_type == grid_types[1]:
	steps = np.ceil(max_angle / res)
	k = np.arange(-steps, steps)
	theta = np.arctan(max_distance) / steps
	i = np.tan(k * theta)

With 5 degree resolution, depending on the rounding the returned value from ceil changed by a whole value, ie. for 9dp the value into ceil was 9+delta (9.000000003023132) -> 10, but with 12dp rounding it was 9-delta -> 9. I implemented a round to 2 dp before ceil to correct for this case, but another case where we should be careful with rounding. The test was also corrected for this.

It comes down to whether everyone is okay with ceil(9.000000003023132) being 9 rather than 10 due to the 2dp rounding.

Perhaps this will also have changed the notebooks? I will check.

[hakonanes]

Perhaps this will also have changed the notebooks? I will check.

They look good to me! Only changes I saw was which of the diamond lattice vectors on the equator were considered in the upper hemisphere.

I implemented a round to 2 dp before ceil to correct for this case, but another case where we should be careful with rounding.

Good catch. The S2 sampling is written by @din14970 in diffsims, which I took into orix (with his permission) to sample orientations for the IPF color keys. I don't consider it a nasty bug, though, as long as the resolution scales more or less as expected with increasing parameter. Haven't actually checked this thoroughly, which we perhaps should do and make into a brief user guide.

[hakonanes]

Sorry, should made this #281 (comment) here instead.

[harripj]

Actually, I find the assumption that the user wants the vectors within the fundamental sector not strong enough to be used in Miller.unique(). We should instead fix Miller.in_fundamental_sector() (inherited from Vector3d) so that it works, so that the unique ones within the sector can be found by projection followed by rejection.

Am I correct in thinking this can now be done as of #288?
And I think you are right about unique, if I am following your thought path then for a cubic material (-1, 0, 1) should be considered unique from (1, 0, 1) in unique, but if we projected to the fundamental sector they would both be considered the same?

If you want the direct or reciprocal lattice vectors within the fundamental sector in one call, we should make a convenience method Miller.from_fundamental_sector() calling these different methods. Will that meet your needs @harripj?

It's not necessarily on my cards at the moment, so happy to hold off. I think this can be done explicitly using a combination of in_fundamental_sector() followed by unique (as per your suggestion)? Maybe best to let the user do that explicitly if that is what they want to do?

[hakonanes]

Am I correct in thinking this can now be done as of #288

Yes, with that PR and the rounding in Miller.unique() from this one

>>> from orix.crystal_map import Phase
>>> from orix.vector import Miller
>>> diamond = Phase(space_group=227)
>>> h = Miller.from_highest_indices(phase=diamond, uvw=[10, 10, 10])
>>> h2 = h.in_fundamental_sector()
>>> h3 = h2.unique(use_symmetry=False)
>>> h4 = h.unique(use_symmetry=True)
>>> print(h.size, h2.size, h3.size, h4.size)
9260 9260 285 285
>>> fig = h3.scatter(return_figure=True, s=5)
>>> h4.scatter(figure=fig, c="g", s=5)
>>> fig.axes[0].plot(diamond.point_group.fundamental_sector.edges, color="r", zorder=0)

[hakonanes]

Good catch, important fix! Two minor suggestions.

[pc494]

I'm happy with this, for normal rounding 12/14 decimals place seems fine, ceil bits can be annoying. In general, given how annoying machine specific bugs are I would tend to round too much (and get an error because we rounded - explicit) than too little (and get an error because of machine precision - implicit).

Approve, will leave @hakonanes to merge.

[pc494]

Turns out this was needed for the next thing I had to look at so I'm merging myself.