Use phase class2

[CSSFrancis]

Description of the change

This simplified the commit history for #201 and focuses the scope of the PR. Hopefully that makes it (slightly) more manageable when checking/merging.

Progress of the PR

• [ ] Docstrings for all functions
• [x] Unit tests with pytest for all lines
• [x] Clean code style by running black

Minimal example of the bug fix or new feature

>>> from diffsims.utils import sim_utils
>>> # Your new feature...

For reviewers

• [ ] The PR title is short, concise, and will make sense 1 year later.
• [ ] New functions are imported in corresponding __init__.py.
• [ ] New features, API changes, and deprecations are mentioned in the unreleased section in CHANGELOG.rst.
• [ ] Contributor(s) are listed correctly in credits in diffsims/release_info.py and in .zenodo.json.

[CSSFrancis]

@viljarjf, can you review this? It seems like we should get this in, and then we can work on getting the rest of diffsims ready for a release.

[CSSFrancis]

@viljarjf This should be better now based on your comments. Thanks, let me know what you think and maybe we can merge this?

[hakonanes]

I'd like to review it given it touches on ReciprocalLatticeVector, which is heavily used in kikuchipy workflows. I'm not so sure about the intensity parameter.

[hakonanes]

I can review within the next few days.

[hakonanes]

So, sorry about the hold-up, @CSSFrancis.

Couple of questions right away:

1. Could you explain the reason for the intensity attribute in ReciprocalLatticeVector? If you're interested in the kinematical scattering intensity I = |F|^2, this can be calculated from F = calculate_structure_factor(). The new attribute lacks a docstring. To me, just looking at the new additions and tests, the attribute looks like a placeholder for any array you want, attached to the ReciprocalLatticeVector class.
2. Why cannot the new simulation classes be part of the current simulation module, diffsims.sims? Having two modules for diffraction simulations is confusing.

[CSSFrancis]

So, sorry about the hold-up, @CSSFrancis.

Couple of questions right away:

1. Could you explain the reason for the intensity attribute in ReciprocalLatticeVector? If you're interested in the kinematical scattering intensity I = |F|^2, this can be calculated from F = calculate_structure_factor(). The new attribute lacks a docstring. To me, just looking at the new additions and tests, the attribute looks like a placeholder for any array you want, attached to the ReciprocalLatticeVector class.

That's a good question; it's mostly that we really want to calculate things once and then not do it again, or it's going to really slow down the orientation mapping. We can add a function calculate_scattering_intensity to the DiffractionVectors class? I can refactor that. That data also contains information about the intersection of the vectors and the Ewald sphere. Maybe this should be a separate object but it needs to be stored somewhere.

2. Why cannot the new simulation classes be part of the current simulation module, diffsims.sims? Having two modules for diffraction simulations is confusing.

Shorting names is fairly un-pythonic. A 1-1 mapping of Signal1D -> Simulation1D and Signal2D --> Simulation2D makes more sense in the long term. I want to make this change long-term in diffsims and don't want to break the pyxem API when that happens. I'm open to suggestions.

[hakonanes]

we really want to calculate things once and then not do it again

I totally understand that. That's why e.g. the structure_factor attribute carries over when slicing ReciprocalLatticeVector already.

We can add a function calculate_scattering_intensity to the DiffractionVectors class?

Hm, where is this class? We have the calculate_structure_factor() method to calculate kinematical scattering intensities for ReciprocalLatticeVector. This has served us well in kikuchipy for calculating kinematical EBSD patterns. We should ideally have one way of doing things. Having attributes for the structure factor F, calculated from calculate_structure_factor(), and then for an intensity I = |F|^2, which may not be derived from F, is confusing (to me).

That data also contains information about the intersection of the vectors and the Ewald sphere. Maybe this should be a separate object but it needs to be stored somewhere.

But, the end result is ultimately structure factors? Or the intensities directly? And what is the current route for calculating these structure factors (steps and input parameters)? It may be best to include that route into the ReciprocalLatticeVector class.

I want to make this change long-term in diffsims and don't want to break the pyxem API when that happens. I'm open to suggestions.

This is fine, and something I agree with in general. Can I suggest to present this idea in an issue for a discussion, albeit a brief one?

[CSSFrancis]

I totally understand that. That's why e.g. the structure_factor attribute carries over when slicing ReciprocalLatticeVector already.

We can add a function calculate_scattering_intensity to the DiffractionVectors class?

Hm, where is this class? We have the calculate_structure_factor() method to calculate kinematical scattering intensities for ReciprocalLatticeVector. This has served us well in kikuchipy for calculating kinematical EBSD patterns. We should ideally have one way of doing things. Having attributes for the structure factor F, calculated from calculate_structure_factor(), and then for an intensity I = |F|^2, which may not be derived from F, is confusing (to me).

I misspoke sorry :) Trying to deal with too many things at one time... I meant ReciprocalLatticeVector (DiffractionVectors is a similar class in pyxem).

Maybe it's best if we think about the entire workflow and determine which parts should be methods for ReciprocalLatticeVector

We want to:

1 - Copy the ReciprocalLatticeVector class 2 - Rotate the ReciprocalLatticeVector (to do this I added the rotate_from_matrix function 3 - Determine which ReciprocalLatticeVector vectors intersect the Ewald Sphere 4 - Calculate the diffraction intensities. This includes the Debye Waller factors and the shape function. 5 - Remove ReciprocalLatticeVector vectors below a cutoff intensity

For step 4, we can derive the intensities from the structure factor. We just need to make sure we are only calculating this once; otherwise, it might slow down the template matching. We must also add Debye Waller factors and the shape function. Having these factors in the ReciprocalLatticeVector class is a bit weird as they are more closely related to the SimulationGenerator.

That data also contains information about the intersection of the vectors and the Ewald sphere. Maybe this should be a separate object but it needs to be stored somewhere.

But, the end result is ultimately structure factors? Or the intensities directly? And what is the current route for calculating these structure factors (steps and input parameters)? It may be best to include that route into the ReciprocalLatticeVector class.

See above. The result that we care about if the diffraction intensity for comparison. I don't think we want those steps in the ReciprocalLatticeVector as it requires that the ReciprocalLatticeVector object has information about the Ewald sphere which seems like it should be part of the SimulationGenerator. I could be convinced otherwise.

I guess the question is: If you wanted to plot the diffraction spots and the Kikuchi lines, how would that look/ how is that handled? It might be nice to do something like: s.plot(diffraction=True, kikuchi_lines=True)

I want to make this change long-term in diffsims and don't want to break the pyxem API when that happens. I'm open to suggestions.

This is fine, and something I agree with in general. Can I suggest to present this idea in an issue for a discussion, albeit a brief one?

Sounds good :)

[hakonanes]

I don't think we want those steps in the ReciprocalLatticeVector as it requires that the ReciprocalLatticeVector object has information about the Ewald sphere which seems like it should be part of the SimulationGenerator.

I agree! It sounds like we want the generator to use ReciprocalLatticeVector (composition), rather than adding to the vector class. I'm of the opinion that no microscope- or geometry-specific parameters should be part of this class.

The only method that takes a microscope parameter, an accelerating voltage, is calculate_theta() for populating the theta attribute, twice the Bragg angle. This should arguably not be an attribute of the class, because information of the voltage is lost upon return. And having a ReciprocalLatticeVector.voltage attribute makes no sense. But, it's been a part of the API for many years, so it should be kept.

We must also add Debye Waller factors

These are already part of the diffpy.structure atoms API in the Atom.Bisoequiv attribute. Can you use this? The factors are used when calculating the atomic scattering factors, f (used in ReciprocalLatticeVector.calculate_structure_factor()):

https://github.com/pyxem/diffsims/blob/bb359b9a6326f22ab414665b005d3e04e5d56340/diffsims/structure_factor/atomic_scattering_factor.py#L62-L64

It might be nice to do something like: s.plot(diffraction=True, kikuchi_lines=True)

What is s here? What information about the microscope, phase, detector, and scattering vectors does it have?

We have plotting of Kikuchi lines in kikuchipy (see this tutorial). It may be that the simulation generator should be upstreamed to diffsims, but that would also require upstreaming EBSDDetector, which powers most of kikuchipy... Perhaps just taking necessary code from kikuchipy is fine for the Kikuchi lines. I don't know.

[CSSFrancis]

I agree! It sounds like we want the generator to use ReciprocalLatticeVector (composition), rather than adding to the vector class. I'm of the opinion that no microscope- or geometry-specific parameters should be part of this class.

The only method that takes a microscope parameter, an accelerating voltage, is calculate_theta() for populating the theta attribute, twice the Bragg angle. This should arguably not be an attribute of the class, because information of the voltage is lost upon return. And having a ReciprocalLatticeVector.voltage attribute makes no sense. But, it's been a part of the API for many years, so it should be kept.

@hakonanes Okay, so what do we do about the diffraction intensity, which depends on the microscope parameters? I can separate them if you like and make that an attribute of Simulation2D. Still, then there is always the possibility that you have more intensity values than ReciprocalLatticeVectors if someone slices either of them.

We must also add Debye Waller factors

These are already part of the diffpy.structure atoms API in the Atom.Bisoequiv attribute. Can you use this? The factors are used when calculating the atomic scattering factors, f (used in ReciprocalLatticeVector.calculate_structure_factor()):

https://github.com/pyxem/diffsims/blob/bb359b9a6326f22ab414665b005d3e04e5d56340/diffsims/structure_factor/atomic_scattering_factor.py#L62-L64

I'm not really trying to change everything all at once... Can we leave that to another PR? This PR is just supposed to use the orix phase class and ReciprocalLatticeVector so we should focus on making that one change and then clean up any other part that we feel isn't working. This PR is holding up a lot and we need to figure out how to get this in.

It might be nice to do something like: s.plot(diffraction=True, kikuchi_lines=True)

What is s here? What information about the microscope, phase, detector, and scattering vectors does it have?

That would be your Simulation2D object. It might make more sense to calculate the two separately and add them to the same plot. It might also be nice to have one SimulationGenerator object that can do both kikuchi and vector diffraction.

We have plotting of Kikuchi lines in kikuchipy (see this tutorial). It may be that the simulation generator should be upstreamed to diffsims, but that would also require upstreaming EBSDDetector, which powers most of kikuchipy... Perhaps just taking necessary code from kikuchipy is fine for the Kikuchi lines. I don't know.

For now maybe we can just have a tutorial that makes a combined plot is someone is interested.

[CSSFrancis]

@hakonanes I guess the question is what needs to be changed for this to be merged? Most of this is just copied from the old function. It's not really anything new, just makes it significantly faster and (somewhat) cleaner.

[hakonanes]

I can separate them if you like and make that an attribute of Simulation2D

That sounds like the best approach. I hope you can find a clean way of combining the intensity and vectors. As I've argued above, I don't think the intensity should ~not~ be part of the reciprocal lattice vector class. I might be convinced otherwise by arguments relating to the nature of the vectors and the intensities.

Still, then there is always the possibility that you have more intensity values than ReciprocalLatticeVectors if someone slices either of them.

Can you add a slicer to the simulator instead? You can then slice the intensity array and vectors accordingly.

Can we leave that to another PR?

Absolutely! I just pointed out that you may not need to have separate arrays for the Debye-Waller factors, and just use diffpy.structure instead.

This PR is just supposed to use the orix phase class and ReciprocalLatticeVector so we should focus on making that one change

I realize that, and I'm sorry for holding this up. But I'm not a fan of adding use-case-specific attributes and function parameters to a general class that is instrumental to workflows in downstream packages, like e.g. kikuchipy.

I guess the question is what needs to be changed for this to be merged?

I suggest to remove the new additions to ReciprocalLatticeVector and instead use the vectors alongside the intensities array in a simulator. I can review after that is done.

[CSSFrancis]

@hakonanes, how do you feel about a subclass of ReciporicalLatticeVector, something like DiffractingVector, which specifically handles ReciporicalLatticeVectors that intersect with the Ewald Sphere?

[hakonanes]

That doesn't sound too bad. If we subclass, and don't use the structure factor calculations already there, these should be overridden to return not implemented, though.

[CSSFrancis]

@hakonanes Can you look over this when you get the chance?

[hakonanes]

Yes, sorry, will do within the next few days!

[CSSFrancis]

Yes, sorry, will do within the next few days!

@hakonanes Any chance we can get this done this week? It would be good to have this done and release a new pyxem version before the Diamond workshop on the 13th.

[CSSFrancis]

@hakonanes Okay so I think we have to sort out some things.

1. Can the RLV be ragged? --> Probably not unless the Miller and Vector3D class supports it.
2. Should we allow the RLV to be rotated? I think the answer to this is yes. But this requires that we add a .rotation attribute to the RLV so that the .hkl value is computed like: https://github.com/CSSFrancis/diffsims/blob/9bf940441bff960c8b11642b26cf79d66f13c8a0/diffsims/simulations/simulation2d.py#L556C1-L560C18 --> Lets just keep this so we can't rotate it just to get this in.
3. Should the RLV have extra information like the intensity? --> No.
4. How do we store only the selected (ragged) diffracting RLV. --> This can just be an array of diffracting vectors.

[CSSFrancis]

Okay, for the sake of getting this in and because we really, really need to get this out before the Diamond workshop in one week, I'm just going to make the DiffractingVector class private. We are getting super stuck trying to make diffsims perfect rather than just making it better... At the end of the day it's a 0.6.0 version, it's not a major release, and we aren't breaking anything.

We can make this a beta feature or something like that if it makes you feel better.

[viljarjf]

If the existence of rotate_from_matrix-method is what is blocking this, and it is only used in SimulationGenerator.get_intersecting_reflections, then how about this instead:

    def get_intersecting_reflections(
        self,
        recip: ReciprocalLatticeVector,
        rot: Rotation,
        wavelength: float,
        max_excitation_error: float,
        shape_factor_width: float = None,
        with_direct_beam: bool = True,
    ):
        ...
        initial_hkl = recip.hkl
        # To ease calculations, we transform the reciprocal lattice vectors rather than the electron beam.
        # Therefore, we need to invert the rotation, as it normally transforms the beam (lab2crystal).
        rotated_vectors = (~rot * recip.to_miller()).data

        if with_direct_beam:
            rotated_vectors = np.vstack([rotated_vectors, [0, 0, 0]])
            initial_hkl = np.vstack([initial_hkl, [0, 0, 0]])

        # Identify the excitation errors of all points (distance from point to Ewald sphere)
        ...

The behaviour is the exact same as before. I can commit this change if you want @CSSFrancis?

[CSSFrancis]

@viljarjf sure you can make a PR to my branch and I can merge it if you want :)

I think that is the right way to do things as the hkl complicates things.

I still think the DiffractingVector class being private might be a good thing. It's only really used for simulation and we need to figure out the best way to handle the hkl. Based on the timeline minimal changes but still being able to do everything we want is probably the best way forward.

[CSSFrancis]

Also just wanted to say thank you @viljarjf and @hakonanes for your help and comments. Sorry if I've been pushing this a bit too much :)

[CSSFrancis]

Okay, so I made a couple of changes:

1. Fixed the failing test. I forgot that the dir is set to the top-level pytest function.
2. Removed the rotate_from_matrix method and used the Rotation class base on @viljarjf
3. Make the DiffractingVector class private. (we could probably make it public.)
4. Added a rotation attribute to the DiffractingVector class to handle hkl.

[CSSFrancis]

@hakonanes I think lots of the confusion is related to not properly handling the rotation with the RLV.

We need to track the rotations applied to the RLV to recover the hkl. This is possible with https://github.com/pyxem/orix/pull/499 and adding the rotation property to the RLV.

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

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

R1 = Rotation.random(2)
random1 = R1 * g
R2 = Rotation.random(2)
random2 =R2*random1

assert np.allclose(g.hkl, random1.hkl) # Previously this would fail
assert np.allclose(g.hkl, random2.hkl) # Previously this would fail

I can push this change if you want but it requires https://github.com/pyxem/orix/pull/499.

I think this helps the logic quite a bit but does change the ReciprocalLatticeVector class a little bit.

[hakonanes]

I'm just going to make the DiffractingVector class private

Perfect. Then I've got no objections to how that class is used internally in diffsims.

[CSSFrancis]

@hakonanes I made some changes based on @viljarjf's suggestion in #211

[CSSFrancis]

Based on:

• https://github.com/pyxem/orix/issues/501

• 211

I just made everything part of the private DiffractingVectors class. That way, we can merge things and start working on the next part, the Simulation Part of this code is much more straight forward and I'm much more confident on that implementation.

@hakonanes can you take another pass through and then we can think about merging?

[CSSFrancis]

Just a few suggestions left. After those are addressed, I'm OK with the parts of this PR that I've reviewed. Will have a final look after that.

Done!

[CSSFrancis]

@hakonanes Did you want to do another pass through or do you want to merge?

[hakonanes]

Fingers crossed the simulations work as you intended in the wild next week!

[CSSFrancis]

Fingers crossed the simulations work as you intended in the wild next week!

At the very least it's an improvement :) Thanks for your help here!