From 1fc6a876a9372afa9035c0894dc7a393311dc657 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 12 Apr 2024 17:35:13 -0500 Subject: [PATCH 01/36] New Feature: Add Simulation classes using Phase Class --- diffsims/generators/simulation_generator.py | 418 +++++++++++++ diffsims/simulations/__init__.py | 20 + diffsims/simulations/simulation1d.py | 88 +++ diffsims/simulations/simulation2d.py | 614 ++++++++++++++++++++ 4 files changed, 1140 insertions(+) create mode 100644 diffsims/generators/simulation_generator.py create mode 100644 diffsims/simulations/__init__.py create mode 100644 diffsims/simulations/simulation1d.py create mode 100644 diffsims/simulations/simulation2d.py diff --git a/diffsims/generators/simulation_generator.py b/diffsims/generators/simulation_generator.py new file mode 100644 index 00000000..c88c9ee3 --- /dev/null +++ b/diffsims/generators/simulation_generator.py @@ -0,0 +1,418 @@ +# -*- coding: utf-8 -*- +# Copyright 2017-2024 The diffsims developers +# +# This file is part of diffsims. +# +# diffsims is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# diffsims is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with diffsims. If not, see . + +from typing import Union, Sequence +import numpy as np + +from orix.quaternion import Rotation +from orix.crystal_map import Phase + +from diffsims.crystallography import ReciprocalLatticeVector +from diffsims.utils.shape_factor_models import ( + linear, + atanc, + lorentzian, + sinc, + sin2c, + lorentzian_precession, + _shape_factor_precession, +) + +from diffsims.utils.sim_utils import ( + get_electron_wavelength, + get_kinematical_intensities, + is_lattice_hexagonal, + get_points_in_sphere, + get_intensities_params, +) + +_shape_factor_model_mapping = { + "linear": linear, + "atanc": atanc, + "sinc": sinc, + "sin2c": sin2c, + "lorentzian": lorentzian, +} + +from diffsims.simulations import Simulation1D, Simulation2D + + +class SimulationGenerator: + """ + A class for generating kinematic diffraction simulations. + """ + + def __repr__(self): + return ( + f"SimulationGenerator(accelerating_voltage={self.accelerating_voltage}, " + f"scattering_params={self.scattering_params}, " + f"approximate_precession={self.approximate_precession})" + ) + + def __init__( + self, + accelerating_voltage: float = 200, + scattering_params: str = "lobato", + precession_angle: float = 0, + shape_factor_model: str = "lorentzian", + approximate_precession: bool = True, + minimum_intensity: float = 1e-20, + **kwargs, + ): + """ + Parameters + ---------- + accelerating_voltage + The accelerating voltage of the electrons in keV. + scattering_params + The scattering parameters to use. One of 'lobato', 'xtables' + precession_angle + The precession angle in degrees. If 0, no precession is applied. + shape_factor_model + The shape factor model to use. One of 'linear', 'atanc', 'sinc', 'sin2c', 'lorentzian' + approximate_precession + If True, the precession is approximated by a Lorentzian function. + minimum_intensity + The minimum intensity of a reflection to be included in the profile. + kwargs + Keyword arguments to pass to the shape factor model. + + """ + self.accelerating_voltage = accelerating_voltage + self.precession_angle = np.abs(precession_angle) + self.approximate_precession = approximate_precession + if isinstance(shape_factor_model, str): + if shape_factor_model in _shape_factor_model_mapping.keys(): + self.shape_factor_model = _shape_factor_model_mapping[ + shape_factor_model + ] + else: + raise NotImplementedError( + f"{shape_factor_model} is not a recognized shape factor " + f"model, choose from: {_shape_factor_model_mapping.keys()} " + f"or provide your own function." + ) + else: + self.shape_factor_model = shape_factor_model + self.minimum_intensity = minimum_intensity + self.shape_factor_kwargs = kwargs + if scattering_params in ["lobato", "xtables", None]: + self.scattering_params = scattering_params + else: + raise NotImplementedError( + "The scattering parameters `{}` is not implemented. " + "See documentation for available " + "implementations.".format(scattering_params) + ) + + @property + def wavelength(self): + return get_electron_wavelength(self.accelerating_voltage) + + def calculate_ed_data( + self, + phase: Union[Phase, Sequence[Phase]], + rotation: Union[Rotation, Sequence[Rotation]] = Rotation.from_euler( + (0, 0, 0), degrees=True + ), + reciprocal_radius: float = 1.0, + with_direct_beam: bool = True, + max_excitation_error: float = 1e-2, + shape_factor_width: float = None, + debye_waller_factors: dict = None, + ): + """Calculates the diffraction pattern for one or more phases given a list + of rotations for each phase. + + Parameters + ---------- + phase: + The phase(s) for which to derive the diffraction pattern. + reciprocal_radius + The maximum radius of the sphere of reciprocal space to + sample, in reciprocal Angstroms. + rotation + The Rotation object(s) to apply to the structure and then + calculate the diffraction pattern. + with_direct_beam + If True, the direct beam is included in the simulated + diffraction pattern. If False, it is not. + max_excitation_error + The cut-off for geometric excitation error in the z-direction + in units of reciprocal Angstroms. Spots with a larger distance + from the Ewald sphere are removed from the pattern. + Related to the extinction distance and roughly equal to 1/thickness. + shape_factor_width + Determines the width of the reciprocal rel-rod, for fine-grained + control. If not set will be set equal to max_excitation_error. + debye_waller_factors + Maps element names to their temperature-dependent Debye-Waller factors. + + Returns + ------- + diffsims.sims.diffraction_simulation.DiffractionSimulation + The data associated with this structure and diffraction setup. + """ + if isinstance(phase, Phase): + phase = [phase] + if isinstance(rotation, Rotation): + rotation = [rotation] + if debye_waller_factors is None: + debye_waller_factors = {} + # Specify variables used in calculation + wavelength = self.wavelength + + # Rotate using all the rotations in the list + vectors = [] + for p, rotate in zip(phase, rotation): + recip = ReciprocalLatticeVector.from_min_dspacing( + p, + min_dspacing=1 / reciprocal_radius, + include_zero_beam=with_direct_beam, + ) + phase_vectors = [] + for rot in rotate.to_matrix(): + # Calculate the reciprocal lattice vectors that intersect the Ewald sphere. + ( + intersected_vectors, + hkl, + shape_factor, + ) = self.get_intersecting_reflections( + recip, + rot, + wavelength, + max_excitation_error, + shape_factor_width=shape_factor_width, + with_direct_beam=with_direct_beam, + ) + + # Calculate diffracted intensities based on a kinematic model. + intensities = get_kinematical_intensities( + p.structure, + hkl, + intersected_vectors.gspacing, + prefactor=shape_factor, + scattering_params=self.scattering_params, + debye_waller_factors=debye_waller_factors, + ) + + # Threshold peaks included in simulation as factor of zero beam intensity. + peak_mask = intensities > np.max(intensities) * self.minimum_intensity + intensities = intensities[peak_mask] + intersected_vectors = intersected_vectors[peak_mask] + intersected_vectors.intensity = intensities + phase_vectors.append(intersected_vectors) + vectors.append(phase_vectors) + + if len(phase) == 1: + vectors = vectors[0] + phase = phase[0] + rotation = rotation[0] + if rotation.size == 1: + vectors = vectors[0] + + # Create a simulation object + sim = Simulation2D( + phases=phase, + coordinates=vectors, + rotations=rotation, + simulation_generator=self, + reciprocal_radius=reciprocal_radius, + ) + return sim + + def calculate_diffraction1d( + self, + phase: Phase, + reciprocal_radius: float = 1.0, + minimum_intensity: float = 1e-3, + debye_waller_factors: dict = None, + ): + """Calculates the 1-D profile of the diffraction pattern for one phases. + + This is useful for plotting the diffracting reflections for some phases. + + Parameters + ---------- + phase: + The phase for which to derive the diffraction pattern. + reciprocal_radius + The maximum radius of the sphere of reciprocal space to + sample, in reciprocal Angstroms. + minimum_intensity + The minimum intensity of a reflection to be included in the profile. + debye_waller_factors + Maps element names to their temperature-dependent Debye-Waller factors. + """ + latt = phase.structure.lattice + + # Obtain crystallographic reciprocal lattice points within range + recip_latt = latt.reciprocal() + spot_indices, _, spot_distances = get_points_in_sphere( + recip_latt, reciprocal_radius + ) + + ##spot_indicies is a numpy.array of the hkls allowed in the recip radius + g_indices, multiplicities, g_hkls = get_intensities_params( + recip_latt, reciprocal_radius + ) + + i_hkl = get_kinematical_intensities( + phase.structure, + g_indices, + np.asarray(g_hkls), + prefactor=multiplicities, + scattering_params=self.scattering_params, + debye_waller_factors=debye_waller_factors, + ) + + if is_lattice_hexagonal(latt): + # Use Miller-Bravais indices for hexagonal lattices. + g_indices = np.array( + [ + g_indices[:, 0], + g_indices[:, 1], + g_indices[:, 0] - g_indices[:, 1], + g_indices[:, 2], + ] + ).T + + hkls_labels = ["".join([str(int(x)) for x in xs]) for xs in g_indices] + + peaks = {} + for l, i, g in zip(hkls_labels, i_hkl, g_hkls): + peaks[l] = [i, g] + + # Scale intensities so that the max intensity is 100. + + max_intensity = max([v[0] for v in peaks.values()]) + reciporical_spacing = [] + intensities = [] + hkls = [] + for k in peaks.keys(): + v = peaks[k] + if v[0] / max_intensity * 100 > minimum_intensity and (k != "000"): + reciporical_spacing.append(v[1]) + intensities.append(v[0]) + hkls.append(k) + + intensities = np.asarray(intensities) / max(intensities) * 100 + + return Simulation1D( + phase=phase, + reciprocal_spacing=reciporical_spacing, + intensities=intensities, + hkl=hkls, + reciprocal_radius=reciprocal_radius, + wavelength=self.wavelength, + ) + + def get_intersecting_reflections( + self, + recip: ReciprocalLatticeVector, + rot: np.ndarray, + wavelength: float, + max_excitation_error: float, + shape_factor_width: float = None, + with_direct_beam: bool = True, + ): + """Calculates the reciprocal lattice vectors that intersect the Ewald sphere. + + Parameters + ---------- + recip + The reciprocal lattice vectors to rotate. + rot + The rotation matrix to apply to the reciprocal lattice vectors. + wavelength + The wavelength of the electrons in Angstroms. + max_excitation_error + The cut-off for geometric excitation error in the z-direction + in units of reciprocal Angstroms. Spots with a larger distance + from the Ewald sphere are removed from the pattern. + Related to the extinction distance and roungly equal to 1/thickness. + shape_factor_width + Determines the width of the reciprocal rel-rod, for fine-grained + control. If not set will be set equal to max_excitation_error. + """ + initial_hkl = recip.hkl + rotated_vectors = recip.rotate_from_matrix(rot) + + if with_direct_beam: + rotated_vectors = np.vstack([rotated_vectors.data, [0, 0, 0]]) + initial_hkl = np.vstack([initial_hkl, [0, 0, 0]]) + else: + rotated_vectors = rotated_vectors.data + # Identify the excitation errors of all points (distance from point to Ewald sphere) + r_sphere = 1 / wavelength + r_spot = np.sqrt(np.sum(np.square(rotated_vectors[:, :2]), axis=1)) + z_spot = rotated_vectors[:, 2] + + z_sphere = -np.sqrt(r_sphere**2 - r_spot**2) + r_sphere + excitation_error = z_sphere - z_spot + + # determine the pre-selection reflections + if self.precession_angle == 0: + intersection = np.abs(excitation_error) < max_excitation_error + else: + # only consider points that intersect the ewald sphere at some point + # the center point of the sphere + P_z = r_sphere * np.cos(np.deg2rad(self.precession_angle)) + P_t = r_sphere * np.sin(np.deg2rad(self.precession_angle)) + # the extremes of the ewald sphere + z_surf_up = P_z - np.sqrt(r_sphere**2 - (r_spot + P_t) ** 2) + z_surf_do = P_z - np.sqrt(r_sphere**2 - (r_spot - P_t) ** 2) + intersection = (z_spot - max_excitation_error <= z_surf_up) & ( + z_spot + max_excitation_error >= z_surf_do + ) + + # select these reflections + intersected_vectors = rotated_vectors[intersection] + intersected_vectors = ReciprocalLatticeVector( + phase=recip.phase, xyz=intersected_vectors + ) + excitation_error = excitation_error[intersection] + r_spot = r_spot[intersection] + hkl = initial_hkl[intersection] + + if shape_factor_width is None: + shape_factor_width = max_excitation_error + # select and evaluate shape factor model + if self.precession_angle == 0: + # calculate shape factor + shape_factor = self.shape_factor_model( + excitation_error, shape_factor_width, **self.shape_factor_kwargs + ) + else: + if self.approximate_precession: + shape_factor = lorentzian_precession( + excitation_error, + shape_factor_width, + r_spot, + np.deg2rad(self.precession_angle), + ) + else: + shape_factor = _shape_factor_precession( + excitation_error, + r_spot, + np.deg2rad(self.precession_angle), + self.shape_factor_model, + shape_factor_width, + **self.shape_factor_kwargs, + ) + return intersected_vectors, hkl, shape_factor diff --git a/diffsims/simulations/__init__.py b/diffsims/simulations/__init__.py new file mode 100644 index 00000000..05864a60 --- /dev/null +++ b/diffsims/simulations/__init__.py @@ -0,0 +1,20 @@ +# -*- coding: utf-8 -*- +# Copyright 2017-2023 The diffsims developers +# +# This file is part of diffsims. +# +# diffsims is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# diffsims is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with diffsims. If not, see . + +from diffsims.simulations.simulation2d import Simulation2D +from diffsims.simulations.simulation1d import Simulation1D diff --git a/diffsims/simulations/simulation1d.py b/diffsims/simulations/simulation1d.py new file mode 100644 index 00000000..49391eff --- /dev/null +++ b/diffsims/simulations/simulation1d.py @@ -0,0 +1,88 @@ +# -*- coding: utf-8 -*- +# Copyright 2017-2024 The diffsims developers +# +# This file is part of diffsims. +# +# diffsims is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# diffsims is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with diffsims. If not, see . + +from typing import TYPE_CHECKING + +import numpy as np +import matplotlib.pyplot as plt +from orix.crystal_map import Phase +from diffsims.utils.sim_utils import get_electron_wavelength + +# to avoid circular imports +if TYPE_CHECKING: # pragma: no cover + from diffsims.generators.simulation_generator import SimulationGenerator + + +class Simulation1D: + """Holds the result of a 1D simulation for some phase""" + + def __init__( + self, + phase: Phase, + reciprocal_spacing: np.ndarray, + intensities: np.ndarray, + hkl: np.ndarray, + reciprocal_radius: float, + wavelength: float, + ): + """Initializes the DiffractionSimulation object with data values for + the coordinates, indices, intensities, calibration and offset. + + Parameters + ---------- + phase + The phase of the simulation + reciprocal_spacing + The spacing of the reciprocal lattice vectors in A^-1 + intensities + The intensities of the diffraction spots + hkl + The hkl indices of the diffraction spots + reciprocal_radius + The radius which the reciprocal lattice spacings are plotted out to + wavelength + The wavelength of the beam in A^-1 + """ + self.phase = phase + self.reciprocal_spacing = reciprocal_spacing + self.intensities = intensities + self.hkl = hkl + self.reciprocal_radius = reciprocal_radius + self.wavelength = wavelength + + def __repr__(self): + return f"Simulation1D(name: {self.phase.name}, wavelength: {self.wavelength})" + + @property + def theta(self): + return np.arctan2(np.array(self.reciprocal_spacing), 1 / self.wavelength) + + def plot(self, ax=None, annotate_peaks=False, fontsize=12, with_labels=True): + """Plots the 1D diffraction pattern,""" + if ax is None: + fig, ax = plt.subplots(1, 1) + for g, i, hkls in zip(self.reciprocal_spacing, self.intensities, self.hkl): + label = hkls + ax.plot([g, g], [0, i], color="k", linewidth=3, label=label) + if annotate_peaks: + ax.annotate(label, xy=[g, i], xytext=[g, i], fontsize=fontsize) + + if with_labels: + ax.set_xlabel("A ($^{-1}$)") + ax.set_ylabel("Intensities (scaled)") + return ax diff --git a/diffsims/simulations/simulation2d.py b/diffsims/simulations/simulation2d.py new file mode 100644 index 00000000..ed6bcdf0 --- /dev/null +++ b/diffsims/simulations/simulation2d.py @@ -0,0 +1,614 @@ +# -*- coding: utf-8 -*- +# Copyright 2017-2024 The diffsims developers +# +# This file is part of diffsims. +# +# diffsims is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# diffsims is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with diffsims. If not, see . + +from typing import Union, Sequence, TYPE_CHECKING, Any +import copy + +import numpy as np +import matplotlib.pyplot as plt +import math +from orix.crystal_map import Phase +from orix.quaternion import Rotation +from orix.vector import Vector3d + +from diffsims.crystallography.reciprocal_lattice_vector import ReciprocalLatticeVector +from diffsims.pattern.detector_functions import add_shot_and_point_spread + +# to avoid circular imports +if TYPE_CHECKING: # pragma: no cover + from diffsims.generators.simulation_generator import SimulationGenerator + + +class PhaseGetter: + """A class for getting the phases of a simulation library. + + Parameters + ---------- + phases : Sequence[Phase] + The phases in the library. + """ + + def __init__(self, simulation): + self.simulation = simulation + + def __getitem__(self, item): + all_phases = self.simulation.phases + if isinstance(all_phases, Phase): + raise ValueError("Only one phase in the simulation") + elif isinstance(item, str): + ind = [phase.name for phase in all_phases].index(item) + elif isinstance(item, (int, slice)): + ind = item + else: + raise ValueError("Item must be a string or integer") + new_coords = self.simulation.coordinates[ind] + new_rotations = self.simulation.rotations[ind] + new_phases = all_phases[ind] + return Simulation2D( + phases=new_phases, + coordinates=new_coords, + rotations=new_rotations, + simulation_generator=self.simulation.simulation_generator, + ) + + +class RotationGetter: + """A class for getting a Rotation of a simulation library. + + Parameters + ---------- + phases : Sequence[Phase] + The phases in the library. + """ + + def __init__(self, simulation): + self.simulation = simulation + + def __getitem__(self, item): + all_phases = self.simulation.phases + if self.simulation.current_size == 1: + raise ValueError("Only one rotation in the simulation") + elif isinstance(all_phases, Phase): # only one phase in the simulation + coords = self.simulation.coordinates[item] + phases = self.simulation.phases + rotations = self.simulation.rotations[item] + else: # multiple phases in the simulation + coords = [c[item] for c in self.simulation.coordinates] + phases = self.simulation.phases + rotations = [rot[item] for rot in self.simulation.rotations] + return Simulation2D( + phases=phases, + coordinates=coords, + rotations=rotations, + simulation_generator=self.simulation.simulation_generator, + ) + + +class Simulation2D: + """Holds the result of a kinematic diffraction simulation for some phase + and rotation. This class is iterable and can be used to iterate through + simulations of different phases and rotations. + """ + + def __init__( + self, + phases: Sequence[Phase], + coordinates: Union[ + Sequence[ReciprocalLatticeVector], + Sequence[Sequence[ReciprocalLatticeVector]], + ], + rotations: Union[Rotation, Sequence[Rotation]], + simulation_generator: "SimulationGenerator", + reciprocal_radius: float = 1.0, + ): + """Initializes the DiffractionSimulation object with data values for + the coordinates, indices, intensities, calibration and offset. + + Parameters + ---------- + coordinates + The list of ReciprocalLatticeVector objects for each phase and rotation. If there + are multiple phases, then this should be a list of lists of ReciprocalLatticeVector objects. + If there is only one phase, then this should be a list of ReciprocalLatticeVector objects. + rotations + The list of Rotation objects for each phase. If there are multiple phases, then this should + be a list of Rotation objects. If there is only one phase, then this should be a single + Rotation object. + phases + The list of Phase objects for each phase. If there is only one phase, then this should be + a single Phase object. + simulation_generator + The SimulationGenerator object used to generate the diffraction patterns. + + """ + # Basic data + if isinstance(rotations, Rotation) and rotations.size == 1: + if not isinstance(coordinates, ReciprocalLatticeVector): + raise ValueError( + "If there is only one rotation, then the coordinates must be a ReciprocalLatticeVector object" + ) + elif isinstance(rotations, Rotation): + coordinates = np.array(coordinates, dtype=object) + if coordinates.size != rotations.size: + raise ValueError( + f"The number of rotations: {rotations.size} must match the number of " + f"coordinates {coordinates.size}" + ) + else: # iterable of Rotation + rotations = np.array(rotations, dtype=object) + coordinates = np.array(coordinates, dtype=object) + phases = np.array(phases) + if rotations.size != phases.size: + raise ValueError( + f"The number of rotations: {rotations.size} must match the number of " + f"phases {phases.size}" + ) + + for r, c in zip(rotations, coordinates): + if isinstance(c, ReciprocalLatticeVector): + c = np.array( + [ + c, + ] + ) + if r.size != len(c): + raise ValueError( + f"The number of rotations: {r.size} must match the number of " + f"coordinates {c.shape[0]}" + ) + self.phases = phases + self.rotations = rotations + self.coordinates = coordinates + self.simulation_generator = simulation_generator + + # for interactive plotting and iterating through the Simulations + self.phase_index = 0 + self.rotation_index = 0 + self._rot_plot = None + self._diff_plot = None + self.reciporical_radius = reciprocal_radius + + # for slicing a simulation + self.iphase = PhaseGetter(self) + self.irot = RotationGetter(self) + + def get_simulation(self, item): + """Return the rotation and the phase index of the simulation""" + if self.has_multiple_phases: + cumsum = np.cumsum(self._num_rotations()) + ind = np.searchsorted(cumsum, item, side="right") + cumsum = np.insert(cumsum, 0, 0) + num_rot = cumsum[ind] + if self.has_multiple_rotations[ind]: + return ( + self.rotations[ind][item - num_rot], + ind, + self.coordinates[ind][item - num_rot], + ) + else: + return self.rotations[ind], ind, self.coordinates[ind] + elif self.has_multiple_rotations: + return self.rotations[item], 0, self.coordinates[item] + else: + return self.rotations[item], 0, self.coordinates + + def _num_rotations(self): + if self.has_multiple_phases: + return [r.size for r in self.rotations] + else: + return self.rotations.size + + def __iter__(self): + return self + + def __next__(self): + if self.phase_index == self.num_phases: + self.phase_index = 0 + raise StopIteration + else: + if self.has_multiple_phases: + coords = self.coordinates[self.phase_index] + else: + coords = self.coordinates + if self.has_multiple_rotations: + coords = coords[self.rotation_index] + else: + coords = coords + if self.rotation_index + 1 == self.current_size: + self.rotation_index = 0 + self.phase_index += 1 + else: + self.rotation_index += 1 + return coords + + @property + def current_size(self): + """Returns the number of rotations in the current phase""" + if self.has_multiple_phases: + return self.rotations[self.phase_index].size + else: + return self.rotations.size + + def deepcopy(self): + return copy.deepcopy(self) + + def _get_transformed_coordinates( + self, + angle: float, + center: Sequence = (0, 0), + mirrored: bool = False, + units: str = "real", + calibration: float = None, + ): + """Translate, rotate or mirror the pattern spot coordinates""" + + coords = self.get_current_coordinates() + + if units != "real": + center = np.array(center) + coords.data = coords.data / calibration + transformed_coords = coords + cx, cy = center + x = transformed_coords.data[:, 0] + y = transformed_coords.data[:, 1] + mirrored_factor = -1 if mirrored else 1 + theta = mirrored_factor * np.arctan2(y, x) + np.deg2rad(angle) + rd = np.sqrt(x**2 + y**2) + transformed_coords[:, 0] = rd * np.cos(theta) + cx + transformed_coords[:, 1] = rd * np.sin(theta) + cy + return transformed_coords + + @property + def current_phase(self): + if self.has_multiple_phases: + return self.phases[self.phase_index] + else: + return self.phases + + def rotate_shift_coordinates( + self, angle: float, center: Sequence = (0, 0), mirrored: bool = False + ): + """Rotate, flip or shift patterns in-plane + + Parameters + ---------- + angle + In plane rotation angle in degrees + center + Center coordinate of the patterns + mirrored + Mirror across the x-axis + """ + coords_new = self._get_transformed_coordinates( + angle, center, mirrored, units="real" + ) + return coords_new + + def polar_flatten_simulations(self, radial_axes=None, azimuthal_axes=None): + """Flattens the simulations into polar coordinates for use in template matching. + The resulting arrays are of shape (n_simulations, n_spots) where n_spots is the + maximum number of spots in any simulation. + + + Returns + ------- + r_templates, theta_templates, intensities_templates + """ + + flattened_vectors = [sim for sim in self] + max_num_spots = max([v.size for v in flattened_vectors]) + + r_templates = np.zeros((len(flattened_vectors), max_num_spots)) + 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() + if radial_axes is not None and azimuthal_axes is not None: + r = get_closest(radial_axes, r) + t = get_closest(azimuthal_axes, t) + r = r[r < len(radial_axes)] + t = t[t < len(azimuthal_axes)] + r_templates[i, : len(r)] = r + theta_templates[i, : len(t)] = t + intensities_templates[i, : len(v.intensity)] = v.intensity + if radial_axes is not None and azimuthal_axes is not None: + r_templates = np.array(r_templates, dtype=int) + theta_templates = np.array(theta_templates, dtype=int) + + return r_templates, theta_templates, intensities_templates + + def get_diffraction_pattern( + self, + shape=None, + sigma=10, + direct_beam_position=None, + in_plane_angle=0, + calibration=0.01, + mirrored=False, + ): + """Returns the diffraction data as a numpy array with + two-dimensional Gaussians representing each diffracted peak. Should only + be used for qualitative work. + + Parameters + ---------- + shape : tuple of ints + The size of a side length (in pixels) + sigma : float + Standard deviation of the Gaussian function to be plotted (in pixels). + direct_beam_position: 2-tuple of ints, optional + The (x,y) coordinate in pixels of the direct beam. Defaults to + the center of the image. + in_plane_angle: float, optional + In plane rotation of the pattern in degrees + mirrored: bool, optional + Whether the pattern should be flipped over the x-axis, + corresponding to the inverted orientation + + Returns + ------- + diffraction-pattern : numpy.array + The simulated electron diffraction pattern, normalised. + + Notes + ----- + If don't know the exact calibration of your diffraction signal using 1e-2 + produces reasonably good patterns when the lattice parameters are on + the order of 0.5nm and the default size and sigma are used. + """ + if direct_beam_position is None: + direct_beam_position = (shape[1] // 2, shape[0] // 2) + transformed = self._get_transformed_coordinates( + in_plane_angle, + direct_beam_position, + mirrored, + units="pixel", + calibration=calibration, + ) + in_frame = ( + (transformed.data[:, 0] >= 0) + & (transformed.data[:, 0] < shape[1]) + & (transformed.data[:, 1] >= 0) + & (transformed.data[:, 1] < shape[0]) + ) + spot_coords = transformed.data[in_frame].astype(int) + + spot_intens = transformed.intensity[in_frame] + pattern = np.zeros(shape) + # checks that we have some spots + if spot_intens.shape[0] == 0: + return pattern + else: + pattern[spot_coords[:, 0], spot_coords[:, 1]] = spot_intens + pattern = add_shot_and_point_spread(pattern.T, sigma, shot_noise=False) + return np.divide(pattern, np.max(pattern)) + + @property + def num_phases(self): + """Returns the number of phases in the simulation""" + if hasattr(self.phases, "__len__"): + return len(self.phases) + else: + return 1 + + @property + def has_multiple_phases(self): + """Returns True if the simulation has multiple phases""" + return self.num_phases > 1 + + @property + def has_multiple_rotations(self): + """Returns True if the simulation has multiple rotations""" + if isinstance(self.rotations, Rotation): + return self.rotations.size > 1 + else: + return [r.size > 1 for r in self.rotations] + + def get_current_coordinates(self): + """Returns the coordinates of the current phase and rotation""" + if self.has_multiple_phases: + return copy.deepcopy( + self.coordinates[self.phase_index][self.rotation_index] + ) + elif not self.has_multiple_phases and self.has_multiple_rotations: + return copy.deepcopy(self.coordinates[self.rotation_index]) + else: + return copy.deepcopy(self.coordinates) + + def get_current_rotation(self): + """Returns the matrix for the current matrix""" + if self.has_multiple_phases: + return copy.deepcopy( + self.rotations[self.phase_index].to_matrix()[self.rotation_index] + ) + else: + return copy.deepcopy(self.rotations.to_matrix()[self.rotation_index]) + + def plot_rotations(self, beam_direction: Vector3d = Vector3d.zvector()): + """Plots the rotations of the current phase in stereographic projection""" + if self.has_multiple_phases: + rots = self.rotations[self.phase_index] + else: + rots = self.rotations + vect_rot = rots * beam_direction + facecolor = ["k"] * rots.size + facecolor[self.rotation_index] = "r" # highlight the current rotation + fig = vect_rot.scatter( + grid=True, + facecolor=facecolor, + return_figure=True, + ) + pointer = vect_rot[self.rotation_index] + _plot = fig.axes[0] + _plot.scatter(pointer.data[0][0], pointer.data[0][1], color="r") + _plot = fig.axes[0] + _plot.set_title("Rotations" + self.current_phase.name) + + def plot( + self, + size_factor=1, + direct_beam_position=None, + in_plane_angle=0, + mirrored=False, + units="real", + show_labels=False, + label_offset=(0, 0), + label_formatting=None, + min_label_intensity=0.1, + include_direct_beam=True, + calibration=0.1, + ax=None, + **kwargs, + ): + """A quick-plot function for a simulation of spots + + Parameters + ---------- + size_factor : float, optional + linear spot size scaling, default to 1 + direct_beam_position: 2-tuple of ints, optional + The (x,y) coordinate in pixels of the direct beam. Defaults to + the center of the image. + in_plane_angle: float, optional + In plane rotation of the pattern in degrees + mirrored: bool, optional + Whether the pattern should be flipped over the x-axis, + corresponding to the inverted orientation + units : str, optional + 'real' or 'pixel', only changes scalebars, falls back on 'real', the default + show_labels : bool, optional + draw the miller indices near the spots + label_offset : 2-tuple, optional + the relative location of the spot labels. Does nothing if `show_labels` + is False. + label_formatting : dict, optional + keyword arguments passed to `ax.text` for drawing the labels. Does + nothing if `show_labels` is False. + min_label_intensity : float, optional + minimum intensity for a spot to be labelled + include_direct_beam : bool, optional + whether to include the direct beam in the plot + ax : matplotlib Axes, optional + axes on which to draw the pattern. If `None`, a new axis is created + **kwargs : + passed to ax.scatter() method + + Returns + ------- + ax,sp + + Notes + ----- + spot size scales with the square root of the intensity. + """ + if label_formatting is None: + label_formatting = {} + if direct_beam_position is None: + direct_beam_position = (0, 0) + if ax is None: + _, ax = plt.subplots() + ax.set_aspect("equal") + coords = self._get_transformed_coordinates( + in_plane_angle, + direct_beam_position, + mirrored, + units=units, + calibration=calibration, + ) + if include_direct_beam: + spots = coords.data[:, :2] + spots = np.concatenate((spots, np.array([direct_beam_position]))) + intensity = np.concatenate((coords.intensity, np.array([1]))) + else: + spots = coords.data[:, :2] + intensity = coords.intensity + sp = ax.scatter( + spots[:, 0], + spots[:, 1], + s=size_factor * np.sqrt(intensity), + **kwargs, + ) + ax.set_xlim(-self.reciporical_radius, self.reciporical_radius) + ax.set_ylim(-self.reciporical_radius, self.reciporical_radius) + + if show_labels: + millers = np.round( + np.matmul( + np.matmul(coords.data, self.get_current_rotation().T), + coords.phase.structure.lattice.base.T, + ) + ).astype(np.int16) + # only label the points inside the axes + xlim = ax.get_xlim() + ylim = ax.get_ylim() + condition = ( + (coords.data[:, 0] > min(xlim)) + & (coords.data[:, 0] < max(xlim)) + & (coords.data[:, 1] > min(ylim)) + & (coords.data[:, 1] < max(ylim)) + ) + millers = millers[condition] + coords = coords.data[condition] + # default alignment options + if ( + "ha" not in label_offset + and "horizontalalignment" not in label_formatting + ): + label_formatting["ha"] = "center" + if "va" not in label_offset and "verticalalignment" not in label_formatting: + label_formatting["va"] = "center" + for miller, coordinate, inten in zip(millers, coords, intensity): + if np.isnan(inten) or inten > min_label_intensity: + label = "(" + for index in miller: + if index < 0: + label += r"$\bar{" + str(abs(index)) + r"}$" + else: + label += str(abs(index)) + label += " " + label = label[:-1] + ")" + ax.text( + coordinate[0] + label_offset[0], + coordinate[1] + label_offset[1], + label, + **label_formatting, + ) + if units == "real": + ax.set_xlabel(r"$\AA^{-1}$") + ax.set_ylabel(r"$\AA^{-1}$") + else: + ax.set_xlabel("pixels") + ax.set_ylabel("pixels") + return ax, sp + + +def get_closest(array, values): + # make sure array is a numpy array + array = np.array(array) + + # get insert positions + idxs = np.searchsorted(array, values, side="left") + + # find indexes where previous index is closer + prev_idx_is_less = (idxs == len(array)) | ( + np.fabs(values - array[np.maximum(idxs - 1, 0)]) + < np.fabs(values - array[np.minimum(idxs, len(array) - 1)]) + ) + idxs[prev_idx_is_less] -= 1 + + return idxs From 82827218b29296f143b0de0df9c7e7529fc9ef54 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 12 Apr 2024 17:39:42 -0500 Subject: [PATCH 02/36] New Feature: Add intensity to reciprocal lattice vector --- .../reciprocal_lattice_vector.py | 40 ++++++++++++++++++- 1 file changed, 38 insertions(+), 2 deletions(-) diff --git a/diffsims/crystallography/reciprocal_lattice_vector.py b/diffsims/crystallography/reciprocal_lattice_vector.py index 8ce16ba7..0fcf4d62 100644 --- a/diffsims/crystallography/reciprocal_lattice_vector.py +++ b/diffsims/crystallography/reciprocal_lattice_vector.py @@ -18,6 +18,7 @@ from collections import defaultdict from copy import deepcopy +from typing import Tuple from diffpy.structure.symmetryutilities import expandPosition from diffpy.structure import Structure @@ -122,10 +123,11 @@ def __init__(self, phase, xyz=None, hkl=None, hkil=None): self._theta = np.full(self.shape, np.nan) self._structure_factor = np.full(self.shape, np.nan, dtype="complex128") + self._intensity = np.full(self.shape, np.nan) def __getitem__(self, key): - miller_new = self.to_miller().__getitem__(key) - rlv_new = self.from_miller(miller_new) + new_data = self.data[key] + rlv_new = self.__class__(self.phase, xyz=new_data) if np.isnan(self.structure_factor).all(): rlv_new._structure_factor = np.full( @@ -139,6 +141,18 @@ def __getitem__(self, key): else: rlv_new._theta = self.theta[key] + if np.isnan(self.intensity).all(): + rlv_new._intensity = np.full(rlv_new.shape, np.nan) + else: + slic = self.intensity[key] + if not hasattr(slic, "__len__"): + slic = np.array( + [ + slic, + ] + ) + rlv_new._intensity = slic + return rlv_new def __repr__(self): @@ -502,6 +516,28 @@ def scattering_parameter(self): return 0.5 * self.gspacing + @property + def intensity(self): + return self._intensity + + @intensity.setter + def intensity(self, value): + if not hasattr(value, "__len__"): + value = np.array( + [ + value, + ] + * self.size + ) + if len(value) != self.size: + raise ValueError("Length of intensity array must match number of vectors") + self._intensity = np.array(value) + + def rotate_from_matrix(self, rotation_matrix): + return ReciprocalLatticeVector( + phase=self.phase, xyz=np.matmul(rotation_matrix.T, self.data.T).T + ) + @property def structure_factor(self): r"""Kinematical structure factors :math:`F`. From bfe0505123e7efe12aeb75eb3a9a1abcf89f3ce5 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 12 Apr 2024 17:43:42 -0500 Subject: [PATCH 03/36] Testing: Add tests for adding intensity of reciprocal_lattive_vector.py --- .../test_reciprocal_lattice_vector.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py b/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py index d8717107..41f29fb3 100644 --- a/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py +++ b/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py @@ -108,6 +108,25 @@ def test_repr(self, ferrite_phase): "[[ 1. 1. 0.]", ] + def test_add_intensity(self, ferrite_phase): + rlv = ReciprocalLatticeVector.from_min_dspacing(ferrite_phase, 1.5) + rlv.intensity = 1 + assert isinstance(rlv.intensity, np.ndarray) + assert np.allclose(rlv.intensity, np.ones(rlv.size)) + + def test_add_intensity_error(self, ferrite_phase): + rlv = ReciprocalLatticeVector.from_min_dspacing(ferrite_phase, 1.5) + with pytest.raises(ValueError): + rlv.intensity = [0, 1] + + @pytest.mark.parametrize("degrees", [True, False]) + def test_to_polar(self, ferrite_phase, degrees): + rlv = ReciprocalLatticeVector.from_min_dspacing(ferrite_phase, 1.5) + r, theta, z = rlv.to_polar(degrees=degrees) + assert r.shape == (rlv.size,) + assert theta.shape == (rlv.size,) + assert z.shape == (rlv.size,) + def test_get_item(self, ferrite_phase): """Indexing gives desired vectors and properties carry over.""" rlv = ReciprocalLatticeVector.from_min_dspacing(ferrite_phase, 1.5) From ac97ee635d46f0ef77bc536ba7f6bed810727b82 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 12 Apr 2024 17:51:14 -0500 Subject: [PATCH 04/36] New Feature: Added zerobeam to reciprocal lattice vector --- diffsims/crystallography/reciprocal_lattice_vector.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/diffsims/crystallography/reciprocal_lattice_vector.py b/diffsims/crystallography/reciprocal_lattice_vector.py index 0fcf4d62..e3de582a 100644 --- a/diffsims/crystallography/reciprocal_lattice_vector.py +++ b/diffsims/crystallography/reciprocal_lattice_vector.py @@ -1106,7 +1106,7 @@ def from_highest_hkl(cls, phase, hkl): return cls(phase, hkl=idx).unique() @classmethod - def from_min_dspacing(cls, phase, min_dspacing=0.7): + def from_min_dspacing(cls, phase, min_dspacing=0.7, include_zero_beam=False): """Create a set of unique reciprocal lattice vectors with a a direct space interplanar spacing greater than a lower threshold. @@ -1119,6 +1119,8 @@ def from_min_dspacing(cls, phase, min_dspacing=0.7): Smallest interplanar spacing to consider. Default is 0.7, in the unit used to define the lattice parameters in ``phase``, which is assumed to be Ångström. + include_zero_beam: bool + If ``True``, include the zero beam in the set of vectors. Examples -------- @@ -1164,6 +1166,8 @@ def from_min_dspacing(cls, phase, min_dspacing=0.7): dspacing = 1 / phase.structure.lattice.rnorm(hkl) idx = dspacing >= min_dspacing hkl = hkl[idx] + if include_zero_beam: + hkl = np.vstack((hkl, np.zeros(3, dtype=int))) return cls(phase, hkl=hkl).unique() @classmethod From b0b67fa901b0eab0c25cb361aaddc7931f629595 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 12 Apr 2024 17:51:41 -0500 Subject: [PATCH 05/36] New Feature: Moved Shape factor calculation to shape_factor_models.py --- diffsims/utils/shape_factor_models.py | 51 +++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/diffsims/utils/shape_factor_models.py b/diffsims/utils/shape_factor_models.py index b8dc7c20..6aab0902 100644 --- a/diffsims/utils/shape_factor_models.py +++ b/diffsims/utils/shape_factor_models.py @@ -17,6 +17,7 @@ # along with diffsims. If not, see . import numpy as np +from scipy.integrate import quad __all__ = [ @@ -217,3 +218,53 @@ def lorentzian_precession( z = np.sqrt(u**2 + 4 * sigma**2 * excitation_error**2) fac = (sigma / np.pi) * np.sqrt(2 * (u + z) / z**2) return fac + + +def _shape_factor_precession( + excitation_error, r_spot, phi, shape_function, max_excitation, **kwargs +): + """ + The rel-rod shape factors for reflections taking into account + precession + + Parameters + ---------- + excitation_error : np.ndarray (N,) + An array of excitation errors + r_spot : np.ndarray (N,) + An array representing the distance of spots from the z-axis in A^-1 + phi : float + The precession angle in radians + shape_function : callable + A function that describes the influence from the rel-rods. Should be + in the form func(excitation_error: np.ndarray, max_excitation: float, + **kwargs) + max_excitation : float + Parameter to describe the "extent" of the rel-rods. + + Other parameters + ---------------- + ** kwargs: passed directly to shape_function + + Notes + ----- + * We calculate excitation_error as z_spot - z_sphere so that it is + negative when the spot is outside the ewald sphere and positive when inside + conform W&C chapter 12, section 12.6 + * We assume that the sample is a thin infinitely wide slab perpendicular + to the optical axis, so that the shape factor function only depends on the + distance from each spot to the Ewald sphere parallel to the optical axis. + """ + shf = np.zeros(excitation_error.shape) + # loop over all spots + 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] + return shf From 6ec15b8514a96179bbec7f1c8979a5bc5fcaf31b Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 12 Apr 2024 17:52:04 -0500 Subject: [PATCH 06/36] Testing: Added testing for simulation_generator. --- .../simulations/test_simulation_generator.py | 252 ++++++++++++++++++ 1 file changed, 252 insertions(+) create mode 100644 diffsims/tests/simulations/test_simulation_generator.py diff --git a/diffsims/tests/simulations/test_simulation_generator.py b/diffsims/tests/simulations/test_simulation_generator.py new file mode 100644 index 00000000..06c6ddea --- /dev/null +++ b/diffsims/tests/simulations/test_simulation_generator.py @@ -0,0 +1,252 @@ +# -*- coding: utf-8 -*- +# Copyright 2017-2023 The diffsims developers +# +# This file is part of diffsims. +# +# diffsims is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# diffsims is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with diffsims. If not, see . + +import numpy as np +import pytest + +import diffpy.structure +from orix.crystal_map import Phase +from orix.quaternion import Rotation + +from diffsims.generators.simulation_generator import SimulationGenerator +from diffsims.utils.shape_factor_models import ( + linear, + binary, + sin2c, + atanc, + lorentzian, + _shape_factor_precession, +) +from diffsims.simulations import Simulation1D +from diffsims.utils.sim_utils import is_lattice_hexagonal + + +@pytest.fixture(params=[(300)]) +def diffraction_calculator(request): + return SimulationGenerator(request.param) + + +@pytest.fixture(scope="module") +def diffraction_calculator_precession_full(): + return SimulationGenerator(300, precession_angle=0.5, approximate_precession=False) + + +@pytest.fixture(scope="module") +def diffraction_calculator_precession_simple(): + return SimulationGenerator(300, precession_angle=0.5, approximate_precession=True) + + +def local_excite(excitation_error, maximum_excitation_error, t): + return (np.sin(t) * excitation_error) / maximum_excitation_error + + +@pytest.fixture(scope="module") +def diffraction_calculator_custom(): + return SimulationGenerator(300, shape_factor_model=local_excite, t=0.2) + + +def make_phase(lattice_parameter=None): + """ + We construct an Fd-3m silicon (with lattice parameter 5.431 as a default) + """ + if lattice_parameter is not None: + a = lattice_parameter + else: + a = 5.431 + latt = diffpy.structure.lattice.Lattice(a, a, a, 90, 90, 90) + # TODO - Make this construction with internal diffpy syntax + atom_list = [] + for coords in [[0, 0, 0], [0.5, 0, 0.5], [0, 0.5, 0.5], [0.5, 0.5, 0]]: + x, y, z = coords[0], coords[1], coords[2] + atom_list.append( + diffpy.structure.atom.Atom(atype="Si", xyz=[x, y, z], lattice=latt) + ) # Motif part A + atom_list.append( + diffpy.structure.atom.Atom( + atype="Si", xyz=[x + 0.25, y + 0.25, z + 0.25], lattice=latt + ) + ) # Motif part B + struct = diffpy.structure.Structure(atoms=atom_list, lattice=latt) + p = Phase(structure=struct, space_group=227) + return p + + +@pytest.fixture() +def local_structure(): + return make_phase() + + +@pytest.mark.parametrize("model", [binary, linear, atanc, sin2c, lorentzian]) +def test_shape_factor_precession(model): + excitation = np.array([-0.1, 0.1]) + r = np.array([1, 5]) + _ = _shape_factor_precession(excitation, r, 0.5, model, 0.1) + + +def test_linear_shape_factor(): + excitation = np.array([-2, -1, -0.5, 0, 0.5, 1, 2]) + totest = linear(excitation, 1) + np.testing.assert_allclose(totest, np.array([0, 0, 0.5, 1, 0.5, 0, 0])) + np.testing.assert_allclose(linear(0.5, 1), 0.5) + + +@pytest.mark.parametrize( + "model, expected", + [("linear", linear), ("lorentzian", lorentzian), (binary, binary)], +) +def test_diffraction_generator_init(model, expected): + generator = SimulationGenerator(300, shape_factor_model=model) + assert generator.shape_factor_model == expected + + +class TestDiffractionCalculator: + def test_init(self, diffraction_calculator: SimulationGenerator): + assert diffraction_calculator.scattering_params == "lobato" + assert diffraction_calculator.precession_angle == 0 + assert diffraction_calculator.shape_factor_model == lorentzian + assert diffraction_calculator.approximate_precession == True + assert diffraction_calculator.minimum_intensity == 1e-20 + + def test_matching_results( + self, diffraction_calculator: SimulationGenerator, local_structure + ): + diffraction = diffraction_calculator.calculate_ed_data( + local_structure, reciprocal_radius=5.0 + ) + assert diffraction.coordinates.size == 69 + + def test_precession_simple( + self, diffraction_calculator_precession_simple, local_structure + ): + diffraction = diffraction_calculator_precession_simple.calculate_ed_data( + local_structure, + reciprocal_radius=5.0, + ) + assert diffraction.coordinates.size == 249 + + def test_precession_full( + self, diffraction_calculator_precession_full, local_structure + ): + diffraction = diffraction_calculator_precession_full.calculate_ed_data( + local_structure, + reciprocal_radius=5.0, + ) + assert diffraction.coordinates.size == 249 + + def test_custom_shape_func(self, diffraction_calculator_custom, local_structure): + diffraction = diffraction_calculator_custom.calculate_ed_data( + local_structure, + reciprocal_radius=5.0, + ) + assert diffraction.coordinates.size == 52 + + def test_appropriate_scaling(self, diffraction_calculator: SimulationGenerator): + """Tests that doubling the unit cell halves the pattern spacing.""" + silicon = make_phase(5) + big_silicon = make_phase(10) + diffraction = diffraction_calculator.calculate_ed_data( + phase=silicon, reciprocal_radius=5.0 + ) + big_diffraction = diffraction_calculator.calculate_ed_data( + phase=big_silicon, reciprocal_radius=5.0 + ) + indices = [tuple(i) for i in diffraction.coordinates.hkl] + big_indices = [tuple(i) for i in big_diffraction.coordinates.hkl] + assert (2, 2, 0) in indices + assert (2, 2, 0) in big_indices + coordinates = diffraction.coordinates[indices.index((2, 2, 0))] + big_coordinates = big_diffraction.coordinates[big_indices.index((2, 2, 0))] + assert np.allclose(coordinates.data, big_coordinates.data * 2) + + def test_appropriate_intensities(self, diffraction_calculator, local_structure): + """Tests the central beam is strongest.""" + diffraction = diffraction_calculator.calculate_ed_data( + local_structure, reciprocal_radius=0.5, with_direct_beam=False + ) # direct beam doesn't work + indices = [tuple(np.round(i).astype(int)) for i in diffraction.coordinates.hkl] + central_beam = indices.index((0, 1, 0)) + + smaller = np.greater_equal( + diffraction.coordinates.intensity[central_beam], + diffraction.coordinates.intensity, + ) + assert np.all(smaller) + + def test_shape_factor_strings(self, diffraction_calculator, local_structure): + _ = diffraction_calculator.calculate_ed_data( + local_structure, + ) + + def test_shape_factor_custom(self, diffraction_calculator, local_structure): + t1 = diffraction_calculator.calculate_ed_data( + local_structure, max_excitation_error=0.02 + ) + t2 = diffraction_calculator.calculate_ed_data( + local_structure, max_excitation_error=0.4 + ) + # softly makes sure the two sims are different + assert np.sum(t1.coordinates.intensity) != np.sum(t2.coordinates.intensity) + + @pytest.mark.parametrize("is_hex", [True, False]) + def test_simulate_1d(self, is_hex): + generator = SimulationGenerator(300) + phase = make_phase() + if is_hex: + phase.structure.lattice.a = phase.structure.lattice.b + phase.structure.lattice.alpha = 90 + phase.structure.lattice.beta = 90 + phase.structure.lattice.gamma = 120 + assert is_lattice_hexagonal(phase.structure.lattice) + else: + assert not is_lattice_hexagonal(phase.structure.lattice) + sim = generator.calculate_diffraction1d(phase, 0.5) + assert isinstance(sim, Simulation1D) + + assert len(sim.intensities) == len(sim.reciprocal_spacing) + assert len(sim.intensities) == len(sim.hkl) + for h in sim.hkl: + h = h.replace("-", "") + if is_hex: + assert len(h) == 4 + else: + assert len(h) == 3 + + +def test_multiphase_multirotation_simulation(): + generator = SimulationGenerator(300) + silicon = make_phase(5) + big_silicon = make_phase(10) + rot = Rotation.from_euler([[0, 0, 0], [0.1, 0.1, 0.1]]) + rot2 = Rotation.from_euler([[0, 0, 0], [0.1, 0.1, 0.1], [0.2, 0.2, 0.2]]) + sim = generator.calculate_ed_data([silicon, big_silicon], rotation=[rot, rot2]) + + +@pytest.mark.parametrize("scattering_param", ["lobato", "xtables"]) +def test_param_check(scattering_param): + generator = SimulationGenerator(300, scattering_params=scattering_param) + + +@pytest.mark.xfail(raises=NotImplementedError) +def test_invalid_scattering_params(): + scattering_param = "_empty" + generator = SimulationGenerator(300, scattering_params=scattering_param) + + +@pytest.mark.xfail(faises=NotImplementedError) +def test_invalid_shape_model(): + generator = SimulationGenerator(300, shape_factor_model="dracula") From f3bd483ea80b2f1dc488f73d7fea33f5c30ac0a1 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 12 Apr 2024 18:07:56 -0500 Subject: [PATCH 07/36] Testing: Testing the zero Beam is most intense. --- diffsims/tests/simulations/test_simulation_generator.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/diffsims/tests/simulations/test_simulation_generator.py b/diffsims/tests/simulations/test_simulation_generator.py index 06c6ddea..5a385813 100644 --- a/diffsims/tests/simulations/test_simulation_generator.py +++ b/diffsims/tests/simulations/test_simulation_generator.py @@ -176,10 +176,10 @@ def test_appropriate_scaling(self, diffraction_calculator: SimulationGenerator): def test_appropriate_intensities(self, diffraction_calculator, local_structure): """Tests the central beam is strongest.""" diffraction = diffraction_calculator.calculate_ed_data( - local_structure, reciprocal_radius=0.5, with_direct_beam=False + local_structure, reciprocal_radius=0.5, with_direct_beam=True ) # direct beam doesn't work indices = [tuple(np.round(i).astype(int)) for i in diffraction.coordinates.hkl] - central_beam = indices.index((0, 1, 0)) + central_beam = indices.index((0, 0, 0)) smaller = np.greater_equal( diffraction.coordinates.intensity[central_beam], From d7c86ec1853fcd0abcbe9288a0f405e0f4b736cf Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 17 Apr 2024 08:24:38 -0500 Subject: [PATCH 08/36] Min Version: Bump orix min version --- .github/workflows/build.yml | 2 +- setup.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 2c9b596c..86dd70a4 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -46,7 +46,7 @@ jobs: - os: ubuntu-latest python-version: 3.8 OLDEST_SUPPORTED_VERSION: true - DEPENDENCIES: diffpy.structure==3.0.2 matplotlib==3.5 numpy==1.17.3 orix==0.9.0 scipy==1.8 tqdm==4.9 + DEPENDENCIES: diffpy.structure==3.0.2 matplotlib==3.5 numpy==1.17.3 orix==0.11.0 scipy==1.8 tqdm==4.9 LABEL: -oldest steps: - uses: actions/checkout@v4 diff --git a/setup.py b/setup.py index 1b9de6c6..0a35c2d1 100644 --- a/setup.py +++ b/setup.py @@ -78,7 +78,7 @@ "matplotlib >= 3.3", "numba", "numpy >= 1.17.3", - "orix >= 0.9", + "orix >= 0.11", "psutil", "scipy >= 1.8", "tqdm >= 4.9", From ef03b393b05563a3214b84b7a383fbfdbfe6e06a Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 17 Apr 2024 08:39:46 -0500 Subject: [PATCH 09/36] Packaging: Include function and classes in init file --- diffsims/generators/__init__.py | 2 ++ diffsims/generators/simulation_generator.py | 2 +- diffsims/simulations/__init__.py | 5 +++++ 3 files changed, 8 insertions(+), 1 deletion(-) diff --git a/diffsims/generators/__init__.py b/diffsims/generators/__init__.py index 56b17d09..ab61bf8a 100644 --- a/diffsims/generators/__init__.py +++ b/diffsims/generators/__init__.py @@ -26,6 +26,7 @@ rotation_list_generators, sphere_mesh_generators, zap_map_generator, + simulation_generator, ) __all__ = [ @@ -33,5 +34,6 @@ "library_generator", "rotation_list_generators", "sphere_mesh_generators", + "simulation_generator", "zap_map_generator", ] diff --git a/diffsims/generators/simulation_generator.py b/diffsims/generators/simulation_generator.py index c88c9ee3..a43442a8 100644 --- a/diffsims/generators/simulation_generator.py +++ b/diffsims/generators/simulation_generator.py @@ -124,7 +124,7 @@ def __init__( def wavelength(self): return get_electron_wavelength(self.accelerating_voltage) - def calculate_ed_data( + def calculate_diffraction2d( self, phase: Union[Phase, Sequence[Phase]], rotation: Union[Rotation, Sequence[Rotation]] = Rotation.from_euler( diff --git a/diffsims/simulations/__init__.py b/diffsims/simulations/__init__.py index 05864a60..2612e7d1 100644 --- a/diffsims/simulations/__init__.py +++ b/diffsims/simulations/__init__.py @@ -18,3 +18,8 @@ from diffsims.simulations.simulation2d import Simulation2D from diffsims.simulations.simulation1d import Simulation1D + +__all__ = [ + "Simulation1D", + "Simulation2D", +] From d1839a532b7cc3cac7edb5f9a889d1b7120723b3 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 17 Apr 2024 08:42:52 -0500 Subject: [PATCH 10/36] Testing: Clean up tests due to changed function name --- .../simulations/test_simulation_generator.py | 24 ++++++++++--------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/diffsims/tests/simulations/test_simulation_generator.py b/diffsims/tests/simulations/test_simulation_generator.py index 5a385813..bfe88654 100644 --- a/diffsims/tests/simulations/test_simulation_generator.py +++ b/diffsims/tests/simulations/test_simulation_generator.py @@ -125,7 +125,7 @@ def test_init(self, diffraction_calculator: SimulationGenerator): def test_matching_results( self, diffraction_calculator: SimulationGenerator, local_structure ): - diffraction = diffraction_calculator.calculate_ed_data( + diffraction = diffraction_calculator.calculate_diffraction2d( local_structure, reciprocal_radius=5.0 ) assert diffraction.coordinates.size == 69 @@ -133,7 +133,7 @@ def test_matching_results( def test_precession_simple( self, diffraction_calculator_precession_simple, local_structure ): - diffraction = diffraction_calculator_precession_simple.calculate_ed_data( + diffraction = diffraction_calculator_precession_simple.calculate_diffraction2d( local_structure, reciprocal_radius=5.0, ) @@ -142,14 +142,14 @@ def test_precession_simple( def test_precession_full( self, diffraction_calculator_precession_full, local_structure ): - diffraction = diffraction_calculator_precession_full.calculate_ed_data( + diffraction = diffraction_calculator_precession_full.calculate_diffraction2d( local_structure, reciprocal_radius=5.0, ) assert diffraction.coordinates.size == 249 def test_custom_shape_func(self, diffraction_calculator_custom, local_structure): - diffraction = diffraction_calculator_custom.calculate_ed_data( + diffraction = diffraction_calculator_custom.calculate_diffraction2d( local_structure, reciprocal_radius=5.0, ) @@ -159,10 +159,10 @@ def test_appropriate_scaling(self, diffraction_calculator: SimulationGenerator): """Tests that doubling the unit cell halves the pattern spacing.""" silicon = make_phase(5) big_silicon = make_phase(10) - diffraction = diffraction_calculator.calculate_ed_data( + diffraction = diffraction_calculator.calculate_diffraction2d( phase=silicon, reciprocal_radius=5.0 ) - big_diffraction = diffraction_calculator.calculate_ed_data( + big_diffraction = diffraction_calculator.calculate_diffraction2d( phase=big_silicon, reciprocal_radius=5.0 ) indices = [tuple(i) for i in diffraction.coordinates.hkl] @@ -175,7 +175,7 @@ def test_appropriate_scaling(self, diffraction_calculator: SimulationGenerator): def test_appropriate_intensities(self, diffraction_calculator, local_structure): """Tests the central beam is strongest.""" - diffraction = diffraction_calculator.calculate_ed_data( + diffraction = diffraction_calculator.calculate_diffraction2d( local_structure, reciprocal_radius=0.5, with_direct_beam=True ) # direct beam doesn't work indices = [tuple(np.round(i).astype(int)) for i in diffraction.coordinates.hkl] @@ -188,15 +188,15 @@ def test_appropriate_intensities(self, diffraction_calculator, local_structure): assert np.all(smaller) def test_shape_factor_strings(self, diffraction_calculator, local_structure): - _ = diffraction_calculator.calculate_ed_data( + _ = diffraction_calculator.calculate_diffraction2d( local_structure, ) def test_shape_factor_custom(self, diffraction_calculator, local_structure): - t1 = diffraction_calculator.calculate_ed_data( + t1 = diffraction_calculator.calculate_diffraction2d( local_structure, max_excitation_error=0.02 ) - t2 = diffraction_calculator.calculate_ed_data( + t2 = diffraction_calculator.calculate_diffraction2d( local_structure, max_excitation_error=0.4 ) # softly makes sure the two sims are different @@ -233,7 +233,9 @@ def test_multiphase_multirotation_simulation(): big_silicon = make_phase(10) rot = Rotation.from_euler([[0, 0, 0], [0.1, 0.1, 0.1]]) rot2 = Rotation.from_euler([[0, 0, 0], [0.1, 0.1, 0.1], [0.2, 0.2, 0.2]]) - sim = generator.calculate_ed_data([silicon, big_silicon], rotation=[rot, rot2]) + sim = generator.calculate_diffraction2d( + [silicon, big_silicon], rotation=[rot, rot2] + ) @pytest.mark.parametrize("scattering_param", ["lobato", "xtables"]) From 3567308e9f5edf067a48f51c6774dc3cd14e3533 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 17 Apr 2024 09:15:35 -0500 Subject: [PATCH 11/36] Testing: Add coverage for simulation classes --- .../test_simulation_generator.py | 10 +- diffsims/tests/simulations/__init__.py | 17 + .../tests/simulations/test_simulations1d.py | 69 +++ .../tests/simulations/test_simulations2d.py | 456 ++++++++++++++++++ 4 files changed, 551 insertions(+), 1 deletion(-) rename diffsims/tests/{simulations => generators}/test_simulation_generator.py (95%) create mode 100644 diffsims/tests/simulations/__init__.py create mode 100644 diffsims/tests/simulations/test_simulations1d.py create mode 100644 diffsims/tests/simulations/test_simulations2d.py diff --git a/diffsims/tests/simulations/test_simulation_generator.py b/diffsims/tests/generators/test_simulation_generator.py similarity index 95% rename from diffsims/tests/simulations/test_simulation_generator.py rename to diffsims/tests/generators/test_simulation_generator.py index bfe88654..2295a5e7 100644 --- a/diffsims/tests/simulations/test_simulation_generator.py +++ b/diffsims/tests/generators/test_simulation_generator.py @@ -177,7 +177,7 @@ def test_appropriate_intensities(self, diffraction_calculator, local_structure): """Tests the central beam is strongest.""" diffraction = diffraction_calculator.calculate_diffraction2d( local_structure, reciprocal_radius=0.5, with_direct_beam=True - ) # direct beam doesn't work + ) indices = [tuple(np.round(i).astype(int)) for i in diffraction.coordinates.hkl] central_beam = indices.index((0, 0, 0)) @@ -187,6 +187,14 @@ def test_appropriate_intensities(self, diffraction_calculator, local_structure): ) assert np.all(smaller) + def test_direct_beam(self, diffraction_calculator, local_structure): + diffraction = diffraction_calculator.calculate_diffraction2d( + local_structure, reciprocal_radius=0.5, with_direct_beam=False + ) + indices = [tuple(np.round(i).astype(int)) for i in diffraction.coordinates.hkl] + with pytest.raises(ValueError): + indices.index((0, 0, 0)) + def test_shape_factor_strings(self, diffraction_calculator, local_structure): _ = diffraction_calculator.calculate_diffraction2d( local_structure, diff --git a/diffsims/tests/simulations/__init__.py b/diffsims/tests/simulations/__init__.py new file mode 100644 index 00000000..8ff2bc07 --- /dev/null +++ b/diffsims/tests/simulations/__init__.py @@ -0,0 +1,17 @@ +# -*- coding: utf-8 -*- +# Copyright 2017-2024 The diffsims developers +# +# This file is part of diffsims. +# +# diffsims is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# diffsims is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with diffsims. If not, see . diff --git a/diffsims/tests/simulations/test_simulations1d.py b/diffsims/tests/simulations/test_simulations1d.py new file mode 100644 index 00000000..4900e664 --- /dev/null +++ b/diffsims/tests/simulations/test_simulations1d.py @@ -0,0 +1,69 @@ +# -*- coding: utf-8 -*- +# Copyright 2017-2024 The diffsims developers +# +# This file is part of diffsims. +# +# diffsims is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# diffsims is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with diffsims. If not, see . +import matplotlib.pyplot as plt +import pytest + +from orix.crystal_map import Phase +import numpy as np + +from diffsims.tests.generators.test_simulation_generator import make_phase +from diffsims.simulations import Simulation1D + + +class TestSingleSimulation: + @pytest.fixture + def simulation1d(self): + al_phase = make_phase() + al_phase.name = "Al" + hkls = np.array(["100", "110", "111"]) + magnitudes = np.array([1, 2, 3]) + inten = np.array([1, 2, 3]) + recip = 4.0 + + return Simulation1D( + phase=al_phase, + hkl=hkls, + reciprocal_spacing=magnitudes, + intensities=inten, + reciprocal_radius=recip, + wavelength=0.025, + ) + + def test_init(self, simulation1d): + assert isinstance(simulation1d, Simulation1D) + assert isinstance(simulation1d.phase, Phase) + assert isinstance(simulation1d.hkl, np.ndarray) + assert isinstance(simulation1d.reciprocal_spacing, np.ndarray) + assert isinstance(simulation1d.intensities, np.ndarray) + assert isinstance(simulation1d.reciprocal_radius, float) + + @pytest.mark.parametrize("annotate", [True, False]) + @pytest.mark.parametrize("ax", [None, "new"]) + @pytest.mark.parametrize("with_labels", [True, False]) + def test_plot(self, simulation1d, annotate, ax, with_labels): + if ax == "new": + fig, ax = plt.subplots() + fig = simulation1d.plot(annotate_peaks=annotate, ax=ax, with_labels=with_labels) + + def test_repr(self, simulation1d): + assert simulation1d.__repr__() == "Simulation1D(name: Al, wavelength: 0.025)" + + def test_theta(self, simulation1d): + np.testing.assert_almost_equal( + simulation1d.theta, np.array([0.02499479, 0.0499584, 0.07485985]) + ) diff --git a/diffsims/tests/simulations/test_simulations2d.py b/diffsims/tests/simulations/test_simulations2d.py new file mode 100644 index 00000000..15fa6790 --- /dev/null +++ b/diffsims/tests/simulations/test_simulations2d.py @@ -0,0 +1,456 @@ +# -*- coding: utf-8 -*- +# Copyright 2017-2024 The diffsims developers +# +# This file is part of diffsims. +# +# diffsims is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# diffsims is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with diffsims. If not, see . + +import numpy as np +import pytest + +from diffpy.structure import Structure, Atom, Lattice +from orix.crystal_map import Phase +from orix.quaternion import Rotation + +from diffsims.simulations import Simulation2D +from diffsims.generators.simulation_generator import SimulationGenerator +from diffsims.crystallography.reciprocal_lattice_vector import ReciprocalLatticeVector + + +@pytest.fixture(scope="module") +def al_phase(): + p = Phase( + name="al", + space_group=225, + structure=Structure( + atoms=[Atom("al", [0, 0, 0])], + lattice=Lattice(0.405, 0.405, 0.405, 90, 90, 90), + ), + ) + return p + + +class TestSingleSimulation: + @pytest.fixture + def single_simulation(self, al_phase): + gen = SimulationGenerator(accelerating_voltage=200) + rot = Rotation.from_axes_angles([1, 0, 0], 45, degrees=True) + coords = ReciprocalLatticeVector(phase=al_phase, xyz=[[1, 0, 0]]) + sim = Simulation2D( + phases=al_phase, simulation_generator=gen, coordinates=coords, rotations=rot + ) + return sim + + def test_init(self, single_simulation): + assert isinstance(single_simulation, Simulation2D) + assert isinstance(single_simulation.phases, Phase) + assert isinstance(single_simulation.simulation_generator, SimulationGenerator) + assert isinstance(single_simulation.rotations, Rotation) + + def test_get_simulation(self, single_simulation): + rotation, phase, coords = single_simulation.get_simulation(0) + assert isinstance(rotation, Rotation) + assert phase == 0 + + def test_iphase(self, single_simulation): + with pytest.raises(ValueError): + single_simulation.iphase[0] + + def test_irot(self, single_simulation): + with pytest.raises(ValueError): + single_simulation.irot[0] + + def test_iter(self, single_simulation): + count = 0 + for sim in single_simulation: + count += 1 + assert isinstance(sim, ReciprocalLatticeVector) + assert count == 1 + + def test_plot(self, single_simulation): + single_simulation.plot() + + def test_num_rotations(self, single_simulation): + assert single_simulation._num_rotations() == 1 + + def test_polar_flatten(self, single_simulation): + ( + r_templates, + theta_templates, + intensities_templates, + ) = single_simulation.polar_flatten_simulations() + assert r_templates.shape == (1, 1) + assert theta_templates.shape == (1, 1) + assert intensities_templates.shape == (1, 1) + + def test_polar_flatten_axes(self, single_simulation): + radial_axes = np.linspace(0, 1, 10) + theta_axes = np.linspace(0, 2 * np.pi, 10) + ( + r_templates, + theta_templates, + intensities_templates, + ) = single_simulation.polar_flatten_simulations( + radial_axes=radial_axes, azimuthal_axes=theta_axes + ) + assert r_templates.shape == (1, 1) + assert theta_templates.shape == (1, 1) + assert intensities_templates.shape == (1, 1) + + def test_deepcopy(self, single_simulation): + copied = single_simulation.deepcopy() + assert copied is not single_simulation + + +class TestSimulationInitFailures: + def test_different_size(self, al_phase): + gen = SimulationGenerator(accelerating_voltage=200) + rot = Rotation.from_axes_angles([1, 0, 0], 45, degrees=True) + coords = ReciprocalLatticeVector(phase=al_phase, xyz=[[1, 0, 0], [1, 1, 1]]) + with pytest.raises(ValueError): + sim = Simulation2D( + phases=al_phase, + simulation_generator=gen, + coordinates=[coords, coords], + rotations=rot, + ) + + def test_different_size2(self, al_phase): + gen = SimulationGenerator(accelerating_voltage=200) + rot = Rotation.from_axes_angles([1, 0, 0], (0, 45), degrees=True) + coords = ReciprocalLatticeVector(phase=al_phase, xyz=[[1, 0, 0], [1, 1, 1]]) + with pytest.raises(ValueError): + sim = Simulation2D( + phases=al_phase, + simulation_generator=gen, + coordinates=[coords, coords, coords], + rotations=rot, + ) + + def test_different_size_multiphase(self, al_phase): + gen = SimulationGenerator(accelerating_voltage=200) + rot = Rotation.from_axes_angles([1, 0, 0], 45, degrees=True) + coords = ReciprocalLatticeVector(phase=al_phase, xyz=[[1, 0, 0], [1, 1, 1]]) + with pytest.raises(ValueError): + sim = Simulation2D( + phases=[al_phase, al_phase], + simulation_generator=gen, + coordinates=[[coords, coords], [coords, coords]], + rotations=[rot, rot], + ) + + def test_different_num_phase(self, al_phase): + gen = SimulationGenerator(accelerating_voltage=200) + rot = Rotation.from_axes_angles([1, 0, 0], 45, degrees=True) + coords = ReciprocalLatticeVector(phase=al_phase, xyz=[[1, 0, 0], [1, 1, 1]]) + with pytest.raises(ValueError): + sim = Simulation2D( + phases=[al_phase, al_phase], + simulation_generator=gen, + coordinates=[[coords, coords], [coords, coords], [coords, coords]], + rotations=[rot, rot], + ) + + def test_different_num_phase_and_rot(self, al_phase): + gen = SimulationGenerator(accelerating_voltage=200) + rot = Rotation.from_axes_angles([1, 0, 0], 45, degrees=True) + coords = ReciprocalLatticeVector(phase=al_phase, xyz=[[1, 0, 0], [1, 1, 1]]) + with pytest.raises(ValueError): + sim = Simulation2D( + phases=[al_phase, al_phase], + simulation_generator=gen, + coordinates=[[coords, coords], [coords, coords], [coords, coords]], + rotations=[rot, rot, rot], + ) + + +class TestSinglePhaseMultiSimulation: + @pytest.fixture + def al_phase(self): + p = Phase( + name="al", + space_group=225, + structure=Structure( + atoms=[Atom("al", [0, 0, 0])], + lattice=Lattice(0.405, 0.405, 0.405, 90, 90, 90), + ), + ) + return p + + @pytest.fixture + def multi_simulation(self, al_phase): + gen = SimulationGenerator(accelerating_voltage=200) + rot = Rotation.from_axes_angles([1, 0, 0], (0, 15, 30, 45), degrees=True) + coords = ReciprocalLatticeVector( + phase=al_phase, xyz=[[1, 0, 0], [0, 1, 0], [1, 1, 0], [1, 1, 1]] + ) + + vectors = [coords, coords, coords, coords] + + sim = Simulation2D( + phases=al_phase, + simulation_generator=gen, + coordinates=vectors, + rotations=rot, + ) + return sim + + def test_get_simulation(self, multi_simulation): + for i in range(4): + rotation, phase, coords = multi_simulation.get_simulation(i) + assert isinstance(rotation, Rotation) + assert phase == 0 + + def test_get_current_rotation(self, multi_simulation): + rot = multi_simulation.get_current_rotation() + np.testing.assert_array_equal(rot, multi_simulation.rotations[0].to_matrix()[0]) + + def test_init(self, multi_simulation): + assert isinstance(multi_simulation, Simulation2D) + assert isinstance(multi_simulation.phases, Phase) + assert isinstance(multi_simulation.simulation_generator, SimulationGenerator) + assert isinstance(multi_simulation.rotations, Rotation) + assert isinstance(multi_simulation.coordinates, np.ndarray) + + def test_iphase(self, multi_simulation): + with pytest.raises(ValueError): + multi_simulation.iphase[0] + + def test_irot(self, multi_simulation): + sliced_sim = multi_simulation.irot[0] + assert isinstance(sliced_sim, Simulation2D) + assert isinstance(sliced_sim.phases, Phase) + assert sliced_sim.rotations.size == 1 + assert sliced_sim.coordinates.size == 4 + + def test_irot_slice(self, multi_simulation): + sliced_sim = multi_simulation.irot[0:2] + assert isinstance(sliced_sim, Simulation2D) + assert isinstance(sliced_sim.phases, Phase) + assert sliced_sim.rotations.size == 2 + assert sliced_sim.coordinates.size == 2 + + def test_plot(self, multi_simulation): + multi_simulation.plot() + + def test_plot_rotation(self, multi_simulation): + multi_simulation.plot_rotations() + + def test_iter(self, multi_simulation): + multi_simulation.phase_index = 0 + multi_simulation.rotation_index = 0 + count = 0 + for sim in multi_simulation: + count += 1 + assert isinstance(sim, ReciprocalLatticeVector) + assert count == 4 + + def test_polar_flatten(self, multi_simulation): + ( + r_templates, + theta_templates, + intensities_templates, + ) = multi_simulation.polar_flatten_simulations() + assert r_templates.shape == (4, 4) + assert theta_templates.shape == (4, 4) + assert intensities_templates.shape == (4, 4) + + +class TestMultiPhaseMultiSimulation: + @pytest.fixture + def al_phase(self): + p = Phase( + name="al", + space_group=225, + structure=Structure( + atoms=[Atom("al", [0, 0, 0])], + lattice=Lattice(0.405, 0.405, 0.405, 90, 90, 90), + ), + ) + return p + + @pytest.fixture + def multi_simulation(self, al_phase): + gen = SimulationGenerator(accelerating_voltage=200) + rot = Rotation.from_axes_angles([1, 0, 0], (0, 15, 30, 45), degrees=True) + rot2 = rot + coords = ReciprocalLatticeVector( + phase=al_phase, + xyz=[ + [1, 0, 0], + [0, -0.3, 0], + [1 / 0.405, 1 / -0.405, 0], + [0.1, -0.1, -0.3], + ], + ) + coords.intensity = 1 + vectors = [coords, coords, coords, coords] + al_phase2 = al_phase.deepcopy() + al_phase2.name = "al2" + sim = Simulation2D( + phases=[al_phase, al_phase2], + simulation_generator=gen, + coordinates=[vectors, vectors], + rotations=[rot, rot2], + ) + return sim + + def test_init(self, multi_simulation): + assert isinstance(multi_simulation, Simulation2D) + assert isinstance(multi_simulation.phases, np.ndarray) + assert isinstance(multi_simulation.simulation_generator, SimulationGenerator) + assert isinstance(multi_simulation.rotations, np.ndarray) + assert isinstance(multi_simulation.coordinates, np.ndarray) + + def test_get_simulation(self, multi_simulation): + for i in range(4): + rotation, phase, coords = multi_simulation.get_simulation(i) + assert isinstance(rotation, Rotation) + assert phase == 0 + for i in range(4, 8): + rotation, phase, coords = multi_simulation.get_simulation(i) + assert isinstance(rotation, Rotation) + assert phase == 1 + + def test_iphase(self, multi_simulation): + phase_slic = multi_simulation.iphase[0] + assert isinstance(phase_slic, Simulation2D) + assert isinstance(phase_slic.phases, Phase) + assert phase_slic.rotations.size == 4 + + def test_iphase_str(self, multi_simulation): + phase_slic = multi_simulation.iphase["al"] + assert isinstance(phase_slic, Simulation2D) + assert isinstance(phase_slic.phases, Phase) + assert phase_slic.rotations.size == 4 + assert phase_slic.phases.name == "al" + + def test_iphase_error(self, multi_simulation): + with pytest.raises(ValueError): + phase_slic = multi_simulation.iphase[3.1] + + def test_irot(self, multi_simulation): + sliced_sim = multi_simulation.irot[0] + assert isinstance(sliced_sim, Simulation2D) + assert isinstance(sliced_sim.phases, np.ndarray) + assert sliced_sim.rotations.size == 2 + + def test_irot_slice(self, multi_simulation): + sliced_sim = multi_simulation.irot[0:2] + assert isinstance(sliced_sim, Simulation2D) + assert isinstance(sliced_sim.phases, np.ndarray) + assert sliced_sim.rotations.size == 2 + + @pytest.mark.parametrize("show_labels", [True, False]) + @pytest.mark.parametrize("units", ["real", "pixel"]) + @pytest.mark.parametrize("include_zero_beam", [True, False]) + def test_plot(self, multi_simulation, show_labels, units, include_zero_beam): + multi_simulation.phase_index = 0 + multi_simulation.rotation_index = 0 + multi_simulation.reciporical_radius = 2 + multi_simulation.coordinates[0][0].intensity = np.nan + multi_simulation.plot( + show_labels=show_labels, + units=units, + min_label_intensity=0.0, + include_direct_beam=include_zero_beam, + calibration=0.1, + ) + + def test_plot_rotation(self, multi_simulation): + multi_simulation.plot_rotations() + + def test_iter(self, multi_simulation): + multi_simulation.phase_index = 0 + multi_simulation.rotation_index = 0 + count = 0 + for sim in multi_simulation: + count += 1 + assert isinstance(sim, ReciprocalLatticeVector) + assert count == 8 + + def test_get_diffraction_pattern(self, multi_simulation): + # No diffraction spots in this pattern + pat = multi_simulation.get_diffraction_pattern( + shape=(50, 50), calibration=0.001 + ) + assert pat.shape == (50, 50) + assert np.max(pat.data) == 0 + + def test_get_diffraction_pattern2(self, multi_simulation): + pat = multi_simulation.get_diffraction_pattern( + shape=(512, 512), calibration=0.01 + ) + assert pat.shape == (512, 512) + assert np.max(pat.data) == 1 + + def test_polar_flatten(self, multi_simulation): + ( + r_templates, + theta_templates, + intensities_templates, + ) = multi_simulation.polar_flatten_simulations() + assert r_templates.shape == (8, 4) + assert theta_templates.shape == (8, 4) + assert intensities_templates.shape == (8, 4) + + def test_rotate_shift_coords(self, multi_simulation): + rot = multi_simulation.rotate_shift_coordinates(angle=0.1) + assert isinstance(rot, ReciprocalLatticeVector) + + +class TestMultiPhaseSingleSimulation: + @pytest.fixture + def al_phase(self): + p = Phase( + name="al", + space_group=225, + structure=Structure( + atoms=[Atom("al", [0, 0, 0])], + lattice=Lattice(0.405, 0.405, 0.405, 90, 90, 90), + ), + ) + return p + + @pytest.fixture + def multi_simulation(self, al_phase): + gen = SimulationGenerator(accelerating_voltage=200) + rot = Rotation.from_axes_angles([1, 0, 0], (0,), degrees=True) + rot2 = rot + coords = ReciprocalLatticeVector( + phase=al_phase, + xyz=[ + [1, 0, 0], + [0, -0.3, 0], + [1 / 0.405, 1 / -0.405, 0], + [0.1, -0.1, -0.3], + ], + ) + coords.intensity = 1 + vectors = coords + al_phase2 = al_phase.deepcopy() + al_phase2.name = "al2" + sim = Simulation2D( + phases=[al_phase, al_phase2], + simulation_generator=gen, + coordinates=[vectors, vectors], + rotations=[rot, rot2], + ) + return sim + + def test_get_simulation(self, multi_simulation): + for i in range(2): + rotation, phase, coords = multi_simulation.get_simulation(i) + assert isinstance(rotation, Rotation) + assert phase == i From 66f01330a7c9aa2ffe8fa8977850698d5d3a3615 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 17 Apr 2024 12:56:45 -0500 Subject: [PATCH 12/36] Documentation: Fixed Broken documentation --- diffsims/generators/simulation_generator.py | 4 ++++ diffsims/simulations/__init__.py | 2 ++ diffsims/simulations/simulation2d.py | 5 +++++ doc/reference/index.rst | 1 + 4 files changed, 12 insertions(+) diff --git a/diffsims/generators/simulation_generator.py b/diffsims/generators/simulation_generator.py index a43442a8..81aa532c 100644 --- a/diffsims/generators/simulation_generator.py +++ b/diffsims/generators/simulation_generator.py @@ -16,6 +16,8 @@ # You should have received a copy of the GNU General Public License # along with diffsims. If not, see . +"""Kinematic Diffraction Simulation Generator.""" + from typing import Union, Sequence import numpy as np @@ -51,6 +53,8 @@ from diffsims.simulations import Simulation1D, Simulation2D +__all__ = ["SimulationGenerator"] + class SimulationGenerator: """ diff --git a/diffsims/simulations/__init__.py b/diffsims/simulations/__init__.py index 2612e7d1..c3909008 100644 --- a/diffsims/simulations/__init__.py +++ b/diffsims/simulations/__init__.py @@ -16,6 +16,8 @@ # You should have received a copy of the GNU General Public License # along with diffsims. If not, see . +"""Kinematic Diffraction Simulation Results.""" + from diffsims.simulations.simulation2d import Simulation2D from diffsims.simulations.simulation1d import Simulation1D diff --git a/diffsims/simulations/simulation2d.py b/diffsims/simulations/simulation2d.py index ed6bcdf0..ade9dae3 100644 --- a/diffsims/simulations/simulation2d.py +++ b/diffsims/simulations/simulation2d.py @@ -33,6 +33,11 @@ if TYPE_CHECKING: # pragma: no cover from diffsims.generators.simulation_generator import SimulationGenerator +__all__ = [ + "Simulation2D", + "get_closest", +] + class PhaseGetter: """A class for getting the phases of a simulation library. diff --git a/doc/reference/index.rst b/doc/reference/index.rst index 91fbe73e..6bfe0428 100644 --- a/doc/reference/index.rst +++ b/doc/reference/index.rst @@ -29,5 +29,6 @@ the `demos `_. libraries pattern sims + simulations structure_factor utils From 3d9c662f315124a1b9493873ff219d195b83fc0b Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Sun, 21 Apr 2024 09:22:17 -0500 Subject: [PATCH 13/36] Refactor: Fixes per @viljarjf --- .github/workflows/build.yml | 2 +- diffsims/generators/simulation_generator.py | 20 ++++++++++++------- diffsims/simulations/__init__.py | 2 +- diffsims/simulations/simulation2d.py | 14 ++++++------- .../generators/test_simulation_generator.py | 12 ++++++++++- .../tests/simulations/test_simulations2d.py | 2 +- setup.py | 2 +- 7 files changed, 35 insertions(+), 19 deletions(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 86dd70a4..cfdd841c 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -46,7 +46,7 @@ jobs: - os: ubuntu-latest python-version: 3.8 OLDEST_SUPPORTED_VERSION: true - DEPENDENCIES: diffpy.structure==3.0.2 matplotlib==3.5 numpy==1.17.3 orix==0.11.0 scipy==1.8 tqdm==4.9 + DEPENDENCIES: diffpy.structure==3.0.2 matplotlib==3.5 numpy==1.17.3 orix==0.12.1 scipy==1.8 tqdm==4.9 LABEL: -oldest steps: - uses: actions/checkout@v4 diff --git a/diffsims/generators/simulation_generator.py b/diffsims/generators/simulation_generator.py index 81aa532c..7279e93b 100644 --- a/diffsims/generators/simulation_generator.py +++ b/diffsims/generators/simulation_generator.py @@ -176,6 +176,12 @@ def calculate_diffraction2d( phase = [phase] if isinstance(rotation, Rotation): rotation = [rotation] + if len(phase) != len(rotation): + raise ValueError( + "The number of phases and rotations must be equal. " + f"Got {len(phase)} phases and {len(rotation)} rotations." + ) + if debye_waller_factors is None: debye_waller_factors = {} # Specify variables used in calculation @@ -298,22 +304,22 @@ def calculate_diffraction1d( hkls_labels = ["".join([str(int(x)) for x in xs]) for xs in g_indices] - peaks = {} + peaks = [] for l, i, g in zip(hkls_labels, i_hkl, g_hkls): - peaks[l] = [i, g] + peaks.append((l, [i, g])) # Scale intensities so that the max intensity is 100. - max_intensity = max([v[0] for v in peaks.values()]) + max_intensity = max([v[1][0] for v in peaks]) reciporical_spacing = [] intensities = [] hkls = [] - for k in peaks.keys(): - v = peaks[k] - if v[0] / max_intensity * 100 > minimum_intensity and (k != "000"): + for p in peaks: + label, v = p # (label, [intensity,g]) + if v[0] / max_intensity * 100 > minimum_intensity and (label != "000"): reciporical_spacing.append(v[1]) intensities.append(v[0]) - hkls.append(k) + hkls.append(label) intensities = np.asarray(intensities) / max(intensities) * 100 diff --git a/diffsims/simulations/__init__.py b/diffsims/simulations/__init__.py index c3909008..91ea1d10 100644 --- a/diffsims/simulations/__init__.py +++ b/diffsims/simulations/__init__.py @@ -1,5 +1,5 @@ # -*- coding: utf-8 -*- -# Copyright 2017-2023 The diffsims developers +# Copyright 2017-2024 The diffsims developers # # This file is part of diffsims. # diff --git a/diffsims/simulations/simulation2d.py b/diffsims/simulations/simulation2d.py index ade9dae3..f943a9d1 100644 --- a/diffsims/simulations/simulation2d.py +++ b/diffsims/simulations/simulation2d.py @@ -44,8 +44,8 @@ class PhaseGetter: Parameters ---------- - phases : Sequence[Phase] - The phases in the library. + simulation : Simulation2D + The simulation to get from. """ def __init__(self, simulation): @@ -77,8 +77,8 @@ class RotationGetter: Parameters ---------- - phases : Sequence[Phase] - The phases in the library. + simulation : Simulation2D + The simulation to get from. """ def __init__(self, simulation): @@ -435,8 +435,8 @@ def get_current_coordinates(self): else: return copy.deepcopy(self.coordinates) - def get_current_rotation(self): - """Returns the matrix for the current matrix""" + def get_current_rotation_matrix(self): + """Returns the current rotation matrix based on the phase and rotation index""" if self.has_multiple_phases: return copy.deepcopy( self.rotations[self.phase_index].to_matrix()[self.rotation_index] @@ -554,7 +554,7 @@ def plot( if show_labels: millers = np.round( np.matmul( - np.matmul(coords.data, self.get_current_rotation().T), + np.matmul(coords.data, self.get_current_rotation_matrix().T), coords.phase.structure.lattice.base.T, ) ).astype(np.int16) diff --git a/diffsims/tests/generators/test_simulation_generator.py b/diffsims/tests/generators/test_simulation_generator.py index 2295a5e7..b9fc2efa 100644 --- a/diffsims/tests/generators/test_simulation_generator.py +++ b/diffsims/tests/generators/test_simulation_generator.py @@ -1,5 +1,5 @@ # -*- coding: utf-8 -*- -# Copyright 2017-2023 The diffsims developers +# Copyright 2017-2024 The diffsims developers # # This file is part of diffsims. # @@ -246,6 +246,16 @@ def test_multiphase_multirotation_simulation(): ) +def test_multiphase_multirotation_simulation_error(): + generator = SimulationGenerator(300) + silicon = make_phase(5) + big_silicon = make_phase(10) + rot = Rotation.from_euler([[0, 0, 0], [0.1, 0.1, 0.1]]) + rot2 = Rotation.from_euler([[0, 0, 0], [0.1, 0.1, 0.1], [0.2, 0.2, 0.2]]) + with pytest.raises(ValueError): + sim = generator.calculate_diffraction2d([silicon, big_silicon], rotation=[rot]) + + @pytest.mark.parametrize("scattering_param", ["lobato", "xtables"]) def test_param_check(scattering_param): generator = SimulationGenerator(300, scattering_params=scattering_param) diff --git a/diffsims/tests/simulations/test_simulations2d.py b/diffsims/tests/simulations/test_simulations2d.py index 15fa6790..f04dca46 100644 --- a/diffsims/tests/simulations/test_simulations2d.py +++ b/diffsims/tests/simulations/test_simulations2d.py @@ -213,7 +213,7 @@ def test_get_simulation(self, multi_simulation): assert phase == 0 def test_get_current_rotation(self, multi_simulation): - rot = multi_simulation.get_current_rotation() + rot = multi_simulation.get_current_rotation_matrix() np.testing.assert_array_equal(rot, multi_simulation.rotations[0].to_matrix()[0]) def test_init(self, multi_simulation): diff --git a/setup.py b/setup.py index 0a35c2d1..1ba7c559 100644 --- a/setup.py +++ b/setup.py @@ -78,7 +78,7 @@ "matplotlib >= 3.3", "numba", "numpy >= 1.17.3", - "orix >= 0.11", + "orix >= 0.12.1", "psutil", "scipy >= 1.8", "tqdm >= 4.9", From 99327ade6f8857a4066bee9ac932472469294a40 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Thu, 25 Apr 2024 14:32:31 -0500 Subject: [PATCH 14/36] Refactor: Clean up intensity in rlv class --- diffsims/crystallography/__init__.py | 2 + .../crystallography/diffracting_vector.py | 111 ++++++++++++++++++ .../reciprocal_lattice_vector.py | 30 ----- diffsims/generators/simulation_generator.py | 4 +- diffsims/simulations/simulation2d.py | 1 + .../test_diffracting_vector.py | 37 ++++++ .../test_reciprocal_lattice_vector.py | 11 -- .../tests/simulations/test_simulations2d.py | 28 ++--- 8 files changed, 167 insertions(+), 57 deletions(-) create mode 100644 diffsims/crystallography/diffracting_vector.py create mode 100644 diffsims/tests/crystallography/test_diffracting_vector.py diff --git a/diffsims/crystallography/__init__.py b/diffsims/crystallography/__init__.py index 736a8e3d..ff42077a 100644 --- a/diffsims/crystallography/__init__.py +++ b/diffsims/crystallography/__init__.py @@ -27,6 +27,7 @@ get_hkl, ) from diffsims.crystallography.reciprocal_lattice_vector import ReciprocalLatticeVector +from diffsims.crystallography.diffracting_vector import DiffractingVector __all__ = [ "get_equivalent_hkl", @@ -34,4 +35,5 @@ "get_hkl", "ReciprocalLatticePoint", "ReciprocalLatticeVector", + "DiffractingVector", ] diff --git a/diffsims/crystallography/diffracting_vector.py b/diffsims/crystallography/diffracting_vector.py new file mode 100644 index 00000000..c2087b35 --- /dev/null +++ b/diffsims/crystallography/diffracting_vector.py @@ -0,0 +1,111 @@ +# -*- coding: utf-8 -*- +# Copyright 2017-2024 The diffsims developers +# +# This file is part of diffsims. +# +# diffsims is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# diffsims is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with diffsims. If not, see . + +from diffsims.crystallography import ReciprocalLatticeVector +import numpy as np + + +class DiffractingVector(ReciprocalLatticeVector): + r"""Reciprocal lattice vectors :math:`(hkl)` for use in electron + diffraction analysis and simulation. + + All lengths are assumed to be given in Å or inverse Å. + + This extends the :class:`ReciprocalLatticeVector` class. Diffracting Vectors + focus on the subset of reciporical lattice vectors that are relevant for + electron diffraction based on the intersection of the Ewald sphere with the + reciprocal lattice. + + Parameters + ---------- + phase : orix.crystal_map.Phase + A phase with a crystal lattice and symmetry. + xyz : numpy.ndarray, list, or tuple, optional + Cartesian coordinates of indices of reciprocal lattice vector(s) + ``hkl``. Default is ``None``. This, ``hkl``, or ``hkil`` is + required. + hkl : numpy.ndarray, list, or tuple, optional + Indices of reciprocal lattice vector(s). Default is ``None``. + This, ``xyz``, or ``hkil`` is required. + hkil : numpy.ndarray, list, or tuple, optional + Indices of reciprocal lattice vector(s), often preferred over + ``hkl`` in trigonal and hexagonal lattices. Default is ``None``. + This, ``xyz``, or ``hkl`` is required. + + Examples + -------- + >>> from diffpy.structure import Atom, Lattice, Structure + >>> from orix.crystal_map import Phase + >>> from diffsims.crystallography import DiffractingVector + >>> phase = Phase( + ... "al", + ... space_group=225, + ... structure=Structure( + ... lattice=Lattice(4.04, 4.04, 4.04, 90, 90, 90), + ... atoms=[Atom("Al", [0, 0, 1])], + ... ), + ... ) + >>> rlv = DiffractingVector(phase, hkl=[[1, 1, 1], [2, 0, 0]]) + >>> rlv + ReciprocalLatticeVector (2,), al (m-3m) + [[1. 1. 1.] + [2. 0. 0.]] + + """ + + def __init__(self, phase, xyz=None, hkl=None, hkil=None, intensity=None): + super().__init__(phase, xyz=xyz, hkl=hkl, hkil=hkil) + if intensity is None: + self._intensity = np.full(self.shape, np.nan) + elif len(intensity) != self.size: + raise ValueError("Length of intensity array must match number of vectors") + else: + self._intensity = np.array(intensity) + + def __getitem__(self, key): + dv_new = super().__getitem__(key) + if np.isnan(self.intensity).all(): + dv_new._intensity = np.full(dv_new.shape, np.nan) + else: + slic = self.intensity[key] + if not hasattr(slic, "__len__"): + slic = np.array( + [ + slic, + ] + ) + dv_new._intensity = slic + + return dv_new + + @property + def intensity(self): + return self._intensity + + @intensity.setter + def intensity(self, value): + if not hasattr(value, "__len__"): + value = np.array( + [ + value, + ] + * self.size + ) + if len(value) != self.size: + raise ValueError("Length of intensity array must match number of vectors") + self._intensity = np.array(value) diff --git a/diffsims/crystallography/reciprocal_lattice_vector.py b/diffsims/crystallography/reciprocal_lattice_vector.py index e3de582a..59360244 100644 --- a/diffsims/crystallography/reciprocal_lattice_vector.py +++ b/diffsims/crystallography/reciprocal_lattice_vector.py @@ -123,7 +123,6 @@ def __init__(self, phase, xyz=None, hkl=None, hkil=None): self._theta = np.full(self.shape, np.nan) self._structure_factor = np.full(self.shape, np.nan, dtype="complex128") - self._intensity = np.full(self.shape, np.nan) def __getitem__(self, key): new_data = self.data[key] @@ -141,18 +140,6 @@ def __getitem__(self, key): else: rlv_new._theta = self.theta[key] - if np.isnan(self.intensity).all(): - rlv_new._intensity = np.full(rlv_new.shape, np.nan) - else: - slic = self.intensity[key] - if not hasattr(slic, "__len__"): - slic = np.array( - [ - slic, - ] - ) - rlv_new._intensity = slic - return rlv_new def __repr__(self): @@ -516,23 +503,6 @@ def scattering_parameter(self): return 0.5 * self.gspacing - @property - def intensity(self): - return self._intensity - - @intensity.setter - def intensity(self, value): - if not hasattr(value, "__len__"): - value = np.array( - [ - value, - ] - * self.size - ) - if len(value) != self.size: - raise ValueError("Length of intensity array must match number of vectors") - self._intensity = np.array(value) - def rotate_from_matrix(self, rotation_matrix): return ReciprocalLatticeVector( phase=self.phase, xyz=np.matmul(rotation_matrix.T, self.data.T).T diff --git a/diffsims/generators/simulation_generator.py b/diffsims/generators/simulation_generator.py index 7279e93b..498cda92 100644 --- a/diffsims/generators/simulation_generator.py +++ b/diffsims/generators/simulation_generator.py @@ -24,7 +24,7 @@ from orix.quaternion import Rotation from orix.crystal_map import Phase -from diffsims.crystallography import ReciprocalLatticeVector +from diffsims.crystallography import ReciprocalLatticeVector, DiffractingVector from diffsims.utils.shape_factor_models import ( linear, atanc, @@ -393,7 +393,7 @@ def get_intersecting_reflections( # select these reflections intersected_vectors = rotated_vectors[intersection] - intersected_vectors = ReciprocalLatticeVector( + intersected_vectors = DiffractingVector( phase=recip.phase, xyz=intersected_vectors ) excitation_error = excitation_error[intersection] diff --git a/diffsims/simulations/simulation2d.py b/diffsims/simulations/simulation2d.py index f943a9d1..7df1c953 100644 --- a/diffsims/simulations/simulation2d.py +++ b/diffsims/simulations/simulation2d.py @@ -114,6 +114,7 @@ def __init__( self, phases: Sequence[Phase], coordinates: Union[ + ReciprocalLatticeVector, Sequence[ReciprocalLatticeVector], Sequence[Sequence[ReciprocalLatticeVector]], ], diff --git a/diffsims/tests/crystallography/test_diffracting_vector.py b/diffsims/tests/crystallography/test_diffracting_vector.py new file mode 100644 index 00000000..adce93db --- /dev/null +++ b/diffsims/tests/crystallography/test_diffracting_vector.py @@ -0,0 +1,37 @@ +from diffsims.crystallography import DiffractingVector +import pytest +import numpy as np + + +class TestDiffractingVector: + def test_init(self, ferrite_phase): + rlv = DiffractingVector( + ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]], intensity=[1, 2] + ) + assert rlv.phase == ferrite_phase + assert rlv.shape == (2,) + assert rlv.hkl.shape == (2, 3) + assert np.allclose(rlv.hkl, [[1, 1, 1], [2, 0, 0]]) + assert np.allclose(rlv.intensity, [1, 2]) + + def test_init_wrong_intensity_length(self, ferrite_phase): + with pytest.raises(ValueError): + DiffractingVector(ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]], intensity=[1]) + + def test_add_intensity(self, ferrite_phase): + rlv = DiffractingVector.from_min_dspacing(ferrite_phase, 1.5) + rlv.intensity = 1 + assert isinstance(rlv.intensity, np.ndarray) + assert np.allclose(rlv.intensity, np.ones(rlv.size)) + + def test_add_intensity_error(self, ferrite_phase): + rlv = DiffractingVector.from_min_dspacing(ferrite_phase, 1.5) + with pytest.raises(ValueError): + rlv.intensity = [0, 1] + + def test_slicing(self, ferrite_phase): + rlv = DiffractingVector.from_min_dspacing(ferrite_phase, 1.5) + rlv.intensity = 1 + rlv_slice = rlv[0:3] + assert rlv_slice.size == 3 + assert np.allclose(rlv_slice.intensity, np.ones(3)) diff --git a/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py b/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py index 41f29fb3..975a382c 100644 --- a/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py +++ b/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py @@ -108,17 +108,6 @@ def test_repr(self, ferrite_phase): "[[ 1. 1. 0.]", ] - def test_add_intensity(self, ferrite_phase): - rlv = ReciprocalLatticeVector.from_min_dspacing(ferrite_phase, 1.5) - rlv.intensity = 1 - assert isinstance(rlv.intensity, np.ndarray) - assert np.allclose(rlv.intensity, np.ones(rlv.size)) - - def test_add_intensity_error(self, ferrite_phase): - rlv = ReciprocalLatticeVector.from_min_dspacing(ferrite_phase, 1.5) - with pytest.raises(ValueError): - rlv.intensity = [0, 1] - @pytest.mark.parametrize("degrees", [True, False]) def test_to_polar(self, ferrite_phase, degrees): rlv = ReciprocalLatticeVector.from_min_dspacing(ferrite_phase, 1.5) diff --git a/diffsims/tests/simulations/test_simulations2d.py b/diffsims/tests/simulations/test_simulations2d.py index f04dca46..f2757319 100644 --- a/diffsims/tests/simulations/test_simulations2d.py +++ b/diffsims/tests/simulations/test_simulations2d.py @@ -25,7 +25,7 @@ from diffsims.simulations import Simulation2D from diffsims.generators.simulation_generator import SimulationGenerator -from diffsims.crystallography.reciprocal_lattice_vector import ReciprocalLatticeVector +from diffsims.crystallography import DiffractingVector @pytest.fixture(scope="module") @@ -46,7 +46,7 @@ class TestSingleSimulation: def single_simulation(self, al_phase): gen = SimulationGenerator(accelerating_voltage=200) rot = Rotation.from_axes_angles([1, 0, 0], 45, degrees=True) - coords = ReciprocalLatticeVector(phase=al_phase, xyz=[[1, 0, 0]]) + coords = DiffractingVector(phase=al_phase, xyz=[[1, 0, 0]]) sim = Simulation2D( phases=al_phase, simulation_generator=gen, coordinates=coords, rotations=rot ) @@ -75,7 +75,7 @@ def test_iter(self, single_simulation): count = 0 for sim in single_simulation: count += 1 - assert isinstance(sim, ReciprocalLatticeVector) + assert isinstance(sim, DiffractingVector) assert count == 1 def test_plot(self, single_simulation): @@ -117,7 +117,7 @@ class TestSimulationInitFailures: def test_different_size(self, al_phase): gen = SimulationGenerator(accelerating_voltage=200) rot = Rotation.from_axes_angles([1, 0, 0], 45, degrees=True) - coords = ReciprocalLatticeVector(phase=al_phase, xyz=[[1, 0, 0], [1, 1, 1]]) + coords = DiffractingVector(phase=al_phase, xyz=[[1, 0, 0], [1, 1, 1]]) with pytest.raises(ValueError): sim = Simulation2D( phases=al_phase, @@ -129,7 +129,7 @@ def test_different_size(self, al_phase): def test_different_size2(self, al_phase): gen = SimulationGenerator(accelerating_voltage=200) rot = Rotation.from_axes_angles([1, 0, 0], (0, 45), degrees=True) - coords = ReciprocalLatticeVector(phase=al_phase, xyz=[[1, 0, 0], [1, 1, 1]]) + coords = DiffractingVector(phase=al_phase, xyz=[[1, 0, 0], [1, 1, 1]]) with pytest.raises(ValueError): sim = Simulation2D( phases=al_phase, @@ -141,7 +141,7 @@ def test_different_size2(self, al_phase): def test_different_size_multiphase(self, al_phase): gen = SimulationGenerator(accelerating_voltage=200) rot = Rotation.from_axes_angles([1, 0, 0], 45, degrees=True) - coords = ReciprocalLatticeVector(phase=al_phase, xyz=[[1, 0, 0], [1, 1, 1]]) + coords = DiffractingVector(phase=al_phase, xyz=[[1, 0, 0], [1, 1, 1]]) with pytest.raises(ValueError): sim = Simulation2D( phases=[al_phase, al_phase], @@ -153,7 +153,7 @@ def test_different_size_multiphase(self, al_phase): def test_different_num_phase(self, al_phase): gen = SimulationGenerator(accelerating_voltage=200) rot = Rotation.from_axes_angles([1, 0, 0], 45, degrees=True) - coords = ReciprocalLatticeVector(phase=al_phase, xyz=[[1, 0, 0], [1, 1, 1]]) + coords = DiffractingVector(phase=al_phase, xyz=[[1, 0, 0], [1, 1, 1]]) with pytest.raises(ValueError): sim = Simulation2D( phases=[al_phase, al_phase], @@ -165,7 +165,7 @@ def test_different_num_phase(self, al_phase): def test_different_num_phase_and_rot(self, al_phase): gen = SimulationGenerator(accelerating_voltage=200) rot = Rotation.from_axes_angles([1, 0, 0], 45, degrees=True) - coords = ReciprocalLatticeVector(phase=al_phase, xyz=[[1, 0, 0], [1, 1, 1]]) + coords = DiffractingVector(phase=al_phase, xyz=[[1, 0, 0], [1, 1, 1]]) with pytest.raises(ValueError): sim = Simulation2D( phases=[al_phase, al_phase], @@ -192,7 +192,7 @@ def al_phase(self): def multi_simulation(self, al_phase): gen = SimulationGenerator(accelerating_voltage=200) rot = Rotation.from_axes_angles([1, 0, 0], (0, 15, 30, 45), degrees=True) - coords = ReciprocalLatticeVector( + coords = DiffractingVector( phase=al_phase, xyz=[[1, 0, 0], [0, 1, 0], [1, 1, 0], [1, 1, 1]] ) @@ -253,7 +253,7 @@ def test_iter(self, multi_simulation): count = 0 for sim in multi_simulation: count += 1 - assert isinstance(sim, ReciprocalLatticeVector) + assert isinstance(sim, DiffractingVector) assert count == 4 def test_polar_flatten(self, multi_simulation): @@ -285,7 +285,7 @@ def multi_simulation(self, al_phase): gen = SimulationGenerator(accelerating_voltage=200) rot = Rotation.from_axes_angles([1, 0, 0], (0, 15, 30, 45), degrees=True) rot2 = rot - coords = ReciprocalLatticeVector( + coords = DiffractingVector( phase=al_phase, xyz=[ [1, 0, 0], @@ -377,7 +377,7 @@ def test_iter(self, multi_simulation): count = 0 for sim in multi_simulation: count += 1 - assert isinstance(sim, ReciprocalLatticeVector) + assert isinstance(sim, DiffractingVector) assert count == 8 def test_get_diffraction_pattern(self, multi_simulation): @@ -407,7 +407,7 @@ def test_polar_flatten(self, multi_simulation): def test_rotate_shift_coords(self, multi_simulation): rot = multi_simulation.rotate_shift_coordinates(angle=0.1) - assert isinstance(rot, ReciprocalLatticeVector) + assert isinstance(rot, DiffractingVector) class TestMultiPhaseSingleSimulation: @@ -428,7 +428,7 @@ def multi_simulation(self, al_phase): gen = SimulationGenerator(accelerating_voltage=200) rot = Rotation.from_axes_angles([1, 0, 0], (0,), degrees=True) rot2 = rot - coords = ReciprocalLatticeVector( + coords = DiffractingVector( phase=al_phase, xyz=[ [1, 0, 0], From 2ad456763db3dcc04b0dd08acfeadf835f9ac968 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Thu, 25 Apr 2024 14:45:33 -0500 Subject: [PATCH 15/36] Refactor: Remove calculate structure factor --- diffsims/crystallography/diffracting_vector.py | 5 +++++ diffsims/tests/crystallography/test_diffracting_vector.py | 5 +++++ 2 files changed, 10 insertions(+) diff --git a/diffsims/crystallography/diffracting_vector.py b/diffsims/crystallography/diffracting_vector.py index c2087b35..53ffb959 100644 --- a/diffsims/crystallography/diffracting_vector.py +++ b/diffsims/crystallography/diffracting_vector.py @@ -109,3 +109,8 @@ def intensity(self, value): if len(value) != self.size: raise ValueError("Length of intensity array must match number of vectors") self._intensity = np.array(value) + + def calculate_structure_factor(self): + raise NotImplementedError( + "Structure factor calculation not implemented for DiffractingVector" + ) diff --git a/diffsims/tests/crystallography/test_diffracting_vector.py b/diffsims/tests/crystallography/test_diffracting_vector.py index adce93db..adf2c4dd 100644 --- a/diffsims/tests/crystallography/test_diffracting_vector.py +++ b/diffsims/tests/crystallography/test_diffracting_vector.py @@ -35,3 +35,8 @@ def test_slicing(self, ferrite_phase): rlv_slice = rlv[0:3] assert rlv_slice.size == 3 assert np.allclose(rlv_slice.intensity, np.ones(3)) + + def test_structure_factor(self, ferrite_phase): + rlv = DiffractingVector.from_min_dspacing(ferrite_phase, 1.5) + with pytest.raises(NotImplementedError): + rlv.calculate_structure_factor() From 273ffccbd2657a935e251623c187dc1ff5c31288 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 3 May 2024 10:00:41 -0500 Subject: [PATCH 16/36] Refactor: Changes from @hakonanes --- .../crystallography/reciprocal_lattice_vector.py | 16 ++++++++++------ diffsims/generators/simulation_generator.py | 2 +- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/diffsims/crystallography/reciprocal_lattice_vector.py b/diffsims/crystallography/reciprocal_lattice_vector.py index 59360244..e377cdee 100644 --- a/diffsims/crystallography/reciprocal_lattice_vector.py +++ b/diffsims/crystallography/reciprocal_lattice_vector.py @@ -505,7 +505,7 @@ def scattering_parameter(self): def rotate_from_matrix(self, rotation_matrix): return ReciprocalLatticeVector( - phase=self.phase, xyz=np.matmul(rotation_matrix.T, self.data.T).T + phase=self.phase, xyz=np.matmul(self.data, rotation_matrix.T) ) @property @@ -1029,7 +1029,7 @@ def symmetrise(self, return_multiplicity=False, return_index=False): return new_out @classmethod - def from_highest_hkl(cls, phase, hkl): + def from_highest_hkl(cls, phase, hkl, include_zero_vector=False): """Create a set of unique reciprocal lattice vectors from three highest indices. @@ -1039,6 +1039,8 @@ def from_highest_hkl(cls, phase, hkl): A phase with a crystal lattice and symmetry. hkl : numpy.ndarray, list, or tuple Three highest reciprocal lattice vector indices. + include_zero_vector : bool + If ``True``, include the zero vector (000) in the set of vectors. Examples -------- @@ -1073,10 +1075,12 @@ def from_highest_hkl(cls, phase, hkl): """ idx = _get_indices_from_highest(highest_indices=hkl) + if include_zero_vector: + idx = np.vstack((idx, np.zeros(3, dtype=int))) return cls(phase, hkl=idx).unique() @classmethod - def from_min_dspacing(cls, phase, min_dspacing=0.7, include_zero_beam=False): + def from_min_dspacing(cls, phase, min_dspacing=0.7, include_zero_vector=False): """Create a set of unique reciprocal lattice vectors with a a direct space interplanar spacing greater than a lower threshold. @@ -1089,8 +1093,8 @@ def from_min_dspacing(cls, phase, min_dspacing=0.7, include_zero_beam=False): Smallest interplanar spacing to consider. Default is 0.7, in the unit used to define the lattice parameters in ``phase``, which is assumed to be Ångström. - include_zero_beam: bool - If ``True``, include the zero beam in the set of vectors. + include_zero_vector: bool + If ``True``, include the zero vector (000) in the set of vectors. Examples -------- @@ -1136,7 +1140,7 @@ def from_min_dspacing(cls, phase, min_dspacing=0.7, include_zero_beam=False): dspacing = 1 / phase.structure.lattice.rnorm(hkl) idx = dspacing >= min_dspacing hkl = hkl[idx] - if include_zero_beam: + if include_zero_vector: hkl = np.vstack((hkl, np.zeros(3, dtype=int))) return cls(phase, hkl=hkl).unique() diff --git a/diffsims/generators/simulation_generator.py b/diffsims/generators/simulation_generator.py index 498cda92..0d7c3e2a 100644 --- a/diffsims/generators/simulation_generator.py +++ b/diffsims/generators/simulation_generator.py @@ -193,7 +193,7 @@ def calculate_diffraction2d( recip = ReciprocalLatticeVector.from_min_dspacing( p, min_dspacing=1 / reciprocal_radius, - include_zero_beam=with_direct_beam, + include_zero_vector=with_direct_beam, ) phase_vectors = [] for rot in rotate.to_matrix(): From f32bfdaae5b8bdec92a2df49dec6210fd239e798 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 3 May 2024 10:04:37 -0500 Subject: [PATCH 17/36] Refactor: Improve Documentation --- diffsims/crystallography/diffracting_vector.py | 3 ++- diffsims/crystallography/reciprocal_lattice_vector.py | 10 ++++++++++ .../crystallography/test_reciprocal_lattice_vector.py | 8 -------- 3 files changed, 12 insertions(+), 9 deletions(-) diff --git a/diffsims/crystallography/diffracting_vector.py b/diffsims/crystallography/diffracting_vector.py index 53ffb959..41418af8 100644 --- a/diffsims/crystallography/diffracting_vector.py +++ b/diffsims/crystallography/diffracting_vector.py @@ -112,5 +112,6 @@ def intensity(self, value): def calculate_structure_factor(self): raise NotImplementedError( - "Structure factor calculation not implemented for DiffractingVector" + "Structure factor calculation not implemented for DiffractionVector. " + "Use ReciprocalLatticeVector instead." ) diff --git a/diffsims/crystallography/reciprocal_lattice_vector.py b/diffsims/crystallography/reciprocal_lattice_vector.py index e377cdee..ee937175 100644 --- a/diffsims/crystallography/reciprocal_lattice_vector.py +++ b/diffsims/crystallography/reciprocal_lattice_vector.py @@ -504,6 +504,16 @@ def scattering_parameter(self): return 0.5 * self.gspacing def rotate_from_matrix(self, rotation_matrix): + """Rotate the reciprocal lattice vectors using a (3x3) rotation matrix. + + This method creates a new instance of :class:`ReciprocalLatticeVector` + with the rotated vectors. + + Parameters + ---------- + rotation_matrix : numpy.ndarray + (3, 3) rotation matrix. + """ return ReciprocalLatticeVector( phase=self.phase, xyz=np.matmul(self.data, rotation_matrix.T) ) diff --git a/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py b/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py index 975a382c..d8717107 100644 --- a/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py +++ b/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py @@ -108,14 +108,6 @@ def test_repr(self, ferrite_phase): "[[ 1. 1. 0.]", ] - @pytest.mark.parametrize("degrees", [True, False]) - def test_to_polar(self, ferrite_phase, degrees): - rlv = ReciprocalLatticeVector.from_min_dspacing(ferrite_phase, 1.5) - r, theta, z = rlv.to_polar(degrees=degrees) - assert r.shape == (rlv.size,) - assert theta.shape == (rlv.size,) - assert z.shape == (rlv.size,) - def test_get_item(self, ferrite_phase): """Indexing gives desired vectors and properties carry over.""" rlv = ReciprocalLatticeVector.from_min_dspacing(ferrite_phase, 1.5) From a8c68a37c34af4e98106f48825041acdff498433 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 3 May 2024 10:20:05 -0500 Subject: [PATCH 18/36] Testing: Test simulation continuity --- .../reciprocal_lattice_vector.py | 2 +- diffsims/tests/generators/old_simulation.npy | Bin 0 -> 131200 bytes .../generators/test_simulation_generator.py | 65 ++++++++++++++++++ 3 files changed, 66 insertions(+), 1 deletion(-) create mode 100644 diffsims/tests/generators/old_simulation.npy diff --git a/diffsims/crystallography/reciprocal_lattice_vector.py b/diffsims/crystallography/reciprocal_lattice_vector.py index ee937175..d52f4e24 100644 --- a/diffsims/crystallography/reciprocal_lattice_vector.py +++ b/diffsims/crystallography/reciprocal_lattice_vector.py @@ -515,7 +515,7 @@ def rotate_from_matrix(self, rotation_matrix): (3, 3) rotation matrix. """ return ReciprocalLatticeVector( - phase=self.phase, xyz=np.matmul(self.data, rotation_matrix.T) + phase=self.phase, xyz=np.matmul(rotation_matrix.T, self.data.T).T ) @property diff --git a/diffsims/tests/generators/old_simulation.npy b/diffsims/tests/generators/old_simulation.npy new file mode 100644 index 0000000000000000000000000000000000000000..90008a11246364ee9551b03727e77bd87df0868d GIT binary patch literal 131200 zcmeFabyQW~7WQukcB5jUf`MXVu@=~YfeHq8i-BS{7Ir6g3yOgq*xfA_c6WZoe)-)! zKl?WZW1Ml$d<JJKo{_eO!R zxqNQtiXF>3yX1E67!cSeuuYFv0fFs1{D0@GwdvNoL+tt9o!j*65E~aLUZSi^?%4k; z=#qPA?tkllcJcoH@2>)X75J;bUj_ax1s=RQd3Su?sqO<^t`3aIdcnPkoqf6oKNG42 zWyX!4a;}KleeuNOW*cg$z9+UF=(WDN(*2TM?mVKc8H1BgJ$&%nIytYOE~8xvUU7f$ zd;7GaH!`YE0hfQ|e^^s(niGF_C$HYB_AAGr`qieZFS{${EjxX=>QH}a=7ImLF=KG@ zsfQ1KYxC(Rx>OQf(u*#oMVFSMOFGddp6EghPCoVU!EbFIfAPNx{7D4{iN8C}tP|?7 zI<2~zWJ?kM@G8prPWI!+F7;5s2`-hbId`H;QMuudvuhTqnZ=(r{NlL6jKRsL9zOVO z{Tm)i-@opEs6_FpKFW7p`RYzbC#Y-t5_Eg;X^SfU<3hpQ(HB*leE)>iFZoavb1hju z{P;671}C3-_~5rTA5YQ6MRb`Yy6lNIba4`0R#q}}@fKZ(!O5o{KKQN8<1hYKfj_AL z>otAr_WSszsp|%+8RF;m1&7Wqo9Cda>KU@WX_reX%&W+jM2Dly7@U0S;e+4Szk&6d zzU}b5TmD0f52zy2=i!f+OzT!OiMF4)dVzOT7X9?#)%Z`3&4OljKRsL9zOW3&8L*;vRBr7`uNK;Lzk#>hAuv$OKH)C z7@U0S;e+4WJpST;75I}1?2&ye*YT|9^mWwFM0F?JPNH4I)Utjv?6p_^Pw9GW&t=Bo zoLUb4UEI)0t3_w@C{?-Ne;nO{dc-|P2i%*9rE?A~a<>f_p*F*y0u!w0{$ z`LKU?mg{`h^Xt<0V?~#np@uGZ@)^1igOg7^eDGVF$6x&aEd`dY3?KjO%QE+-9o?24 zoAtmwW%HB;`=!sKrfvRoA%96v)pBs_8>dgTR3VoSw)CFVQ)S3eE@cb9zGe(gKK1ay zZ>wrGzjmure{*w{u}WBqsPIKunh_EnLyP zkiOE`#qoNo;${p^KK1ayZ|i(w|Y@6EnwLJM?<&{*?C)LWI-m8Pz8{`sW;5@>L;nBI0l6L^->ert_^r*ShUijDbV($-j1yg& zh%QM)moJG8U5LTSryf4|t<3|RK7HQF@J$Wz&Bw`xZ{C(Od=n(T`Buj84KXlL*iZRo0@y?&`Q`R%Wl^vJFg*SfwW*@%|qFi~?M@H<0zdhxC1m9%DD3?RQUg*pjf?!ht_kiv(*+mfes>w_B3D-@~MXpep~0m{+#_f`*-&9?CYI&)ZAipW4+pmv`H2`p3_e-FNu8=-E@B?pWtl z!i>Smryf4|ZJmA4=E~Qi7WY!e+^e(-FS}V~Z(!$=b?FmTuj$b)kD}A+%SRh8E>xnZ z{?c;A(^};{bp6C@rU$*RV#eU)Qx6~fw$6w9^W3lJ{yq2e>1)>aKXw1VgXpqWbcyjc zbSW*mtP)+KMVDft3o$tP)WZkAt-9=%{r!8n&SyQRuRX*!oy0e##W%#@!~UM@eAajRe&>~;FWwANy3MJ#FH5bftck(Nryf4|ZJphchdX|zO*B%~ zJ2k6Z(`AQLjtlSCrFZt+L*Cl)ld*zKJ<$_=XsqeCpwYpYi@4|BM3uCC;>6 z9eY1N(YFIh9WrlM6C6^043FL4-$*}h`q9|y{5>9rr@V+=?*~^I)o(y--*?Jgaq+&k zHO&~DeCpwY-`3gZxRCMU%>rZ9hpRq!>^>e>eqX)(???Vr_g=Xz*M*(*$I}ig?5|hW zoA!_GQnrzgP8GCn#H(OGGX^K0didbCbv~79uMJDL(O$hP?$vx+a8otfuTrfy5mQx_ zs>g;LUAjYE?eO}1pwm_Lre(LPbN4(^wQhPW-sSt+jKRsL9zOW3&F9B?LzfGpOL5WV zfanq>y37z=o(wm1AqFR(didbCHV^Qx#SPy~`D*wkCcELAZsMD!;+qWOn|rwo-w=b7 zPd$9_Gv43hpHblYk5p$;Pw%Pvbcjy1DNlsz5mGj0L)FJ>=4dZndq*1WzP9@9YwcWg z=ILwn~oA{V9IQi7W2fwZJ`LV#QNuHv`RF>W^ z-nZS;Qx!~|I$P+mmFnBbq;H;lJ*P%ixc9s6?w9IENADGZXMZZM2S+b^R!FFc!O5o{ zKKQN8XTIptS9Iwhy6hBPs)#OaMVAYr3o$tP)WZkAwRwP#5Z`1I-~60s_~z9p!#8Wf z4BsRb-xLww5QCFXJ$&#p-rwV&QQ%qLz^;yCTd6@wrVJ}qXQ}dx9M$Bq`xW&pYo7s+ ziu_XNyH4yLkubX+S9aCcLdlEh^a*zs35hIW#^B^r42z)4jiyEP{tL}AWtyb*G@m>v2KZGmV{r1RhYx;R z=d<-d&Xl>gv{RiAWE+{i)k^hah|linwQs28ZMwegbSr`G@in+!q5HY@)H)N|2IVWE zw>_=1sGz@_8H1BgJ$&$6o6lC!rKIT6UUa!Bx?~kyriw0sq6;xN`P9P)zqNUQdx~%J ziEpm%F?>@)eB&X$dE;jIWn*;=Th`A@GdY`<8I+&=M;!^TT$)`}L*y5IV)9v_?8clujlIdl}rPjTwWJ&*uPq@LQWti0CpwbZH^Fv=m*oi7xM&8M+J_?2Ql`L7!_Uc>Fbft2CTd8JmfBPd<%A0CipMve5_$Sndvt1gtuuUHQ!EQjJ zsT)e@VIP~!@^vU{#^B^r4ax)6htPd$9_ zTbl>?1o2Hj@y*s_hHu=(H_gO1uiF{ESs=b41}C3-_~2)}zsEnL!0W5dtIB
#RP z*S6~TWhztGeQoO;yQ-d#`4BR;LOku$z1)gnX>#cE%V$Lwtmvu-uJk$n*twJ$gOg7^ zeDK>k`+U{1tXy5PkLp}DYk@?MwyL*9H@ypa_gvi>wmMgUZw7rnadc>{da>X8JJGVz zo%xk?RPV?HuM^ZTV{r1RhYx;R=hNXz#H~`5B2_+bw+GM4{#0D&v)a%K!pKK1ayZ|m%r zWpyu7;^#nBCA!?5O!s!HRlA;e4Bq%oy?NEsb-_J*J)z7*=Y3Vm=HB?K$NEnUpKX2)fRj%> zeDGVFkB{i`<+hS|+-$dmwe3Mdq zQ(Js9O?*>Od_xRQKK1ay&v<{2e@20%rMI=(ld6p>>fSBQqUFm}k&3CZ9({aO&2hUJ z^P@pLy}!-#t3%7>(Akb&8ds-_tF~|O&%z{COPMh^`P9P)zpb-h+ob!i_)Yt${!YzT zuCm*vV$$~rou2!Jy0&Cvr{aS$=%#z)fAk6{u1_a++gx>TCEajA_bUxuYnU-O`P9P) zzpe9OeP+F8{boI9eP_MrbAT9}eCpwY-`acvM3-`+OStH=Safj`U1p0esYDlIaPp~# z4}NR&0N-EA@Xf=DhHuVSFnr@JzVQ^_%sgWF=AgIX8)9(ssfQ1K#`}BxGYYgn^rgW5 zR1Hb*OA_G_3n_IclkH^nFA(S?6JJUwrf zlNp1PPd$9_+dBK24|kL*<=RQj2za}3!Phk^@A<$%$&N;=v9Br&Z&*Bu-tU*b@%ZSx z`bJpihUuD>)_1BM>oQ{+>SnQ~UpwqZ_a!-d_xRQKK1ay zZ|i*i=KXsL^gfnzOP@wf)X}m>8s07tqB7*FKV+8sdG&Z__|aswzNqev=7pWRo>5m^ z7u$86opiP6I!;Mmx|lII`P9P)zpb-#TkzqgjtWq7l74@=u-_&XynE`|O1&PbHT&K# znl&i3_Pmkgly8l~y1~}PH9j_Q*NYZ}?k)YQk{N@OPd$9_+d3ca&vU=l(wH z`=7f1e{a5_OJ~uguISQGbouAHp-X(xrL*Wl3{F1v@WF4ZE*Zr)ZN)bsISt=T5Z`1M z-#l7t_~y!O!#Bj>cdN^eitm%_&Q9xBQ1}C3-_~5s7cHfKaYwv3}PQ|a5 z=u?^7CsgjXTO&K)`K6BENSh^dLT8^ZM%2!>6x)J)^l9gOg7^ zeDK>kAD++iyq@Rxd=79uz;!+Qdtz|%sfQ1KYxCJ7x*Vu#=yI-&q00o(<<=rYm+qp= zi(!T?#Ngyp4TM zwduVw#Eikoryf4|ZPhEvcj)8cB}S>kTN?GSzj;_?%=U3}=hZPP;zPR{HCj9Bvi1ej zehT!^dpDH)YPZ2#C#cars`^D=GX^K0didbCbw1m?q7U6z6rwiHTb*v`(3@&}OOKgV z0@LbtZ(CmT9aC0!ED-N|ijO|}%M8yY!{>I=$v#Kb8z0!ijKRsL9zOVOollar&bjNR zE2*-kJ|BH#M{gC`^hNoTgF}^PldJJYdR$Ty-`74r()XRpma)#G0aN4aDY`g{E-s?WF41M@CqoxvaPp~#4}NR&0G}+r=`X(7meuf0 z7xB&4EQW6eif`6cH+(}3PCoVU!OwWsj<;VMeHDJ+=&J^=jlSx;(deropN+n9l)l>j z&gd&*aPp~#4}NR&u!hs8d4`@Vm?C|3wfStvVnM5hsABQz#Q)Z8l}ek-t>M}nhgFgq zjbB}xctvFm9$K&d=$mE?PCoVU!EbFI;H~W&IIK7^PW2h~J#qYHCshLX-cJr)`mIV7 zY4Azya@IM!Hm=}Tv8uiu^{h$W%f5Q@>xeEn`nE7*aPp~#4}Qk8cASuAS>`jbzxxfV zJNa{wZ1+`yiJR*;ub4?!|1mYHeVEoKZq!`dp>cB^QSWr~TK#(H@a#K3mJaQ0#^B^r z4-1G6pH|fPU<0~4zx$a{4=B&Trn=#@W2k{LtIQi7W2Y=jgEyeU{zad6nMV~SH zYT9qFfA!UfoJL>usBZLCqcKKb5rdOYJ$&$68sFCQtk?8ySWCaP-%muSEYjybo@;gl zM1D~bg$6%x7?oUCeeKs|M8C9}7@U0S;e+4S`IISgqtS=6Fqnt zq6TFjKegz&n<_egi=8c2rqwg26kC|=LRoz%W!It4O8e?_dEP&2ez%h@vDAB8jq^Rs z7@U0S;e+4S`PA#69&Qf|QvQcFyVpx_P$f!sZQahCpVg;IJ7Y=?&86>UtD8N_wW1!F zARyPV+I6-2{ekU{O=@Vy;N(*eANT!cEdK7+ zizc5@`PR3qD#N@~MXp zerxjp|Mu)m)%d}oYQnhk1G6@NZuWhfOu>iy4Xm!4x?X5fIbjF=WW}ZPM`!obU3M=y z^&in(rQ{q~l=uPJ?%TKa17T%)i0 zxEp;%3{F1v@WF3sd|S^ak^Mc_@r`7?pTYbJ52>+gz1Ay+*q<(+#9MzJ60mbijE@;F z7My(Q;e+4S`Hb-S=k)9?T~(RF*S2+A6|UYj>Nn% zqNB&B7*$TUpSswoUq25s1}C3-_~5s7KHui;7%-^ZXyyDO;|sc@r}RC!c!w;J0-?cfUWFu=3+k)wo*f z95uHjG1q(g`sMA>y`0XE)XQ{-6usSM=^7XG(2iN=n=v@~To1qpzqR?C7hS$SH*~r8 z(au9HKvYVhUzo6s%1%Uwaq^9 zx$5;{;hLj%8MI4*x80|uD6V~4#9N=fK_z{-x^t>;E32C^IQi7W2S4LkJEn}-HruJz z6ty=r`Q_OS&Z{q1&)=*cm{1?|dB3awFBd(jwePUZpK9t+etYY>Klalt)5JUN7O$-t zgOg7^eDGVF2e?P`i6MFCMX9!BJG$oTXm4KEKbG~LzCXOVk>9r6tMq{G)%Fxm5w3~B z$!DJcAN-7G>3C_Ap-YhHa!GUAOELHHy379ly?L)pG~g&+23=W&w5Ya&z|3GTAM*V%+CRE@~MXpep}~raDZdv>P-WbU9Q}^ z}yyu2q@sJT3g984O8x|rZ11Io^e5wI%W(`KK1ay zZ|i(kd5oCS;?WG{)ywWn?eUk@$1x3FWUQ1-M>&q2Ff>yMJ>4U5srbd}$A14PbbY3d zZFKDeAA5G`)6tB<$)_GZ_-&ofv$K2OHSc{##mIi2{k^|j=d<3^_X`hQIJ_m@w$)_GZ_^r(Yyhop8ffH44RpsGL$H`B&sECFk_V>0u zQ+@iscyqpMdTqZyVbN+si|J+;bDVPA<*6^fo%L|u-D+kGPCoVU!OwWsj?Fff&$oQ! zWHq-*ohg+^pHsQE^*&MdLjwKv>ALhow2Myc+#;c4&6+ywtWVZqMf`MN^{JCqjc;Sd z;N(*eANRHdd~XJde7$oF*y0u!v{a(Svoq3E-gftZ(9vr4vQ{VoeW+2 zi7wuv3o$tP)WZkArSapQ&mq2f|G@A~P4SJB_~!CH!#8QfH-X|CVsP@QhY$X^<64SY zpINV2zgf>&-&ybJdtz|%sfQ1KOXJ&mp8Yxdb@uP<=h@%?r}qDuceIdY_tubv-PQoBDs`^l$uIWq<)pL+P< zw^g6W5W5fGqJmVF+LKx)?RrqHaerI*z|zku>CAZf$H!iy)XDtuWb5dPI>++@3A~cj z(^vfNm(QEHkr{)NPd$9_+d3ci=j_+nzq6lbe@~yY-V=k9Pd$9_Tbs`c(WQ{+5-hsB zn_=jZM|64p%+O_<=t2xmKK1ayZ*3mnJGUFY=_$TxCcX(4--L*74m>w}bMBGh8)9(s zsfQ1K#pd|z`P9P)zqNT-!?{1t{d(@-b3dQ^`>gL==M#gI zPd$9_Tbl>?+bLeR9VYmzj2o+-IM-{PiU_Fk^kee-s$EIfb!S)UV#BhYP=kKDonP|mr*bMc zplDQGCp~rU*6>eXD(fr{GDY_5=cChf9#%Q$7C$ovC!c!w;I}ppaPH4@zn=T|+|Q@4 z|MdO;7oy8;(WRg0l2>$@EV^V8T~dlJ#Ngyp4d_xRQKK1ayZ*3mn+@I%uJ@@aqpHE-2zH^;V3{F1v@WIb`|L&3J^E|KT z`906`dA`qeJ^OoNaPp~#4}M$cx2*HpHA%*eQc=r)-56iuh+5k!b+_^vzp87`_78K~ z?5LmHt?phi(nIH(bT!-gAKv=Jqch{)mT7Fp;N(*eAN;n?r~RC27ltfdrlQIZ+1jf7 z12ww9+nuGCXV$w*ylK<_u7_UrP}}_|)l!$9=DOfqg`WDBee$<`KK3zVaPp~#4}M$c z!}EEb*Yo_I&jHqRuIt&~6N8gaJ$&$6n@=av<(B9&UUbePuZ1K&JiH2_)if@R)$)_GZ_!-aIk>~R~ujlza&+~b{&-4EO z*7JY9pU3z0=v#CB{ttiuuU~~r6LuV*r3#$=vL(;Zt7>}HLG^ZRu+wQ{hOH{}sHA?D zAxFuA=e%{Ek?tuf6>jmn$^Cio|1E3v-Vyi^qWETPI)v8He+z|sfQ1K#+)&<@3JlSN%;N(*eAN(=CkMY&X?5D)IZtARewhLN*>UMw`gOg7^eDE`# zwWEFBc{3MR*r~Rzo!h!dny;$lsMeh~<#5%*8h&3?sYGMFbgk=w;KO}%icD+kR`Cwj z4JP@=A8~e)8H1BgJ$&$6n+N#09}N$TT2)kK8ntLq+K2#EvBDweui-1y+WG$7Uiq9= z)z4jO)T+`8b=#%it~TXXBZJ zE4SCL``rfW@wEq0&H-bUeS+;5^neI8EaR9(BdVQOLykYqarWFj<@GXPaiNIEW(-a~ z_3*)O>wKpBE%zOEYL@b^Q+RK@psT9xqHOVMZMM^svR+Pp`AJFrB=J!92^YMzd%P^6 zC-b+{XKHPDRAzf;GX^K0didbCbv{j|IJNG4X0w_x@$%9!JwK{{7W`c2owtAv?{jYJ zlwaQZ+mMvW24(N5x8L%pk+}W{y)1f5=-bA@W(-a~_3*)O>wMDBuIoK2T^F?=d9~gx zUWKXVD@tW8Y=2KJbhw@`{_mu^|Hdb4&kV}1vlZ$4+#|D_u6)$v@}AS~W(-a~_3*)O zZ9XkTm${FRk@~MXperxjpFCf18D;N(*eAN-7G?RY!f=&SE#jlL=^eKkA2(N{Aw7=88gn9)~ZFO0q-1}C3- z_~5rT4{LafRiE6v!>Xw#GQXDoS<~dqx>wbbKHpI3c-o7x@42j2vgr83cQWer1MXe$ zd7I6Q!O5o{KKQN81H7u=w4BBEk5Umi6Fr*M^@z%__vo9{eqUAd_w`?VeB!7(XC6~0 z=%t5Vd+th$g{}>Br;P`a2MuU!#^B^r4WD55qRUs&<(25NQ*@~;x)6htPd$9_TN*#^`G?I6-_#c0%#JjC^YfkIn~dU{ z`r?}p8w}qNgOg7^eDKE|*HXN7rO{Vqq_0MIG5RW0`rJYKDy)RjR}cFeeMJmTKK1ay zZ)tp6&ksnp>-gbEJI+c7fRxx96 z@~P)~0DfEN(_wk`wuv_csB{gh?5y2plS&%XXz3~6N9t*&aSL-aOs##67pjsle_|*RaPp~#4}M$cbFkv0IjuvcDDU~HhK*WzUX9=Ut9H><3AMA= z^ul#x|9ii+>jZv}ajT`bU+-J6RMF;oYlWENbzIw;F*y0u!w0{u^YL%}C0n;r$JGK^ z@9E>OTgHi|BrL>FRk@~MXperxjpzqris&5;I%Z^FJ9zG)`Di8^NZrkD6;TQ$Qs z#Ngyp48s?@SIc6IzA7PoMGQ_p_3*)OZ64O}sjk6I2dF9L zb^dz!T)HpwJF_jNPsSHD_3qPOyAQtU7yEtJvAe2$ZDYpZ;_e{B8kdUk3lo%+J=gq1J3V4$V69~1x|lII`P9P)zqNUQ*G-eHM8BhtR6n`SXFcC4eILJ3-Vx=` zEYqnL!UHnBC zVsP@QhYx;Bebec#Chm#W&wy8NML~C!c!w;Ey}5rFdOYqp!G* z&n|tHO!~@B`s(>9qpuQ4UlD_oPd$9_TN>Zi^X1x1Tlmytskxu$I-m7^kM#Z8!i`%d zeb-$tJo};F4Da5W7@U0SxnBUkt@ByiDf8IoV@9c}CF}P);d(?>wy!<$T-vW{Z}85E zA+hiM{h99i$B(-_bnRW)9yg1<|6ahQ+0i2<8=EmW`P9P)zpe9m;A;++M|ATo=Tke6Kc$Jm$)_GZ_^r+7qUh2` zbjc&S)D&HIi!KX9mlmQ6F*y0u!w0{$d4RVR-~3a;@XdoPhHrX{Z}y095{qxn7B+lC z3{F1v@WIb`){gsH8-2xop6mK0vfk6@Ckh#TWhZ?_3{F1v@WF3w9@cQz_vH#DkNxjj za9@x8J=gh@<#TC)%>8Z)l)9R* zO{$`$;_1mP67N4!pn#sY|5X~-rB(H(hXuL?sK)xyoBP3Yf44AWaPp~#4}Qk8cHBAe zQK?0T*Qv?d8l{f6{gvu+dXioF%X#$W@tX?GnORFGANX~0_EVkoq|vWK&P^Yv>wjK< zv`5^)@#;p)^pZ(*84xb{{KqVHhuLRV^pDhbBiy$dQ?p)yTA4vhwm!co#LM9%H-1( z>P^fzFiRz!IP#cB(!hp#M!#Y%3+Gs;I>XB2~4)@iClU6#OH1iNM1}C3-_~5s7KCI8I z*R0>H=dACn_k0czgOg7^eDGVFPY2N@o9L2Sbg3@7EEHYJiZ1Cy7h-VosfQ1KYx4kK zFTS~&*YHgq@y(uS!#AJn8ot@R$M8*Y@eMIJ`P9P)KjT?DvOcq3vwpLlv%a(5)Az*S zk?t@M+n-a6h4e=qL?zGe(gKK1ayZ%r5Q&c837?y_=?^4S$S z{zBk0752vK@;JYo`qKJ4A!C+R(_b2fcKPDjUSI2^T@xPbt84cR8ge4*5HkiRpL+P< zXFO|1_UG)^*}t=&XMazh|Ec|dW6`Cu=rU7uIV8F~6=l(v|16=15gOg7^eDGV-3%sZEh-dBmW~%C~a&~W%=!%M;en+wwH-1AQ2o9`_R->gVy_$D%s;hXoV4c~l-Fnp6md_xRQKK1ayA9q|! zG56=WU(fw}?&ouVpZop)t^5BxpXYf!&+qADp6|0SV1G{xPCoVU!EZ^Yxaa+{dFIJK zWra$(y3GFQvX9g_$12w{Ib_ud7v1U-_h`IlDx z%@~|~vtIb%k2|iVIK#G?^U97sr3!D!o@qzd6#8qQ-qQ;w@zlTS1t-1Lri*rZ@Gfu9 zox%F*vnS6&yUx>zgL|xRI(ms2gOg7^eDGTuKkj*+&-1*V=lAq6*8^PF|EcHyNkx~O zqRTtc<-O?QAiA6sUH%zk=t2xmKK1ayANTLK6sHs4Y+q&grhxb+a-rdy8seJ^9Sz@% z72gnplTST-@LL)`?s=Zi^Sqwt_k0fUe4p$3|JL*WVsAg)D1AJG+EMXg=)K;JRILQA z@tb!Tuacggac0?@&1y{gvPoNixS&1_{BZbpnFlKGu8_3JM?N)UaPp~#5B|9NSc)G7 zC-3mdAxxz#PswZO_*2q#pFg^T&au?1|Eb*lb?0Ig zs`}g*YR2H?Qx6~fmd1~Je$S>`D?i#@RFg+0OS5-LIz6QLu+k+TRMo>PuHV?Oeox(^ z@>u8Hr6%je7q87J+J1??{wdemot`Vr7@U0S;e$W!xR&B7pY|>-zNd%!TxsmV!N)hM zOCJu^c(nVW8a+ScX}EtXz5c`dEgu>c(#?9>ziLvaqAs&w>CX<8W1o{s3{F1v@WF3s ze9mtcT}F#8y+oJ(qDwu|Wvb|MPIMs#C!c!w;AcE*$I;@Ou!4qfKDRb}Q$T#PX|my) zqe~3m{Nrl)h8UcD>fwXm+C0D)yfOO9LHcS!e50=ho-+C>vGkR{^i|+aqpyg;$)_GZ z_!;ltJ@)@R*>!fCZ0gPPyHy`o_ft(n{E&6iu%1ZV{r1RhYx;R=eINdgGVR(uTU%UuyLEywy*iaMIH+g&oh?p^AR}Ytg&YFMV|K z9PKkss@2?#!O5o{KKQN8XQSwHU375~T?&dWbw!t@qRTGPg&3TC>fwXm+C0E}h;Q!Z zF?^F&d^1LT(?@(08D{t^)R(pR^d8+|pXj?q_3s~CM1 z^1`z~mHf zl%Ln=2A>|;X_u+rQssH(ptq!-T(9~XM>7T|pL+PA&(0Y}}VUj@!k_l7-r@gu`! zm7;j`xUQR$>UvR;X`Y4_*E!cbJ+)v_J$px!V){#ga_w%XZ=$Q!Z4)(PQeQo<`iUui1;^^y z**}D~dosxedRmF=&OR#SH$4tQx6~fmd3aB{F5nX_U7u?Uv1bsWMjt? zJC&F8eZAMYS3K^JSvR?}-DzIRQaaDhos(Bhs;bX*Dq3&<`Pya-PCoVU!Eft))`#0= zUl;q|A1M*_viIR7dsWj-S62G9|De8|+q*HQYj)lDwfE;UvF{(xU8a4lUu|mXfv=i+ zHHv*NF)=v#)WZkAt@9aB(Ld(Q&Xp=tg6GebTKp z+8k`7i!NVNXIQ4*vG4!x+{$&?05b+BpL+PPk;KPnNU%NV)eIInN z-B*u;qxA^;`GwQYTCNA4EAr;umo56h4t+1+=q@t`C!fy&_~5rTpPr&iO3@{+=#o%$ zNiVuMi7p*P7h-VosfQ1KYx4lFE55mK&G3zb_-2XtCZG7`L}|k}9^xBfaPp~#4}Qk8 zb}TD>^)1rqtM1bGp|6a-Ivj2E)rcKNU;XZA^c68U`P9P)zqNT-!yAllGWO&9S?WxN z3Z2e{UNzTy`g&x)QmvYP@YX4B7wVe3Ry%$5^s%oFk2{+&IQd);zz4szd4Rjyhi5Bs zDoE8ovCZ{Z^Mk5I_>g$1rhZmU4js{r)92Plu6^!z@MuNdX-b&uvsZQX#Q_Un9Qf4G zjKRsL9zOUP&)V_!si_mbv|Oithwt0@vd$}YJx%BD3uffeLpG;gs!G<buh?4sxjyWX9m+Qx6~f*5(0zGi8x{^ExNgjr{MwNtV8fx!%*)c}qrJ zsk7{W9#DAB!rON*>SN7+gv?AAWyav-b3Fhb{ETPmcu92mM|8O?y7UoU%8M@JMVDKm z3o$tP)WZkArSapQ&nUjhFTN=tz9}xg8798TF1}eNz99xDpL+P5@$uL?_F z&6M??zW&wM=&PZ(jJ|S^z9I%EpL+PcfQca+WL91vh%0fwbVO$1q_I2+s=%^$)_GZ_-&m}@;=8JJWsPn zRd#IeQ~1Ys^=S3o)i+xe*GFzOv7hm@sSfS;sBfNP{dESnh4JdW8mGqwmdRQ?(=;;% zC!c!w;J0-?iBnCsYkfPTo*m&mdDsje^SYk(p1yD3P&nAL%tM`TwS9E|)6X<9IQi@o z;Dg`VeBz5P6-1XiqD!RcvP^VIDY}dlU5LTSryf4|t<3{`uK1>p_@;sQ=9#_Wn?&N9 zmExO;;u~Uc@~MXpe#WzQJndxk)s$C8U-`&&{zzHx>2qJ{t4`O9z9I%EpL+PA%J zc-QEKj!#pqSM}t&pZz`8`B&s~iN0T*@I~`#RYsVf1K{LS4cln=v@~ z)WZiq<5@eNaxC(yV8??h^!)kdX@14mo*UZV^KM*5Cw#rJz=}1o_aTeuV#zNK)r($5 zKF*S4n)b8r+T3x&EHeftpL+PA&(HAmtv^>~#{d&qvD{XN(DtoQk)?~S^5ni!mX z>fwW*@hlyWh%U=Tm-3>^uUCdHUtA1b=87(jL>FRk@~MXpeoN!WJ%7KE;Tu=+O=|JY zaPiG(@y%}W&9AkFZ-~Lkryf4|Nna6zlTST- z@LL++*7K~-tk%Xt)qC}TS(Pg~oQe1SoDZ11WU5LTSryf4| zZPg{O_$Gz;#!-Cp=855(E3FOR1c`5EiEoI($)_GZ_-&mJ>oe;$>o@B;>pScHe{20` zf6jiL{X6@4_V@JppWgr92s&76@q`hoiqFYPhw>g$Cy%U&mnZli^}136-*KmM>qgVP z-i}D4bg8Ig3)VW<*FS%iU(xSIBQpjkpL+PH#Sh$KN~JOBai3yFb_1F8b({q;+%k9IVfrJ~CqE))8h5PCoVU!Eft)*q^gs zXaCNAp8Y+a1FZMN;N(*eAN?ax0-{T#=yFVSAqFR(didbCHV^Pk z;+t&Zo7v)E&VHT!JNtR|_k0ep-V=k9Pd$9_ zTbqY9ocr_Kujl?f_w)H2_)p#czjws3X_s>&Rp#ym?dJs^R!O&cG<}pmMonEdDRN*u zM}6(s;OWJ4c<7)y6}Lrud+U164O@CQ^EG2|@~MXpervjb=e*S9V};of>QSS#A(2Ht zsN{*_&ptHMDfT|Zxcw1V>g%5F^YU-W(?cg+HLg#u;z7F7!-1zWJ{fJs;N(*eAN-7G z?a2Lk?$>kwp8NUqHP-`N=M#gIPd$9_TboA~(dD-25>Iq-7hT$mF7ZT{@}dhdIQi7W z2fwv>fJcaL9(6N(b8@TUn>ON`aPdte@l7W24KXwq*J1`p~VV8-C&Qx6~fjA!i_lJsWF zq$5tMu2pWO-q^}ce_iD|;&3Yuefd;S>H}Rm>v|dT+-`bvjIOjQVXsd^=4!W^wTD-) zzu1hy$)_GZ_^r(YoaggAujlzaeaw2!zJUEbF*y0u!v{a(SvpP>U0RAR??jj5qDvLg zrHts}EV>YblTST-@LL)`?)j(ghHrwzHzDGi%Ho^6;+vx4n;haBVsP@QhY$X^<64S& zKF{-dp5OC4pXd8L?`MBc3{F1v@WF3sd|S_7D$)8zI`CD< zCC~JowM*Z9ubbT+sC%bc-Z|;!ATtIhpL+Ps>>}#Vb%=cralydEHIh5%#_hdV z{_5ndb=Ls*Z662fA?ohoU*+fOxgpVAe|!nmVfC{X-0iT*jKRsL9zOVOoln=tiOLP^ zHdf^sIxyr**%RuLd=Dt0d>?@8_&M@?fctvH;N(*eANTbl>?>7Isf-jy(XQ(t`3Ong&9eDk!A;hP%b8)9(ssfQ1K#4zgv8+FQx6~f*5+XiZ&IT|{8O)=yJz#cTFm)EA+@2x zrbKgBwN^9w&f9hC;y875^W?=7+*T{+^i%5PyuU+D+jD2KB^7@U0S;e+4WJis>w zHs4k@=Vq07x?_zFksnn2-WMu_R(96AqKEam)Yw}On3y4^cgg_$T&Jxze`Aoom$k~x zC*#JLF*y0u!v{a(Sv&TeSaa&^A@|jjz;p?8CP#fX()V=sKE8U2|Kb{1_Xg?yU9#rt z`(=^N-mmhlz1cSD$@{vrJ#k{I8H1BgJ$&$6n+N!rA$#3-t*)l%)6MdIz?Jg7AolM` zP>^Ay}7@U0S;e+4O__m&J-M`kMW*HqmbsPnx)TX25t`cDi_KK1ayZ|i)j zetUVhSI%21tVsCb5-A*X?s-GJUAuVeuqA1N27ex?)3`NkH@o6oUGPSm0TEwAb+uV% zU9aTWWX9m+oAts6zpeB6U&$*cx-68x!(Mb*EV?8TU2=#nT67@>C!c!w;J5Yf?IXTf zD84x^z8NgOsV=^GJ=O3{%saz3#Ngyp4GRFC-MiWsOsx}VTKUSq zbZ+f)F4^1CFY}o(IQi7`IRL+<@#CH^U-qHflu7f|^3Tos`}Db`QfzSTFeXJhJ^%Hc zwpY89(@xuWc+cAJt5e@@_Q0oiXMH8BU6|*T05b+BpL+PpgwEZpyVuI%ihB?bxocmxIgb;|He}TU4lq&VJb~W2a*E z%ov<}>bV|(-_rPT&%4xpJ@(h;QR>R%x;J849#M5NCyVG<`>WbjS>3FD%Te!Hwt8-= z`yTpxw#lR0v2EN zA1|O6ZU1HOI*j8xpC!uo&(&oe^tlP+6WYJ*YR2H?Q_sEteoN!WJ-;PbZA!a$ zv~msWek(`dQMK!z(XUc0`=*9&_}L@oMqd4-Va)bqM?H1*3e$FX&Dv0(bUXjWV|5cV z1}C3-_~4H_uBCWfqJGXz?6#>Nx$^!Fm!wxn)l<-i*P?ryf4|Esg)bdVYlH@>z7rExNQ9U5bk?g+-SPq6;xN z`P9P)zpa1&Ao0zoP{TJp#W&T(H`~QG8N@eZ#W%#@sK^`5?7 zQQ7D#FX^j9(pSXbM5_)>R z)IX|xuCLcSxD;|I-&W@;I`C1-;7(=?PCoVU!5?>AOL6p|98KSsJfcp|Yt+U2Q6k;y z%|Fp!-<8vYM$KQo-M5_%?;kIGLcJ0Ca?;Id@-CR6%N1{Zu4=-0W(-a~_3*)OY5f1y z^Vda}qoPY=(WQ^*vPN_{E4q{rU5LTSryf4|ZT zx!SMkb;lte7Mn3R`P9P)zoqg2SI<`yU4BOyx_FB&^F)`#qRT_k#anbC1}C3-_~5to z?|<;n@XhBnhHqYuF?{19zBw$usUyDWB)%aAC!c!w;J0-??WM1{ugCu0O|J73%6flP z`uc_@nNFS;{@w=ufRNPYKm>rzK&&eRnv z)b%V-CB%%u$)_GZ_~VZIzaBRgT{4I+DMgn+(Pf6{azk_}DY_7YlTST-@Y||G7x7I3 z@l9#*O&ak{a`DZ4@y+LghHr?$$)_GZ_-&mJ>oe;$>o@B;>pScHe{20`f6jiL{X6@4 z_V@p({eR)t8#{mQIZZ{q%@!~$-bMBH?4vRL9TMx?Ze#mL4RqD(a{EVJaq-eIJ9n3z zP@|=O^nhe)sL$I@8an-12*`NO;%dh@7u0- zj^?fO#$N;PmQ6oIr?`}N!n?v#^!L0ShE^Cb(~QB%ryf4|ZJp2mO5Rb?<&NmGM0BYt zx~vdg;)yPwL>FRk@~MXpep~At!-Tq`~>&t-^^w_GMSG=9xM0*rjF?>#Ef9)}`U&L&?US=u${@sV}-b5M64BE@ecQNumoeIQi7W2fwXnLvy0)I)Z!araPp~#4}M$c!~J>g*K_}#`}uqha6Rzfy8qAfd7jtv z{GQJNp6|0S_|wn-7p1$@B6F>9rPt3Yp2qXF3fPf*?bL#fIyk!2>AD$fYu9FXgBrPa z)~ou?idda$kZx5mPyMIogUlG5eCpwY-`0P}#_3)qIbHvjdhlsL|H|4y@9(*{)5!7O zdQ1BR4c6x$q_;h4+hSLzx%%#s2A9(m2-Anh`D8v&VUrnylTST-@Y_0{|CPMjqRSQ0 zWvA$}Ms!&#x?C1rB1IQsaPp~#4}M$!-rC}um*IwQDu{2&iEmysHGJbIz9}lcAqFR( zdidbCbv`_w=XpKP@A(|y`99D4+20d`lTST-@Y^~cOL>KNb!t`SUJjLbe&Pek9=B9I z-u85=Q*@?!ek9MVG$E1dMDjtA^+Rr`%}-Jv3Lo`Sm6_x_veLo#W(-a~_3*)OY5ch7 zUpySy`$F!`s&|)%p&9mnP%ar8kN8o|Szo`_<#e59-a4a0Q~$?l1N8oEttbDC2-0bN z_8tB{af}&*lTST-@W&n3Qarx#h^7^W-&Z>l1wa3u#ZkACF{Dz{{Cm^;y%rTV!!{} zpsr`~k+Hx3DcC#g=Y-hb|NJ=7dPZ36?|(*&b1qOm_WQrY;N(*eAN;n?C#uckQEy^@ z|5Ks!tK!RJfB(}Wna9#HoxSyCuOr)%eHy6WWo=&6~zy<$RjhRoByzmEO= z4>36T)WZkAt@EMp`MnOm-{JQ>{Jw|Z`~0cj|C|zCTtydO(dD@4@>q1)CAy>#U5LTS zryf4|ZPmqDd{b6@Q&W7CNPP1s(C|%!_@=)2h8UcD>fwXm*7;OwYV_5POh#XIm%a*5 zVDwcH>8qMQjK2EnZ}b&0IQi7W2fwZJv6NR@-v3=t-v511-cwvz-Urq|-v8|=@Bdye z?-19}{_rE9Qcbv!Ncbwzocbu)|cboz8JI-G6J5FM7@~MXp{ ziurvmzt`pWyZoM)-}nCg{jc%6XnxPj?|b>Z@1Oqt?*sY$?_lxG9Py2Z_$HV5rl9!7 zUVK9gPCoVU!EebYanIM4_kVAa_kY)tzG^SO^P4E|1ND;k0y;`x5rdOYJ$&%T9ryn> zj+XZnua@_JdrIF|m-m0ym-l~XllOrhk@uGogOg7^eDK>kuZQv;hph7c4?B7P#|(Lo z!*zL&!$^6LLj!q_12H)H)WZkAt@Ftx@40v^@42v(_gqw!_gney=p%i&9o} zAqFR(didbCRhMw_&0Tr_$6fwXm*7+2dzWOM?_iZPA zzg*t`9Vzes9xU$z{YT#aO$<&x_3*)O>wGNb)spvrAC=$vvz~iP-w&1de_r*iThr$oPv_kR(ClTST-@W&n3 zQv5^SbG1g^bCq7+b2U`nb5&a2b9GAIb9G+cb43hJKK1ayZ)yCv=h@$LozMF`=zHGp z#rwYg)ce1lh%VDamqemVXVJw`blEPt91>lK!O5o{KKSGQ{g&dD;+s_Bn=#^>-Qt^= zX~uiw5{Pg1$@{;E!O5o{KKLz-ANTx5dH;7_=_}TA`udB!|GS^O|2w1f6)`yZ)WZk= z|7q{7qq5q%w|}wo*nyqcB7z{sT8JH}h+<%1ccPeBcr5J3#1>nyu#YV&c6Tc(sMz^B zo;^SN55|vi-_Ljae7?{74*%L?-q)PhnqKExd+oJ1XYcsq{=0E^dH;7gdH*-p^{nsA z^Jns&;==O&Z(?xrsfQ1KYx`;_?{RdI_c*qc_c&IS_kZS<_kW&~_c)%A_c#)RlTST- z@LStY2YJtBJ$cXN5qZz$TzSuB6?xC)V|mY|r@ZHq7@U0S;e+4We*RnD|JhJfwXm+J276`@hG@`@gxaXMJa$H0RptHP#Ngyp z43 z-tWbAKI=X2`!?$T|9}7YLD6N2=;9>0e3AEm$CfhQ|2<1|X(+l7gOg7^eDKHr_a_ve zm-m0?72kY$YP>hTxA?|ie6wDB6D{xmCI%;;didZ^X#DuUPm=e4cb4mZ_V--pv)(iB z^UJ&<1}C3-_~4H}?!OzeKC@o4ezTslzO&vl?}@?5ryf4|t?lcue8*v#e8*vue8-`Z zeE*|?e8=Ive8=H{e8+(poP6rxgWuYI(#v-)>d1F4p2>GEPRe&KHpq7_a>#csCdqd$ zh{4IH9zOW3?dQK`{r}=;d>6%0bh$3NY!_V?iZ01SmjKa)7@U0S;e+2=U4q0nUBox; z;+wPb{f`)X44uI9*ht{%vDu86_Oryf4|A z$@jmaM3)|-i=*f=T6DQ0x)6htPd$9_C!|CC-&2ZjY{fUW;+rDko9^P9q~e>I;u~Uc z@~MXp{`li26th2Pzs~-h{XF}7<~i#*V{Le9x2bd-A={-~Ihh?$2|-p8NOQ&u3n5Re8)LlzT@0TzT->`PCoVU!5@FzgyQe=o$JB!o$I0Uo$F8Xo$Ffio$Fxv&hC!c!w;EzA^@KgE7dyPWFL^Sd=3*YZ(btrw4t ziJh)WH~91-%q>E>)ksm|(1I(fM%>3Phw433-aQ|rF5?qx#^B^r4|D|zup}6$xXYQKCs|S{cgHtxqFZ1ts18{ z_W7}?OSAbp((mxv^4*u3F*y0u!w0{${dm|r@9leLyc)Q3|LQ7>PpJF#3nzA&{zHvD zTd3Bl6LxyWh0x=9J5|@Oek^`{>P1tXDo@8Olk2uJV{r1RhYx;B{j3&UqC^)v(WQXs zQd4x9FS=|JU5LTSryf4|E%gCDP<->GfZ-ck@y%rM%@Fa;Ve!po@eMIJ`P9P)KjT?C zu6t+9t6Va#uCy}dReuj-UIkY-=GDx%#=Pp&!5#XFe>y)xYXy>)G;a3?^~Kt2RMu<@A#mG9#whc zX0Jwlb=zv=i|%%sqC?JKnOmmgQhohhzID5*t~6tC@~MXpeoK9TyS>{N;;^T;dgDI+ zpJB&0t8;G;)qb@5p&Gp)`00l3>GZm{Z^GX;DWRM9&GE9iM`i7_FywouYVKwXPCoVU z!OwUJ9oLI4BSe>OqDx=VrMl=cR&+Tgx)6htPd$9_Cp3Qi-=~Og!o)XUd=1~&i*I&_ zZ%!{Ye3M9gLkvzn_3*(Tf82!P=`qH<${_Pr5p^`Ot8tB3yMtCMs&b*NwB{yMV5`tIe`ATtIhpL+P7Au@=OXj$zI@fRM{rziZUFoRX#XYC0=$EeZKVMwmQtvL_spsx> zJ+;rieg)Q!>~F^4dF|RHRGUgRAIQi7W2fw90EaAOUMjSuRijC^&aiv-ZV>fGX^K0^8kGCTj~S6bxMp zUY%;j;N(*eAN-d30Po`aF=tPw<7!^+YlUrorZDGy;OqNtjh#E|tvO%+uz3)u-&~GZ zcDKzO9drLq8NT_n)$mPT@eMIJ`P9P)fBbP1iepk3^J(|iKKQNeCp0*H^NNK7RrTtdYI~nL zthUcD*D34lPilngp2SNV7u4aw>4G9UyXlepdX|3@;;Ex&bg%Du(A$i`$)_GZ_^s`y z+zxm0!(?g4x>pk=Or&poT6;ChIsTUS#TC_m8 zUf%x7rAk|NnK3x|To1qpzomW}iY}Ezmx-dwZqem|=yFhWDJ{AXgOg7^eDGW91KeMH zlTLioUwq>sz8NpRX&}DYD83;EC!c!w;AcEb$KW!?yc(6+m{+Xl%>DITLkHah=b4zgp(~oRZDjqKEwv~|3o`~MpL+P*PMqRDlxr z-w$n;S2t~c>e0w0HTBxN;aA;tM_seVk;t5f2kL?;+>fWoHr$NC$)_GZ_$~DTexhmS zi_N`0sBEPNK6lIJVqVv?-lvgyKVzBu!MW3p=;q;P|0+5$QWJxd&prV@_!%#uV;j*W zv*?mcba4}1=7=ukM3+>e3o$tP)WZjVLgUB(eS`SsQX#`P`NcQ8#W(L77`}-R-;@>K z5QCFXJ$&%TA2*?RyRR{?;*J^fitBpTcjozy62`pxaoU(y#Ngyp4U%p zNcClpty{rxMs&S$-ECiDKv*vdLM-+4_8PCoVU z!EdRb*`iBV(WSBI5+=Gt4>ok^DY{e=U5LTSryf4|E%gDu{b%ie`NmOvGgo|5MtoCU zd=o6b=`Fq?1}C3-_~2(eOUJRVjd`_D_WS2$f6sOPXIbx=_g{_}^NJXpeCpwY-%=lz zaMowmYu0bpbJlm(d(H#I;N(*eAN-d30H4<-+xS+02Py~0M*b)4kEl{P>IFq*_^j5A z+Z7bNxuAZP@mgG5gqyAyk@IoOFP=JMsg_5NI5sn5aPp~#4}Qk8bi7+|Ns(_oSE(^G zvlU2H?1?IOId}QkF}bvN4X4qq$5hwHqT44s>eF7Y4RP`7v&T>O>bSIOi$O!o7@U0S z;e+2&AKUuXZ$exCjPe`^2lc&kXj@_T-&`uF`+OHSRUCJn6B!e-oab+_*R z6rYl3)?@ZLJu8;dL63iN@pZQb?%M5DK(4AGwagfteCpwY-&(zvw`?^r*|0h4#EJ9i zr#N3xHYrBAzn`5(f86%f{>(Hdohowu_fuUOYtKsCqFSBkpogyPUpM81Ze|QlKK1ay zZ*4#9&)KiDe`i0>{{DCG|AR!At)k20ekLLOYPgTAP(v1jQc~;E~`%$mVswBF6{TU@2 z)Ge)>d-VVM$)&D7d~INHr!uYd{z{)5JnY+>F*y0u!w0{${cwMt`}N$v=YBrt0j>vr z*Zu#QqDukMrHAP9xt^g*A<<=$=n^ct5QCFXJ$&$6tING*hHp-bZ^FbkUBx#Kj~Tuh zA->r!z99xDpL+P=RCmr&UHR9IQi7W2fwxbB$UVVd7jtv{GRgw z*8^PF|L*7id3>hlEqfqPjnA9((d-^a)V;k&WB&5_tjfP>6dm`ZpstbaZ;w$g-1NAp zORX2#H`Wz4A51lBa5FOoC!c!w;7>@8_`esVsP@QhY$XQ#^?8UISpNki!N_+8@gO?ZRoPi&d_Cn=#ooxAqFR(didaH zJWI#emWFRyh;LSjZ&HbG@``VK#5ak>H^kuNQx6~fmihqa`8?0-d4A74=6ZnZ`v0xx z|1CCNO=atrNey$|=hN$MWApETxE|;t&js1fb6<}boP5qB@WGERmX4i$mV1wgoUOD+ z$$g1NT~^r_=S);LEVb@ndok6;Cyx41vVba+&Uxx=iL$RfS)_vwuDj`x^Ny}&3{F1v z@WF4X5AdE-?R>kQ4pXaxE{6QwH%?_(^nJZ&pLfT(R#e|YLN0cleX&DB?*yy;du zS)`)Swda4UFVd~96ovLE7_-hvft8G59%;#dE~7=W(-a~_3**Zc$SV1h33y%QfZfp zS~t(Pbo$TAE6}&=mOS=)j@Q@4?&X^4?d$9hjypU+=g+#XL3Ph@y6@!fHlxo>He+z| zsfQ1KOMQSJ{pNKra8((VCUEiM3|sps=Sqi)e%`P`g)Hda^JUXBN=IF2($+m%T`c7v z(a!aodLB0XV4!niO$<&x_3**ZcnKY23mUpi6kQUDF3m)j0MX@f7DJbUq6;xN`P9P) ze?sHO|J_S`^Q^hyoA)gZ-z*c~yy;~4W|8=&wD^V?oP6rxgFpVb3B`*;jCr+iu`#bg zrWo^T>0M)9*~q-w6J*S*ZZfZk!O5o{KKK(F-`d}+JjlMJX2awvT|iXvG?{9td1pG8 z9kpt>%9Y5&=1a>}>PKD|uXTA2tGHUtUS0{hq!Nt_XxM1%4KoHOpL+P1VfPZu55 z)r`T(ryf4|t?g%IShHP&PJ}7

QZ_l#4aS1*lq&s$8Nolvyg$zRVO3bx#K)HX+7 z9g=r`%Z~L&>tw?kM!I>8Gh=Y_sfQ1KYx}X;Rmq|Jh`MURvB70Vb_!Cyhr+5fO#F`; zHlgAW+g8t2otl5;sS}=DH>uFsK40e?`my`2PmTfk%ov<}>fwXmQa|-Ymvo{_SJ7pO z=<=b3p-W%U<-|}!7h-VosfQ1KOMQUn65m9MZ*q%oBE&Zf#5ZTeH&`7`v00&y=7jtka^|$w=u7X!O5o{KKL#5VF^zjlV$nZ z@Ji;qpWnu3-PaRa)pzMfRb8{Q&)$z}eu-fZatEf;&aZr$kM_@?iNVRIp7Q|wmiho6 z(skXM2Y&x=BOuIpM-SnY{ zI`ucFHoC@i`-M@J`s&v?QpFC48(_xZwjd?{3PCoVU!Jp9h*8cu&vyH=+)g4uC)_dvi z3Y>3yvy3Jf_`?(wQb8FZN}i_ zQx6~f*7ox~@3C>I8jM%55#`H{usxxseRny#^yPO|-gR)9s}1b*;(6OPy!%*9+das- z*WbUXF4T2IwY=dzW(-a~_3*)OZ9gxbovvXsZmqfAGmpKpjyvo>q?WF4f3A78B%SoN z6&KDPo#U^Y?_L@iT_nJa!O5qd>jC&J^>ajY*(SQw6kUE)Fm&lDy1cz_=yIovp$joM z`P9P)zokCF2a9il#Wx+rH*d-rzS$_gNiDvaExsWJC!c!w;AcEb$96KWSihOaPZ<-KB+F{GrDUqZnY!Op$3As` zeW_1UeQ{Zh#AWIg)IXjNTD>6DQRjVp=3P@?S2G4DpL+P}i}ML+aPp~#5B~V$CKPvD9OV&Uaurl^6Lch_d8MHZDM`>>H1$^{Vt`S+qF(ouudI4 z`Ak#W5v6@}+gj5mubSA-jKRsL9zOW3?Z=^2P;jC7R~7qvuJf#S`|ee` zYp zj~*DlaTVX>6yKcKZ}=vu_@;ySh8UcD>fwW*@hlzF%e-Pge@(9QgJivDo>!N7HM5K{ zuZY3Pryf4|E%jjuZ~f*=JHM5|=6yZ(_gv?*-dB@(KX^gQ>FtL0*2LiCQ_p=uliyM! zg0CH1aPR6ZgVpQ2`L*MO2$eCQV2=;svC7eNN~Jo>v+EwK%0E6+w7h;fwyIa!eKqyA zF}1QRZ0=#k;N(*eAN*FflW4lz=&7wA%~a0)Qh%&B@uIr^cjM?R?kV-*g5xIzWG$yh zx+Qb6DceZT-o18X*3Rv;cNF*a%E|to>--9` z-ZSrK9=dR3dyX*kJOEBU_3*)Osh{DZO99bkzv!}Abm=9!j22yNMHgal@~MXpeoK9T z*A?GP72nJsZusVUBf~eg;+r|*n`+`4VsP@QhYxBghd@|RXxlN5vakG7{ zpy#Uf+eK@Rrp~N$6^rdPJ&l8QZ=GmkrpE4ieXXMDzO1Ze#^B^r4{W&0ev zx>MD{wW%)7X?#{aynOaXqy9f)yqRSxB<&5YuN_1%~y3`O| zoJAL6aPp~#4}MF1fUgqYvE1}C3-_~5tH2Y8pJA5P5) z@1b&)yt2LLstxL16aSg9C+@4mZ9i$w{{C)e z3{F1v@WIb`mX7!4?;JeTb*xGgz4cb==p*W5oEnq$(`S{sV5!1MS{2fJwvBpOyGd2u z&G~f6${~%l%fP=%-AL8kjKRsL9zOUj^#RWQoc%icclPt_@0sVnYyW>rbb0^W(B-D+ z@}YpCOQ`7LExP0rU5LTSryf4|`FC47erRI&W=#phH}k|d&BQm~Ul_iz6W?SM-w=b7 zPd$9_Tj~Rx{W<$}_V4WH+21qIS?`I#$)_GZ_!;kidgT5*_v^WT&;5Mn^?&OAf4YYu z^TNFbsG-OGqKhTluIewi(fuFWXqEcXqg;6gX4X-5Rg?VF-9h(Xy>VUpU+)9kKE}^I zS@GIt3{F1v@WF4bUgcBt9x`oFkZNeN+jsGl(<;Tr4jx5!+vp3Kiez)@T|&>AAD3~} z+FE*4rI5T=d$rJBwi%p{cl9-6aPp~#4}NR=;r=}L>$!i={d~>?zx)1wF44tFbosEs z&?Q222@qYniZ0be7h-VosfQ1KYjp{GVE87t_~u%$;hXS_hHplRZ}N(7>WFWM!O5o{ zKKQNehx_x~ujl?f_w$+8-0$Z)pBS8c>fwXm+I|wsv?|9^L*BGp7*oACk7{< zdidZ^X#DuU7cYIi$=i;B>Rw6bTGjqKqVoBLY> zQ_C(&d9H$9oVG{6v+~}0b%8gJTHWoUzc2ONUi)lsGX^K0didZ^X#DuU^L(D?^*q1l zJizq;*Y&^a`Tt$f#Zh$GD!MchT_%Vw(?pj$q6;xN`P9P)fBb)cLh<=hhHu7-Z#MNe zeB&a%DI>n|72oV}G<-u0PCoVU!Jp9h@qg#}JkRTSe$PB+J?D8p`+H(=@~MXp{`lkm zyD@*C$KUHQZ~1#ZuH*kxfB)AYN8{WpPE1hkMtn_Xv+SgLQQ7awp$k7%!P1T2sfeOF zeUE09d@9$_r>{P1Ug)B?9{g%+w>$$|n=v@~)WZkAwL0}mzbxD7i6LrogDLMz=e)1( zObTn%s&ZDH>D#oc9YeLAf4$C@q9nO+$r>`pO+FByP886{YVrCw#!g+x8ntcKL_Z9Lyqm?`_AOx+ZcnE9(S1_aztz1q0Vs@^|?t`OViiud0Yv|F)6U(>`j)n{i`T{?1J7k%N~ z)rJ%M_cmj2@~MXpeoK9TzY8szzd=SvmHe->w~y@fQ!86US3Egvtt#96a-uPA7u4`K z_0Ep*eyvhv@pv?Nn$55Gfd)>w7@AxYgOg7^eDE_~LdX81OIOjQjp$N9bSWvi>=Ipe zh%Ut7vs0<+{fd7J~24?)WZjVLgWA6`umO+LD}A>?5~3IRDaQNb(ng1Zq_vK zVo%kN80Rwss%6xh=WQ$X<7pYa?cws0i5I!)Olt<@T->LI8H1BgJ$&$6tIzq;`#XB4 zo}jMONcztC)(I8eZrk3jcYdg)*E3|#mZYeTS=@77{ov|)Odap$_X~OJq$6jnel@d| z8H1BgJ$&$6+t1R2dDG?(@1VR6<{Xo&?Mn4vc+=f8>RngweS5@qxtUnE`aG^-iTnBW zaF0ptM-?uo7d`b@T)exB8H1BgJ$&$6>Swp;azJ!BExJq+U2cjlJw=yiBMe=L!O5o{ zKKL#50p3x3<08HpD!$n)z9}NU=_bAjTW$D;7@U0S;e(&?{v3ag0&S1w4Ij{?x!O|U zh}W%R!RkxiM#E=UIjb)0+Hf>w-H)nOllh@h*Rts1>nGkW)YVRxx$R+>{6#4<1}C3- z_~5s?`-N|B=&OCy;1pk9Eb`x?+V7rr#=YMo)p!4!#j}V0r88bn9_d}Xq|Uo-N$t4C zRrJt>Yxk9Z>2Aj0- z6;*%ylSJXGfW zm!>(gq#NH^zgqM3+$`t8I?c!2o3HATW(-a~_3*)OsUI8BrIzS&OmwLxx{MZG=7}y7 zMHgal@~MXpeoK9TfBnbs&F#{LZ(bcYe3L?a^Qyn$n_A+VOejQ8jG zdlY!<3)?-v9N?NqzdLeaf=4%d z{loLaqU1H4%ov<}>fwXm+U`d+@Abo`#Q;^;uGPv_skf_tGWA|NBY(8oxO8(DhoPCZ z(>|NH`h6VqiexTfHTJpd0th%VDamlUE4 zF*y0u!w0{mKEOl8H>WEbzS&mE@J(g$O>Xf`5An?!@eMIJ`P9P)KjZy5{vHKltKQl; zDyowzcV=SZ`;%9z(s}o;uuF7H)sGoaX2;#6dhCYm;aF=? zdgw^khz7B>^ibE?TW&pUV8-C&Qx6~f*7nnEgY)$Al@F;8a^27Vp6h(pd*=O`ZR5KK zZJ1-82f)dv9zOUj_0vFfiM?g$(pPj@D!Nn?UD}E+%|#btaPp~#4}MF1fX@)$oI7p! zCOnVfo42nF-;@*I_=|5+if@R)$)_GZ_!;ld@%Jcj?Q+pouJ3)-sA-AX*O<6Wed)fx zy~nZ3>dN15gU44&q#b&>t{9O%k6yQY_U+=8?X~a9rpF%_buwdc@~MXpervn`R3rP! z)s6#H-EWQudUnALFYZu#{JF$cfHbYZ{k;p zYnw4R`P9P)zqS3aKC@o4ezTslzO&wQ9v}uMpL+PX`$1B!Qi^6R^(AvrIM zSk$h7KAn1S(rKH@=>~DlXM5+aV8-C&Qx6~f*7n2xoc%icclPt_@0sVU_r&1jQx6~f zmip16%MQ`SU33{Fx&(5*`r;c8s<|nmO55FjrY2wg z_HI?xaio} zQ0~nal;4WhEqmSks!kl6HE>334xQ9z!$a3ICG@z^oDt*dIhZjx`P9P)zqS2vf1dmG z+`s33KKJ)o-?`2w1}C3-_~5tHPjAtsjp$NUba4`0rim^mS{k~v6tC3N+;C?K5eS;Sk)R^U)S52UVjd)b?ZupQrdU-<~N~%N_X}MEPZuL zH8Tb$pL+PB61 zeX3=N={~8EzNeeJJ#XIBjKRsL9zOW3?I(JnOY;I{%Bp03(Qn%C>8mpR^;gcd$5yHv zV^YLC`5dJ>R=W4I!R{C8c4yBO{m*?iY~K6 z7aP$o5JE7VsP@QhYxRz-F~&tKy~)>!4$c(?ND`cr;Xdt^MzWQX~K-7{j%ukZilBv*Dt3%st5WH zUSCa@$X|KM{`Pgu7@U0S;e+4W?$e&ja{flK@#_BNrgu`u9al9!*YAFR?{{_nrOR?% z(oR1XokjtxH=vQusE^y+MXyUXfSo1Qi1?Rlc|+;CeG;r+^t!O5o{KKL#5 z6M5FqWs&Ie`jer{4AEtS=+a4aSueT}gOg7^eDGW913bpT@Xc89&5K-yZ`zA*JjFMu z#W$Do8NML~C!c!w;Agx)$KRtszl;C0@toXOg|*M4Pv+jLLYH5>KVLsq856H_tgtD) z-t~QDubn=nbl*z>(J@1nE*Ba5!s*v}fEb*7>fwXm+U{pQ-1$92(lN?Ca<*%WWrtL{ zb8psXD)>?LO+C2aNmW3P7(Qz8vb##(D_$=ua9$(5Z?dP~xo%C(7@U0S;e+4Wej@u^ z{8r>)9aUm|#aebpC#ibRdlosgqn{$bUwTwT&n?W|J0 z?S=|rDpljurEEiEc=l;Y!AG`}!$z|Rzsd&FnmLZvR?;;yV+D$C0OO9Qz;K|XdI?epI z4Kls1X~y8>Qx6~fminnGx{MWFB0UXVGKnsoMVGyzOES@g7@U0S;e+2&AK>eU8oqhp zYWOC%_-0R-;TuQs4fDLT_=XsqeCpwYpYi@2e~$uP!rq-L;#gG`o!Z2%!O79e)#YW` zb6IyN=icMP&iLO``Q$vn`d)5nsx~8k?F)#($)_GZ_^s{UdF`p`*}t__v$IwX4*m82 z0hfB2zDV?hOKQQDeP0i|`lOd-)Z6FZtQIQi7W2fwxb zlzmxnRHK?RRMv0jis!$5UQLwso_W3P>!gzd7I^5myZ(>~Hg)`X5`3{F1v z@WF4XpLc5vUHXeINko?xqDzqI;xD=s7F~$J$)_GZ_$~DTURr$fVv6CLSZBjGKYAIy z*)8XRHF7-=v&8TXF*y0u!v{a({W<;~1upxfE$*KwyPCSicIShe&D7&=DgE<@Ofk;` ztmj*m5O?D@Jw5Eb?&qu;r&9CT~qI@C3;*^wKJ?; z70~a8I#xKxHAlBxdPm(WL!Ewo4>--1=4odp9o;*~K@+;N(*eAN-8>=lFXR zDAn0z*|FIVsub+H#p#E0A@e-I`ptDb`}s<8eLxIOKJy+v_^s{UYs4;(t)cCe_k{_= zm)2jQGNq~Lwg1CaHS2Pb<6gIIbdSdKyZe95r!y=%vOaBD2c0B&(3%ZpoXr@VeCpwY z-`ak3qmXPvKCDqG<+^{n?C-hGXT4|MH%>mQrRT8W=6L{|eCpwY-%>yBqDu|YbOPSDg&s-0focrm~y|qUCzvp6mNY4wZywk#NnFoJWN9?liUOqR6?wUF4q7}VL=;h5y6}*jAFY*}uEWeZu+T8)9(ssfQ1KYy0_=_y2|hJXI}sA`~M|Gm#v~pZPDe`GeZ~Vbv@C=UUVS_ zC!c!w;I~$nKY4!?_#F!5O`Rj-gYQYy{N&-KyKktjUj8_~W!m~7=6Rr;ta}-%25)T9 z?Skr1zjS!g!&enCIP;!*_~5s;cb?Dlyq@RxJkRI({_lSNpH_4kC%Rk`UFyhrfc5;R zo1ses(S;bCeCpwY-&$S%IHKg(mS+v*JFd4Tns>-cT%Uf<5F zS*(1-|HRRzIZSJ+O3L2>?U26*$Rf`Lxj(>tJz{Y3 zsfQ1KOZ`+9UF<{`<}K@YZrK+Ei!OCV7h-VosfQ1KOMU#o|0wV~6&U~W zA%7Rd{+;`T+|MTlC!c!w;J4IIPtm1{oCjF1ONuV+-&2Y%p9dPc5QCFXJ$&$6>f;ao RM}a>I{88YK0{`C%{2!HUjTZm_ literal 0 HcmV?d00001 diff --git a/diffsims/tests/generators/test_simulation_generator.py b/diffsims/tests/generators/test_simulation_generator.py index b9fc2efa..d5de57c9 100644 --- a/diffsims/tests/generators/test_simulation_generator.py +++ b/diffsims/tests/generators/test_simulation_generator.py @@ -270,3 +270,68 @@ def test_invalid_scattering_params(): @pytest.mark.xfail(faises=NotImplementedError) def test_invalid_shape_model(): generator = SimulationGenerator(300, shape_factor_model="dracula") + + +def test_same_simulation_results(): + # This test is to ensure that the new SimulationGenerator produces the same + # results as the old DiffractionLibraryGenerator. Based on comments from + # viljarjf in https://github.com/pyxem/diffsims/pull/201 + # 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 + reciprocal_radius = 1.6768 + with_direct_beam = False + max_excitation_error = 0.1 + accelerating_voltage = 200 + shape = (128, 128) + sigma = 1.4 + generator_kwargs = { + "accelerating_voltage": accelerating_voltage, + "scattering_params": "lobato", + "precession_angle": 0, + "shape_factor_model": "lorentzian", + "approximate_precession": True, + "minimum_intensity": 1e-20, + } + # The euler angles are different, as orix.Phase enforces x||a, z||c* + # euler_angles_old = np.array([[0, 90, 120]]) + euler_angles_new = np.array([[0, 90, 90]]) + + # Old way. For creating the old data. + # struct_library = StructureLibrary(["Graphite"], [structure_matrix], [euler_angles_old]) + # diff_gen = DiffractionGenerator(**generator_kwargs) + # lib_gen = DiffractionLibraryGenerator(diff_gen) + # diff_lib = lib_gen.get_diffraction_library(struct_library, + # calibration=calibration, + # reciprocal_radius=reciprocal_radius, + # with_direct_beam=with_direct_beam, + # max_excitation_error=max_excitation_error, + # half_shape=shape[0] // 2, + # ) + # old_data = diff_lib["Graphite"]["simulations"][0].get_diffraction_pattern(shape=shape, sigma=sigma) + + # New + p = Phase("Graphite", point_group="6/mmm", structure=structure_matrix) + gen = SimulationGenerator(**generator_kwargs) + rot = Rotation.from_euler(euler_angles_new, degrees=True) + sim = gen.calculate_diffraction2d( + phase=p, + rotation=rot, + reciprocal_radius=reciprocal_radius, + max_excitation_error=max_excitation_error, + with_direct_beam=with_direct_beam, + ) + new_data = sim.get_diffraction_pattern( + shape=shape, sigma=sigma, calibration=calibration + ) + old_data = np.load( + "old_simulation.npy", + ) + np.testing.assert_allclose(new_data, old_data, atol=1e-8) From 3a27a7b7370bd353d39165d6e95059b564c03d95 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 3 May 2024 10:49:38 -0500 Subject: [PATCH 19/36] Documentation: Add intensity parameter --- diffsims/crystallography/diffracting_vector.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/diffsims/crystallography/diffracting_vector.py b/diffsims/crystallography/diffracting_vector.py index 41418af8..cb04e3f8 100644 --- a/diffsims/crystallography/diffracting_vector.py +++ b/diffsims/crystallography/diffracting_vector.py @@ -46,6 +46,9 @@ class DiffractingVector(ReciprocalLatticeVector): Indices of reciprocal lattice vector(s), often preferred over ``hkl`` in trigonal and hexagonal lattices. Default is ``None``. This, ``xyz``, or ``hkl`` is required. + intensity : numpy.ndarray, list, or tuple, optional + Intensity of the diffraction vector(s). Default is ``None``. + Examples -------- From 9bf940441bff960c8b11642b26cf79d66f13c8a0 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Fri, 3 May 2024 11:04:25 -0500 Subject: [PATCH 20/36] Manifest: Add numpy old file --- MANIFEST.in | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/MANIFEST.in b/MANIFEST.in index c7068aa3..d848c43c 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -7,4 +7,4 @@ include readthedocs.yaml include setup.cfg include setup.py -recursive-include doc Makefile make.bat *.rst *.py *.png +recursive-include doc Makefile make.bat *.rst *.py *.png *.npy From 9c216f820288258127152c02fe7872dbec860f34 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Sun, 5 May 2024 07:59:47 -0500 Subject: [PATCH 21/36] Refactor: Refactor DiffractingVector code --- diffsims/crystallography/__init__.py | 2 +- ...cting_vector.py => _diffracting_vector.py} | 26 ++++++++++++++++++- .../reciprocal_lattice_vector.py | 19 +++----------- diffsims/generators/simulation_generator.py | 10 +++---- .../test_diffracting_vector.py | 2 +- .../generators/test_simulation_generator.py | 8 +++--- 6 files changed, 41 insertions(+), 26 deletions(-) rename diffsims/crystallography/{diffracting_vector.py => _diffracting_vector.py} (81%) diff --git a/diffsims/crystallography/__init__.py b/diffsims/crystallography/__init__.py index ff42077a..e232e95d 100644 --- a/diffsims/crystallography/__init__.py +++ b/diffsims/crystallography/__init__.py @@ -27,7 +27,7 @@ get_hkl, ) from diffsims.crystallography.reciprocal_lattice_vector import ReciprocalLatticeVector -from diffsims.crystallography.diffracting_vector import DiffractingVector +from diffsims.crystallography._diffracting_vector import DiffractingVector __all__ = [ "get_equivalent_hkl", diff --git a/diffsims/crystallography/diffracting_vector.py b/diffsims/crystallography/_diffracting_vector.py similarity index 81% rename from diffsims/crystallography/diffracting_vector.py rename to diffsims/crystallography/_diffracting_vector.py index cb04e3f8..e2ad96b6 100644 --- a/diffsims/crystallography/diffracting_vector.py +++ b/diffsims/crystallography/_diffracting_vector.py @@ -18,6 +18,7 @@ from diffsims.crystallography import ReciprocalLatticeVector import numpy as np +from orix.vector.miller import _transform_space class DiffractingVector(ReciprocalLatticeVector): @@ -31,6 +32,9 @@ class DiffractingVector(ReciprocalLatticeVector): electron diffraction based on the intersection of the Ewald sphere with the reciprocal lattice. + 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. + Parameters ---------- phase : orix.crystal_map.Phase @@ -48,6 +52,10 @@ class DiffractingVector(ReciprocalLatticeVector): This, ``xyz``, or ``hkl`` is required. intensity : numpy.ndarray, list, or tuple, optional Intensity of the diffraction vector(s). Default is ``None``. + rotation : 3x3 numpy.ndarray, list, or tuple, optional + Rotation matrix previously applied to the reciprocal lattice vector(s) and the + lattice of the phase. Default is ``None`` which corresponds to the + identity matrix. Examples @@ -71,7 +79,9 @@ class DiffractingVector(ReciprocalLatticeVector): """ - def __init__(self, phase, xyz=None, hkl=None, hkil=None, intensity=None): + def __init__( + self, phase, xyz=None, hkl=None, hkil=None, intensity=None, rotation=None + ): super().__init__(phase, xyz=xyz, hkl=hkl, hkil=hkil) if intensity is None: self._intensity = np.full(self.shape, np.nan) @@ -79,6 +89,7 @@ def __init__(self, phase, xyz=None, hkl=None, hkil=None, intensity=None): raise ValueError("Length of intensity array must match number of vectors") else: self._intensity = np.array(intensity) + self._rotation = rotation def __getitem__(self, key): dv_new = super().__getitem__(key) @@ -113,6 +124,19 @@ def intensity(self, value): raise ValueError("Length of intensity array must match number of vectors") self._intensity = np.array(value) + @property + def lattice_aligned_data(self): + if self._rotation is None: + return self.data + else: + return np.matmul(self.data, self._rotation.T) + + @property + def hkl(self): + return _transform_space( + self.lattice_aligned_data, "c", "r", self.phase.structure.lattice + ) + def calculate_structure_factor(self): raise NotImplementedError( "Structure factor calculation not implemented for DiffractionVector. " diff --git a/diffsims/crystallography/reciprocal_lattice_vector.py b/diffsims/crystallography/reciprocal_lattice_vector.py index d52f4e24..222a9a25 100644 --- a/diffsims/crystallography/reciprocal_lattice_vector.py +++ b/diffsims/crystallography/reciprocal_lattice_vector.py @@ -78,6 +78,10 @@ class ReciprocalLatticeVector(Vector3d): Indices of reciprocal lattice vector(s), often preferred over ``hkl`` in trigonal and hexagonal lattices. Default is ``None``. This, ``xyz``, or ``hkl`` is required. + lattice_rotation : 3x3 numpy.ndarray, list, or tuple, optional + Rotation matrix applied to the reciprocal lattice vector(s) and the + lattice of the phase. Default is ``None`` which corresponds to the + identity matrix. Examples -------- @@ -503,21 +507,6 @@ def scattering_parameter(self): return 0.5 * self.gspacing - def rotate_from_matrix(self, rotation_matrix): - """Rotate the reciprocal lattice vectors using a (3x3) rotation matrix. - - This method creates a new instance of :class:`ReciprocalLatticeVector` - with the rotated vectors. - - Parameters - ---------- - rotation_matrix : numpy.ndarray - (3, 3) rotation matrix. - """ - return ReciprocalLatticeVector( - phase=self.phase, xyz=np.matmul(rotation_matrix.T, self.data.T).T - ) - @property def structure_factor(self): r"""Kinematical structure factors :math:`F`. diff --git a/diffsims/generators/simulation_generator.py b/diffsims/generators/simulation_generator.py index 0d7c3e2a..136be512 100644 --- a/diffsims/generators/simulation_generator.py +++ b/diffsims/generators/simulation_generator.py @@ -361,13 +361,11 @@ def get_intersecting_reflections( control. If not set will be set equal to max_excitation_error. """ initial_hkl = recip.hkl - rotated_vectors = recip.rotate_from_matrix(rot) + rotated_vectors = np.matmul(rot.T, recip.data.T).T if with_direct_beam: - rotated_vectors = np.vstack([rotated_vectors.data, [0, 0, 0]]) + rotated_vectors = np.vstack([rotated_vectors, [0, 0, 0]]) initial_hkl = np.vstack([initial_hkl, [0, 0, 0]]) - else: - rotated_vectors = rotated_vectors.data # Identify the excitation errors of all points (distance from point to Ewald sphere) r_sphere = 1 / wavelength r_spot = np.sqrt(np.sum(np.square(rotated_vectors[:, :2]), axis=1)) @@ -394,7 +392,9 @@ def get_intersecting_reflections( # select these reflections intersected_vectors = rotated_vectors[intersection] intersected_vectors = DiffractingVector( - phase=recip.phase, xyz=intersected_vectors + phase=recip.phase, + xyz=intersected_vectors, + rotation=rot, ) excitation_error = excitation_error[intersection] r_spot = r_spot[intersection] diff --git a/diffsims/tests/crystallography/test_diffracting_vector.py b/diffsims/tests/crystallography/test_diffracting_vector.py index adf2c4dd..29c810ea 100644 --- a/diffsims/tests/crystallography/test_diffracting_vector.py +++ b/diffsims/tests/crystallography/test_diffracting_vector.py @@ -1,4 +1,4 @@ -from diffsims.crystallography import DiffractingVector +from diffsims.crystallography import DiffractingVector, ReciprocalLatticeVector import pytest import numpy as np diff --git a/diffsims/tests/generators/test_simulation_generator.py b/diffsims/tests/generators/test_simulation_generator.py index d5de57c9..df16fbd7 100644 --- a/diffsims/tests/generators/test_simulation_generator.py +++ b/diffsims/tests/generators/test_simulation_generator.py @@ -18,6 +18,7 @@ import numpy as np import pytest +from pathlib import Path import diffpy.structure from orix.crystal_map import Phase @@ -35,6 +36,9 @@ from diffsims.simulations import Simulation1D from diffsims.utils.sim_utils import is_lattice_hexagonal +TEST_DATA_DIR = Path(__file__).parent +FILE1 = TEST_DATA_DIR / "old_simulation.npy" + @pytest.fixture(params=[(300)]) def diffraction_calculator(request): @@ -331,7 +335,5 @@ def test_same_simulation_results(): new_data = sim.get_diffraction_pattern( shape=shape, sigma=sigma, calibration=calibration ) - old_data = np.load( - "old_simulation.npy", - ) + old_data = np.load(FILE1) np.testing.assert_allclose(new_data, old_data, atol=1e-8) From 2ada5679c51eba0b13f14a4b485ce0b16079bf6d Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Sun, 5 May 2024 08:12:32 -0500 Subject: [PATCH 22/36] Refactor: Refactor Rotation implementation --- diffsims/crystallography/_diffracting_vector.py | 4 ++-- diffsims/generators/simulation_generator.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/diffsims/crystallography/_diffracting_vector.py b/diffsims/crystallography/_diffracting_vector.py index e2ad96b6..e2ee4250 100644 --- a/diffsims/crystallography/_diffracting_vector.py +++ b/diffsims/crystallography/_diffracting_vector.py @@ -52,7 +52,7 @@ class DiffractingVector(ReciprocalLatticeVector): This, ``xyz``, or ``hkl`` is required. intensity : numpy.ndarray, list, or tuple, optional Intensity of the diffraction vector(s). Default is ``None``. - rotation : 3x3 numpy.ndarray, list, or tuple, optional + rotation : orix.quaternion.Rotation, optional Rotation matrix previously applied to the reciprocal lattice vector(s) and the lattice of the phase. Default is ``None`` which corresponds to the identity matrix. @@ -129,7 +129,7 @@ def lattice_aligned_data(self): if self._rotation is None: return self.data else: - return np.matmul(self.data, self._rotation.T) + return np.matmul(self.data, self._rotation.to_matrix()[0]) @property def hkl(self): diff --git a/diffsims/generators/simulation_generator.py b/diffsims/generators/simulation_generator.py index 136be512..eefd02b3 100644 --- a/diffsims/generators/simulation_generator.py +++ b/diffsims/generators/simulation_generator.py @@ -196,7 +196,7 @@ def calculate_diffraction2d( include_zero_vector=with_direct_beam, ) phase_vectors = [] - for rot in rotate.to_matrix(): + for rot in rotate: # Calculate the reciprocal lattice vectors that intersect the Ewald sphere. ( intersected_vectors, @@ -361,7 +361,7 @@ def get_intersecting_reflections( control. If not set will be set equal to max_excitation_error. """ initial_hkl = recip.hkl - rotated_vectors = np.matmul(rot.T, recip.data.T).T + rotated_vectors = (~rot * recip.to_miller()).data if with_direct_beam: rotated_vectors = np.vstack([rotated_vectors, [0, 0, 0]]) From b598a5a40420461c31d442aa61b31575d3347cde Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Sun, 5 May 2024 08:20:15 -0500 Subject: [PATCH 23/36] BugFix: Fix HKL computation --- diffsims/crystallography/_diffracting_vector.py | 2 +- .../tests/crystallography/test_diffracting_vector.py | 9 +++++++++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/diffsims/crystallography/_diffracting_vector.py b/diffsims/crystallography/_diffracting_vector.py index e2ee4250..c29b32e5 100644 --- a/diffsims/crystallography/_diffracting_vector.py +++ b/diffsims/crystallography/_diffracting_vector.py @@ -129,7 +129,7 @@ def lattice_aligned_data(self): if self._rotation is None: return self.data else: - return np.matmul(self.data, self._rotation.to_matrix()[0]) + return np.matmul(self.data, self._rotation.to_matrix()[0].T) @property def hkl(self): diff --git a/diffsims/tests/crystallography/test_diffracting_vector.py b/diffsims/tests/crystallography/test_diffracting_vector.py index 29c810ea..87e6adb9 100644 --- a/diffsims/tests/crystallography/test_diffracting_vector.py +++ b/diffsims/tests/crystallography/test_diffracting_vector.py @@ -1,4 +1,6 @@ from diffsims.crystallography import DiffractingVector, ReciprocalLatticeVector +from orix.quaternion import Rotation + import pytest import numpy as np @@ -40,3 +42,10 @@ def test_structure_factor(self, ferrite_phase): rlv = DiffractingVector.from_min_dspacing(ferrite_phase, 1.5) with pytest.raises(NotImplementedError): rlv.calculate_structure_factor() + + def test_hl(self, ferrite_phase): + rlv = ReciprocalLatticeVector(ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]]) + rot = Rotation.from_euler([90, 90, 0], degrees=True) + rotated_vectors = (~rot * rlv.to_miller()).data + dv = DiffractingVector(ferrite_phase, xyz=rotated_vectors, rotation=rot) + assert np.allclose(rlv.hkl, dv.hkl) From e0815d54936c63cef40781900107384fee196f3c Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Sun, 5 May 2024 12:02:57 -0500 Subject: [PATCH 24/36] Refactor: Fix MANIFEST.in --- MANIFEST.in | 2 +- setup.cfg | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/MANIFEST.in b/MANIFEST.in index d848c43c..c7068aa3 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -7,4 +7,4 @@ include readthedocs.yaml include setup.cfg include setup.py -recursive-include doc Makefile make.bat *.rst *.py *.png *.npy +recursive-include doc Makefile make.bat *.rst *.py *.png diff --git a/setup.cfg b/setup.cfg index 46204e6d..b804aa93 100644 --- a/setup.cfg +++ b/setup.cfg @@ -28,4 +28,5 @@ known_excludes = .git/** doc/build/** htmlcov/** - *.code-workspace \ No newline at end of file + *.code-workspace + **/*.npy \ No newline at end of file From dff72bbb3957432721d3827bab71af0b7dac55ef Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Mon, 6 May 2024 10:34:44 -0500 Subject: [PATCH 25/36] Plotting: Add interactive plotting --- diffsims/simulations/simulation2d.py | 217 +++++++++++++++++++++------ 1 file changed, 174 insertions(+), 43 deletions(-) diff --git a/diffsims/simulations/simulation2d.py b/diffsims/simulations/simulation2d.py index 7df1c953..dfafc440 100644 --- a/diffsims/simulations/simulation2d.py +++ b/diffsims/simulations/simulation2d.py @@ -21,7 +21,7 @@ import numpy as np import matplotlib.pyplot as plt -import math +from matplotlib.widgets import Slider from orix.crystal_map import Phase from orix.quaternion import Rotation from orix.vector import Vector3d @@ -192,6 +192,8 @@ def __init__( # for slicing a simulation self.iphase = PhaseGetter(self) self.irot = RotationGetter(self) + self._rotation_slider = None + self._phase_slider = None def get_simulation(self, item): """Return the rotation and the phase index of the simulation""" @@ -465,6 +467,60 @@ def plot_rotations(self, beam_direction: Vector3d = Vector3d.zvector()): _plot = fig.axes[0] _plot.set_title("Rotations" + self.current_phase.name) + def _get_spots( + self, + in_plane_angle, + direct_beam_position, + mirrored, + units, + calibration, + include_direct_beam, + ): + """Returns the spots of the current phase and rotation for plotting""" + coords = self._get_transformed_coordinates( + in_plane_angle, + direct_beam_position, + mirrored, + units=units, + calibration=calibration, + ) + if include_direct_beam: + spots = coords.data[:, :2] + spots = np.concatenate((spots, np.array([direct_beam_position]))) + intensity = np.concatenate((coords.intensity, np.array([1]))) + else: + spots = coords.data[:, :2] + intensity = coords.intensity + return spots, intensity, coords + + def _get_labels(self, coords, intensity, min_label_intensity, xlim, ylim): + condition = ( + (coords.data[:, 0] > min(xlim)) + & (coords.data[:, 0] < max(xlim)) + & (coords.data[:, 1] > min(ylim)) + & (coords.data[:, 1] < max(ylim)) + ) + in_range_coords = coords.data[condition] + millers = np.round( + np.matmul( + np.matmul(in_range_coords, self.get_current_rotation_matrix().T), + coords.phase.structure.lattice.base.T, + ) + ).astype(np.int16) + labels = [] + for miller, coordinate, inten in zip(millers, in_range_coords, intensity): + if np.isnan(inten) or inten > min_label_intensity: + label = "(" + for index in miller: + if index < 0: + label += r"$\bar{" + str(abs(index)) + r"}$" + else: + label += str(abs(index)) + label += " " + label = label[:-1] + ")" + labels.append((coordinate, label)) + return labels + def plot( self, size_factor=1, @@ -479,6 +535,7 @@ def plot( include_direct_beam=True, calibration=0.1, ax=None, + interactive=False, **kwargs, ): """A quick-plot function for a simulation of spots @@ -522,27 +579,23 @@ def plot( ----- spot size scales with the square root of the intensity. """ + if label_formatting is None: label_formatting = {} if direct_beam_position is None: direct_beam_position = (0, 0) if ax is None: - _, ax = plt.subplots() + fig, ax = plt.subplots() ax.set_aspect("equal") - coords = self._get_transformed_coordinates( - in_plane_angle, - direct_beam_position, - mirrored, + + spots, intensity, coords = self._get_spots( + in_plane_angle=in_plane_angle, + direct_beam_position=direct_beam_position, + mirrored=mirrored, units=units, calibration=calibration, + include_direct_beam=include_direct_beam, ) - if include_direct_beam: - spots = coords.data[:, :2] - spots = np.concatenate((spots, np.array([direct_beam_position]))) - intensity = np.concatenate((coords.intensity, np.array([1]))) - else: - spots = coords.data[:, :2] - intensity = coords.intensity sp = ax.scatter( spots[:, 0], spots[:, 1], @@ -551,25 +604,13 @@ def plot( ) ax.set_xlim(-self.reciporical_radius, self.reciporical_radius) ax.set_ylim(-self.reciporical_radius, self.reciporical_radius) - + texts = [] if show_labels: - millers = np.round( - np.matmul( - np.matmul(coords.data, self.get_current_rotation_matrix().T), - coords.phase.structure.lattice.base.T, - ) - ).astype(np.int16) - # only label the points inside the axes xlim = ax.get_xlim() ylim = ax.get_ylim() - condition = ( - (coords.data[:, 0] > min(xlim)) - & (coords.data[:, 0] < max(xlim)) - & (coords.data[:, 1] > min(ylim)) - & (coords.data[:, 1] < max(ylim)) + labels = self._get_labels( + coords, intensity, min_label_intensity, xlim, ylim ) - millers = millers[condition] - coords = coords.data[condition] # default alignment options if ( "ha" not in label_offset @@ -578,28 +619,118 @@ def plot( label_formatting["ha"] = "center" if "va" not in label_offset and "verticalalignment" not in label_formatting: label_formatting["va"] = "center" - for miller, coordinate, inten in zip(millers, coords, intensity): - if np.isnan(inten) or inten > min_label_intensity: - label = "(" - for index in miller: - if index < 0: - label += r"$\bar{" + str(abs(index)) + r"}$" - else: - label += str(abs(index)) - label += " " - label = label[:-1] + ")" + for coordinate, label in labels: + texts.append( ax.text( coordinate[0] + label_offset[0], coordinate[1] + label_offset[1], label, **label_formatting, ) - if units == "real": - ax.set_xlabel(r"$\AA^{-1}$") - ax.set_ylabel(r"$\AA^{-1}$") + ) + if units == "real": + ax.set_xlabel(r"$\AA^{-1}$") + ax.set_ylabel(r"$\AA^{-1}$") + else: + ax.set_xlabel("pixels") + ax.set_ylabel("pixels") + if interactive and self.has_multiple_rotations or self.has_multiple_phases: + axrot = fig.add_axes([0.5, 0.05, 0.4, 0.03]) + axphase = fig.add_axes([0.1, 0.05, 0.2, 0.03]) + + fig.subplots_adjust(left=0.25, bottom=0.25) + if self.has_multiple_rotations and self.has_multiple_phases: + max_rot = np.max([r.size for r in self.rotations]) + rotation_slider = Slider( + ax=axrot, + label="Rotation", + valmin=0, + valmax=max_rot - 1, + valinit=self.rotation_index, + valstep=1, + orientation="horizontal", + ) + phase_slider = Slider( + ax=axphase, + label="Phase ", + valmin=0, + valmax=self.phases.size - 1, + valinit=self.phase_index, + valstep=1, + orientation="horizontal", + ) + elif self.has_multiple_rotations: + rotation_slider = Slider( + ax=axrot, + label="Rotation", + valmin=0, + valmax=self.rotations.size - 1, + valinit=self.rotation_index, + valstep=1, + orientation="horizontal", + ) + phase_slider = None else: - ax.set_xlabel("pixels") - ax.set_ylabel("pixels") + rotation_slider = None + phase_slider = Slider( + ax=axrot, + label="Phase", + valmin=0, + valmax=self.phases.size - 1, + valinit=self.phase_index, + valstep=1, + orientation="horizontal", + ) + self._rotation_slider = rotation_slider + self._phase_slider = phase_slider + + def update(val): + if self.has_multiple_rotations and self.has_multiple_phases: + self.rotation_index = int(rotation_slider.val) + self.phase_index = int(phase_slider.val) + self._rotation_slider.valmax = ( + self.rotations[self.phase_index].size - 1 + ) + elif self.has_multiple_rotations: + self.rotation_index = int(rotation_slider.val) + else: + self.phase_index = int(phase_slider.val) + spots, intensity, coords = self._get_spots( + in_plane_angle, + direct_beam_position, + mirrored, + units, + calibration, + include_direct_beam, + ) + sp.set( + offsets=spots, + sizes=size_factor * np.sqrt(intensity), + ) + for t in texts: + t.remove() + texts.clear() + if show_labels: + xlim = ax.get_xlim() + ylim = ax.get_ylim() + labels = self._get_labels( + coords, intensity, min_label_intensity, xlim, ylim + ) + for coordinate, label in labels: + texts.append( + ax.text( + coordinate[0] + label_offset[0], + coordinate[1] + label_offset[1], + label, + **label_formatting, + ) + ) + fig.canvas.draw_idle() + + if self._rotation_slider is not None: + self._rotation_slider.on_changed(update) + if self._phase_slider is not None: + self._phase_slider.on_changed(update) return ax, sp From 8b92dce82a42919632cd48a1269dfa78efe49509 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Mon, 6 May 2024 10:57:19 -0500 Subject: [PATCH 26/36] Testing: no cover interactive plotting --- diffsims/simulations/simulation2d.py | 20 ++++++------------- .../tests/simulations/test_simulations2d.py | 4 +++- 2 files changed, 9 insertions(+), 15 deletions(-) diff --git a/diffsims/simulations/simulation2d.py b/diffsims/simulations/simulation2d.py index dfafc440..2758053d 100644 --- a/diffsims/simulations/simulation2d.py +++ b/diffsims/simulations/simulation2d.py @@ -634,12 +634,14 @@ def plot( else: ax.set_xlabel("pixels") ax.set_ylabel("pixels") - if interactive and self.has_multiple_rotations or self.has_multiple_phases: + if ( + interactive and self.has_multiple_rotations or self.has_multiple_phases + ): # pragma: no cover axrot = fig.add_axes([0.5, 0.05, 0.4, 0.03]) axphase = fig.add_axes([0.1, 0.05, 0.2, 0.03]) fig.subplots_adjust(left=0.25, bottom=0.25) - if self.has_multiple_rotations and self.has_multiple_phases: + if self.has_multiple_phases: max_rot = np.max([r.size for r in self.rotations]) rotation_slider = Slider( ax=axrot, @@ -659,7 +661,7 @@ def plot( valstep=1, orientation="horizontal", ) - elif self.has_multiple_rotations: + else: # self.has_multiple_rotations: rotation_slider = Slider( ax=axrot, label="Rotation", @@ -670,17 +672,6 @@ def plot( orientation="horizontal", ) phase_slider = None - else: - rotation_slider = None - phase_slider = Slider( - ax=axrot, - label="Phase", - valmin=0, - valmax=self.phases.size - 1, - valinit=self.phase_index, - valstep=1, - orientation="horizontal", - ) self._rotation_slider = rotation_slider self._phase_slider = phase_slider @@ -717,6 +708,7 @@ def update(val): coords, intensity, min_label_intensity, xlim, ylim ) for coordinate, label in labels: + # this could be faster using a TextCollection when available in matplotlib texts.append( ax.text( coordinate[0] + label_offset[0], diff --git a/diffsims/tests/simulations/test_simulations2d.py b/diffsims/tests/simulations/test_simulations2d.py index f2757319..0f4a99ff 100644 --- a/diffsims/tests/simulations/test_simulations2d.py +++ b/diffsims/tests/simulations/test_simulations2d.py @@ -193,7 +193,9 @@ def multi_simulation(self, al_phase): gen = SimulationGenerator(accelerating_voltage=200) rot = Rotation.from_axes_angles([1, 0, 0], (0, 15, 30, 45), degrees=True) coords = DiffractingVector( - phase=al_phase, xyz=[[1, 0, 0], [0, 1, 0], [1, 1, 0], [1, 1, 1]] + phase=al_phase, + xyz=[[1, 0, 0], [0, 1, 0], [1, 1, 0], [1, 1, 1]], + intensity=[1, 2, 3, 4], ) vectors = [coords, coords, coords, coords] From f3d9b957122cf41fa3776577d28f61fcacd59ab1 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Mon, 6 May 2024 11:00:24 -0500 Subject: [PATCH 27/36] Documentation: Add docstrings for interactive plotting --- diffsims/simulations/simulation2d.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/diffsims/simulations/simulation2d.py b/diffsims/simulations/simulation2d.py index 2758053d..44993ef8 100644 --- a/diffsims/simulations/simulation2d.py +++ b/diffsims/simulations/simulation2d.py @@ -568,6 +568,9 @@ def plot( whether to include the direct beam in the plot ax : matplotlib Axes, optional axes on which to draw the pattern. If `None`, a new axis is created + interactive : bool, optional + Whether to add sliders for selecting the rotation and phase. This + is an experimental feature and will evolve/change in the future. **kwargs : passed to ax.scatter() method From 82b33339508d5ca22afe179be8b9e4916fa5403e Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Mon, 6 May 2024 12:56:57 -0500 Subject: [PATCH 28/36] Documentation: Remove lattice rotation --- diffsims/crystallography/reciprocal_lattice_vector.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/diffsims/crystallography/reciprocal_lattice_vector.py b/diffsims/crystallography/reciprocal_lattice_vector.py index 222a9a25..b424c0d7 100644 --- a/diffsims/crystallography/reciprocal_lattice_vector.py +++ b/diffsims/crystallography/reciprocal_lattice_vector.py @@ -78,10 +78,6 @@ class ReciprocalLatticeVector(Vector3d): Indices of reciprocal lattice vector(s), often preferred over ``hkl`` in trigonal and hexagonal lattices. Default is ``None``. This, ``xyz``, or ``hkl`` is required. - lattice_rotation : 3x3 numpy.ndarray, list, or tuple, optional - Rotation matrix applied to the reciprocal lattice vector(s) and the - lattice of the phase. Default is ``None`` which corresponds to the - identity matrix. Examples -------- From c39e35775e39c550aa78396d82b8c8f3717de091 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Tue, 7 May 2024 09:07:38 -0500 Subject: [PATCH 29/36] Bugfix: Fix polar conversion --- .../crystallography/_diffracting_vector.py | 9 ++++++++ diffsims/simulations/simulation2d.py | 22 +++++++++---------- .../test_diffracting_vector.py | 6 +++++ 3 files changed, 26 insertions(+), 11 deletions(-) diff --git a/diffsims/crystallography/_diffracting_vector.py b/diffsims/crystallography/_diffracting_vector.py index c29b32e5..46c412ab 100644 --- a/diffsims/crystallography/_diffracting_vector.py +++ b/diffsims/crystallography/_diffracting_vector.py @@ -142,3 +142,12 @@ def calculate_structure_factor(self): "Structure factor calculation not implemented for DiffractionVector. " "Use ReciprocalLatticeVector instead." ) + + 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) + theta = np.arctan2( + self.data[:, 1], + self.data[:, 0], + ) + return r, theta diff --git a/diffsims/simulations/simulation2d.py b/diffsims/simulations/simulation2d.py index 44993ef8..9852c5c4 100644 --- a/diffsims/simulations/simulation2d.py +++ b/diffsims/simulations/simulation2d.py @@ -26,7 +26,7 @@ from orix.quaternion import Rotation from orix.vector import Vector3d -from diffsims.crystallography.reciprocal_lattice_vector import ReciprocalLatticeVector +from diffsims.crystallography._diffracting_vector import DiffractingVector from diffsims.pattern.detector_functions import add_shot_and_point_spread # to avoid circular imports @@ -114,9 +114,9 @@ def __init__( self, phases: Sequence[Phase], coordinates: Union[ - ReciprocalLatticeVector, - Sequence[ReciprocalLatticeVector], - Sequence[Sequence[ReciprocalLatticeVector]], + DiffractingVector, + Sequence[DiffractingVector], + Sequence[Sequence[DiffractingVector]], ], rotations: Union[Rotation, Sequence[Rotation]], simulation_generator: "SimulationGenerator", @@ -128,9 +128,9 @@ def __init__( Parameters ---------- coordinates - The list of ReciprocalLatticeVector objects for each phase and rotation. If there - are multiple phases, then this should be a list of lists of ReciprocalLatticeVector objects. - If there is only one phase, then this should be a list of ReciprocalLatticeVector objects. + The list of DiffractingVector objects for each phase and rotation. If there + are multiple phases, then this should be a list of lists of DiffractingVector objects. + If there is only one phase, then this should be a list of DiffractingVector objects. rotations The list of Rotation objects for each phase. If there are multiple phases, then this should be a list of Rotation objects. If there is only one phase, then this should be a single @@ -144,9 +144,9 @@ def __init__( """ # Basic data if isinstance(rotations, Rotation) and rotations.size == 1: - if not isinstance(coordinates, ReciprocalLatticeVector): + if not isinstance(coordinates, DiffractingVector): raise ValueError( - "If there is only one rotation, then the coordinates must be a ReciprocalLatticeVector object" + "If there is only one rotation, then the coordinates must be a DiffractingVector object" ) elif isinstance(rotations, Rotation): coordinates = np.array(coordinates, dtype=object) @@ -166,7 +166,7 @@ def __init__( ) for r, c in zip(rotations, coordinates): - if isinstance(c, ReciprocalLatticeVector): + if isinstance(c, DiffractingVector): c = np.array( [ c, @@ -325,7 +325,7 @@ def polar_flatten_simulations(self, radial_axes=None, azimuthal_axes=None): 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() + r, t = v.to_flat_polar() if radial_axes is not None and azimuthal_axes is not None: r = get_closest(radial_axes, r) t = get_closest(azimuthal_axes, t) diff --git a/diffsims/tests/crystallography/test_diffracting_vector.py b/diffsims/tests/crystallography/test_diffracting_vector.py index 87e6adb9..c08c01e7 100644 --- a/diffsims/tests/crystallography/test_diffracting_vector.py +++ b/diffsims/tests/crystallography/test_diffracting_vector.py @@ -49,3 +49,9 @@ def test_hl(self, ferrite_phase): rotated_vectors = (~rot * rlv.to_miller()).data dv = DiffractingVector(ferrite_phase, xyz=rotated_vectors, rotation=rot) assert np.allclose(rlv.hkl, dv.hkl) + + def test_flat_polar(self, ferrite_phase): + dv = DiffractingVector(ferrite_phase, xyz=[[1, 1, 1], [0.5, -0.5, 0]]) + r, t = dv.to_flat_polar() + assert np.allclose(r, [np.sqrt(2), 0.70710678]) + assert np.allclose(t, [np.pi / 4, -np.pi / 4]) From 05b373691bc4651d4497d681322748575633b3ca Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Tue, 7 May 2024 21:33:01 -0500 Subject: [PATCH 30/36] Bugfix: Fix Rotation with RLV --- .../crystallography/_diffracting_vector.py | 17 +----- .../reciprocal_lattice_vector.py | 52 ++++++++++++++++--- diffsims/generators/simulation_generator.py | 5 +- diffsims/simulations/simulation2d.py | 17 +++--- .../test_diffracting_vector.py | 7 ++- .../test_reciprocal_lattice_vector.py | 44 ++++++++++++++-- .../tests/simulations/test_simulations2d.py | 4 ++ 7 files changed, 106 insertions(+), 40 deletions(-) diff --git a/diffsims/crystallography/_diffracting_vector.py b/diffsims/crystallography/_diffracting_vector.py index 46c412ab..1732a297 100644 --- a/diffsims/crystallography/_diffracting_vector.py +++ b/diffsims/crystallography/_diffracting_vector.py @@ -18,7 +18,6 @@ from diffsims.crystallography import ReciprocalLatticeVector import numpy as np -from orix.vector.miller import _transform_space class DiffractingVector(ReciprocalLatticeVector): @@ -82,14 +81,13 @@ class DiffractingVector(ReciprocalLatticeVector): def __init__( self, phase, xyz=None, hkl=None, hkil=None, intensity=None, rotation=None ): - super().__init__(phase, xyz=xyz, hkl=hkl, hkil=hkil) + super().__init__(phase, xyz=xyz, hkl=hkl, hkil=hkil, rotation=rotation) if intensity is None: self._intensity = np.full(self.shape, np.nan) elif len(intensity) != self.size: raise ValueError("Length of intensity array must match number of vectors") else: self._intensity = np.array(intensity) - self._rotation = rotation def __getitem__(self, key): dv_new = super().__getitem__(key) @@ -124,19 +122,6 @@ def intensity(self, value): raise ValueError("Length of intensity array must match number of vectors") self._intensity = np.array(value) - @property - def lattice_aligned_data(self): - if self._rotation is None: - return self.data - else: - return np.matmul(self.data, self._rotation.to_matrix()[0].T) - - @property - def hkl(self): - return _transform_space( - self.lattice_aligned_data, "c", "r", self.phase.structure.lattice - ) - def calculate_structure_factor(self): raise NotImplementedError( "Structure factor calculation not implemented for DiffractionVector. " diff --git a/diffsims/crystallography/reciprocal_lattice_vector.py b/diffsims/crystallography/reciprocal_lattice_vector.py index b424c0d7..7a734ae9 100644 --- a/diffsims/crystallography/reciprocal_lattice_vector.py +++ b/diffsims/crystallography/reciprocal_lattice_vector.py @@ -25,6 +25,7 @@ import numba as nb import numpy as np from orix.vector import Miller, Vector3d +from orix.quaternion import Rotation from orix.vector.miller import ( _check_hkil, _get_highest_hkl, @@ -78,6 +79,8 @@ class ReciprocalLatticeVector(Vector3d): Indices of reciprocal lattice vector(s), often preferred over ``hkl`` in trigonal and hexagonal lattices. Default is ``None``. This, ``xyz``, or ``hkl`` is required. + rotation : orix.quaternion.Rotation, optional + Rotation to apply to the vectors. Default is ``None``. Examples -------- @@ -100,7 +103,7 @@ class ReciprocalLatticeVector(Vector3d): """ - def __init__(self, phase, xyz=None, hkl=None, hkil=None): + def __init__(self, phase, xyz=None, hkl=None, hkil=None, rotation=None): self.phase = phase self._raise_if_no_point_group() @@ -120,13 +123,24 @@ def __init__(self, phase, xyz=None, hkl=None, hkil=None): self._coordinate_format = "hkl" xyz = _transform_space(hkl, "r", "c", phase.structure.lattice) super().__init__(xyz) - + if rotation is None: + self._rotation = Rotation.from_euler([0, 0, 0], degrees=True) + elif rotation.size == 1 or self.size == rotation.size: + self._rotation = rotation + else: + raise ValueError( + "Rotation must be a single Rotation or have the same size as the vectors" + ) self._theta = np.full(self.shape, np.nan) self._structure_factor = np.full(self.shape, np.nan, dtype="complex128") def __getitem__(self, key): new_data = self.data[key] - rlv_new = self.__class__(self.phase, xyz=new_data) + if self.rotation.size == 1: + rotation = self.rotation + else: + rotation = self.rotation[key] + rlv_new = self.__class__(self.phase, xyz=new_data, rotation=rotation) if np.isnan(self.structure_factor).all(): rlv_new._structure_factor = np.full( @@ -151,6 +165,28 @@ def __repr__(self): phase_name = self.phase.name return f"{name} {shape}, {phase_name} ({symmetry})\n" f"{data}" + def get_kwargs(self, new_data, *args, **kwargs): + if "rotation" in kwargs: + mult_rotation = kwargs.pop("rotation") + new_rotation = mult_rotation * self.rotation + else: + new_rotation = self.rotation + return dict(phase=self.phase, xyz=new_data, rotation=new_rotation) + + @property + def rotation(self): + """ + Rotation applied to the vectors, for calculating hkl + """ + return self._rotation + + @property + def unrotated_data(self): + """ + Un-rotated data for calculating hkl + """ + return (~self.rotation * self).data + @property def hkl(self): """Miller indices. @@ -173,7 +209,9 @@ def hkl(self): """ - return _transform_space(self.data, "c", "r", self.phase.structure.lattice) + return _transform_space( + self.unrotated_data, "c", "r", self.phase.structure.lattice + ) @property def hkil(self): @@ -1070,9 +1108,11 @@ def from_highest_hkl(cls, phase, hkl, include_zero_vector=False): """ idx = _get_indices_from_highest(highest_indices=hkl) + new = cls(phase, hkl=idx).unique() if include_zero_vector: - idx = np.vstack((idx, np.zeros(3, dtype=int))) - return cls(phase, hkl=idx).unique() + new_data = np.vstack((new.hkl, np.zeros(3, dtype=int))) + new = ReciprocalLatticeVector(phase, hkl=new_data) + return new @classmethod def from_min_dspacing(cls, phase, min_dspacing=0.7, include_zero_vector=False): diff --git a/diffsims/generators/simulation_generator.py b/diffsims/generators/simulation_generator.py index eefd02b3..3902edb2 100644 --- a/diffsims/generators/simulation_generator.py +++ b/diffsims/generators/simulation_generator.py @@ -361,7 +361,8 @@ def get_intersecting_reflections( control. If not set will be set equal to max_excitation_error. """ initial_hkl = recip.hkl - rotated_vectors = (~rot * recip.to_miller()).data + + rotated_vectors = (~rot * recip).data if with_direct_beam: rotated_vectors = np.vstack([rotated_vectors, [0, 0, 0]]) @@ -394,7 +395,7 @@ def get_intersecting_reflections( intersected_vectors = DiffractingVector( phase=recip.phase, xyz=intersected_vectors, - rotation=rot, + rotation=~rot, ) excitation_error = excitation_error[intersection] r_spot = r_spot[intersection] diff --git a/diffsims/simulations/simulation2d.py b/diffsims/simulations/simulation2d.py index 9852c5c4..c6885a27 100644 --- a/diffsims/simulations/simulation2d.py +++ b/diffsims/simulations/simulation2d.py @@ -253,6 +253,8 @@ def current_size(self): return self.rotations.size def deepcopy(self): + self._rotation_slider = None + self._phase_slider = None return copy.deepcopy(self) def _get_transformed_coordinates( @@ -379,6 +381,8 @@ def get_diffraction_pattern( produces reasonably good patterns when the lattice parameters are on the order of 0.5nm and the default size and sigma are used. """ + if shape is None: + shape = (256, 256) if direct_beam_position is None: direct_beam_position = (shape[1] // 2, shape[0] // 2) transformed = self._get_transformed_coordinates( @@ -501,12 +505,7 @@ def _get_labels(self, coords, intensity, min_label_intensity, xlim, ylim): & (coords.data[:, 1] < max(ylim)) ) in_range_coords = coords.data[condition] - millers = np.round( - np.matmul( - np.matmul(in_range_coords, self.get_current_rotation_matrix().T), - coords.phase.structure.lattice.base.T, - ) - ).astype(np.int16) + millers = coords.hkl labels = [] for miller, coordinate, inten in zip(millers, in_range_coords, intensity): if np.isnan(inten) or inten > min_label_intensity: @@ -640,11 +639,10 @@ def plot( if ( interactive and self.has_multiple_rotations or self.has_multiple_phases ): # pragma: no cover - axrot = fig.add_axes([0.5, 0.05, 0.4, 0.03]) - axphase = fig.add_axes([0.1, 0.05, 0.2, 0.03]) - fig.subplots_adjust(left=0.25, bottom=0.25) if self.has_multiple_phases: + axphase = fig.add_axes([0.1, 0.05, 0.2, 0.03]) + axrot = fig.add_axes([0.5, 0.05, 0.4, 0.03]) max_rot = np.max([r.size for r in self.rotations]) rotation_slider = Slider( ax=axrot, @@ -665,6 +663,7 @@ def plot( orientation="horizontal", ) else: # self.has_multiple_rotations: + axrot = fig.add_axes([0.5, 0.05, 0.4, 0.03]) rotation_slider = Slider( ax=axrot, label="Rotation", diff --git a/diffsims/tests/crystallography/test_diffracting_vector.py b/diffsims/tests/crystallography/test_diffracting_vector.py index c08c01e7..c7782008 100644 --- a/diffsims/tests/crystallography/test_diffracting_vector.py +++ b/diffsims/tests/crystallography/test_diffracting_vector.py @@ -43,12 +43,11 @@ def test_structure_factor(self, ferrite_phase): with pytest.raises(NotImplementedError): rlv.calculate_structure_factor() - def test_hl(self, ferrite_phase): + def test_hkl(self, ferrite_phase): rlv = ReciprocalLatticeVector(ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]]) rot = Rotation.from_euler([90, 90, 0], degrees=True) - rotated_vectors = (~rot * rlv.to_miller()).data - dv = DiffractingVector(ferrite_phase, xyz=rotated_vectors, rotation=rot) - assert np.allclose(rlv.hkl, dv.hkl) + rotated = rot * rlv + assert np.allclose(rlv.hkl, rotated.hkl) def test_flat_polar(self, ferrite_phase): dv = DiffractingVector(ferrite_phase, xyz=[[1, 1, 1], [0.5, -0.5, 0]]) diff --git a/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py b/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py index d8717107..3ba938e5 100644 --- a/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py +++ b/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py @@ -20,6 +20,7 @@ import numpy as np from orix.crystal_map import Phase from orix.vector import Miller, Vector3d +from orix.quaternion import Rotation import pytest from diffsims.crystallography import ReciprocalLatticeVector @@ -78,6 +79,7 @@ def test_init_from_min_dspacing(self, ferrite_phase, d, desired_size): rlv = ReciprocalLatticeVector.from_min_dspacing(ferrite_phase, d) assert rlv.size == desired_size + @pytest.mark.parametrize("with_zb", [True, False]) @pytest.mark.parametrize( "hkl, desired_highest_hkl, desired_lowest_hkl, desired_size", [ @@ -93,12 +95,33 @@ def test_init_from_highest_hkl( desired_highest_hkl, desired_lowest_hkl, desired_size, + with_zb, ): """Class method gives desired number of vectors and indices.""" - rlv = ReciprocalLatticeVector.from_highest_hkl(silicon_carbide_phase, hkl) - assert np.allclose(rlv[0].hkl, desired_highest_hkl) - assert np.allclose(rlv[-1].hkl, desired_lowest_hkl) + rlv = ReciprocalLatticeVector.from_highest_hkl( + silicon_carbide_phase, hkl, include_zero_vector=with_zb + ) + if with_zb: + desired_lowest_hkl = [0, 0, 0] + desired_size += 1 assert rlv.size == desired_size + assert np.allclose(rlv[-1].hkl, desired_lowest_hkl) + assert np.allclose(rlv[0].hkl, desired_highest_hkl) + + def test_rotate(self, ferrite_phase): + rlv = ReciprocalLatticeVector(ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]]) + rot = Rotation.from_euler([90, 90, 0], degrees=True) + rotated = rot * rlv + assert np.allclose(rlv.hkl, rotated.hkl) + + def test_2rotate(self, ferrite_phase): + rlv = ReciprocalLatticeVector(ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]]) + rot = Rotation.from_euler([90, 90, 0], degrees=True) + rot2 = Rotation.from_euler([90, 45, 0], degrees=True) + rotated = rot * rlv + rotated2 = rot2 * rotated + assert np.allclose(rlv.hkl, rotated.hkl) + assert np.allclose(rlv.hkl, rotated2.hkl) def test_repr(self, ferrite_phase): """String representation gives desired (start of) string.""" @@ -129,6 +152,14 @@ def test_get_item(self, ferrite_phase): rlv[20:23].phase.structure.lattice.abcABG(), ) + def test_get_item_multi_rotation(self, ferrite_phase): + """Indexing gives desired vectors and properties carry over.""" + rlv = ReciprocalLatticeVector(ferrite_phase, [1, 1, 1]) + rot = Rotation.from_euler([[90, 90, 0], [35, 25, 0]], degrees=True) + rlv = rot * rlv + assert rlv[0].rotation.size == 1 + assert rlv[:2].rotation.size == 2 + def test_get_hkil(self, silicon_carbide_phase): """Miller index properties give desired values.""" rlv = ReciprocalLatticeVector.from_min_dspacing( @@ -255,6 +286,13 @@ def test_get_allowed_without_space_group_raises(self): with pytest.raises(ValueError, match=f"The phase {phase} must have a"): _ = rlv.allowed + def test_init_with_wrong_rotation_size_raises(self, nickel_phase): + rot = Rotation.from_euler([[90, 90, 0], [90, 90, 0], [90, 90, 0]], degrees=True) + with pytest.raises(ValueError, match="Rotation must be a single "): + rlv_ni = ReciprocalLatticeVector( + nickel_phase, hkl=[[1, 1, 1], [2, 0, 0]], rotation=rot + ) + def test_multiplicity(self, nickel_phase, silicon_carbide_phase): """Correct vector multiplicity for cubic and hexagonal phases.""" rlv_ni = ReciprocalLatticeVector(nickel_phase, hkl=[[1, 1, 1], [2, 0, 0]]) diff --git a/diffsims/tests/simulations/test_simulations2d.py b/diffsims/tests/simulations/test_simulations2d.py index 0f4a99ff..9d1586ea 100644 --- a/diffsims/tests/simulations/test_simulations2d.py +++ b/diffsims/tests/simulations/test_simulations2d.py @@ -342,6 +342,10 @@ def test_iphase_error(self, multi_simulation): with pytest.raises(ValueError): phase_slic = multi_simulation.iphase[3.1] + def test_get_current_rotation(self, multi_simulation): + rot = multi_simulation.get_current_rotation_matrix() + np.testing.assert_array_equal(rot, multi_simulation.rotations[0].to_matrix()[0]) + def test_irot(self, multi_simulation): sliced_sim = multi_simulation.irot[0] assert isinstance(sliced_sim, Simulation2D) From bfbd299915565d9f259294f8adc14dcf733fe41a Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 8 May 2024 10:38:09 -0500 Subject: [PATCH 31/36] Revert "Bugfix: Fix Rotation with RLV" This reverts commit 05b373691bc4651d4497d681322748575633b3ca. --- .../crystallography/_diffracting_vector.py | 17 +++++- .../reciprocal_lattice_vector.py | 52 +++---------------- diffsims/generators/simulation_generator.py | 5 +- diffsims/simulations/simulation2d.py | 17 +++--- .../test_diffracting_vector.py | 7 +-- .../test_reciprocal_lattice_vector.py | 44 ++-------------- .../tests/simulations/test_simulations2d.py | 4 -- 7 files changed, 40 insertions(+), 106 deletions(-) diff --git a/diffsims/crystallography/_diffracting_vector.py b/diffsims/crystallography/_diffracting_vector.py index 1732a297..46c412ab 100644 --- a/diffsims/crystallography/_diffracting_vector.py +++ b/diffsims/crystallography/_diffracting_vector.py @@ -18,6 +18,7 @@ from diffsims.crystallography import ReciprocalLatticeVector import numpy as np +from orix.vector.miller import _transform_space class DiffractingVector(ReciprocalLatticeVector): @@ -81,13 +82,14 @@ class DiffractingVector(ReciprocalLatticeVector): def __init__( self, phase, xyz=None, hkl=None, hkil=None, intensity=None, rotation=None ): - super().__init__(phase, xyz=xyz, hkl=hkl, hkil=hkil, rotation=rotation) + super().__init__(phase, xyz=xyz, hkl=hkl, hkil=hkil) if intensity is None: self._intensity = np.full(self.shape, np.nan) elif len(intensity) != self.size: raise ValueError("Length of intensity array must match number of vectors") else: self._intensity = np.array(intensity) + self._rotation = rotation def __getitem__(self, key): dv_new = super().__getitem__(key) @@ -122,6 +124,19 @@ def intensity(self, value): raise ValueError("Length of intensity array must match number of vectors") self._intensity = np.array(value) + @property + def lattice_aligned_data(self): + if self._rotation is None: + return self.data + else: + return np.matmul(self.data, self._rotation.to_matrix()[0].T) + + @property + def hkl(self): + return _transform_space( + self.lattice_aligned_data, "c", "r", self.phase.structure.lattice + ) + def calculate_structure_factor(self): raise NotImplementedError( "Structure factor calculation not implemented for DiffractionVector. " diff --git a/diffsims/crystallography/reciprocal_lattice_vector.py b/diffsims/crystallography/reciprocal_lattice_vector.py index 7a734ae9..b424c0d7 100644 --- a/diffsims/crystallography/reciprocal_lattice_vector.py +++ b/diffsims/crystallography/reciprocal_lattice_vector.py @@ -25,7 +25,6 @@ import numba as nb import numpy as np from orix.vector import Miller, Vector3d -from orix.quaternion import Rotation from orix.vector.miller import ( _check_hkil, _get_highest_hkl, @@ -79,8 +78,6 @@ class ReciprocalLatticeVector(Vector3d): Indices of reciprocal lattice vector(s), often preferred over ``hkl`` in trigonal and hexagonal lattices. Default is ``None``. This, ``xyz``, or ``hkl`` is required. - rotation : orix.quaternion.Rotation, optional - Rotation to apply to the vectors. Default is ``None``. Examples -------- @@ -103,7 +100,7 @@ class ReciprocalLatticeVector(Vector3d): """ - def __init__(self, phase, xyz=None, hkl=None, hkil=None, rotation=None): + def __init__(self, phase, xyz=None, hkl=None, hkil=None): self.phase = phase self._raise_if_no_point_group() @@ -123,24 +120,13 @@ def __init__(self, phase, xyz=None, hkl=None, hkil=None, rotation=None): self._coordinate_format = "hkl" xyz = _transform_space(hkl, "r", "c", phase.structure.lattice) super().__init__(xyz) - if rotation is None: - self._rotation = Rotation.from_euler([0, 0, 0], degrees=True) - elif rotation.size == 1 or self.size == rotation.size: - self._rotation = rotation - else: - raise ValueError( - "Rotation must be a single Rotation or have the same size as the vectors" - ) + self._theta = np.full(self.shape, np.nan) self._structure_factor = np.full(self.shape, np.nan, dtype="complex128") def __getitem__(self, key): new_data = self.data[key] - if self.rotation.size == 1: - rotation = self.rotation - else: - rotation = self.rotation[key] - rlv_new = self.__class__(self.phase, xyz=new_data, rotation=rotation) + rlv_new = self.__class__(self.phase, xyz=new_data) if np.isnan(self.structure_factor).all(): rlv_new._structure_factor = np.full( @@ -165,28 +151,6 @@ def __repr__(self): phase_name = self.phase.name return f"{name} {shape}, {phase_name} ({symmetry})\n" f"{data}" - def get_kwargs(self, new_data, *args, **kwargs): - if "rotation" in kwargs: - mult_rotation = kwargs.pop("rotation") - new_rotation = mult_rotation * self.rotation - else: - new_rotation = self.rotation - return dict(phase=self.phase, xyz=new_data, rotation=new_rotation) - - @property - def rotation(self): - """ - Rotation applied to the vectors, for calculating hkl - """ - return self._rotation - - @property - def unrotated_data(self): - """ - Un-rotated data for calculating hkl - """ - return (~self.rotation * self).data - @property def hkl(self): """Miller indices. @@ -209,9 +173,7 @@ def hkl(self): """ - return _transform_space( - self.unrotated_data, "c", "r", self.phase.structure.lattice - ) + return _transform_space(self.data, "c", "r", self.phase.structure.lattice) @property def hkil(self): @@ -1108,11 +1070,9 @@ def from_highest_hkl(cls, phase, hkl, include_zero_vector=False): """ idx = _get_indices_from_highest(highest_indices=hkl) - new = cls(phase, hkl=idx).unique() if include_zero_vector: - new_data = np.vstack((new.hkl, np.zeros(3, dtype=int))) - new = ReciprocalLatticeVector(phase, hkl=new_data) - return new + idx = np.vstack((idx, np.zeros(3, dtype=int))) + return cls(phase, hkl=idx).unique() @classmethod def from_min_dspacing(cls, phase, min_dspacing=0.7, include_zero_vector=False): diff --git a/diffsims/generators/simulation_generator.py b/diffsims/generators/simulation_generator.py index 3902edb2..eefd02b3 100644 --- a/diffsims/generators/simulation_generator.py +++ b/diffsims/generators/simulation_generator.py @@ -361,8 +361,7 @@ def get_intersecting_reflections( control. If not set will be set equal to max_excitation_error. """ initial_hkl = recip.hkl - - rotated_vectors = (~rot * recip).data + rotated_vectors = (~rot * recip.to_miller()).data if with_direct_beam: rotated_vectors = np.vstack([rotated_vectors, [0, 0, 0]]) @@ -395,7 +394,7 @@ def get_intersecting_reflections( intersected_vectors = DiffractingVector( phase=recip.phase, xyz=intersected_vectors, - rotation=~rot, + rotation=rot, ) excitation_error = excitation_error[intersection] r_spot = r_spot[intersection] diff --git a/diffsims/simulations/simulation2d.py b/diffsims/simulations/simulation2d.py index c6885a27..9852c5c4 100644 --- a/diffsims/simulations/simulation2d.py +++ b/diffsims/simulations/simulation2d.py @@ -253,8 +253,6 @@ def current_size(self): return self.rotations.size def deepcopy(self): - self._rotation_slider = None - self._phase_slider = None return copy.deepcopy(self) def _get_transformed_coordinates( @@ -381,8 +379,6 @@ def get_diffraction_pattern( produces reasonably good patterns when the lattice parameters are on the order of 0.5nm and the default size and sigma are used. """ - if shape is None: - shape = (256, 256) if direct_beam_position is None: direct_beam_position = (shape[1] // 2, shape[0] // 2) transformed = self._get_transformed_coordinates( @@ -505,7 +501,12 @@ def _get_labels(self, coords, intensity, min_label_intensity, xlim, ylim): & (coords.data[:, 1] < max(ylim)) ) in_range_coords = coords.data[condition] - millers = coords.hkl + millers = np.round( + np.matmul( + np.matmul(in_range_coords, self.get_current_rotation_matrix().T), + coords.phase.structure.lattice.base.T, + ) + ).astype(np.int16) labels = [] for miller, coordinate, inten in zip(millers, in_range_coords, intensity): if np.isnan(inten) or inten > min_label_intensity: @@ -639,10 +640,11 @@ def plot( if ( interactive and self.has_multiple_rotations or self.has_multiple_phases ): # pragma: no cover + axrot = fig.add_axes([0.5, 0.05, 0.4, 0.03]) + axphase = fig.add_axes([0.1, 0.05, 0.2, 0.03]) + fig.subplots_adjust(left=0.25, bottom=0.25) if self.has_multiple_phases: - axphase = fig.add_axes([0.1, 0.05, 0.2, 0.03]) - axrot = fig.add_axes([0.5, 0.05, 0.4, 0.03]) max_rot = np.max([r.size for r in self.rotations]) rotation_slider = Slider( ax=axrot, @@ -663,7 +665,6 @@ def plot( orientation="horizontal", ) else: # self.has_multiple_rotations: - axrot = fig.add_axes([0.5, 0.05, 0.4, 0.03]) rotation_slider = Slider( ax=axrot, label="Rotation", diff --git a/diffsims/tests/crystallography/test_diffracting_vector.py b/diffsims/tests/crystallography/test_diffracting_vector.py index c7782008..c08c01e7 100644 --- a/diffsims/tests/crystallography/test_diffracting_vector.py +++ b/diffsims/tests/crystallography/test_diffracting_vector.py @@ -43,11 +43,12 @@ def test_structure_factor(self, ferrite_phase): with pytest.raises(NotImplementedError): rlv.calculate_structure_factor() - def test_hkl(self, ferrite_phase): + def test_hl(self, ferrite_phase): rlv = ReciprocalLatticeVector(ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]]) rot = Rotation.from_euler([90, 90, 0], degrees=True) - rotated = rot * rlv - assert np.allclose(rlv.hkl, rotated.hkl) + rotated_vectors = (~rot * rlv.to_miller()).data + dv = DiffractingVector(ferrite_phase, xyz=rotated_vectors, rotation=rot) + assert np.allclose(rlv.hkl, dv.hkl) def test_flat_polar(self, ferrite_phase): dv = DiffractingVector(ferrite_phase, xyz=[[1, 1, 1], [0.5, -0.5, 0]]) diff --git a/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py b/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py index 3ba938e5..d8717107 100644 --- a/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py +++ b/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py @@ -20,7 +20,6 @@ import numpy as np from orix.crystal_map import Phase from orix.vector import Miller, Vector3d -from orix.quaternion import Rotation import pytest from diffsims.crystallography import ReciprocalLatticeVector @@ -79,7 +78,6 @@ def test_init_from_min_dspacing(self, ferrite_phase, d, desired_size): rlv = ReciprocalLatticeVector.from_min_dspacing(ferrite_phase, d) assert rlv.size == desired_size - @pytest.mark.parametrize("with_zb", [True, False]) @pytest.mark.parametrize( "hkl, desired_highest_hkl, desired_lowest_hkl, desired_size", [ @@ -95,33 +93,12 @@ def test_init_from_highest_hkl( desired_highest_hkl, desired_lowest_hkl, desired_size, - with_zb, ): """Class method gives desired number of vectors and indices.""" - rlv = ReciprocalLatticeVector.from_highest_hkl( - silicon_carbide_phase, hkl, include_zero_vector=with_zb - ) - if with_zb: - desired_lowest_hkl = [0, 0, 0] - desired_size += 1 - assert rlv.size == desired_size - assert np.allclose(rlv[-1].hkl, desired_lowest_hkl) + rlv = ReciprocalLatticeVector.from_highest_hkl(silicon_carbide_phase, hkl) assert np.allclose(rlv[0].hkl, desired_highest_hkl) - - def test_rotate(self, ferrite_phase): - rlv = ReciprocalLatticeVector(ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]]) - rot = Rotation.from_euler([90, 90, 0], degrees=True) - rotated = rot * rlv - assert np.allclose(rlv.hkl, rotated.hkl) - - def test_2rotate(self, ferrite_phase): - rlv = ReciprocalLatticeVector(ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]]) - rot = Rotation.from_euler([90, 90, 0], degrees=True) - rot2 = Rotation.from_euler([90, 45, 0], degrees=True) - rotated = rot * rlv - rotated2 = rot2 * rotated - assert np.allclose(rlv.hkl, rotated.hkl) - assert np.allclose(rlv.hkl, rotated2.hkl) + assert np.allclose(rlv[-1].hkl, desired_lowest_hkl) + assert rlv.size == desired_size def test_repr(self, ferrite_phase): """String representation gives desired (start of) string.""" @@ -152,14 +129,6 @@ def test_get_item(self, ferrite_phase): rlv[20:23].phase.structure.lattice.abcABG(), ) - def test_get_item_multi_rotation(self, ferrite_phase): - """Indexing gives desired vectors and properties carry over.""" - rlv = ReciprocalLatticeVector(ferrite_phase, [1, 1, 1]) - rot = Rotation.from_euler([[90, 90, 0], [35, 25, 0]], degrees=True) - rlv = rot * rlv - assert rlv[0].rotation.size == 1 - assert rlv[:2].rotation.size == 2 - def test_get_hkil(self, silicon_carbide_phase): """Miller index properties give desired values.""" rlv = ReciprocalLatticeVector.from_min_dspacing( @@ -286,13 +255,6 @@ def test_get_allowed_without_space_group_raises(self): with pytest.raises(ValueError, match=f"The phase {phase} must have a"): _ = rlv.allowed - def test_init_with_wrong_rotation_size_raises(self, nickel_phase): - rot = Rotation.from_euler([[90, 90, 0], [90, 90, 0], [90, 90, 0]], degrees=True) - with pytest.raises(ValueError, match="Rotation must be a single "): - rlv_ni = ReciprocalLatticeVector( - nickel_phase, hkl=[[1, 1, 1], [2, 0, 0]], rotation=rot - ) - def test_multiplicity(self, nickel_phase, silicon_carbide_phase): """Correct vector multiplicity for cubic and hexagonal phases.""" rlv_ni = ReciprocalLatticeVector(nickel_phase, hkl=[[1, 1, 1], [2, 0, 0]]) diff --git a/diffsims/tests/simulations/test_simulations2d.py b/diffsims/tests/simulations/test_simulations2d.py index 9d1586ea..0f4a99ff 100644 --- a/diffsims/tests/simulations/test_simulations2d.py +++ b/diffsims/tests/simulations/test_simulations2d.py @@ -342,10 +342,6 @@ def test_iphase_error(self, multi_simulation): with pytest.raises(ValueError): phase_slic = multi_simulation.iphase[3.1] - def test_get_current_rotation(self, multi_simulation): - rot = multi_simulation.get_current_rotation_matrix() - np.testing.assert_array_equal(rot, multi_simulation.rotations[0].to_matrix()[0]) - def test_irot(self, multi_simulation): sliced_sim = multi_simulation.irot[0] assert isinstance(sliced_sim, Simulation2D) From 96a48ec3a532a24dcd52574f2842a2aa9e18e6d7 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 8 May 2024 11:48:11 -0500 Subject: [PATCH 32/36] Testing: Improve testing --- diffsims/crystallography/__init__.py | 2 - .../crystallography/_diffracting_vector.py | 18 +----- .../reciprocal_lattice_vector.py | 57 ++++++++++++++++++- diffsims/generators/simulation_generator.py | 13 +++-- .../test_diffracting_vector.py | 9 ++- .../test_reciprocal_lattice_vector.py | 45 ++++++++++++++- .../tests/simulations/test_simulations2d.py | 2 +- 7 files changed, 113 insertions(+), 33 deletions(-) diff --git a/diffsims/crystallography/__init__.py b/diffsims/crystallography/__init__.py index e232e95d..736a8e3d 100644 --- a/diffsims/crystallography/__init__.py +++ b/diffsims/crystallography/__init__.py @@ -27,7 +27,6 @@ get_hkl, ) from diffsims.crystallography.reciprocal_lattice_vector import ReciprocalLatticeVector -from diffsims.crystallography._diffracting_vector import DiffractingVector __all__ = [ "get_equivalent_hkl", @@ -35,5 +34,4 @@ "get_hkl", "ReciprocalLatticePoint", "ReciprocalLatticeVector", - "DiffractingVector", ] diff --git a/diffsims/crystallography/_diffracting_vector.py b/diffsims/crystallography/_diffracting_vector.py index 46c412ab..f5e08fa2 100644 --- a/diffsims/crystallography/_diffracting_vector.py +++ b/diffsims/crystallography/_diffracting_vector.py @@ -79,9 +79,7 @@ class DiffractingVector(ReciprocalLatticeVector): """ - def __init__( - self, phase, xyz=None, hkl=None, hkil=None, intensity=None, rotation=None - ): + def __init__(self, phase, xyz=None, hkl=None, hkil=None, intensity=None): super().__init__(phase, xyz=xyz, hkl=hkl, hkil=hkil) if intensity is None: self._intensity = np.full(self.shape, np.nan) @@ -89,7 +87,6 @@ def __init__( raise ValueError("Length of intensity array must match number of vectors") else: self._intensity = np.array(intensity) - self._rotation = rotation def __getitem__(self, key): dv_new = super().__getitem__(key) @@ -124,19 +121,6 @@ def intensity(self, value): raise ValueError("Length of intensity array must match number of vectors") self._intensity = np.array(value) - @property - def lattice_aligned_data(self): - if self._rotation is None: - return self.data - else: - return np.matmul(self.data, self._rotation.to_matrix()[0].T) - - @property - def hkl(self): - return _transform_space( - self.lattice_aligned_data, "c", "r", self.phase.structure.lattice - ) - def calculate_structure_factor(self): raise NotImplementedError( "Structure factor calculation not implemented for DiffractionVector. " diff --git a/diffsims/crystallography/reciprocal_lattice_vector.py b/diffsims/crystallography/reciprocal_lattice_vector.py index b424c0d7..de9a2d5a 100644 --- a/diffsims/crystallography/reciprocal_lattice_vector.py +++ b/diffsims/crystallography/reciprocal_lattice_vector.py @@ -25,6 +25,7 @@ import numba as nb import numpy as np from orix.vector import Miller, Vector3d +from orix.quaternion import Rotation from orix.vector.miller import ( _check_hkil, _get_highest_hkl, @@ -78,6 +79,8 @@ class ReciprocalLatticeVector(Vector3d): Indices of reciprocal lattice vector(s), often preferred over ``hkl`` in trigonal and hexagonal lattices. Default is ``None``. This, ``xyz``, or ``hkl`` is required. + rotation : orix.quaternion.Rotation, optional + Rotation to apply to the vectors. Default is ``None``. Examples -------- @@ -120,7 +123,6 @@ def __init__(self, phase, xyz=None, hkl=None, hkil=None): self._coordinate_format = "hkl" xyz = _transform_space(hkl, "r", "c", phase.structure.lattice) super().__init__(xyz) - self._theta = np.full(self.shape, np.nan) self._structure_factor = np.full(self.shape, np.nan, dtype="complex128") @@ -151,6 +153,53 @@ def __repr__(self): phase_name = self.phase.name return f"{name} {shape}, {phase_name} ({symmetry})\n" f"{data}" + @property + def basis_rotation(self): + """ + Returns the lattice basis rotation. + """ + return Rotation.from_matrix(self.phase.structure.lattice.baserot) + + @basis_rotation.setter + def basis_rotation(self, value): + """ + Sets the lattice basis rotation. + """ + if isinstance(value, Rotation): + if value.size != 1: + raise ValueError("Rotation must be a single rotation") + matrix = value.to_matrix().squeeze() + elif isinstance(value, np.ndarray): + if value.shape != (3, 3): + raise ValueError("Rotation matrix must be 3x3") + matrix = value + else: + raise ValueError("Rotation must be a Rotation or a 3x3 numpy array") + self.phase.structure.lattice.setLatPar(baserot=matrix) + + def rotate_with_basis(self, 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. + """ + + if rotation.size != 1: + raise ValueError("Rotation must be a single rotation") + # rotate basis + new_phase = self.phase.deepcopy() + br = new_phase.structure.lattice.baserot + # In case the base rotation is set already + new_br = br @ rotation.to_matrix().squeeze() + new_phase.structure.lattice.setLatPar(baserot=new_br) + # rotate vectors + vecs = ~rotation * self.to_miller() + return ReciprocalLatticeVector(new_phase, xyz=vecs.data) + @property def hkl(self): """Miller indices. @@ -1070,9 +1119,11 @@ def from_highest_hkl(cls, phase, hkl, include_zero_vector=False): """ idx = _get_indices_from_highest(highest_indices=hkl) + new = cls(phase, hkl=idx).unique() if include_zero_vector: - idx = np.vstack((idx, np.zeros(3, dtype=int))) - return cls(phase, hkl=idx).unique() + new_data = np.vstack((new.hkl, np.zeros(3, dtype=int))) + new = ReciprocalLatticeVector(phase, hkl=new_data) + return new @classmethod def from_min_dspacing(cls, phase, min_dspacing=0.7, include_zero_vector=False): diff --git a/diffsims/generators/simulation_generator.py b/diffsims/generators/simulation_generator.py index eefd02b3..956cdf99 100644 --- a/diffsims/generators/simulation_generator.py +++ b/diffsims/generators/simulation_generator.py @@ -24,7 +24,8 @@ from orix.quaternion import Rotation from orix.crystal_map import Phase -from diffsims.crystallography import ReciprocalLatticeVector, DiffractingVector +from diffsims.crystallography import ReciprocalLatticeVector +from diffsims.crystallography._diffracting_vector import DiffractingVector from diffsims.utils.shape_factor_models import ( linear, atanc, @@ -361,10 +362,11 @@ def get_intersecting_reflections( control. If not set will be set equal to max_excitation_error. """ initial_hkl = recip.hkl - rotated_vectors = (~rot * recip.to_miller()).data - + rotated_vectors = recip.rotate_with_basis(rotation=rot) + rotated_phase = rotated_vectors.phase + rotated_vectors = rotated_vectors.data if with_direct_beam: - rotated_vectors = np.vstack([rotated_vectors, [0, 0, 0]]) + rotated_vectors = np.vstack([rotated_vectors.data, [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) r_sphere = 1 / wavelength @@ -392,9 +394,8 @@ def get_intersecting_reflections( # select these reflections intersected_vectors = rotated_vectors[intersection] intersected_vectors = DiffractingVector( - phase=recip.phase, + phase=rotated_phase, xyz=intersected_vectors, - rotation=rot, ) excitation_error = excitation_error[intersection] r_spot = r_spot[intersection] diff --git a/diffsims/tests/crystallography/test_diffracting_vector.py b/diffsims/tests/crystallography/test_diffracting_vector.py index c08c01e7..e30b570b 100644 --- a/diffsims/tests/crystallography/test_diffracting_vector.py +++ b/diffsims/tests/crystallography/test_diffracting_vector.py @@ -1,4 +1,5 @@ -from diffsims.crystallography import DiffractingVector, ReciprocalLatticeVector +from diffsims.crystallography import ReciprocalLatticeVector +from diffsims.crystallography._diffracting_vector import DiffractingVector from orix.quaternion import Rotation import pytest @@ -43,11 +44,13 @@ def test_structure_factor(self, ferrite_phase): with pytest.raises(NotImplementedError): rlv.calculate_structure_factor() - def test_hl(self, ferrite_phase): + def test_hkl(self, ferrite_phase): rlv = ReciprocalLatticeVector(ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]]) rot = Rotation.from_euler([90, 90, 0], degrees=True) rotated_vectors = (~rot * rlv.to_miller()).data - dv = DiffractingVector(ferrite_phase, xyz=rotated_vectors, rotation=rot) + ferrite_phase2 = ferrite_phase.deepcopy() + ferrite_phase2.structure.lattice.setLatPar(baserot=rot.to_matrix()[0]) + dv = DiffractingVector(ferrite_phase2, xyz=rotated_vectors) assert np.allclose(rlv.hkl, dv.hkl) def test_flat_polar(self, ferrite_phase): diff --git a/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py b/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py index d8717107..bb1067bc 100644 --- a/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py +++ b/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py @@ -20,6 +20,8 @@ import numpy as np from orix.crystal_map import Phase from orix.vector import Miller, Vector3d +from orix.quaternion import Rotation + import pytest from diffsims.crystallography import ReciprocalLatticeVector @@ -78,6 +80,7 @@ def test_init_from_min_dspacing(self, ferrite_phase, d, desired_size): rlv = ReciprocalLatticeVector.from_min_dspacing(ferrite_phase, d) assert rlv.size == desired_size + @pytest.mark.parametrize("include_zero_beam", [True, False]) @pytest.mark.parametrize( "hkl, desired_highest_hkl, desired_lowest_hkl, desired_size", [ @@ -93,9 +96,15 @@ def test_init_from_highest_hkl( desired_highest_hkl, desired_lowest_hkl, desired_size, + include_zero_beam, ): """Class method gives desired number of vectors and indices.""" - rlv = ReciprocalLatticeVector.from_highest_hkl(silicon_carbide_phase, hkl) + rlv = ReciprocalLatticeVector.from_highest_hkl( + silicon_carbide_phase, hkl, include_zero_vector=include_zero_beam + ) + if include_zero_beam: + desired_size += 1 + desired_lowest_hkl = [0, 0, 0] assert np.allclose(rlv[0].hkl, desired_highest_hkl) assert np.allclose(rlv[-1].hkl, desired_lowest_hkl) assert rlv.size == desired_size @@ -139,6 +148,40 @@ def test_get_hkil(self, silicon_carbide_phase): assert np.allclose(rlv.i, [0, 0, 0, 0, 0, 0]) assert np.allclose(rlv.l, [3, 2, 1, -1, -2, -3]) + def test_get_lattice_basis_rotation(self, ferrite_phase): + """Rotation matrix to align the lattice basis with the Cartesian + basis is correct. + """ + rlv = ReciprocalLatticeVector(ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]]) + rot = rlv.basis_rotation + assert np.allclose(rot.to_matrix(), np.eye(3)) + + def test_set_lattice_basis_rotation(self, ferrite_phase): + """Setting the lattice basis rotation works.""" + rlv = ReciprocalLatticeVector(ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]]) + rot = Rotation.from_euler([90, 90, 0], degrees=True) + rlv.basis_rotation = rot + assert np.allclose(rot.to_matrix(), rlv.basis_rotation.to_matrix()) + rlv.basis_rotation = rot.to_matrix()[0] + assert np.allclose(rot.to_matrix(), rlv.basis_rotation.to_matrix()) + + def test_set_lattice_basis_rotation_raises(self, ferrite_phase): + """Setting the lattice basis rotation works.""" + rlv = ReciprocalLatticeVector(ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]]) + rot = Rotation.from_euler([[90, 90, 0], [90, 90, 1]], degrees=True) + with pytest.raises(ValueError): + rlv.basis_rotation = rot + with pytest.raises(ValueError): + rlv.basis_rotation = rot.to_matrix() + with pytest.raises(ValueError): + rlv.basis_rotation = [1, 2, 3] + + def test_rotation_with_basis_raises(self, ferrite_phase): + rlv = ReciprocalLatticeVector(ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]]) + rot = Rotation.from_euler([[90, 90, 0], [90, 90, 1]], degrees=True) + with pytest.raises(ValueError): + rlv.rotate_with_basis(rotation=rot) + def test_gspacing_dspacing_scattering_parameter(self, ferrite_phase): """Length and scattering parameter properties give desired values. diff --git a/diffsims/tests/simulations/test_simulations2d.py b/diffsims/tests/simulations/test_simulations2d.py index 0f4a99ff..43d5f419 100644 --- a/diffsims/tests/simulations/test_simulations2d.py +++ b/diffsims/tests/simulations/test_simulations2d.py @@ -25,7 +25,7 @@ from diffsims.simulations import Simulation2D from diffsims.generators.simulation_generator import SimulationGenerator -from diffsims.crystallography import DiffractingVector +from diffsims.crystallography._diffracting_vector import DiffractingVector @pytest.fixture(scope="module") From 1f2150db388047464bf0ac1f37b3ca4baa823f3b Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 8 May 2024 12:56:54 -0500 Subject: [PATCH 33/36] Refactor: Improve based on @hakonanes and @viljarjf's suggestions --- .../crystallography/_diffracting_vector.py | 14 +++++---- .../reciprocal_lattice_vector.py | 24 +++----------- .../test_diffracting_vector.py | 7 +++++ .../test_reciprocal_lattice_vector.py | 31 ++++++------------- 4 files changed, 28 insertions(+), 48 deletions(-) diff --git a/diffsims/crystallography/_diffracting_vector.py b/diffsims/crystallography/_diffracting_vector.py index f5e08fa2..4ede0284 100644 --- a/diffsims/crystallography/_diffracting_vector.py +++ b/diffsims/crystallography/_diffracting_vector.py @@ -27,13 +27,14 @@ class DiffractingVector(ReciprocalLatticeVector): All lengths are assumed to be given in Å or inverse Å. - 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 electron diffraction based on the intersection of the Ewald sphere with the reciprocal lattice. 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. + :class:`~diffsims.simulations.DiffractionSimulation` class. It is not (currently) + intended to be used directly by the user. Parameters ---------- @@ -129,9 +130,10 @@ def calculate_structure_factor(self): 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) + flat_self = self.flatten() + r = np.linalg.norm(flat_self.data[:, :2], axis=1) theta = np.arctan2( - self.data[:, 1], - self.data[:, 0], + flat_self.data[:, 1], + flat_self.data[:, 0], ) return r, theta diff --git a/diffsims/crystallography/reciprocal_lattice_vector.py b/diffsims/crystallography/reciprocal_lattice_vector.py index de9a2d5a..e0492de2 100644 --- a/diffsims/crystallography/reciprocal_lattice_vector.py +++ b/diffsims/crystallography/reciprocal_lattice_vector.py @@ -18,7 +18,6 @@ from collections import defaultdict from copy import deepcopy -from typing import Tuple from diffpy.structure.symmetryutilities import expandPosition from diffpy.structure import Structure @@ -160,23 +159,6 @@ def basis_rotation(self): """ return Rotation.from_matrix(self.phase.structure.lattice.baserot) - @basis_rotation.setter - def basis_rotation(self, value): - """ - Sets the lattice basis rotation. - """ - if isinstance(value, Rotation): - if value.size != 1: - raise ValueError("Rotation must be a single rotation") - matrix = value.to_matrix().squeeze() - elif isinstance(value, np.ndarray): - if value.shape != (3, 3): - raise ValueError("Rotation matrix must be 3x3") - matrix = value - else: - raise ValueError("Rotation must be a Rotation or a 3x3 numpy array") - self.phase.structure.lattice.setLatPar(baserot=matrix) - def rotate_with_basis(self, rotation): """Rotate both vectors and the basis with a given `Rotation`. This differs from simply multiplying with a `Rotation`, @@ -1186,9 +1168,11 @@ def from_min_dspacing(cls, phase, min_dspacing=0.7, include_zero_vector=False): dspacing = 1 / phase.structure.lattice.rnorm(hkl) idx = dspacing >= min_dspacing hkl = hkl[idx] + new = cls(phase, hkl=hkl).unique() if include_zero_vector: - hkl = np.vstack((hkl, np.zeros(3, dtype=int))) - return cls(phase, hkl=hkl).unique() + new_data = np.vstack((new.hkl, np.zeros(3, dtype=int))) + new = ReciprocalLatticeVector(phase, hkl=new_data) + return new @classmethod def from_miller(cls, miller): diff --git a/diffsims/tests/crystallography/test_diffracting_vector.py b/diffsims/tests/crystallography/test_diffracting_vector.py index e30b570b..ef3f4d30 100644 --- a/diffsims/tests/crystallography/test_diffracting_vector.py +++ b/diffsims/tests/crystallography/test_diffracting_vector.py @@ -58,3 +58,10 @@ def test_flat_polar(self, ferrite_phase): r, t = dv.to_flat_polar() assert np.allclose(r, [np.sqrt(2), 0.70710678]) assert np.allclose(t, [np.pi / 4, -np.pi / 4]) + dv = DiffractingVector( + ferrite_phase, + xyz=[[[1, 1, 1], [0.5, -0.5, 0]], [[1, 1, 1], [0.5, -0.5, 0]]], + ) + r, t = dv.to_flat_polar() + assert np.allclose(r, [np.sqrt(2), np.sqrt(2), 0.70710678, 0.70710678]) + assert np.allclose(t, [np.pi / 4, np.pi / 4, -np.pi / 4, -np.pi / 4]) diff --git a/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py b/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py index bb1067bc..5e226466 100644 --- a/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py +++ b/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py @@ -74,10 +74,17 @@ def test_init_raises(self, nickel_phase): with pytest.raises(ValueError, match="Exactly one of "): _ = ReciprocalLatticeVector(nickel_phase) + @pytest.mark.parametrize("include_zero_beam", [True, False]) @pytest.mark.parametrize("d, desired_size", [(2, 18), (1, 92), (0.5, 750)]) - def test_init_from_min_dspacing(self, ferrite_phase, d, desired_size): + def test_init_from_min_dspacing( + self, ferrite_phase, d, desired_size, include_zero_beam + ): """Class method gives desired number of vectors.""" - rlv = ReciprocalLatticeVector.from_min_dspacing(ferrite_phase, d) + rlv = ReciprocalLatticeVector.from_min_dspacing( + ferrite_phase, d, include_zero_vector=include_zero_beam + ) + if include_zero_beam: + desired_size += 1 assert rlv.size == desired_size @pytest.mark.parametrize("include_zero_beam", [True, False]) @@ -156,26 +163,6 @@ def test_get_lattice_basis_rotation(self, ferrite_phase): rot = rlv.basis_rotation assert np.allclose(rot.to_matrix(), np.eye(3)) - def test_set_lattice_basis_rotation(self, ferrite_phase): - """Setting the lattice basis rotation works.""" - rlv = ReciprocalLatticeVector(ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]]) - rot = Rotation.from_euler([90, 90, 0], degrees=True) - rlv.basis_rotation = rot - assert np.allclose(rot.to_matrix(), rlv.basis_rotation.to_matrix()) - rlv.basis_rotation = rot.to_matrix()[0] - assert np.allclose(rot.to_matrix(), rlv.basis_rotation.to_matrix()) - - def test_set_lattice_basis_rotation_raises(self, ferrite_phase): - """Setting the lattice basis rotation works.""" - rlv = ReciprocalLatticeVector(ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]]) - rot = Rotation.from_euler([[90, 90, 0], [90, 90, 1]], degrees=True) - with pytest.raises(ValueError): - rlv.basis_rotation = rot - with pytest.raises(ValueError): - rlv.basis_rotation = rot.to_matrix() - with pytest.raises(ValueError): - rlv.basis_rotation = [1, 2, 3] - def test_rotation_with_basis_raises(self, ferrite_phase): rlv = ReciprocalLatticeVector(ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]]) rot = Rotation.from_euler([[90, 90, 0], [90, 90, 1]], degrees=True) From efde31d42eb57cf9d4c4f94bb90835b0286274ee Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 8 May 2024 12:58:59 -0500 Subject: [PATCH 34/36] Refactor: Tests including zero-vector --- diffsims/tests/generators/test_simulation_generator.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/diffsims/tests/generators/test_simulation_generator.py b/diffsims/tests/generators/test_simulation_generator.py index df16fbd7..a40c6f50 100644 --- a/diffsims/tests/generators/test_simulation_generator.py +++ b/diffsims/tests/generators/test_simulation_generator.py @@ -132,7 +132,7 @@ def test_matching_results( diffraction = diffraction_calculator.calculate_diffraction2d( local_structure, reciprocal_radius=5.0 ) - assert diffraction.coordinates.size == 69 + assert diffraction.coordinates.size == 70 def test_precession_simple( self, diffraction_calculator_precession_simple, local_structure @@ -141,7 +141,7 @@ def test_precession_simple( local_structure, reciprocal_radius=5.0, ) - assert diffraction.coordinates.size == 249 + assert diffraction.coordinates.size == 250 def test_precession_full( self, diffraction_calculator_precession_full, local_structure @@ -150,7 +150,7 @@ def test_precession_full( local_structure, reciprocal_radius=5.0, ) - assert diffraction.coordinates.size == 249 + assert diffraction.coordinates.size == 250 def test_custom_shape_func(self, diffraction_calculator_custom, local_structure): diffraction = diffraction_calculator_custom.calculate_diffraction2d( From da45a3270aac8335717d143ed2daaee8d78b0c1d Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Wed, 8 May 2024 16:16:47 -0500 Subject: [PATCH 35/36] Refactor: Remove changes from RLV --- .../crystallography/_diffracting_vector.py | 44 ++++++++++++++++++- .../reciprocal_lattice_vector.py | 39 ++-------------- diffsims/generators/simulation_generator.py | 5 +-- .../test_diffracting_vector.py | 14 ++++++ .../test_reciprocal_lattice_vector.py | 14 ------ 5 files changed, 62 insertions(+), 54 deletions(-) diff --git a/diffsims/crystallography/_diffracting_vector.py b/diffsims/crystallography/_diffracting_vector.py index 4ede0284..fb2a1b13 100644 --- a/diffsims/crystallography/_diffracting_vector.py +++ b/diffsims/crystallography/_diffracting_vector.py @@ -19,6 +19,7 @@ from diffsims.crystallography import ReciprocalLatticeVector import numpy as np from orix.vector.miller import _transform_space +from orix.quaternion import Rotation class DiffractingVector(ReciprocalLatticeVector): @@ -90,7 +91,18 @@ def __init__(self, phase, xyz=None, hkl=None, hkil=None, intensity=None): self._intensity = np.array(intensity) def __getitem__(self, key): - dv_new = super().__getitem__(key) + new_data = self.data[key] + dv_new = self.__class__(self.phase, xyz=new_data) + + if np.isnan(self.structure_factor).all(): + dv_new._structure_factor = np.full(dv_new.shape, np.nan, dtype="complex128") + + else: + dv_new._structure_factor = self.structure_factor[key] + if np.isnan(self.theta).all(): + dv_new._theta = np.full(dv_new.shape, np.nan) + else: + dv_new._theta = self.theta[key] if np.isnan(self.intensity).all(): dv_new._intensity = np.full(dv_new.shape, np.nan) else: @@ -105,6 +117,36 @@ def __getitem__(self, key): return dv_new + @property + def basis_rotation(self): + """ + Returns the lattice basis rotation. + """ + return Rotation.from_matrix(self.phase.structure.lattice.baserot) + + def rotate_with_basis(self, 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. + """ + + if rotation.size != 1: + raise ValueError("Rotation must be a single rotation") + # rotate basis + new_phase = self.phase.deepcopy() + br = new_phase.structure.lattice.baserot + # In case the base rotation is set already + new_br = br @ rotation.to_matrix().squeeze() + new_phase.structure.lattice.setLatPar(baserot=new_br) + # rotate vectors + vecs = ~rotation * self.to_miller() + return ReciprocalLatticeVector(new_phase, xyz=vecs.data) + @property def intensity(self): return self._intensity diff --git a/diffsims/crystallography/reciprocal_lattice_vector.py b/diffsims/crystallography/reciprocal_lattice_vector.py index e0492de2..12f190f2 100644 --- a/diffsims/crystallography/reciprocal_lattice_vector.py +++ b/diffsims/crystallography/reciprocal_lattice_vector.py @@ -24,7 +24,6 @@ import numba as nb import numpy as np from orix.vector import Miller, Vector3d -from orix.quaternion import Rotation from orix.vector.miller import ( _check_hkil, _get_highest_hkl, @@ -78,8 +77,6 @@ class ReciprocalLatticeVector(Vector3d): Indices of reciprocal lattice vector(s), often preferred over ``hkl`` in trigonal and hexagonal lattices. Default is ``None``. This, ``xyz``, or ``hkl`` is required. - rotation : orix.quaternion.Rotation, optional - Rotation to apply to the vectors. Default is ``None``. Examples -------- @@ -126,8 +123,8 @@ def __init__(self, phase, xyz=None, hkl=None, hkil=None): self._structure_factor = np.full(self.shape, np.nan, dtype="complex128") def __getitem__(self, key): - new_data = self.data[key] - rlv_new = self.__class__(self.phase, xyz=new_data) + miller_new = self.to_miller().__getitem__(key) + rlv_new = self.from_miller(miller_new) if np.isnan(self.structure_factor).all(): rlv_new._structure_factor = np.full( @@ -152,36 +149,6 @@ def __repr__(self): phase_name = self.phase.name return f"{name} {shape}, {phase_name} ({symmetry})\n" f"{data}" - @property - def basis_rotation(self): - """ - Returns the lattice basis rotation. - """ - return Rotation.from_matrix(self.phase.structure.lattice.baserot) - - def rotate_with_basis(self, 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. - """ - - if rotation.size != 1: - raise ValueError("Rotation must be a single rotation") - # rotate basis - new_phase = self.phase.deepcopy() - br = new_phase.structure.lattice.baserot - # In case the base rotation is set already - new_br = br @ rotation.to_matrix().squeeze() - new_phase.structure.lattice.setLatPar(baserot=new_br) - # rotate vectors - vecs = ~rotation * self.to_miller() - return ReciprocalLatticeVector(new_phase, xyz=vecs.data) - @property def hkl(self): """Miller indices. @@ -1171,7 +1138,7 @@ def from_min_dspacing(cls, phase, min_dspacing=0.7, include_zero_vector=False): new = cls(phase, hkl=hkl).unique() if include_zero_vector: new_data = np.vstack((new.hkl, np.zeros(3, dtype=int))) - new = ReciprocalLatticeVector(phase, hkl=new_data) + new = cls(phase, hkl=new_data) return new @classmethod diff --git a/diffsims/generators/simulation_generator.py b/diffsims/generators/simulation_generator.py index 956cdf99..0a701e3c 100644 --- a/diffsims/generators/simulation_generator.py +++ b/diffsims/generators/simulation_generator.py @@ -24,7 +24,6 @@ from orix.quaternion import Rotation from orix.crystal_map import Phase -from diffsims.crystallography import ReciprocalLatticeVector from diffsims.crystallography._diffracting_vector import DiffractingVector from diffsims.utils.shape_factor_models import ( linear, @@ -191,7 +190,7 @@ def calculate_diffraction2d( # Rotate using all the rotations in the list vectors = [] for p, rotate in zip(phase, rotation): - recip = ReciprocalLatticeVector.from_min_dspacing( + recip = DiffractingVector.from_min_dspacing( p, min_dspacing=1 / reciprocal_radius, include_zero_vector=with_direct_beam, @@ -335,7 +334,7 @@ def calculate_diffraction1d( def get_intersecting_reflections( self, - recip: ReciprocalLatticeVector, + recip: DiffractingVector, rot: np.ndarray, wavelength: float, max_excitation_error: float, diff --git a/diffsims/tests/crystallography/test_diffracting_vector.py b/diffsims/tests/crystallography/test_diffracting_vector.py index ef3f4d30..bf4ccdf6 100644 --- a/diffsims/tests/crystallography/test_diffracting_vector.py +++ b/diffsims/tests/crystallography/test_diffracting_vector.py @@ -65,3 +65,17 @@ def test_flat_polar(self, ferrite_phase): r, t = dv.to_flat_polar() assert np.allclose(r, [np.sqrt(2), np.sqrt(2), 0.70710678, 0.70710678]) assert np.allclose(t, [np.pi / 4, np.pi / 4, -np.pi / 4, -np.pi / 4]) + + def test_get_lattice_basis_rotation(self, ferrite_phase): + """Rotation matrix to align the lattice basis with the Cartesian + basis is correct. + """ + rlv = DiffractingVector(ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]]) + rot = rlv.basis_rotation + assert np.allclose(rot.to_matrix(), np.eye(3)) + + def test_rotation_with_basis_raises(self, ferrite_phase): + rlv = DiffractingVector(ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]]) + rot = Rotation.from_euler([[90, 90, 0], [90, 90, 1]], degrees=True) + with pytest.raises(ValueError): + rlv.rotate_with_basis(rotation=rot) diff --git a/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py b/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py index 5e226466..fe02cd07 100644 --- a/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py +++ b/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py @@ -155,20 +155,6 @@ def test_get_hkil(self, silicon_carbide_phase): assert np.allclose(rlv.i, [0, 0, 0, 0, 0, 0]) assert np.allclose(rlv.l, [3, 2, 1, -1, -2, -3]) - def test_get_lattice_basis_rotation(self, ferrite_phase): - """Rotation matrix to align the lattice basis with the Cartesian - basis is correct. - """ - rlv = ReciprocalLatticeVector(ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]]) - rot = rlv.basis_rotation - assert np.allclose(rot.to_matrix(), np.eye(3)) - - def test_rotation_with_basis_raises(self, ferrite_phase): - rlv = ReciprocalLatticeVector(ferrite_phase, hkl=[[1, 1, 1], [2, 0, 0]]) - rot = Rotation.from_euler([[90, 90, 0], [90, 90, 1]], degrees=True) - with pytest.raises(ValueError): - rlv.rotate_with_basis(rotation=rot) - def test_gspacing_dspacing_scattering_parameter(self, ferrite_phase): """Length and scattering parameter properties give desired values. From 0e97f2281854f3c78e035248717329ceae5c4c38 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Thu, 9 May 2024 08:35:09 -0500 Subject: [PATCH 36/36] Refactor: suggestions for @hakonanes --- diffsims/crystallography/_diffracting_vector.py | 13 +++++++++++++ .../test_reciprocal_lattice_vector.py | 17 ++++++++--------- 2 files changed, 21 insertions(+), 9 deletions(-) diff --git a/diffsims/crystallography/_diffracting_vector.py b/diffsims/crystallography/_diffracting_vector.py index fb2a1b13..be9a61f3 100644 --- a/diffsims/crystallography/_diffracting_vector.py +++ b/diffsims/crystallography/_diffracting_vector.py @@ -133,6 +133,19 @@ def rotate_with_basis(self, rotation): ---------- rot : orix.quaternion.Rotation A rotation to apply to vectors and the basis. + + Returns + ------- + DiffractingVector + A new DiffractingVector with the rotated vectors and basis. This maintains + the hkl indices of the vectors, but the underlying vector xyz coordinates + are rotated by the given rotation. + + Notes + ----- + Rotating the lattice basis may lead to undefined behavior in orix as it violates + the assumption that the basis is aligned with the crystal axes. Particularly, + applying symmetry operations to the phase may lead to unexpected results. """ if rotation.size != 1: diff --git a/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py b/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py index fe02cd07..bd9ada2f 100644 --- a/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py +++ b/diffsims/tests/crystallography/test_reciprocal_lattice_vector.py @@ -20,7 +20,6 @@ import numpy as np from orix.crystal_map import Phase from orix.vector import Miller, Vector3d -from orix.quaternion import Rotation import pytest @@ -74,20 +73,20 @@ def test_init_raises(self, nickel_phase): with pytest.raises(ValueError, match="Exactly one of "): _ = ReciprocalLatticeVector(nickel_phase) - @pytest.mark.parametrize("include_zero_beam", [True, False]) + @pytest.mark.parametrize("include_zero_vector", [True, False]) @pytest.mark.parametrize("d, desired_size", [(2, 18), (1, 92), (0.5, 750)]) def test_init_from_min_dspacing( - self, ferrite_phase, d, desired_size, include_zero_beam + self, ferrite_phase, d, desired_size, include_zero_vector ): """Class method gives desired number of vectors.""" rlv = ReciprocalLatticeVector.from_min_dspacing( - ferrite_phase, d, include_zero_vector=include_zero_beam + ferrite_phase, d, include_zero_vector=include_zero_vector ) - if include_zero_beam: + if include_zero_vector: desired_size += 1 assert rlv.size == desired_size - @pytest.mark.parametrize("include_zero_beam", [True, False]) + @pytest.mark.parametrize("include_zero_vector", [True, False]) @pytest.mark.parametrize( "hkl, desired_highest_hkl, desired_lowest_hkl, desired_size", [ @@ -103,13 +102,13 @@ def test_init_from_highest_hkl( desired_highest_hkl, desired_lowest_hkl, desired_size, - include_zero_beam, + include_zero_vector, ): """Class method gives desired number of vectors and indices.""" rlv = ReciprocalLatticeVector.from_highest_hkl( - silicon_carbide_phase, hkl, include_zero_vector=include_zero_beam + silicon_carbide_phase, hkl, include_zero_vector=include_zero_vector ) - if include_zero_beam: + if include_zero_vector: desired_size += 1 desired_lowest_hkl = [0, 0, 0] assert np.allclose(rlv[0].hkl, desired_highest_hkl)