Merge pull request #9 from Marcraven/main

Vectorization of xrd_simulator
FABLE-3DXRD · Jun 13, 2024 · 030f530 · 030f530
2 parents e8748d9 + 7acef6c
commit 030f530
Show file tree

Hide file tree

Showing 23 changed files with 1,663 additions and 1,028 deletions.
diff --git a/.gitignore b/.gitignore
@@ -35,7 +35,7 @@ vscode
 .vscode/
 .vscode/*
 */.VSCodeCounter/*
-
+makefile
 *.pyc
 
 # Local History for Visual Studio Code
@@ -123,6 +123,7 @@ venv/
 ENV/
 env.bak/
 venv.bak/
+.micromamba
 
 # Spyder project settings
 .spyderproject
@@ -150,4 +151,3 @@ cspell.json
 tmp_profile_dump
 pypi_release/
 extras/
-
diff --git a/docs/source/examples/example_readme.py b/docs/source/examples/example_readme.py
@@ -96,6 +96,5 @@
 polycrystal.diffract(beam, detector, motion)
 ##example:
 
-import os
 path = os.path.join( os.path.dirname(__file__), "..", 'images', "diffraction_pattern.png")
 fig.savefig(os.path.abspath(path), bbox_inches='tight', transparent=True)
diff --git a/setup.py b/setup.py
@@ -33,5 +33,6 @@
                       "dill",
                       "xfab",
                       "netcdf4",
-                      "h5py"]
+                      "h5py",
+                      "pandas"]
 )
diff --git a/tests/end_to_end_tests/powder.py b/tests/end_to_end_tests/powder.py
@@ -4,6 +4,7 @@
 from xrd_simulator.beam import Beam
 from xrd_simulator.motion import RigidBodyMotion
 from xrd_simulator.templates import get_uniform_powder_sample
+import os
 
 pixel_size = 75.
 detector_size = pixel_size * 1024

diff --git a/tests/end_to_end_tests/saves/fast_polycrystal_from_odf.pc b/tests/end_to_end_tests/saves/fast_polycrystal_from_odf.pc
diff --git a/tests/end_to_end_tests/saves/fast_polycrystal_from_odf.xdmf b/tests/end_to_end_tests/saves/fast_polycrystal_from_odf.xdmf
@@ -1 +1 @@
-<Xdmf Version="3.0"><Domain><Grid Name="Grid"><Geometry GeometryType="XYZ"><DataItem DataType="Float" Dimensions="27 3" Format="HDF" Precision="4">fast_polycrystal_from_odf.h5:/data0</DataItem></Geometry><Topology TopologyType="Tetrahedron" NumberOfElements="55" NodesPerElement="4"><DataItem DataType="Int" Dimensions="55 4" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data1</DataItem></Topology><Attribute Name="$\epsilon_{11}$" AttributeType="Scalar" Center="Cell"><DataItem DataType="Float" Dimensions="55" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data2</DataItem></Attribute><Attribute Name="$\epsilon_{22}$" AttributeType="Scalar" Center="Cell"><DataItem DataType="Float" Dimensions="55" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data3</DataItem></Attribute><Attribute Name="$\epsilon_{33}$" AttributeType="Scalar" Center="Cell"><DataItem DataType="Float" Dimensions="55" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data4</DataItem></Attribute><Attribute Name="$\epsilon_{12}$" AttributeType="Scalar" Center="Cell"><DataItem DataType="Float" Dimensions="55" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data5</DataItem></Attribute><Attribute Name="$\epsilon_{13}$" AttributeType="Scalar" Center="Cell"><DataItem DataType="Float" Dimensions="55" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data6</DataItem></Attribute><Attribute Name="$\epsilon_{23}$" AttributeType="Scalar" Center="Cell"><DataItem DataType="Float" Dimensions="55" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data7</DataItem></Attribute><Attribute Name="$\varphi_{1}$" AttributeType="Scalar" Center="Cell"><DataItem DataType="Float" Dimensions="55" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data8</DataItem></Attribute><Attribute Name="$\Phi$" AttributeType="Scalar" Center="Cell"><DataItem DataType="Float" Dimensions="55" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data9</DataItem></Attribute><Attribute Name="$\varphi_{2}$" AttributeType="Scalar" Center="Cell"><DataItem DataType="Float" Dimensions="55" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data10</DataItem></Attribute><Attribute Name="Phases" AttributeType="Scalar" Center="Cell"><DataItem DataType="Int" Dimensions="55" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data11</DataItem></Attribute></Grid></Domain></Xdmf>
+<Xdmf Version="3.0"><Domain><Grid Name="Grid"><Geometry GeometryType="XYZ"><DataItem DataType="Float" Dimensions="27 3" Format="HDF" Precision="4">fast_polycrystal_from_odf.h5:/data0</DataItem></Geometry><Topology TopologyType="Tetrahedron" NumberOfElements="55" NodesPerElement="4"><DataItem DataType="Int" Dimensions="55 4" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data1</DataItem></Topology><Attribute Name="$\epsilon_{11}$" AttributeType="Scalar" Center="Cell"><DataItem DataType="Float" Dimensions="55" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data2</DataItem></Attribute><Attribute Name="$\epsilon_{22}$" AttributeType="Scalar" Center="Cell"><DataItem DataType="Float" Dimensions="55" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data3</DataItem></Attribute><Attribute Name="$\epsilon_{33}$" AttributeType="Scalar" Center="Cell"><DataItem DataType="Float" Dimensions="55" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data4</DataItem></Attribute><Attribute Name="$\epsilon_{12}$" AttributeType="Scalar" Center="Cell"><DataItem DataType="Float" Dimensions="55" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data5</DataItem></Attribute><Attribute Name="$\epsilon_{13}$" AttributeType="Scalar" Center="Cell"><DataItem DataType="Float" Dimensions="55" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data6</DataItem></Attribute><Attribute Name="$\epsilon_{23}$" AttributeType="Scalar" Center="Cell"><DataItem DataType="Float" Dimensions="55" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data7</DataItem></Attribute><Attribute Name="$\varphi_{1}$" AttributeType="Scalar" Center="Cell"><DataItem DataType="Float" Dimensions="55" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data8</DataItem></Attribute><Attribute Name="$\Phi$" AttributeType="Scalar" Center="Cell"><DataItem DataType="Float" Dimensions="55" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data9</DataItem></Attribute><Attribute Name="$\varphi_{2}$" AttributeType="Scalar" Center="Cell"><DataItem DataType="Float" Dimensions="55" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data10</DataItem></Attribute><Attribute Name="Phases" AttributeType="Scalar" Center="Cell"><DataItem DataType="Int" Dimensions="55" Format="HDF" Precision="8">fast_polycrystal_from_odf.h5:/data11</DataItem></Attribute></Grid></Domain></Xdmf>
diff --git a/tests/test_detector.py b/tests/test_detector.py
@@ -66,6 +66,9 @@ def test_centroid_render(self):
         scattered_wave_vector = 2 * np.pi * scattered_wave_vector / \
             (np.linalg.norm(scattered_wave_vector) * wavelength)
 
+        zd1,yd1 = tuple(self.detector.get_intersection(scattered_wave_vector,verts1.mean(axis=0)[np.newaxis,:])[0])
+        zd2,yd2 = tuple(self.detector.get_intersection(scattered_wave_vector,verts2.mean(axis=0)[np.newaxis,:])[0])
+
         data = os.path.join(
             os.path.join(
                 os.path.dirname(__file__),
@@ -86,7 +89,10 @@ def test_centroid_render(self):
                                time=0,
                                phase=phase,
                                hkl_indx=0,
-                               element_index=0)
+                               element_index=0,
+                               zd=zd1,
+                               yd=yd1)
+
         scattering_unit2 = ScatteringUnit(ch2,
                                scattered_wave_vector=scattered_wave_vector,
                                incident_wave_vector=incident_wave_vector,
@@ -96,7 +102,9 @@ def test_centroid_render(self):
                                time=0,
                                phase=phase,
                                hkl_indx=0,
-                               element_index=0)
+                               element_index=0,
+                               zd=zd2,
+                               yd=yd2)
 
         self.detector.frames.append([scattering_unit1, scattering_unit2])
         diffraction_pattern = self.detector.render(
@@ -162,6 +170,9 @@ def test_centroid_render_with_scintillator(self):
         scattered_wave_vector = 2 * np.pi * scattered_wave_vector / \
             (np.linalg.norm(scattered_wave_vector) * wavelength)
 
+        zd1,yd1 = tuple(self.detector.get_intersection(scattered_wave_vector,verts1.mean(axis=0)[np.newaxis,:])[0])
+        zd2,yd2 = tuple(self.detector.get_intersection(scattered_wave_vector,verts2.mean(axis=0)[np.newaxis,:])[0])
+
         data = os.path.join(
             os.path.join(
                 os.path.dirname(__file__),
@@ -182,7 +193,9 @@ def test_centroid_render_with_scintillator(self):
                                time=0,
                                phase=phase,
                                hkl_indx=0,
-                               element_index=0)
+                               element_index=0,
+                               zd=zd1,
+                               yd=yd1)
         scattering_unit2 = ScatteringUnit(ch2,
                                scattered_wave_vector=scattered_wave_vector,
                                incident_wave_vector=incident_wave_vector,
@@ -192,8 +205,9 @@ def test_centroid_render_with_scintillator(self):
                                time=0,
                                phase=phase,
                                hkl_indx=0,
-                               element_index=0)
-
+                               element_index=0,
+                               zd=zd2,
+                               yd=yd2)
         self.detector.frames.append([scattering_unit1, scattering_unit2])
         self.detector.point_spread_kernel_shape = (3,3)
         diffraction_pattern = self.detector.render(
@@ -208,6 +222,8 @@ def test_centroid_render_with_scintillator(self):
         verts1 = np.array([[0, 0, 0], [0, 0, 1], [0, 1, 0], [1, 0, 0]]) + \
             v * np.sqrt(2) * self.detector_size / 2.  + self.pixel_size_z*0.1 # tetra at detector center perturbed
         ch1 = ConvexHull(verts1)
+        zd1,yd1 = tuple(self.detector.get_intersection(scattered_wave_vector,verts1.mean(axis=0)[np.newaxis,:])[0])
+
         scattering_unit1 = ScatteringUnit(ch1,
                                scattered_wave_vector=scattered_wave_vector,
                                incident_wave_vector=incident_wave_vector,
@@ -217,7 +233,9 @@ def test_centroid_render_with_scintillator(self):
                                time=0,
                                phase=phase,
                                hkl_indx=0,
-                               element_index=0)
+                               element_index=0,
+                               zd=zd1,
+                               yd=yd1)
         self.detector.frames[-1] = [scattering_unit1, scattering_unit2]
         diffraction_pattern_2 = self.detector.render(
             frames_to_render=0,
@@ -232,6 +250,8 @@ def test_centroid_render_with_scintillator(self):
         v = v / np.linalg.norm(v)
         verts1 = np.array([[0, 0, 0], [0, 0, 1], [0, 1, 0], [1, 0, 0]]) - np.array([0, 1.5, 2.0 ])*self.pixel_size_z
         ch1 = ConvexHull(verts1)
+        zd1,yd1 = tuple(self.detector.get_intersection(scattered_wave_vector,verts1.mean(axis=0)[np.newaxis,:])[0])
+
         scattering_unit1 = ScatteringUnit(ch1,
                                scattered_wave_vector=scattered_wave_vector,
                                incident_wave_vector=incident_wave_vector,
@@ -241,7 +261,9 @@ def test_centroid_render_with_scintillator(self):
                                time=0,
                                phase=phase,
                                hkl_indx=0,
-                               element_index=0)
+                               element_index=0,
+                               zd=zd1,
+                               yd=yd1)
         self.detector.frames[-1] = [scattering_unit1, scattering_unit2]
         diffraction_pattern_3 = self.detector.render(
             frames_to_render=0,
@@ -255,6 +277,8 @@ def test_centroid_render_with_scintillator(self):
         v = v / np.linalg.norm(v)
         verts1 = np.array([[0, 0, 0], [0, 0, 1], [0, 1, 0], [1, 0, 0]]) + np.array([0, 246*self.pixel_size_y, 196   *self.pixel_size_z ])
         ch1 = ConvexHull(verts1)
+        zd1,yd1 = tuple(self.detector.get_intersection(scattered_wave_vector,verts1.mean(axis=0)[np.newaxis,:])[0])
+
         scattering_unit1 = ScatteringUnit(ch1,
                                scattered_wave_vector=scattered_wave_vector,
                                incident_wave_vector=incident_wave_vector,
@@ -264,7 +288,9 @@ def test_centroid_render_with_scintillator(self):
                                time=0,
                                phase=phase,
                                hkl_indx=0,
-                               element_index=0)
+                               element_index=0,
+                               zd=zd1,
+                               yd=yd2)
         self.detector.frames[-1] = [scattering_unit1, scattering_unit2]
         diffraction_pattern_3 = self.detector.render(
             frames_to_render=0,
@@ -362,6 +388,8 @@ def test_projection_render(self):
         phase.setup_diffracting_planes(
             wavelength, 0, 20 * np.pi / 180)
 
+        zd1,yd1 = tuple(self.detector.get_intersection(scattered_wave_vector,hull_points.mean(axis=0)[np.newaxis,:])[0])
+
         scattering_unit = ScatteringUnit(sphere_hull,
                               scattered_wave_vector=scattered_wave_vector,
                               incident_wave_vector=incident_wave_vector,
@@ -371,7 +399,9 @@ def test_projection_render(self):
                               time=0,
                               phase=phase,
                               hkl_indx=0,
-                              element_index=0)
+                              element_index=0,
+                              zd=zd1,
+                              yd=yd1)
 
         self.detector.frames.append([scattering_unit])
         diffraction_pattern = self.detector.render(
@@ -429,7 +459,8 @@ def test_get_intersection(self):
         # central normal algined ray
         ray_direction = np.array([2.23, 0., 0.])
         source_point = np.array([0., 0., 0.])
-        z, y = self.detector.get_intersection(ray_direction, source_point)
+
+        z, y = tuple(self.detector.get_intersection(ray_direction, source_point[np.newaxis,:])[0])
         self.assertAlmostEqual(
             z, 0, msg="central detector-normal algined ray does not intersect at 0")
         self.assertAlmostEqual(
@@ -438,7 +469,7 @@ def test_get_intersection(self):
         # translate the ray
         source_point += self.detector.ydhat * self.pixel_size_y
         source_point -= self.detector.zdhat * 2 * self.pixel_size_z
-        z, y = self.detector.get_intersection(ray_direction, source_point)
+        z, y = tuple(self.detector.get_intersection(ray_direction, source_point[np.newaxis,:])[0])
         self.assertAlmostEqual(
             z,
             -2 * self.pixel_size_z,
@@ -452,7 +483,7 @@ def test_get_intersection(self):
         ang = np.arctan(self.pixel_size_y / self.detector_size)
         frac = np.tan(ang) * np.linalg.norm(ray_direction)
         ray_direction += self.detector.ydhat * frac * 3
-        z, y = self.detector.get_intersection(ray_direction, source_point)
+        z, y = tuple(self.detector.get_intersection(ray_direction, source_point[np.newaxis,:])[0])
         self.assertAlmostEqual(
             z, -2 * self.pixel_size_z, msg="translated and tilted ray does not intersect properly")
         self.assertAlmostEqual(
@@ -578,7 +609,7 @@ def test_eta0_render(self):
         phase = Phase(unit_cell, sgname, path_to_cif_file=data)
         phase.setup_diffracting_planes(
             wavelength, 0, 20 * np.pi / 180)
-
+        zd,yd = tuple(self.detector.get_intersection(scattered_wave_vector,verts1.mean(axis=0)[np.newaxis,:])[0])
         scattering_unit = ScatteringUnit(ch,
                                scattered_wave_vector=scattered_wave_vector,
                                incident_wave_vector=incident_wave_vector,
@@ -588,7 +619,9 @@ def test_eta0_render(self):
                                time=0,
                                phase=phase,
                                hkl_indx=0,
-                               element_index=0)
+                               element_index=0,
+                               zd=zd,
+                               yd=yd)
 
         self.detector.frames.append([scattering_unit])