Merge pull request #238 from CSSFrancis/bugfix_simulation

Bugfix: Small Bugfix for polar flatten simulation

Author: CSSFrancis

.github/workflows/build.yml

@@ -46,7 +46,7 @@ jobs:
           - os: ubuntu-latest
             python-version: 3.9
             OLDEST_SUPPORTED_VERSION: true
-            DEPENDENCIES: diffpy.structure==3.0.2 matplotlib==3.5 numpy==1.17.3 orix==0.12.1 scipy==1.8 tqdm==4.61.2
+            DEPENDENCIES: diffpy.structure==3.0.2 matplotlib==3.7 numpy==1.24 orix==0.12.1 scipy==1.10 tqdm==4.61.2
             LABEL: -oldest
     steps:
       - uses: actions/checkout@v4

diffsims/generators/simulation_generator.py

@@ -27,6 +27,7 @@
 from orix.crystal_map import Phase
 
 from diffsims.crystallography._diffracting_vector import DiffractingVector
+from diffsims.crystallography import ReciprocalLatticeVector
 from diffsims.utils.shape_factor_models import (
     linear,
     atanc,
@@ -280,65 +281,37 @@ def calculate_diffraction1d(
         debye_waller_factors
             Maps element names to their temperature-dependent Debye-Waller factors.
         """
-        latt = phase.structure.lattice
+        rpl = ReciprocalLatticeVector.from_min_dspacing(phase)
+        rpl = rpl.unique(use_symmetry=True)
 
-        # 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(
+        intensities = get_kinematical_intensities(
             phase.structure,
-            g_indices,
-            np.asarray(g_hkls),
-            prefactor=multiplicities,
+            rpl.hkl,
+            rpl.gspacing,
+            prefactor=rpl.multiplicity,
             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.append((l, [i, g]))
+        if rpl.has_hexagonal_lattice:
+            hkl = rpl.hkil
+        else:
+            hkl = rpl.hkl
+        g_vectors = rpl.gspacing
 
         # Scale intensities so that the max intensity is 100.
 
-        max_intensity = max([v[1][0] for v in peaks])
-        reciporical_spacing = []
-        intensities = []
-        hkls = []
-        for p in peaks:
-            label, v = p  # (label, [intensity,g])
-            if v[0] / max_intensity * 100 > minimum_intensity and (label != "000"):
-                reciporical_spacing.append(v[1])
-                intensities.append(v[0])
-                hkls.append(label)
-
+        max_intensity = np.max(intensities)
         intensities = np.asarray(intensities) / max(intensities) * 100
+        mask = intensities > minimum_intensity
+        intensities = intensities[mask]
+        g_vectors = g_vectors[mask]
+        hkl = hkl[mask]
 
         return Simulation1D(
             phase=phase,
-            reciprocal_spacing=reciporical_spacing,
+            reciprocal_spacing=g_vectors,
             intensities=intensities,
-            hkl=hkls,
+            hkl=hkl,
             reciprocal_radius=reciprocal_radius,
             wavelength=self.wavelength,
         )

diffsims/simulations/simulation1d.py

@@ -72,15 +72,45 @@ def __repr__(self):
     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,"""
+    def plot(
+        self,
+        ax=None,
+        annotate_peaks=False,
+        fontsize=12,
+        with_labels=True,
+        rotation=90,
+        va="bottom",
+        ha="center",
+        **kwargs,
+    ):
+        """Plots the 1D diffraction pattern,
+
+        Parameters
+        ----------
+        ax : matplotlib.axes.Axes, optional
+            The axes to plot on. If None, a new figure and axes are created.
+        annotate_peaks : bool, optional
+            Whether to annotate the peaks with their hkl indices. Default is False.
+
+        """
         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)
+                label = "[" + "".join([str(int(np.round(x))) for x in label]) + "]"
+
+                ax.annotate(
+                    label,
+                    xy=[g, i],
+                    xytext=[g, i + 4],
+                    fontsize=fontsize,
+                    rotation=rotation,
+                    va=va,
+                    ha=ha,
+                    **kwargs,
+                )
 
             if with_labels:
                 ax.set_xlabel("A ($^{-1}$)")

diffsims/simulations/simulation2d.py

@@ -327,14 +327,25 @@ def polar_flatten_simulations(self, radial_axes=None, azimuthal_axes=None):
         intensities_templates = np.zeros((len(flattened_vectors), max_num_spots))
         for i, v in enumerate(flattened_vectors):
             r, t = v.to_flat_polar()
+            inten = v.intensity
             if radial_axes is not None and azimuthal_axes is not None:
-                r = get_closest(radial_axes, r)
-                t = get_closest(azimuthal_axes, t)
-                r = r[r < len(radial_axes)]
-                t = t[t < len(azimuthal_axes)]
-            r_templates[i, : len(r)] = r
-            theta_templates[i, : len(t)] = t
-            intensities_templates[i, : len(v.intensity)] = v.intensity
+                r = get_closest(
+                    radial_axes, r
+                )  # convert from real to pixel coordinates
+                t = get_closest(
+                    azimuthal_axes, t
+                )  # convert from real to pixel coordinates
+                mask = (r < len(radial_axes) - 1) & (
+                    t < len(azimuthal_axes) - 1
+                )  # combined mask for out-of-bounds indices
+                r = r[mask]  # apply combined mask
+                t = t[mask]  # apply combined mask
+                inten = inten[mask]  # apply combined mask
+            r_templates[i, : len(r)] = (
+                r  # set the r coordinates (len r and t should be the same)
+            )
+            theta_templates[i, : len(r)] = t  # set the theta coordinates
+            intensities_templates[i, : len(inten)] = inten
         if radial_axes is not None and azimuthal_axes is not None:
             r_templates = np.array(r_templates, dtype=int)
             theta_templates = np.array(theta_templates, dtype=int)

diffsims/tests/generators/test_simulation_generator.py

@@ -19,6 +19,7 @@
 import numpy as np
 import pytest
 from pathlib import Path
+import re
 
 import diffpy.structure
 from orix.crystal_map import Phase
@@ -239,7 +240,6 @@ def test_simulate_1d(self, is_hex):
         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:
@@ -360,10 +360,9 @@ def test_calculate_diffraction2d_progressbar_single_phase(capsys):
     sims = gen.calculate_diffraction2d(phase, rots, show_progressbar=True)
 
     captured = capsys.readouterr()
-    expected = "test phase: 100%|██████████| 10/10"  # also some more, but that is compute-time dependent
-    # ignore possible flushing
-    captured = captured.err.split("\r")[-1]
-    assert captured[: len(expected)] == expected
+    # Accept any number of "█" characters
+    expected = "test phase: 100\%\|█+\| 10\/10"
+    assert re.findall(expected, captured.err)
 
 
 def test_calculate_diffraction2d_progressbar_multi_phase(capsys):
@@ -389,15 +388,8 @@ def test_calculate_diffraction2d_progressbar_multi_phase(capsys):
     )
 
     captured = capsys.readouterr()
-    expected1 = "A: 100%|██████████| 10/10 "
-    expected2 = "B: 100%|██████████| 10/10 "
-    # Find the correct output in the stream, i.e. final line containing the name of the phase
-    captured1 = ""
-    captured2 = ""
-    for line in captured.err.split("\r"):
-        if "A" in line:
-            captured1 = line
-        if "B" in line:
-            captured2 = line
-    assert captured1[: len(expected1)] == expected1
-    assert captured2[: len(expected2)] == expected2
+    # Accept any number of "█" characters
+    expected1 = "A: 100\%\|█+\| 10\/10 "
+    expected2 = "B: 100\%\|█+\| 10\/10 "
+    assert re.findall(expected1, captured.err)
+    assert re.findall(expected2, captured.err)

diffsims/tests/simulations/test_simulations1d.py

@@ -30,7 +30,7 @@ class TestSingleSimulation:
     def simulation1d(self):
         al_phase = make_phase()
         al_phase.name = "Al"
-        hkls = np.array(["100", "110", "111"])
+        hkls = np.array([[1, 0, 0], [1, 1, 0], [1, 1, 1]])
         magnitudes = np.array([1, 2, 3])
         inten = np.array([1, 2, 3])
         recip = 4.0

diffsims/tests/simulations/test_simulations2d.py

@@ -46,7 +46,20 @@ class TestSingleSimulation:
     def single_simulation(self, al_phase):
         gen = SimulationGenerator(accelerating_voltage=200)
         rot = Rotation.from_axes_angles([1, 0, 0], 45, degrees=True)
-        coords = DiffractingVector(phase=al_phase, xyz=[[1, 0, 0]])
+        coords = DiffractingVector(
+            phase=al_phase,
+            xyz=[
+                [1, 0, 0],
+                [2, 0, 0],
+                [3, 3, 0],
+                [-4, 0, 0],
+                [-5, 0, 0],
+                [-6, 0, 0],
+                [-7, 0, 0],
+                [-8, 0, 0],
+            ],
+            intensity=[1, 2, 3, 4, 5, 6, 7, 8],
+        )
         sim = Simulation2D(
             phases=al_phase, simulation_generator=gen, coordinates=coords, rotations=rot
         )
@@ -91,12 +104,12 @@ def test_polar_flatten(self, single_simulation):
             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)
+        assert r_templates.shape == (1, 8)
+        assert theta_templates.shape == (1, 8)
+        assert intensities_templates.shape == (1, 8)
 
     def test_polar_flatten_axes(self, single_simulation):
-        radial_axes = np.linspace(0, 1, 10)
+        radial_axes = np.linspace(0, 7, 5)
         theta_axes = np.linspace(0, 2 * np.pi, 10)
         (
             r_templates,
@@ -105,9 +118,13 @@ def test_polar_flatten_axes(self, single_simulation):
         ) = single_simulation.polar_flatten_simulations(
             radial_axes=radial_axes, azimuthal_axes=theta_axes
         )
-        assert r_templates.shape == (1, 1)
-        assert theta_templates.shape == (1, 1)
-        assert intensities_templates.shape == (1, 1)
+        assert r_templates.shape == (1, 8)
+        assert theta_templates.shape == (1, 8)
+        assert intensities_templates.shape == (1, 8)
+        # The last 2 elements should be zero
+        np.testing.assert_array_equal(r_templates[:, 6:], 0)
+        np.testing.assert_array_equal(theta_templates[:, 6:], 0)
+        np.testing.assert_array_equal(intensities_templates[:, 6:], 0)
 
     def test_deepcopy(self, single_simulation):
         copied = single_simulation.deepcopy()