Rotation of ReciprocalLatticeVector throws an error

[hakonanes]

Describe the bug
Rotation of a ReciprocalLatticeVector, g, with a Rotation from orix, R, throws an error. The multiplication uses Quaternion.__mul__(). There, the other being multiplied, here g, fails the test isinstance(other, Miller) and thus phase information is not passed on when creating the rotated ReciprocalLatticeVector in its __init__().

Code in orix: https://github.com/pyxem/orix/blob/a1be2697dc0cf02974b00d06d94272de77efeafa/orix/quaternion/quaternion.py#L238-L243

Minimal example

from orix.crystal_map import Phase
from orix.quaternion import Rotation
from diffsims.crystallography import ReciprocalLatticeVector


phase = Phase(point_group="m-3m")
g = ReciprocalLatticeVector(phase, [1, 1, 1])

R = Rotation.random()
g1 = R * g

throws

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
File /Users/hakon/kode/scratches/ds_scratch.py:5
      2 R = Rotation.random()
      3 om = R.to_matrix()
----> 5 g1 = R * g
      6 #g2 = g.rotate_from_matrix(om)

File /opt/homebrew/Caskroom/miniforge/base/envs/ds/lib/python3.11/site-packages/orix/quaternion/rotation.py:123, in Rotation.__mul__(self, other)
    121     return Q
    122 if isinstance(other, Vector3d):
--> 123     v = Quaternion(self) * other
    124     improper = (self.improper * np.ones(other.shape)).astype(bool)
    125     v[improper] = -v[improper]

File /opt/homebrew/Caskroom/miniforge/base/envs/ds/lib/python3.11/site-packages/orix/quaternion/quaternion.py:243, in Quaternion.__mul__(self, other)
    241         return m
    242     else:
--> 243         return other.__class__(v)
    244 return NotImplemented

File ~/kode/personal/diffsims/diffsims/crystallography/reciprocal_lattice_vector.py:105, in ReciprocalLatticeVector.__init__(self, phase, xyz, hkl, hkil)
    103 def __init__(self, phase, xyz=None, hkl=None, hkil=None):
    104     self.phase = phase
--> 105     self._raise_if_no_point_group()
    107     if np.sum([i is not None for i in [xyz, hkl, hkil]]) != 1:
    108         raise ValueError(\"Exactly one of `xyz`, `hkl`, or `hkil` must be passed\")

File ~/kode/personal/diffsims/diffsims/crystallography/reciprocal_lattice_vector.py:1265, in ReciprocalLatticeVector._raise_if_no_point_group(self)
   1259 def _raise_if_no_point_group(self):
   1260     \"\"\"Raise ValueError if the phase attribute has no point group
   1261     set.
   1262 
   1263     \"\"\"
-> 1265     if self.phase.point_group is None:
   1266         raise ValueError(f\"The phase {self.phase} must have a point group set\")

AttributeError: 'numpy.ndarray' object has no attribute 'point_group'"

The rotated ReciprocalLatticeVector with the phase information intact should be returned.

I'm not 100% sure how to fix this. ReciprocalLatticeVector should not subclass Miller, as the coordinates of the former can only be (hkl). Any input would be appreciated.

[viljarjf]

Phase information kan be kept by changing https://github.com/pyxem/orix/blob/develop/orix/quaternion/quaternion.py#L238 to if hasattr(other, "phase"):, which feels a little hacky.
If we override __neg__ in ReciprocalLatticeVector too, this issue is solved.

[CSSFrancis]

Should this even be possible?

1. The hkl value is not preserved
2. How do you handle multiple orations. Should that return multiple RLV (that is actually not ideal for the simulation purposes as the size of that array becomes large!)

[viljarjf]

I might be misunderstanding what rotating a hkl vector entails. What do you mean by preserving the hkl value? When rotating, should they not change?
If the shapes are compatible, a Rotation object with multiple rotations will distribute as expected, e.g. a RLV with shape (10,) and a Rotation with shape (4, 1) would create a RLV with shape (4, 10). This is how Miller works currently.

[CSSFrancis]

@viljarjf The HKL is dependent on the RLV and the lattice. If we have a set of RLV and we rotate them (and don't rotate the lattice in some way) calculating the hkl no longer works as expected and gives non physical values.

[viljarjf]

This is exactly the same operation as your rotate_from_matrix, apart from the default of Rotation being lab2crystal when rotate_from_matrix implicitly uses crystal2lab. Is that not what you would want?

[CSSFrancis]

import diffpy
import numpy as np
from diffsims.crystallography import ReciprocalLatticeVector
from orix.crystal_map import Phase
from orix.quaternion import Rotation

# Shared parameters
latt = diffpy.structure.lattice.Lattice(2.464, 2.464, 6.711, 90, 90, 120)
atoms = [diffpy.structure.atom.Atom(atype='C', xyz=[0.0, 0.0, 0.25], lattice=latt),
         diffpy.structure.atom.Atom(atype='C', xyz=[0.0, 0.0, 0.75], lattice=latt),
         diffpy.structure.atom.Atom(atype='C', xyz=[1/3, 2/3, 0.25], lattice=latt),
         diffpy.structure.atom.Atom(atype='C', xyz=[2/3, 1/3, 0.75], lattice=latt),]
structure_matrix = diffpy.structure.Structure(atoms=atoms, lattice=latt)
calibration = 0.0262
p = Phase("Graphite", point_group="6/mmm", structure=structure_matrix)
euler_angles_new = np.array([[0, 90, 90]])
rot = Rotation.from_euler(euler_angles_new, degrees=True)
rlv = ReciprocalLatticeVector.from_min_dspacing(phase=p, min_dspacing= 1)


rotated_rlv = rlv.rotate_from_matrix(rot.to_matrix()[0])
rlv.hkl ==rotated_rlv.hkl # False

Maybe this will help show the problem related to rotating the RLV.

[CSSFrancis]

@hakonanes have you ever used the https://www.diffpy.org/diffpy.structure/mod_lattice.html baserot parameter. I think that we could have the rotation also rotate the lattice. That would at least allow us to rotate the RLV and be consistent with the hkl values.

[CSSFrancis]

@hakonanes This is partially fixed in pyxem/orix#499. Proper HKL values requires some more thought and better handling.

[hakonanes]

If we override __neg__ in ReciprocalLatticeVector too, this issue is solved.

@viljarjf, I assume you mean that we can override __mul__()? That is a valid solution, although at the cost of code duplication.

Should this even be possible?

@CSSFrancis, good point! Now that I think about it, I'd say no. ReciprocalLatticeVector was initially made to generate the set of strongly scattering reflections for a given unit cell. The intention was never to rotate them. I did this first only in #205 (comment), but was surprised at the result... I thought we subclassed Vector3d from orix perfectly, which we clearly do not.

---

I see you've suggested a change. May I suggest that we instead overwrite __mul__() as @viljarjf suggest, but throw an error explaining that this is not implemented? If we need to rotate the vectors, we can cast them to regular Miller (hkl) first.

[hakonanes]

have you ever used the https://www.diffpy.org/diffpy.structure/mod_lattice.html baserot parameter.

Nope.

I think that we could have the rotation also rotate the lattice. That would at least allow us to rotate the RLV and be consistent with the hkl values.

Hm... But, then what is the point of rotating the vectors in the first place? When I rotate vectors, I always want new vectors.

[CSSFrancis]

@hakonanes So what I've come to is that it should be possible to rotate the RLV as long as we can track the rotation.

This allows us to have consistent hkl values but change the underlying xyz values.

[hakonanes]

Can this be a feature of the private DiffractingVector class in #205 instead?

[CSSFrancis]

@hakonanes sure it can be... What is the argument for not having it part of the RLV. How does kikchuipy use the RLV for doing simulations? Do you need to do a rotation at some point or how are you handling the geometry considerations?

[hakonanes]

I'm fine with leaving my reported bug a bit longer.

What is the argument for not having it part of the RLV.

I've got to think a bit on it. Quoting python -c "import this":

Now is better than never.
Although never is often better than right now.

How does kikchuipy use the RLV for doing simulations?

We use composition instead of inheritance.

Steps for creating geometrical EBSD simulations:

1. Create a Kikuchi pattern simulator from a set of reciprocal lattice vectors (re-use the phase)
2. Supply a detector with projection parameters and one or more crystal rotations to simulate patterns for
3. Return simulations storing the detector, the rotations, the (un-rotated) reciprocal lattice vectors (that were visible on the detector for any rotation), and the gnomonic coordinates of Kikuchi lines and zone axes per rotation (the actual simulation)

Relevant part of tutorial: https://kikuchipy.org/en/stable/tutorials/geometrical_ebsd_simulations.html#Create-simulations-on-detector.

[CSSFrancis]

I've got to think a bit on it. Quoting python -c "import this":

Now is better than never.
Although never is often better than right now.

That is fine; that is kind of how we coded ourselves into this mess in pyxem :) At some point, we do have to figure out a way to fix it. From a purely pragmatic standpoint, I have to teach people how to do this next week :) and I don't think any of these changes break the API (only fix broken things) so it would be nice to push a fix before then.

I can just move everything to the DiffractingVector class and make it private but pyxem/orix#499 would be a nice addition.

How does kikchuipy use the RLV for doing simulations?

We use composition instead of inheritance.

Steps for creating geometrical EBSD simulations:

1. Create a Kikuchi pattern simulator from a set of reciprocal lattice vectors (re-use the phase)
2. Supply a detector with projection parameters and one or more crystal rotations to simulate patterns for
3. Return simulations storing the detector, the rotations, the (un-rotated) reciprocal lattice vectors (that were visible on the detector for any rotation), and the gnomonic coordinates of Kikuchi lines and zone axes per rotation (the actual simulation)

Relevant part of tutorial: https://kikuchipy.org/en/stable/tutorials/geometrical_ebsd_simulations.html#Create-simulations-on-detector.

Hmmm. Okay, that works as well...

[hakonanes]

But do you need to rotate reciprocal lattice vectors for the tutorial? This "bug" was reported based on our discussions in #205, not by a user in a real use-case.

[CSSFrancis]

But do you need to rotate reciprocal lattice vectors for the tutorial? This "bug" was reported based on our discussions in #205, not by a user in a real use-case.

It would be nice to be able to rotate them. The only reason I didn't rotate them in the simulation at first was to reduce the changes to the RLV.

Essentially what the rotation does is rotates the underlying lattice object (as well as the vectors). I guess the question is: is this only useful for electron diffraction or is this universally useful. In the case of Kikichi diffraction wouldn't this simplify your simulation calculations?

I'm struggling with the concept of kinematic diffraction differences obviously :)

[CSSFrancis]

@hakonanes I've looked at the kikuchipy simulation code and do you end up rotating the hkl vectors rather the xyz vectors?

[viljarjf]

If we override __neg__ in ReciprocalLatticeVector too, this issue is solved.

@viljarjf, I assume you mean that we can override __mul__()? That is a valid solution, although at the cost of code duplication.

I did mean __neg__, as it is Rotation.__mul__ that is called when multiplying with a Rotation on the left. Overriding __mul__ in RLV has no effect, as the error happens before Rotation.__mul__ returns NotImplemented (which would then try to call RLV's __mul__).
Quaternion.__mul__ handles the phase (if RLV inherits from Miller, or if we expand the subclass check to instead look for a phase-member), but the override in Rotation.__mul__ uses other.__neg__, which is inherited from Vector3d.
This does not throw an error for Miller, as it allows no phase to be supplied to the constructor. RLV requires the phase, and as the first argument at that. This is incompatible with a couple methods where Vector3d or its subclasses expects the first argument of the constructor to be xyz. Regardless, this is just implementation details, and this approach seems inelegant in my opinion anyway.

Actually, it seems Miller.__neg__ silently discards the phase, since it is inherited from Vector3d:

from orix.vector import Miller
from orix.crystal_map import Phase

p = Phase(point_group="m-3m")
v = Miller([1, 0, 0], phase=p)

print(v)
print(-v)

>>> Miller (1,), point group m-3m, xyz
[[1 0 0]]
>>> Miller (1,), point group None, xyz
[[-1  0  0]]

[CSSFrancis]

@hakonanes @viljarjf both of these are solved by the PR to orix.

[viljarjf]

@hakonanes have you ever used the https://www.diffpy.org/diffpy.structure/mod_lattice.html baserot parameter. I think that we could have the rotation also rotate the lattice. That would at least allow us to rotate the RLV and be consistent with the hkl values.

This is actually really promising, since what we really want is to rotate the vectors AND the basis. This is different from the normal handling of rotations, which only concerns the former.
How about a new method of RLV for this usecase? Something like

class ReciprocalLatticeVector(Vector3d):
    ...
    def rotate_with_basis(self, rot: Rotation):
        """Rotate both vectors and the basis with a given `Rotation`.
        This differs from simply multiplying with a `Rotation`, 
        as that would NOT update the basis.

        Parameters
        ----------
        rot : orix.quaternion.Rotation
            A rotation to apply to vectors and the basis.
        """
        # rotate basis
        new_phase = self.phase.deepcopy()
        br = new_phase.structure.lattice.baserot
        # In case the base rotation is set already
        new_br = br @ rot.to_matrix().squeeze()
        new_phase.structure.lattice.setLatPar(baserot=new_br)

        # rotate vectors
        vecs = ~rot * self.to_miller()

        self.data = vecs.data
        self.phase = new_phase

This avoids having to store the rotation and do a bunch of conversions, while keeping hkl values intact AND rotating the lattice as desired.

from diffsims.crystallography import DiffractingVector, ReciprocalLatticeVector
from orix.crystal_map import Phase
from orix.quaternion import Rotation
import numpy as np

p = Phase(point_group="6/mmm")
r = Rotation.random()
rlv = ReciprocalLatticeVector(p, hkl=[1, 2, 3])

# Somewhat convoluted initialization process, but it is private after all...
miller = ~r * rlv.to_miller()
dv = DiffractingVector(p, miller.data, rotation=r)

old_data = rlv.data.copy()
rlv.rotate_with_basis(r)

assert np.allclose(dv.hkl, rlv.hkl)
assert np.allclose(dv.data, rlv.data)
assert not np.allclose(rlv.data, old_data)

print(rlv)
>>> ReciprocalLatticeVector (1,),  (6/mmm)
[[1. 2. 3.]]

[CSSFrancis]

@viljarjf this would only work for Rotations with size==(1,) but that isn't really a huge issue as we already are creating seperate RLV for storing the simulation results.

[viljarjf]

True, that was intentional but not communicated

[CSSFrancis]

@hakonanes Do you feel more comfortable with this? This doesn't change anything about the RLV; it only uses already existing attributes in the Lattice object.

[hakonanes]

I've looked at the kikuchipy simulation code and do you end up rotating the hkl vectors rather the xyz vectors?

Yes, we rotate the (hkl) from the reciprocal lattice directly to the EBSD detector using a transformation matrix combining the three transformations detector -> sample, sample -> cartesian crystal, cartesian crystal -> reciprocal crystal.

https://github.com/pyxem/kikuchipy/blob/9b76c653918e1b1ea309a763b982f62d9df17ed5/kikuchipy/simulations/kikuchi_pattern_simulator.py#L290-L310

I did mean neg, as it is Rotation.mul that is called when multiplying with a Rotation on the left.

Of course, thanks for the thorough explanation. And for spotting the bug with the discarded phase in the negated crystal vector. I've opened pyxem/orix#501 to track that issue.

Do you feel more comfortable with this? This doesn't change anything about the RLV; it only uses already existing attributes in the Lattice object.

The symmetry operations in the Symmetry object attached to ReciprocalLatticeVector.phase.point_group assumes the e1 || a, e3 || c*. After rotating the basis, this assumption is invalid, no?

[viljarjf]

The symmetry operations in the Symmetry object attached to ReciprocalLatticeVector.phase.point_group assumes the e1 || a, e3 || c*. After rotating the basis, this assumption is invalid, no?

You're right...
I tried messing around by rotating the symmetry operations but I can't make it work properly

[CSSFrancis]

We could always go back to the idea of saving the Rotation seperately and then applying the inverse rotation when calculating the hkl. It maybe isn't as elegant of a solution but works.

[viljarjf]

I'm not 100% sure how to fix this. ReciprocalLatticeVector should not subclass Miller, as the coordinates of the former can only be (hkl). Any input would be appreciated.

This could be enforced in the constructor of RLV, no? i.e.

class ReciprocalLatticeVector(Miller):

    def __init__(self, xyz=None, hkl=None, hkil=None, phase=None):
        if phase is None:
            raise ValueError("Must provide phase")
        # Temporarily allow xyz
        self._allowed_coordinate_formats = ["hkl", "hkil", "xyz"]
        super().__init__(xyz=xyz, hkl=hkl, hkil=hkil, phase=phase)
        self._allowed_coordinate_formats.pop() # xyz is now illegal
        if self.coordinate_format == "xyz":
            self.coordinate_format = "hkl"

        self._raise_if_no_point_group()
        self._theta = np.full(self.shape, np.nan)
        self._structure_factor = np.full(self.shape, np.nan, dtype="complex128")
    
    @Miller.coordinate_format.setter
    def coordinate_format(self, value):
        """Set the vector coordinate format, either ``"hkl"``, or
        ``"hkil"``.

        """
        if value not in self._allowed_coordinate_formats:
            raise ValueError(f"Available coordinate formats are {self._allowed_coordinate_formats}")
        self._coordinate_format = value

With this (+ explicitly setting phase=phase everywhere, as it is no longer the first parameter but instead a keyword parameter), all tests pass, and ReciprocalLatticeVectors can be rotated. Additionally, a lot of the code duplication between the two classes can be avoided. The above example is a little hacky

Are there any more reasons for wanting to keep the two classes seperate, @hakonanes? I tried looking through previous discussions and did not find much, but I am not familiar with kikuchipy and its use of RLVs.