Getting closer, still some tests left

Author: viljarjf

diffsims/crystallography/reciprocal_lattice_vector.py

@@ -20,7 +20,7 @@
 from copy import deepcopy
 
 from diffpy.structure.symmetryutilities import expandPosition
-from diffpy.structure import Structure
+from diffpy.structure import Structure, Lattice
 import numba as nb
 import numpy as np
 from orix.vector import Miller, Vector3d
@@ -696,33 +696,17 @@ def calculate_structure_factor(
         """
 
         # Compute one structure factor per set {hkl}
-        hkl_sets = self.get_hkl_sets()
-
-        # For each set, get the indices of the first vector in the
-        # present vectors, accounting for potential multiple dimensions
-        # and avoding computing the unique vectors again
-        first_idx = []
-        for arr in list(hkl_sets.values()):
-            i = []
-            for arr_i in arr:
-                i.append(arr_i[0])
-            first_idx.append(i)
-        first_idx_arr = np.array(first_idx).T
-
-        # Get 2D array of unique vectors, one for each set
-        hkl_unique = self.hkl[tuple(first_idx_arr)]
-
+        unique, inds = self.unique(use_symmetry=True, return_inverse=True)
+        hkl_unique = unique.hkl
         structure_factor = _get_kinematical_structure_factor(
             structure=self.phase.structure,
             g_indices=hkl_unique,
             g_hkls_array=self.phase.structure.lattice.rnorm(hkl_unique),
             debye_waller_factors=debye_waller_factors,
             scattering_params=scattering_params,
         )
-
         # Set structure factors of all symmetrically equivalent vectors
-        for i, idx in enumerate(hkl_sets.values()):
-            self._structure_factor[idx] = structure_factor[i]
+        self._structure_factor[:] = structure_factor[inds]
 
     def calculate_theta(self, voltage):
         r"""Populate :attr:`theta` with the Bragg angle :math:`theta_B`
@@ -1517,7 +1501,7 @@ def transpose(self, *axes):
         new._update_shapes()
         return new
 
-    def unique(self, use_symmetry=False, return_index=False):
+    def unique(self, use_symmetry=False, return_index=False, return_inverse=False):
         """The unique vectors.
 
         Parameters
@@ -1528,48 +1512,51 @@ def unique(self, use_symmetry=False, return_index=False):
         return_index : bool, optional
             Whether to return the indices of the (flattened) data where
             the unique entries were found. Default is ``False``.
+        return_inverse: bool, optional
+            Whether to also return the indices to reconstruct the
+            (flattened) data from the unique data.
 
         Returns
         -------
         ReciprocalLatticeVector
             Flattened unique vectors.
         idx : numpy.ndarray
             Indices of the unique data in the (flattened) array.
-
+        inv : numpy.ndarray
+            Indices to reconstruct the original vectors from the unique
         """
 
         # TODO: Reduce floating point precision in orix!
         def miller_unique(miller, use_symmetry=False):
-            v, idx = Vector3d(miller).unique(return_index=True)
+            v, idx, inv = Vector3d(miller).unique(return_index=True, return_inverse=True)
+            idx = idx[::-1]
+            inv = inv[::-1]
 
             if use_symmetry:
                 operations = miller.phase.point_group
                 n_v = v.size
-                v2 = operations.outer(v).flatten().reshape(*(n_v, operations.size))
+                v2 = operations.outer(v[::-1]).flatten().reshape(*(n_v, operations.size))
                 data = v2.data.round(6)  # Reduced precision
-                data_sorted = np.zeros_like(data)
-                for i in range(n_v):
-                    a = data[i]
-                    order = np.lexsort(a.T)  # Sort by column 1, 2, then 3
-                    data_sorted[i] = a[order]
-                _, idx = np.unique(data_sorted, return_index=True, axis=0)
-                v = v[idx[::-1]]
+                _, idx, inv = np.unique(data, return_index=True, return_inverse=True, axis=0)
+                v = v[idx]
 
             m = miller.__class__(xyz=v.data, phase=miller.phase)
             m.coordinate_format = miller.coordinate_format
-            return m, idx
 
-        #        kwargs = dict(use_symmetry=use_symmetry, return_index=True)
-        #        miller, idx = self.to_miller().unique(**kwargs)
-        miller, idx = miller_unique(self.to_miller(), use_symmetry)
-        idx = idx[::-1]
+            return m, idx, inv
+
+        miller, idx, inv = miller_unique(self.to_miller(), use_symmetry)
 
         new_rlv = self.from_miller(miller)
         new_rlv._structure_factor = self.structure_factor.ravel()[idx]
         new_rlv._theta = self.theta.ravel()[idx]
 
-        if return_index:
+        if return_index and return_inverse:
+            return new_rlv, idx, inv
+        elif return_index:
             return new_rlv, idx
+        elif return_inverse:
+            return new_rlv, inv
         else:
             return new_rlv
 
@@ -1623,7 +1610,7 @@ def _atom_eq(atom1, atom2):
 
     return (
         atom1.element == atom2.element
-        and np.allclose(atom1.xyz, atom2.xyz, atol=1e-7)
+        and np.allclose(atom1.xyz_cartn, atom2.xyz_cartn, atol=1e-7)
         and atom1.occupancy == atom2.occupancy
         and np.allclose(atom1.U, atom2.U, atol=1e-7)
         and np.allclose(atom1.Uisoequiv, atom2.Uisoequiv, atol=1e-7)
@@ -1646,9 +1633,17 @@ def _expand_unit_cell(space_group, structure):
     new_structure : diffpy.structure.Structure
 
     """
-
+    # Transform to diffpy axis conventions
+    structure_matrix = structure.lattice.base
+    pos = structure.xyz_cartn
+    structure = Structure(
+        atoms=list(structure),
+        lattice=Lattice(*structure.lattice.abcABG()),
+    )
     new_structure = Structure(lattice=structure.lattice)
+    structure.xyz_cartn = pos
 
+    # Perform symmetry expansion
     for atom in structure:
         equal = []
         for atom2 in new_structure:
@@ -1660,6 +1655,10 @@ def _expand_unit_cell(space_group, structure):
                 new_atom.xyz = new_position
                 new_structure.append(new_atom)
 
+    # Transform back to orix convention
+    old_xyz_cartn = new_structure.xyz_cartn
+    new_structure.lattice.setLatBase(structure_matrix)
+    new_structure.xyz_cartn = old_xyz_cartn
     return new_structure
 
 

diffsims/generators/simulation_generator.py

@@ -133,9 +133,7 @@ def wavelength(self):
     def calculate_diffraction2d(
         self,
         phase: Union[Phase, Sequence[Phase]],
-        rotation: Union[Rotation, Sequence[Rotation]] = Rotation.from_euler(
-            (0, 0, 0), degrees=True
-        ),
+        rotation: Union[Rotation, Sequence[Rotation]] = Rotation.identity(),
         reciprocal_radius: float = 1.0,
         with_direct_beam: bool = True,
         max_excitation_error: float = 1e-2,
@@ -201,12 +199,10 @@ def calculate_diffraction2d(
                 include_zero_vector=with_direct_beam,
             )
             recip.sanitise_phase()
-            # No need to pre-calculate structure factors if only a few rotations are present
-            if rotate.size > 4:
-                recip.calculate_structure_factor(
-                    scattering_params=self.scattering_params,
-                    debye_waller_factors=debye_waller_factors,
-                    )
+            recip.calculate_structure_factor(
+                scattering_params=self.scattering_params,
+                debye_waller_factors=debye_waller_factors,
+            )
             phase_vectors = []
             for rot in rotate:
                 # Calculate the reciprocal lattice vectors that intersect the Ewald sphere.
@@ -227,15 +223,7 @@ def calculate_diffraction2d(
                     excitation_error, max_excitation_error, r_spot, shape_factor_width
                 )
                 # Calculate intensity
-                if rotate.size <= 4:
-                    v = recip[intersection]
-                    v.calculate_structure_factor(
-                        scattering_params=self.scattering_params,
-                        debye_waller_factors=debye_waller_factors,
-                        )
-                    f_hkls = v.structure_factor
-                else:
-                    f_hkls = recip[intersection].structure_factor
+                f_hkls = recip.structure_factor[intersection]
                 intensities = (f_hkls * f_hkls.conjugate()).real * shape_factor
 
                 # Threshold peaks included in simulation as factor of zero beam intensity.

diffsims/tests/generators/test_simulation_generator.py

@@ -196,12 +196,13 @@ def test_appropriate_scaling(self, diffraction_calculator: SimulationGenerator):
         big_diffraction = diffraction_calculator.calculate_diffraction2d(
             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))]
+        indices = [tuple(i.tolist()) for i in diffraction.coordinates.hkl.astype(int)]
+        big_indices = [tuple(i.tolist()) for i in big_diffraction.coordinates.hkl.astype(int)]
+        target = (2, 2, 0)
+        assert target in indices
+        assert target in big_indices
+        coordinates = diffraction.coordinates[indices.index(target)]
+        big_coordinates = big_diffraction.coordinates[big_indices.index(target)]
         assert np.allclose(coordinates.data, big_coordinates.data * 2)
 
     def test_appropriate_intensities(self, diffraction_calculator, local_structure):
@@ -345,7 +346,7 @@ def test_same_simulation_results():
     # old_data = diff_lib["Graphite"]["simulations"][0].get_diffraction_pattern(shape=shape, sigma=sigma)
 
     # New
-    p = Phase("Graphite", structure=structure_matrix, space_group=1) # Space group 1 ensures no symmetry is applied
+    p = Phase("Graphite", structure=structure_matrix, space_group=194)
     gen = SimulationGenerator(**generator_kwargs)
     rot = Rotation.from_euler(euler_angles_new, degrees=True)
     sim = gen.calculate_diffraction2d(