Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Use phase class2 #205

Merged
merged 36 commits into from
May 9, 2024
Merged

Use phase class2 #205

merged 36 commits into from
May 9, 2024

Conversation

CSSFrancis
Copy link
Member

@CSSFrancis CSSFrancis commented Apr 12, 2024

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

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 CSSFrancis mentioned this pull request Apr 17, 2024
7 tasks
@CSSFrancis
Copy link
Member Author

@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.

Copy link

@viljarjf viljarjf left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good! See comments for some minor things. I did not go through all the tests very thoroughly, but they all pass at least.

diffsims/generators/simulation_generator.py Show resolved Hide resolved
diffsims/generators/simulation_generator.py Outdated Show resolved Hide resolved
diffsims/simulations/__init__.py Outdated Show resolved Hide resolved
Comment on lines 47 to 48
phases : Sequence[Phase]
The phases in the library.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
phases : Sequence[Phase]
The phases in the library.
simulation: Simulation2D
The simulation to get from.

Comment on lines 80 to 81
phases : Sequence[Phase]
The phases in the library.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
phases : Sequence[Phase]
The phases in the library.
simulation: Simulation2D
The simulation to get from.

diffsims/simulations/simulation2d.py Outdated Show resolved Hide resolved
return copy.deepcopy(self.coordinates)

def get_current_rotation(self):
"""Returns the matrix for the current matrix"""

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
"""Returns the matrix for the current matrix"""
"""Returns the matrix for the current rotation"""

diffsims/tests/generators/test_simulation_generator.py Outdated Show resolved Hide resolved
Comment on lines +260 to +269
for i, (excitation_error_i, r_spot_i) in enumerate(zip(excitation_error, r_spot)):

def integrand(theta):
# Equation 8 in L.Palatinus et al. Acta Cryst. (2019) B75, 512-522
S_zero = excitation_error_i
variable_term = r_spot_i * (phi) * np.cos(theta)
return shape_function(S_zero + variable_term, max_excitation, **kwargs)

# average factor integrated over the full revolution of the beam
shf[i] = (1 / (2 * np.pi)) * quad(integrand, 0, 2 * np.pi)[0]

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This can be vectorized with quad_vec, but that can wait

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's make an issue for this and we can handle this later.

setup.py Outdated Show resolved Hide resolved
@CSSFrancis
Copy link
Member Author

CSSFrancis commented Apr 21, 2024

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

Copy link

@viljarjf viljarjf left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this can be merged :)

@hakonanes
Copy link
Member

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
Copy link
Member

I can review within the next few days.

@hakonanes
Copy link
Member

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
Copy link
Member Author

CSSFrancis commented Apr 21, 2024

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.

  1. 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
Copy link
Member

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
Copy link
Member Author

CSSFrancis commented Apr 21, 2024

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
Copy link
Member

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()):

# Correct for occupancy and the Debye-Waller factor
dw_factor = np.exp(-atom.Bisoequiv * s2)
f *= atom.occupancy * dw_factor

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
Copy link
Member Author

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()):

# Correct for occupancy and the Debye-Waller factor
dw_factor = np.exp(-atom.Bisoequiv * s2)
f *= atom.occupancy * dw_factor

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
Copy link
Member Author

@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.

theta_templates = np.zeros((len(flattened_vectors), max_num_spots))
intensities_templates = np.zeros((len(flattened_vectors), max_num_spots))
for i, v in enumerate(flattened_vectors):
r, t, _ = v.to_polar()
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am getting some strange results when I try to template match, and I believe this line is the cause. Should it not be t, _, r = v.to_polar()? Seeing as Vector3d.to_polar returns (azimuthal, polar, radius) spherical coordinates. By inspection, the current way yields t-values close to pi/2, which would make sense for the polar angle with a mostly flat Ewald's sphere.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@viljarjf Sorry, I've gotten kind of confused with all of these comments and changes. We don't want to use that function at all... We want the radius of the projection of the vector onto the x,y plane.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@viljarjf This (should) be better

Copy link
Member

@hakonanes hakonanes left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've had a look at the new changes.

Comment on lines 30 to 31
This extends the :class:`ReciprocalLatticeVector` class. Diffracting Vectors
focus on the subset of reciporical lattice vectors that are relevant for
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
This extends the :class:`ReciprocalLatticeVector` class. Diffracting Vectors
focus on the subset of reciporical lattice vectors that are relevant for
This extends the :class:`ReciprocalLatticeVector` class. `DiffractingVector`
focus on the subset of reciprocal lattice vectors that are relevant for

Comment on lines 35 to 36
This class is only used internally to store the DiffractionVectors generated from the
DiffractionSimulation class. It is not (currently) intended to be used directly by the user.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
This class is only used internally to store the DiffractionVectors generated from the
DiffractionSimulation class. It is not (currently) intended to be used directly by the user.
This class is only used internally to store the diffracting vectors generated from a
:class:`~diffsims.simulations.DiffractionSimulation`. It is not intended to be used
directly by the user.

self._intensity = np.array(value)

@property
def lattice_aligned_data(self):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So these vectors are the same as they were before rotated by a rotation. Can't we separately keep track of a set of starting vectors, a rotation, and a set of rotated vectors? Why do these separate things have to be stored in a single class instance?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm.... I think this is related to above. What use is the RLV if it doesn't have information about the Rotation otherwise it is just the same as a Vector3D with confusing hkl values that don't really make physical sense and look correct because they are rounded.


def to_flat_polar(self):
"""Return the vectors in polar coordinates as projected onto the x,y plane"""
r = np.linalg.norm(self.data[:, :2], axis=1)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What if self.data.ndim > 2? The reciprocal lattice vectors, as normal 3D vectors, are allowed to have a "navigation shape" > 1. I suggest to flatten before extracting the x, y coordinates.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good point.

@@ -18,6 +18,7 @@

from collections import defaultdict
from copy import deepcopy
from typing import Tuple
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This line can be dropped.

from orix.quaternion import Rotation
from orix.crystal_map import Phase

from diffsims.crystallography import ReciprocalLatticeVector, DiffractingVector
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can the diffracting vector class be private?

You state in its docstring that it is not meant to be used by end-users. This tells me that the logic should be a simpler class or sets of functions to handle this particular use-case, instead of a public class with many caveats.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If the user wants the original reciprocal lattice vectors, rotations, or rotated vectors, can they instead get ReciprocalLatticeVector or Rotation (preferably the more powerful Orientation)?

@CSSFrancis
Copy link
Member Author

CSSFrancis commented May 7, 2024

@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 pyxem/orix#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 pyxem/orix#499.

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

setup.cfg Show resolved Hide resolved
@hakonanes
Copy link
Member

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
Copy link
Member Author

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

matrix = value
else:
raise ValueError("Rotation must be a Rotation or a 3x3 numpy array")
self.phase.structure.lattice.setLatPar(baserot=matrix)
Copy link

@viljarjf viljarjf May 8, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Phases are passed by refrence, so this would affect all other vectors with the same phase.
This is why I copied the phase in rotate_with_basis.

Suggested change
self.phase.structure.lattice.setLatPar(baserot=matrix)
self.phase = self.phase.deepcopy()
self.phase.structure.lattice.setLatPar(baserot=matrix)

The base rotation is not always the identity as well, especially for hexagonal crystals from cif or initialized from a structure. This comes from the enforced x || a, z || c* convention used in orix (https://orix.readthedocs.io/en/stable/tutorials/crystal_reference_frame.html).
What do you intend the basis_rotation member to be used for?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fair point. I'll probably just remove this then. I added it for completeness sake but didn't realize that the base rot was used for hexagonal phases.

@CSSFrancis
Copy link
Member Author

Based on:

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?

Copy link
Member

@hakonanes hakonanes left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

@CSSFrancis
Copy link
Member Author

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 CSSFrancis mentioned this pull request May 9, 2024
16 tasks
@CSSFrancis
Copy link
Member Author

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

@hakonanes hakonanes merged commit 507dda2 into pyxem:main May 9, 2024
11 checks passed
@hakonanes
Copy link
Member

hakonanes commented May 9, 2024

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

@CSSFrancis
Copy link
Member Author

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!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants