diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index e127d820..620a3784 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -45,7 +45,7 @@ jobs: - os: ubuntu-latest python-version: 3.7 OLDEST_SUPPORTED_VERSION: true - DEPENDENCIES: diffpy.structure==3.0.0 matplotlib==3.3 numpy==1.17 orix==0.9.0 scipy==1.1 tqdm==4.9 + DEPENDENCIES: diffpy.structure==3.0.2 matplotlib==3.3 numpy==1.17 orix==0.11.0 scipy==1.2 tqdm==4.9 LABEL: -oldest steps: - uses: actions/checkout@v3 diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 00000000..d959f632 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,8 @@ +repos: + - repo: https://github.com/psf/black + rev: 23.1.0 + hooks: + - id: black + - id: black-jupyter + files: \.ipynb + args: [--line-length=77] \ No newline at end of file diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 0b15cd2a..b80b651f 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -13,6 +13,13 @@ Unreleased Added ----- - Explicit support for Python 3.11. +- Added deprecation tools for deprecating functions and arguments. +- Added new class :class:`diffsims.generators.simulation_generator.SimulationGenerator` for + generating kinemetic diffraction from a phase and a set of reflections. +- Added new class :class:`diffsims.generators.simulations.Simulation` for storing results + from a :class:`diffsims.generators.simulation_generator.SimulationGenerator`. +- Added new class :class:`diffsims.simulations.Simulation1D` for storing the results of a + 1D simulation. Changed ------- @@ -20,6 +27,12 @@ Changed Deprecated ---------- +- Deprecated :class:`diffsims.generators.diffraction_generator.DiffractionGenerator` in + favor of :class:`diffsims.generators.simulation_generator.SimulationGenerator`. +- Deprecated :class:`diffsims.generators.library_generator.DiffractionLibraryGenerator` in + favor of :class:`diffsims.generators.simulation_generator.SimulationGenerator`. +- Deprecated :class:`diffsims.diffraction_library.DiffractionLibrary` in favor of + :class:`diffsims.generators.simulations.Simulation`. Removed ------- diff --git a/CONTRIBUTING.rst b/CONTRIBUTING.rst index 1f5ab479..f9b8c383 100644 --- a/CONTRIBUTING.rst +++ b/CONTRIBUTING.rst @@ -124,6 +124,27 @@ Useful hints on testing: error-prone. See `pytest documentation for more details `_. + +Deprecations +------------ +'This project attempts to adhere to (pre 1.0.0) semantic versioning and ideally no functionality +should be broken between minor releases. Deprecation warnings are raised to give user a full minor +release cycle to adjust their code before a breaking change is implemented.' + +The decorator should be placed right above the object signature to be deprecated:: + +.. code-block:: python + >>> from diffsims.utils._deprecated import deprecate + >>> @deprecate(since=0.8, removal=0.9, alternative="bar") + >>> def foo(self, n): + >>> return n + 1 + + >>> @property + >>> @deprecate(since=0.9, removal=0.10, alternative="another", is_function=True) + >>> def this_property(self): + >>> return 2 + + Build and write documentation ----------------------------- diff --git a/diffsims/crystallography/reciprocal_lattice_point.py b/diffsims/crystallography/reciprocal_lattice_point.py index 656b5b72..c9a9f658 100644 --- a/diffsims/crystallography/reciprocal_lattice_point.py +++ b/diffsims/crystallography/reciprocal_lattice_point.py @@ -29,7 +29,6 @@ get_refraction_corrected_wavelength, ) - _FLOAT_EPS = np.finfo(float).eps # Used to round values below 1e-16 to zero diff --git a/diffsims/crystallography/reciprocal_lattice_vector.py b/diffsims/crystallography/reciprocal_lattice_vector.py index 8ce16ba7..ef54c9f8 100644 --- a/diffsims/crystallography/reciprocal_lattice_vector.py +++ b/diffsims/crystallography/reciprocal_lattice_vector.py @@ -15,7 +15,7 @@ # # You should have received a copy of the GNU General Public License # along with diffsims. If not, see . - +from typing import Tuple from collections import defaultdict from copy import deepcopy @@ -122,10 +122,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 +140,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): @@ -169,7 +182,6 @@ def hkl(self): >>> rlv.hkl array([[1., 1., 1.], [2., 0., 0.]]) - """ return _transform_space(self.data, "c", "r", self.phase.structure.lattice) @@ -502,6 +514,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, self.data.T).T + ) + @property def structure_factor(self): r"""Kinematical structure factors :math:`F`. @@ -1070,7 +1104,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. @@ -1128,6 +1162,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 @@ -1180,6 +1216,54 @@ def from_miller(cls, miller): ) return cls(miller.phase, **{miller.coordinate_format: miller.coordinates}) + def to_polar( + self, degrees: bool = False + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """Convert the vectors to polar coordinates. + + Parameters + ---------- + degrees : bool, optional + Whether to return angles in degrees (default is False). + + Returns + ------- + r : numpy.ndarray + Length of the vectors. + theta : numpy.ndarray + Polar angle of the vectors. + phi : numpy.ndarray + Azimuthal angle of the vectors. + + Examples + -------- + See :class:`ReciprocalLatticeVector` for the creation of ``rlv`` + + >>> rlv + ReciprocalLatticeVector (2,), al (m-3m) + [[1. 1. 1.] + [2. 0. 0.]] + >>> r, theta, phi = rlv.to_polar() + >>> r + array([1.73205081, 2. ]) + >>> theta + array([0.95531662, 0. ]) + >>> phi + array([0., 0.]) + + """ + + x = self.data[:, 0] + y = self.data[:, 1] + z = self.data[:, 2] + + r = np.sqrt(x**2 + y**2) + theta = np.arctan2(y, x) + + if degrees: + theta = np.rad2deg(theta) + return r, theta, z + def to_miller(self): """Return the vectors as a ``Miller`` instance. 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/diffraction_generator.py b/diffsims/generators/diffraction_generator.py index 69cc11de..85727903 100644 --- a/diffsims/generators/diffraction_generator.py +++ b/diffsims/generators/diffraction_generator.py @@ -41,6 +41,7 @@ lorentzian_precession, ) +from diffsims.utils._deprecated import deprecated __all__ = [ "AtomicDiffractionGenerator", @@ -153,6 +154,9 @@ class DiffractionGenerator(object): `custom shape_factor_model` is used. """ + @deprecated(since="0.6.0", + removal="0.7.0", + alternative="diffsims.generators.simulation_generator.SimulationGenerator") def __init__( self, accelerating_voltage, diff --git a/diffsims/generators/library_generator.py b/diffsims/generators/library_generator.py index 82cefb11..908ac4ce 100644 --- a/diffsims/generators/library_generator.py +++ b/diffsims/generators/library_generator.py @@ -27,7 +27,7 @@ from diffsims.utils.sim_utils import get_points_in_sphere from diffsims.utils.vector_utils import get_angle_cartesian_vec - +from diffsims.utils._deprecated import deprecated __all__ = [ "DiffractionLibraryGenerator", @@ -39,7 +39,7 @@ class DiffractionLibraryGenerator: """Computes a library of electron diffraction patterns for specified atomic structures and orientations. """ - + @deprecated(since="0.6.0", alternative="Diffsims.generators.SimulationGenerator") def __init__(self, electron_diffraction_calculator): """Initialises the generator with a diffraction calculator. diff --git a/diffsims/generators/simulation_generator.py b/diffsims/generators/simulation_generator.py new file mode 100644 index 00000000..51d8dcfa --- /dev/null +++ b/diffsims/generators/simulation_generator.py @@ -0,0 +1,399 @@ +# -*- 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 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, + ): + 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/libraries/__init__.py b/diffsims/libraries/__init__.py index 1db3a521..a08b1968 100644 --- a/diffsims/libraries/__init__.py +++ b/diffsims/libraries/__init__.py @@ -21,7 +21,7 @@ from diffsims.libraries import ( diffraction_library, structure_library, - vector_library, + vector_library ) __all__ = [ diff --git a/diffsims/libraries/diffraction_library.py b/diffsims/libraries/diffraction_library.py index b0c084e1..1929724b 100644 --- a/diffsims/libraries/diffraction_library.py +++ b/diffsims/libraries/diffraction_library.py @@ -20,6 +20,8 @@ import numpy as np +from diffsims.generators.diffraction_generator import DiffractionGenerator +from diffsims.utils._deprecated import deprecated __all__ = [ "DiffractionLibrary", @@ -114,6 +116,12 @@ class DiffractionLibrary(dict): """ + @deprecated( + since="0.6.0", + alternative="diffsims.generators.simulation_generator.SimulationGenerator", + alternative_is_function=False, + removal="0.8.0", + ) def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) self.identifiers = None @@ -127,11 +135,11 @@ def get_library_entry(self, phase=None, angle=None): Parameters ---------- - phase : str + phase : str or int Key for the phase of interest. If unspecified the choice is random. angle : tuple The orientation of interest as a tuple of Euler angles following the - Bunge convention [z, x, z] in degrees. If unspecified the choise is + Bunge convention [z, x, z] in degrees. If unspecified the choice is random (the first hit). Returns diff --git a/diffsims/sims/diffraction_simulation.py b/diffsims/sims/diffraction_simulation.py index 8f9fdb83..04b13445 100644 --- a/diffsims/sims/diffraction_simulation.py +++ b/diffsims/sims/diffraction_simulation.py @@ -22,6 +22,7 @@ from diffsims.pattern.detector_functions import add_shot_and_point_spread from diffsims.utils import mask_utils +from diffsims.utils._deprecated import deprecated __all__ = [ @@ -50,6 +51,12 @@ class DiffractionSimulation: zero in each direction. """ + @deprecated( + since="0.6.0", + alternative="diffsims.simulation.Simulation2D", + alternative_is_function=False, + removal="0.8.0", + ) def __init__( self, coordinates, @@ -210,7 +217,7 @@ def _get_transformed_coordinates( y = coords_transformed[:, 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) + rd = np.sqrt(x**2 + y**2) coords_transformed[:, 0] = rd * np.cos(theta) + cx coords_transformed[:, 1] = rd * np.sin(theta) + cy return coords_transformed 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..2e36bb17 --- /dev/null +++ b/diffsims/simulations/simulation1d.py @@ -0,0 +1,88 @@ +# -*- 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 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..7a9d6494 --- /dev/null +++ b/diffsims/simulations/simulation2d.py @@ -0,0 +1,542 @@ +# -*- 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 typing import Union, Sequence, TYPE_CHECKING, Any +import copy + +import numpy as np +import matplotlib.pyplot as plt +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) + if len(coordinates.shape) != 2: + coordinates = coordinates[:, np.newaxis] + 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 r.size != c.size: + 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 __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): + """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() + r_templates[i, : len(r)] = r + theta_templates[i, : len(t)] = t + intensities_templates[i, : len(v.intensity)] = v.intensity + + 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""" + return self.rotations.size > 1 + + 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 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 = coords.hkl.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 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) diff --git a/diffsims/tests/generators/test_simulation_generator.py b/diffsims/tests/generators/test_simulation_generator.py new file mode 100644 index 00000000..f1116e6b --- /dev/null +++ b/diffsims/tests/generators/test_simulation_generator.py @@ -0,0 +1,242 @@ +# -*- 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 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 + + +@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") diff --git a/diffsims/tests/simulations/__init__.py b/diffsims/tests/simulations/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/diffsims/tests/simulations/test_simulation.py b/diffsims/tests/simulations/test_simulation.py new file mode 100644 index 00000000..4ddcba31 --- /dev/null +++ b/diffsims/tests/simulations/test_simulation.py @@ -0,0 +1,368 @@ +# -*- 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 + +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_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_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_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_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_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) diff --git a/diffsims/tests/simulations/test_simulation1d.py b/diffsims/tests/simulations/test_simulation1d.py new file mode 100644 index 00000000..b2121b20 --- /dev/null +++ b/diffsims/tests/simulations/test_simulation1d.py @@ -0,0 +1,69 @@ +# -*- 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 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/utils/test_deprecation.py b/diffsims/tests/utils/test_deprecation.py new file mode 100644 index 00000000..8bfdfa53 --- /dev/null +++ b/diffsims/tests/utils/test_deprecation.py @@ -0,0 +1,154 @@ +# -*- 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 warnings + +import numpy as np +import pytest + +from diffsims.utils._deprecated import deprecated, deprecated_argument + + +class TestDeprecationWarning: + def test_deprecation_since(self): + """Ensure functions decorated with the custom deprecated + decorator returns desired output, raises a desired warning, and + gets the desired additions to their docstring. + """ + + @deprecated(since=0.7, alternative="bar", removal=0.8) + def foo(n): + """Some docstring.""" + return n + 1 + + with pytest.warns(np.VisibleDeprecationWarning) as record: + assert foo(4) == 5 + desired_msg = ( + "Function `foo()` is deprecated and will be removed in version 0.8. Use " + "`bar()` instead." + ) + assert str(record[0].message) == desired_msg + assert foo.__doc__ == ( + "[*Deprecated*] Some docstring.\n\n" + "Notes\n-----\n" + ".. deprecated:: 0.7\n" + f" {desired_msg}" + ) + + @deprecated(since=1.9) + def foo2(n): + """Another docstring. + Notes + ----- + Some existing notes. + """ + return n + 2 + + with pytest.warns(np.VisibleDeprecationWarning) as record2: + assert foo2(4) == 6 + desired_msg2 = "Function `foo2()` is deprecated." + assert str(record2[0].message) == desired_msg2 + assert foo2.__doc__ == ( + "[*Deprecated*] Another docstring." + "\nNotes\n-----\n" + "Some existing notes.\n\n" + ".. deprecated:: 1.9\n" + f" {desired_msg2}" + ) + + def test_deprecation_no_old_doc(self): + @deprecated(since=0.7, alternative="bar", removal=0.8) + def foo(n): + return n + 1 + + with pytest.warns(np.VisibleDeprecationWarning) as record: + assert foo(4) == 5 + desired_msg = ( + "Function `foo()` is deprecated and will be removed in version 0.8. Use " + "`bar()` instead." + ) + assert str(record[0].message) == desired_msg + assert foo.__doc__ == ( + "[*Deprecated*] \n" + "\nNotes\n-----\n" + ".. deprecated:: 0.7\n" + f" {desired_msg}" + ) + + def test_deprecation_not_function(self): + @deprecated( + since=0.7, alternative="bar", removal=0.8, alternative_is_function=False + ) + def foo(n): + return n + 1 + + with pytest.warns(np.VisibleDeprecationWarning) as record: + assert foo(4) == 5 + desired_msg = ( + "Function `foo()` is deprecated and will be removed in version 0.8. Use " + "`bar` instead." + ) + assert str(record[0].message) == desired_msg + assert foo.__doc__ == ( + "[*Deprecated*] \n" + "\nNotes\n-----\n" + ".. deprecated:: 0.7\n" + f" {desired_msg}" + ) + + +class TestDeprecateArgument: + def test_deprecate_argument(self): + """Functions decorated with the custom `deprecated_argument` + decorator returns desired output and raises a desired warning + only if the argument is passed. + """ + + class Foo: + @deprecated_argument(name="a", since="1.3", removal="1.4") + def bar_arg(self, **kwargs): + return kwargs + + @deprecated_argument(name="a", since="1.3", removal="1.4", alternative="b") + def bar_arg_alt(self, **kwargs): + return kwargs + + my_foo = Foo() + + # Does not warn + with warnings.catch_warnings(): + warnings.simplefilter("error") + assert my_foo.bar_arg(b=1) == {"b": 1} + + # Warns + with pytest.warns(np.VisibleDeprecationWarning) as record2: + assert my_foo.bar_arg(a=2) == {"a": 2} + assert str(record2[0].message) == ( + r"Argument `a` is deprecated and will be removed in version 1.4. " + r"To avoid this warning, please do not use `a`. See the documentation of " + r"`bar_arg()` for more details." + ) + + # Warns with alternative + with pytest.warns(np.VisibleDeprecationWarning) as record3: + assert my_foo.bar_arg_alt(a=3) == {"b": 3} + assert str(record3[0].message) == ( + r"Argument `a` is deprecated and will be removed in version 1.4. " + r"To avoid this warning, please do not use `a`. Use `b` instead. See the " + r"documentation of `bar_arg_alt()` for more details." + ) diff --git a/diffsims/utils/_deprecated.py b/diffsims/utils/_deprecated.py new file mode 100644 index 00000000..a4f1bb18 --- /dev/null +++ b/diffsims/utils/_deprecated.py @@ -0,0 +1,154 @@ +# -*- 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 . + + +"""Helper functions and classes for managing diffsims. +This module and documentation is only relevant for diffsims developers, +not for users. +.. warning: + This module and its submodules are for internal use only. Do not + use them in your own code. We may change the API at any time with no + warning. +""" + +import functools +import inspect +from typing import Callable, Optional, Union +import warnings + +import numpy as np + + +class deprecated: + """Decorator to mark deprecated functions with an informative + warning. + Adapted from + `scikit-image + `_ + and `matplotlib + `_. + """ + + def __init__( + self, + since: Union[str, int, float], + alternative: Optional[str] = None, + alternative_is_function: bool = True, + removal: Union[str, int, float, None] = None, + ): + """Visible deprecation warning. + Parameters + ---------- + since + The release at which this API became deprecated. + alternative + An alternative API that the user may use in place of the + deprecated API. + alternative_is_function + Whether the alternative is a function. Default is ``True``. + removal + The expected removal version. + """ + self.since = since + self.alternative = alternative + self.alternative_is_function = alternative_is_function + self.removal = removal + + def __call__(self, func: Callable): + # Wrap function to raise warning when called, and add warning to + # docstring + if self.alternative is not None: + if self.alternative_is_function: + alt_msg = f" Use `{self.alternative}()` instead." + else: + alt_msg = f" Use `{self.alternative}` instead." + else: + alt_msg = "" + if self.removal is not None: + rm_msg = f" and will be removed in version {self.removal}" + else: + rm_msg = "" + msg = f"Function `{func.__name__}()` is deprecated{rm_msg}.{alt_msg}" + + @functools.wraps(func) + def wrapped(*args, **kwargs): + warnings.simplefilter( + action="always", category=np.VisibleDeprecationWarning, append=True + ) + func_code = func.__code__ + warnings.warn_explicit( + message=msg, + category=np.VisibleDeprecationWarning, + filename=func_code.co_filename, + lineno=func_code.co_firstlineno + 1, + ) + return func(*args, **kwargs) + + # Modify docstring to display deprecation warning + old_doc = inspect.cleandoc(func.__doc__ or "").strip("\n") + notes_header = "\nNotes\n-----" + new_doc = ( + f"[*Deprecated*] {old_doc}\n" + f"{notes_header if notes_header not in old_doc else ''}\n" + f".. deprecated:: {self.since}\n" + f" {msg.strip()}" # Matplotlib uses three spaces + ) + wrapped.__doc__ = new_doc + + return wrapped + + +class deprecated_argument: + """Decorator to remove an argument from a function or method's + signature. + Adapted from `scikit-image + `_. + """ + + def __init__(self, name, since, removal, alternative=None): + self.name = name + self.since = since + self.removal = removal + self.alternative = alternative + + def __call__(self, func): + @functools.wraps(func) + def wrapped(*args, **kwargs): + if self.name in kwargs.keys(): + msg = ( + f"Argument `{self.name}` is deprecated and will be removed in " + f"version {self.removal}. To avoid this warning, please do not use " + f"`{self.name}`. " + ) + if self.alternative is not None: + msg += f"Use `{self.alternative}` instead. " + kwargs[self.alternative] = kwargs.pop(self.name) + msg += f"See the documentation of `{func.__name__}()` for more details." + warnings.simplefilter( + action="always", category=np.VisibleDeprecationWarning + ) + func_code = func.__code__ + warnings.warn_explicit( + message=msg, + category=np.VisibleDeprecationWarning, + filename=func_code.co_filename, + lineno=func_code.co_firstlineno + 1, + ) + return func(*args, **kwargs) + + return wrapped diff --git a/diffsims/utils/shape_factor_models.py b/diffsims/utils/shape_factor_models.py index 97015c4a..aa271ec1 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__ = [ @@ -69,7 +70,7 @@ def linear(excitation_error, max_excitation_error): if isinstance(excitation_error, np.ndarray): sf[sf < 0.0] = 0.0 else: - sf = max(sf, 0.) + sf = max(sf, 0.0) return sf @@ -178,7 +179,7 @@ def lorentzian(excitation_error, max_excitation_error): sigma = np.pi / max_excitation_error fac = ( sigma - / (np.pi * (sigma ** 2 * excitation_error ** 2 + 1)) + / (np.pi * (sigma**2 * excitation_error**2 + 1)) * max_excitation_error ) return fac @@ -217,7 +218,57 @@ def lorentzian_precession( [1] L. Palatinus, P. Brázda, M. Jelínek, J. Hrdá, G. Steciuk, M. Klementová, Specifics of the data processing of precession electron diffraction tomography data and their implementation in the program PETS2.0, Acta Crystallogr. Sect. B Struct. Sci. Cryst. Eng. Mater. 75 (2019) 512–522. doi:10.1107/S2052520619007534. """ sigma = np.pi / max_excitation_error - u = sigma ** 2 * (r_spot ** 2 * precession_angle ** 2 - excitation_error ** 2) + 1 - z = np.sqrt(u ** 2 + 4 * sigma ** 2 * excitation_error ** 2) - fac = (sigma / np.pi) * np.sqrt(2 * (u + z) / z ** 2) + u = sigma**2 * (r_spot**2 * precession_angle**2 - excitation_error**2) + 1 + 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 diff --git a/diffsims/utils/sim_utils.py b/diffsims/utils/sim_utils.py index 09af2e80..b6f5c19e 100644 --- a/diffsims/utils/sim_utils.py +++ b/diffsims/utils/sim_utils.py @@ -15,7 +15,6 @@ # # You should have received a copy of the GNU General Public License # along with diffsims. If not, see . - import collections import math @@ -267,7 +266,7 @@ def _get_kinematical_structure_factor( scattering_params=scattering_params, ) - gspacing_squared = g_hkls_array ** 2 + gspacing_squared = g_hkls_array**2 if scattering_params is not None: atomic_scattering_factor = get_atomic_scattering_factors( @@ -368,7 +367,7 @@ def simulate_kinematic_scattering( accelerating_voltage : float Accelerating voltage in keV. simulation_size : int - Simulation size, n, specifies the n x n array size for + Simulation2D size, n, specifies the n x n array size for the simulation calculation. max_k : float Maximum scattering vector magnitude in reciprocal angstroms. @@ -395,7 +394,7 @@ def simulate_kinematic_scattering( kx, ky = np.meshgrid(l, l) # Convert 2D k-vectors into 3D k-vectors accounting for Ewald sphere. - k = np.array((kx, ky, (wavelength / 2) * (kx ** 2 + ky ** 2))) + k = np.array((kx, ky, (wavelength / 2) * (kx**2 + ky**2))) # Calculate scattering vector squared for each k-vector. gs_sq = np.linalg.norm(k, axis=0) ** 2 @@ -411,7 +410,7 @@ def simulate_kinematic_scattering( elif illumination == "gaussian_probe": for r in atomic_coordinates: probe = (1 / (np.sqrt(2 * np.pi) * sigma)) * np.exp( - (-np.abs(((r[0] ** 2) - (r[1] ** 2)))) / (4 * sigma ** 2) + (-np.abs(((r[0] ** 2) - (r[1] ** 2)))) / (4 * sigma**2) ) scattering = scattering + (probe * fs * np.exp(np.dot(k.T, r) * np.pi * 2j)) else: @@ -503,7 +502,6 @@ def uvtw_to_uvw(uvtw): def get_intensities_params(reciprocal_lattice, reciprocal_radius): - """Calculates the variables needed for get_kinematical_intensities Parameters @@ -574,7 +572,7 @@ def get_holz_angle(electron_wavelength, lattice_parameter): k0 = 1.0 / electron_wavelength kz = 1.0 / lattice_parameter in_root = kz * ((2 * k0) - kz) - sin_angle = (in_root ** 0.5) / k0 + sin_angle = (in_root**0.5) / k0 angle = np.arcsin(sin_angle) return angle @@ -716,7 +714,7 @@ def acceleration_voltage_to_velocity(acceleration_voltage): """ - part1 = (1 + (acceleration_voltage * e) / (m_e * c ** 2)) ** 2 + part1 = (1 + (acceleration_voltage * e) / (m_e * c**2)) ** 2 v = c * (1 - (1 / part1)) ** 0.5 return v @@ -741,7 +739,7 @@ def acceleration_voltage_to_relativistic_mass(acceleration_voltage): """ v = acceleration_voltage_to_velocity(acceleration_voltage) - part1 = 1 - (v ** 2) / (c ** 2) + part1 = 1 - (v**2) / (c**2) mr = m_e / (part1) ** 0.5 return mr @@ -773,7 +771,7 @@ def et_to_beta(et, acceleration_voltage): wavelength = acceleration_voltage_to_wavelength(acceleration_voltage) m = acceleration_voltage_to_relativistic_mass(acceleration_voltage) - beta = e * (wavelength ** 2) * m * et / (h ** 2) + beta = e * (wavelength**2) * m * et / (h**2) return beta @@ -792,7 +790,7 @@ def acceleration_voltage_to_wavelength(acceleration_voltage): """ energy = acceleration_voltage * e - wavelength = h / (2 * m_e * energy * (1 + (energy / (2 * m_e * c ** 2)))) ** 0.5 + wavelength = h / (2 * m_e * energy * (1 + (energy / (2 * m_e * c**2)))) ** 0.5 return wavelength @@ -825,6 +823,6 @@ def diffraction_scattering_angle(acceleration_voltage, lattice_size, miller_inde wavelength = acceleration_voltage_to_wavelength(acceleration_voltage) h, k, l = miller_index a = lattice_size - d = a / (h ** 2 + k ** 2 + l ** 2) ** 0.5 + d = a / (h**2 + k**2 + l**2) ** 0.5 scattering_angle = 2 * np.arcsin(wavelength / (2 * d)) return scattering_angle diff --git a/examples/creating_a_simulation_library/simulating_diffraction_patterns.py b/examples/creating_a_simulation_library/simulating_diffraction_patterns.py new file mode 100644 index 00000000..022e7972 --- /dev/null +++ b/examples/creating_a_simulation_library/simulating_diffraction_patterns.py @@ -0,0 +1,92 @@ +# ===================================================== +# Simulating One Diffraction Pattern for a Single Phase +# ===================================================== + +from orix.crystal_map import Phase +from orix.quaternion import Rotation +from diffpy.structure import Atom, Lattice, Structure +import matplotlib.pyplot as plt + +from diffsims.generators.simulation_generator import SimulationGenerator + +a = 5.431 +latt = Lattice(a, a, a, 90, 90, 90) +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(Atom(atype="Si", xyz=[x, y, z], lattice=latt)) # Motif part A + atom_list.append( + Atom(atype="Si", xyz=[x + 0.25, y + 0.25, z + 0.25], lattice=latt) + ) # Motif part B +struct = Structure(atoms=atom_list, lattice=latt) +p = Phase(structure=struct, space_group=227) + +gen = SimulationGenerator( + accelerating_voltage=200, +) +rot = Rotation.from_axes_angles( + [1, 0, 0], 45, degrees=True +) # 45 degree rotation around x-axis +sim = gen.calculate_ed_data(phase=p, rotation=rot) + +sim.plot(show_labels=True) # plot the first (and only) diffraction pattern + +# %% + +sim.coordinates # coordinates of the first (and only) diffraction pattern + +# %% + + +# =========================================================== +# Simulating Multiple Rotations for a Single Phase +# =========================================================== + +rot = Rotation.from_axes_angles( + [1, 0, 0], (0, 15, 30, 45, 60, 75, 90), degrees=True +) # 45 degree rotation around x-axis +sim = gen.calculate_ed_data(phase=p, rotation=rot) + +sim.plot(show_labels=True) # plot the first diffraction pattern + +# %% + +sim.irot[3].plot(show_labels=True) # plot the fourth(45 degrees) diffraction pattern +# %% + +sim.coordinates # coordinates of all the diffraction patterns + +# ============================================================ +# Simulating Multiple Rotations for Multiple Phases +# ============================================================ + +p2 = p.deepcopy() # copy the phase + +p2.name = "al_2" + +rot = Rotation.from_axes_angles( + [1, 0, 0], (0, 15, 30, 45, 60, 75, 90), degrees=True +) # 45 degree rotation around x-axis +sim = gen.calculate_ed_data(phase=[p, p2], rotation=[rot, rot]) + +sim.plot( + include_direct_beam=True, show_labels=True, min_label_intensity=0.1 +) # plot the first diffraction pattern + +# %% + +sim.iphase["al_2"].irot[3].plot( + show_labels=True, min_label_intensity=0.1 +) # plot the fourth(45 degrees) diffraction pattern + + +# =================================== +# Plotting a Real Diffraction Pattern +# =================================== +dp = sim.get_diffraction_pattern( + shape=(512, 512), + calibration=0.01, +) +plt.figure() +plt.imshow(dp) +# %% diff --git a/setup.cfg b/setup.cfg index 46204e6d..16c923af 100644 --- a/setup.cfg +++ b/setup.cfg @@ -20,6 +20,7 @@ precision = 2 [manifix] known_excludes = + examples/** .* .*/** **/*.nbi diff --git a/setup.py b/setup.py index b830d2dd..eeacbcb6 100644 --- a/setup.py +++ b/setup.py @@ -75,13 +75,13 @@ packages=find_packages(), extras_require=extra_feature_requirements, install_requires=[ - "diffpy.structure >= 3.0.0", # First Python 3 support + "diffpy.structure >= 3.0.2", # First Python 3 support "matplotlib >= 3.3", "numba", "numpy >= 1.17", - "orix >= 0.9", + "orix >= 0.11", "psutil", - "scipy >= 1.1", + "scipy >= 1.2", "tqdm >= 4.9", "transforms3d", ],