From 488dd65470551c811050358db67cfadda665585b Mon Sep 17 00:00:00 2001 From: edgar1993a Date: Thu, 5 Dec 2024 14:43:31 +0100 Subject: [PATCH 1/7] fiberunit numexpr --- src/pyFAI/units.py | 60 +++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 59 insertions(+), 1 deletion(-) diff --git a/src/pyFAI/units.py b/src/pyFAI/units.py index 2a9c5a3ea..eef632438 100644 --- a/src/pyFAI/units.py +++ b/src/pyFAI/units.py @@ -160,7 +160,7 @@ def __init__(self, name, scale=1, label=None, equation=None, formula=None, scale=scale, label=label, equation=equation, - formula=formula, + # formula=formula, center=center, corner=corner, delta=delta, @@ -169,9 +169,48 @@ def __init__(self, name, scale=1, label=None, equation=None, formula=None, positive=positive, period=period, ) + self.formula = formula + self.formula_so1 = formula self._incident_angle = incident_angle self._tilt_angle = tilt_angle self._sample_orientation = sample_orientation + self._update_ne_equation() + + def _update_ne_equation(self): + if (numexpr is not None) and isinstance(self.formula, str): + signature = [(key, numpy.float64) for key in "xyzλπηχ" if key in self.formula] + + if self._sample_orientation == 1: + ... + elif self._sample_orientation == 2: + self.formula = self.formula_so1.replace('x', 'y_') + self.formula = self.formula_so1.replace('y', 'x_') + elif self._sample_orientation == 3: + self.formula = self.formula_so1.replace('x', '(-1)*x_') + elif self._sample_orientation == 4: + self.formula = self.formula_so1.replace('x', '(-1)*y_') + self.formula = self.formula_so1.replace('y', '(-1)*x_') + self.formula = self.formula_so1.replace('x_', 'x') + self.formula = self.formula_so1.replace('y_', 'y') + + ne_formula = numexpr.NumExpr(self.formula, signature) + + def ne_equation(x, y, z=None, wavelength=None, + incident_angle=self._incident_angle, + tilt_angle=self._tilt_angle, + sample_orientation=self._sample_orientation, + ne_formula=ne_formula): + π = numpy.pi + λ = wavelength + η = self._incident_angle + χ = self._tilt_angle + ldict = locals() + args = tuple(ldict[i] for i in ne_formula.input_names) + return ne_formula(*args) + + self.equation = ne_equation + else: + self.equation = self._equation def __repr__(self): return f""" @@ -195,12 +234,15 @@ def sample_orientation(self): def set_incident_angle(self, value:float): self._incident_angle = value + self._update_ne_equation() def set_tilt_angle(self, value:float): self._tilt_angle = value + self._update_ne_equation() def set_sample_orientation(self, value: int): self._sample_orientation = value + self._update_ne_equation() RADIAL_UNITS = {} @@ -618,6 +660,17 @@ def eq_q_total(x, y, z, wavelength, incident_angle=0.0, tilt_angle=0.0, sample_o formula_qx = "4.0e-9*π/λ*sin(arctan2(x, z)/2.0)" formula_qy = "4.0e-9*π/λ*sin(arctan2(y, z)/2.0)" +formula_exit_angle = "arctan2(y,sqrt(z**2+x**2))" +formula_exit_angle_horz = "arctan2(x,z)" +formula_qbeam_lab = f"2.0e-9/λ*π*(cos({formula_exit_angle})*cos({formula_exit_angle_horz}) - 1)" +formula_qhorz_lab = f"2.0e-9/λ*π*cos({formula_exit_angle})*sin({formula_exit_angle_horz})" +formula_qvert_lab = f"2.0e-9/λ*π*sin({formula_exit_angle})" +formula_qbeam_rot = f"cos(η)*({formula_qbeam_lab})+sin(η)*({formula_qvert_lab})" +formula_qhorz_rot = f"cos(χ)*({formula_qhorz_lab})-sin(χ)*sin(η)*({formula_qbeam_lab})+sin(χ)*cos(η)*({formula_qvert_lab})" +formula_qvert_rot = f"-sin(χ)*({formula_qhorz_lab})-cos(χ)*sin(η)*({formula_qbeam_lab})+cos(χ)*cos(η)*({formula_qvert_lab})" +formula_qip = f"sqrt(({formula_qbeam_rot})**2+({formula_qhorz_rot})**2)*((({formula_qhorz_rot} > 0) * 2) - 1)" +formula_qoop = formula_qvert_rot + register_radial_unit("r_mm", center="rArray", delta="deltaR", @@ -810,6 +863,7 @@ def eq_q_total(x, y, z, wavelength, incident_angle=0.0, tilt_angle=0.0, sample_o register_radial_fiber_unit("qxgi_nm^-1", scale=1.0, label=r"Scattering vector $q_x$ ($nm^{-1}$)", + formula=formula_qhorz_rot, equation=eq_qhorz_gi, short_name="qxgi", unit_symbol="nm^{-1}", @@ -818,6 +872,7 @@ def eq_q_total(x, y, z, wavelength, incident_angle=0.0, tilt_angle=0.0, sample_o register_radial_fiber_unit("qygi_nm^-1", scale=1.0, label=r"Scattering vector $q_y$ ($nm^{-1}$)", + formula=formula_qvert_rot, equation=eq_qvert_gi, short_name="qygi", unit_symbol="nm^{-1}", @@ -826,6 +881,7 @@ def eq_q_total(x, y, z, wavelength, incident_angle=0.0, tilt_angle=0.0, sample_o register_radial_fiber_unit("qzgi_nm^-1", scale=1.0, label=r"Scattering vector $q_z$ ($nm^{-1}$)", + formula=formula_qbeam_rot, equation=eq_qbeam_gi, short_name="qzgi", unit_symbol="nm^{-1}", @@ -834,6 +890,7 @@ def eq_q_total(x, y, z, wavelength, incident_angle=0.0, tilt_angle=0.0, sample_o register_radial_fiber_unit("qip_nm^-1", scale=1.0, label=r"Scattering vector $q_{IP}$ ($nm^{-1}$)", + formula=formula_qip, equation=eq_qip, short_name="qip", unit_symbol="nm^{-1}", @@ -842,6 +899,7 @@ def eq_q_total(x, y, z, wavelength, incident_angle=0.0, tilt_angle=0.0, sample_o register_radial_fiber_unit("qoop_nm^-1", scale=1.0, label=r"Scattering vector $q_{OOP}$ ($nm^{-1}$)", + formula=formula_qoop, equation=eq_qoop, short_name="qoop", unit_symbol="nm^{-1}", From 98c72f339e0578e5991dfe08f41d3a58733a9fdf Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 5 Dec 2024 13:46:05 +0000 Subject: [PATCH 2/7] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/pyFAI/units.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/pyFAI/units.py b/src/pyFAI/units.py index eef632438..87ebf3d1f 100644 --- a/src/pyFAI/units.py +++ b/src/pyFAI/units.py @@ -195,10 +195,10 @@ def _update_ne_equation(self): ne_formula = numexpr.NumExpr(self.formula, signature) - def ne_equation(x, y, z=None, wavelength=None, - incident_angle=self._incident_angle, + def ne_equation(x, y, z=None, wavelength=None, + incident_angle=self._incident_angle, tilt_angle=self._tilt_angle, - sample_orientation=self._sample_orientation, + sample_orientation=self._sample_orientation, ne_formula=ne_formula): π = numpy.pi λ = wavelength From 0cc5231acd630bd170d0c8c7bb8660b2267f5a0f Mon Sep 17 00:00:00 2001 From: edgar1993a Date: Thu, 5 Dec 2024 17:30:29 +0100 Subject: [PATCH 3/7] sample orientation --- src/pyFAI/units.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/pyFAI/units.py b/src/pyFAI/units.py index eef632438..cf6fa44ab 100644 --- a/src/pyFAI/units.py +++ b/src/pyFAI/units.py @@ -180,19 +180,20 @@ def _update_ne_equation(self): if (numexpr is not None) and isinstance(self.formula, str): signature = [(key, numpy.float64) for key in "xyzλπηχ" if key in self.formula] + formula_ = self.formula if self._sample_orientation == 1: ... elif self._sample_orientation == 2: - self.formula = self.formula_so1.replace('x', 'y_') - self.formula = self.formula_so1.replace('y', 'x_') + formula_ = self.formula_so1 + formula_ = formula_.replace('x', 'ψ').replace('y', 'ξ') + formula_ = formula_.replace('ψ', 'y').replace('ξ', 'x') elif self._sample_orientation == 3: - self.formula = self.formula_so1.replace('x', '(-1)*x_') + formula_ = self.formula_so1.replace('x', '(-1)*x') elif self._sample_orientation == 4: - self.formula = self.formula_so1.replace('x', '(-1)*y_') - self.formula = self.formula_so1.replace('y', '(-1)*x_') - self.formula = self.formula_so1.replace('x_', 'x') - self.formula = self.formula_so1.replace('y_', 'y') - + formula_ = self.formula_so1 + formula_ = formula_.replace('x', 'ψ').replace('y', 'ξ') + formula_ = formula_.replace('ψ', '(-1)*y').replace('ξ', '(-1)*x') + self.formula = formula_ ne_formula = numexpr.NumExpr(self.formula, signature) def ne_equation(x, y, z=None, wavelength=None, From 3a3dc4ff40ba06856a4e8814f97d26c73e49d818 Mon Sep 17 00:00:00 2001 From: edgar1993a Date: Thu, 5 Dec 2024 17:37:24 +0100 Subject: [PATCH 4/7] correction --- src/pyFAI/units.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pyFAI/units.py b/src/pyFAI/units.py index cf6fa44ab..d78a16df0 100644 --- a/src/pyFAI/units.py +++ b/src/pyFAI/units.py @@ -188,11 +188,11 @@ def _update_ne_equation(self): formula_ = formula_.replace('x', 'ψ').replace('y', 'ξ') formula_ = formula_.replace('ψ', 'y').replace('ξ', 'x') elif self._sample_orientation == 3: - formula_ = self.formula_so1.replace('x', '(-1)*x') + formula_ = self.formula_so1.replace('x', '(-x)') elif self._sample_orientation == 4: formula_ = self.formula_so1 formula_ = formula_.replace('x', 'ψ').replace('y', 'ξ') - formula_ = formula_.replace('ψ', '(-1)*y').replace('ξ', '(-1)*x') + formula_ = formula_.replace('ψ', '(-y)').replace('ξ', '(-x)') self.formula = formula_ ne_formula = numexpr.NumExpr(self.formula, signature) From 2b2517eed04527fd8b2b35bcb025371e1f8ee7d2 Mon Sep 17 00:00:00 2001 From: edgar1993a Date: Thu, 5 Dec 2024 18:30:57 +0100 Subject: [PATCH 5/7] test_numexpr --- src/pyFAI/test/test_fiber_integrator.py | 46 +++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/src/pyFAI/test/test_fiber_integrator.py b/src/pyFAI/test/test_fiber_integrator.py index 013815fcd..4e70340af 100644 --- a/src/pyFAI/test/test_fiber_integrator.py +++ b/src/pyFAI/test/test_fiber_integrator.py @@ -542,6 +542,52 @@ def test_sample_orientation_equivalence(self): self.assertLess((abs(res_so_4.intensity) - abs(res_so_3.intensity)).max(), threshold) + def test_numexpr_equations(self): + incident_angle = 0.2 + tilt_angle = -1.0 + sample_orientation = 3 + + qhorz = get_unit_fiber(name="qxgi_nm^-1", incident_angle=incident_angle, tilt_angle=tilt_angle, sample_orientation=sample_orientation) + qvert = get_unit_fiber(name="qygi_nm^-1", incident_angle=incident_angle, tilt_angle=tilt_angle, sample_orientation=sample_orientation) + qbeam = get_unit_fiber(name="qzgi_nm^-1", incident_angle=incident_angle, tilt_angle=tilt_angle, sample_orientation=sample_orientation) + qip = get_unit_fiber(name="qip_nm^-1", incident_angle=incident_angle, tilt_angle=tilt_angle, sample_orientation=sample_orientation) + qoop = get_unit_fiber(name="qoop_nm^-1", incident_angle=incident_angle, tilt_angle=tilt_angle, sample_orientation=sample_orientation) + + self.fi.reset() + arr_qhorz_fast = self.fi.array_from_unit(unit=qhorz) + arr_qvert_fast = self.fi.array_from_unit(unit=qvert) + arr_qbeam_fast = self.fi.array_from_unit(unit=qbeam) + arr_qip_fast = self.fi.array_from_unit(unit=qip) + arr_qoop_fast = self.fi.array_from_unit(unit=qoop) + res2d_fast = self.fi.integrate2d_grazing_incidence(data=self.data, unit_ip=qip, unit_oop=qoop) + + qhorz.formula = None + qhorz.equation = qhorz._equation + qvert.formula = None + qvert.equation = qvert._equation + qbeam.formula = None + qbeam.equation = qbeam._equation + qip.formula = None + qip.equation = qip._equation + qoop.formula = None + qoop.equation = qoop._equation + + self.fi.reset() + arr_qhorz_slow = self.fi.array_from_unit(unit=qhorz) + arr_qvert_slow = self.fi.array_from_unit(unit=qvert) + arr_qbeam_slow = self.fi.array_from_unit(unit=qbeam) + arr_qip_slow = self.fi.array_from_unit(unit=qip) + arr_qoop_slow = self.fi.array_from_unit(unit=qoop) + res2d_slow = self.fi.integrate2d_grazing_incidence(data=self.data, unit_ip=qip, unit_oop=qoop) + + self.assertAlmostEqual((arr_qhorz_fast - arr_qhorz_slow).max(), 0.0) + self.assertAlmostEqual((arr_qvert_fast - arr_qvert_slow).max(), 0.0) + self.assertAlmostEqual((arr_qbeam_fast - arr_qbeam_slow).max(), 0.0) + self.assertAlmostEqual((arr_qip_fast - arr_qip_slow).max(), 0.0) + self.assertAlmostEqual((arr_qoop_fast - arr_qoop_slow).max(), 0.0) + self.assertAlmostEqual((res2d_fast.intensity - res2d_slow.intensity).max(), 0.0) + self.assertAlmostEqual((res2d_fast.radial - res2d_slow.radial).max(), 0.0) + self.assertAlmostEqual((res2d_fast.azimuthal - res2d_slow.azimuthal).max(), 0.0) def suite(): testsuite = unittest.TestSuite() From 5340608fb747a69b0a108de295b4946e6e747715 Mon Sep 17 00:00:00 2001 From: edgar1993a Date: Thu, 5 Dec 2024 18:52:58 +0100 Subject: [PATCH 6/7] sandbox test --- sandbox/test_grazingincidence.py | 47 ++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) create mode 100644 sandbox/test_grazingincidence.py diff --git a/sandbox/test_grazingincidence.py b/sandbox/test_grazingincidence.py new file mode 100644 index 000000000..8128e2945 --- /dev/null +++ b/sandbox/test_grazingincidence.py @@ -0,0 +1,47 @@ +from pyFAI.units import Unit, get_unit_fiber, to_unit +from pyFAI.integrator.azimuthal import AzimuthalIntegrator +from pyFAI.calibrant import get_calibrant +from pyFAI.integrator.fiber import FiberIntegrator +from pyFAI import detector_factory +import time +from pyFAI.gui.jupyter import subplots, display, plot2d +import matplotlib.pyplot as plt +import numpy +from pyFAI.test.utilstest import UtilsTest +import fabio +from pyFAI import load + +if __name__ == "__main__": + fi_1 = FiberIntegrator(dist=0.1, poni1=0.1, poni2=0.1, detector=detector_factory("Eiger_4M"), wavelength=1e-10) + poni = UtilsTest.getimage("LaB6_5.poni") + fi_2 = load(filename=poni, type_="pyFAI.integrator.fiber.FiberIntegrator") + + cal = get_calibrant("LaB6") + data_1 = cal.fake_calibration_image(ai=fi_1) + data_file = UtilsTest.getimage("Y6.edf") + data_2 = fabio.open(data_file).data + + for fi, data in zip((fi_1, fi_2), (data_1, data_2)): + res2d_1 = fi.integrate2d_grazing_incidence(data=data, incident_angle=0.0, tilt_angle=0.0, sample_orientation=1) + res2d_2 = fi.integrate2d_grazing_incidence(data=data, incident_angle=0.2, tilt_angle=0.0, sample_orientation=1) + res2d_3 = fi.integrate2d_grazing_incidence(data=data, incident_angle=0.0, tilt_angle=-0.2, sample_orientation=1) + res2d_4 = fi.integrate2d_grazing_incidence(data=data, incident_angle=0.2, tilt_angle=-1.54, sample_orientation=1) + + res2d_5 = fi.integrate2d_grazing_incidence(data=data, incident_angle=0.0, tilt_angle=0.0, sample_orientation=2) + res2d_6 = fi.integrate2d_grazing_incidence(data=data, incident_angle=0.0, tilt_angle=0.0, sample_orientation=3) + res2d_7 = fi.integrate2d_grazing_incidence(data=data, incident_angle=0.0, tilt_angle=0.0, sample_orientation=4) + + fig, axes = subplots(ncols=4, nrows=2) + plot2d(res2d_1, ax=axes[0,0]) + plot2d(res2d_2, ax=axes[0,1]) + plot2d(res2d_3, ax=axes[0,2]) + plot2d(res2d_4, ax=axes[0,3]) + plot2d(res2d_5, ax=axes[1,0]) + plot2d(res2d_6, ax=axes[1,1]) + plot2d(res2d_7, ax=axes[1,2]) + for ax in axes.ravel(): + if len(ax.get_images()) == 0: + continue + ax.get_images()[0].set_cmap("viridis") + # ax.get_images()[0].set_clim(0,1500) + plt.show() From 91056ac611c701e367fd57605f0fc41fea31f589 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 5 Dec 2024 17:54:20 +0000 Subject: [PATCH 7/7] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- sandbox/test_grazingincidence.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sandbox/test_grazingincidence.py b/sandbox/test_grazingincidence.py index 8128e2945..46aad7cca 100644 --- a/sandbox/test_grazingincidence.py +++ b/sandbox/test_grazingincidence.py @@ -13,9 +13,9 @@ if __name__ == "__main__": fi_1 = FiberIntegrator(dist=0.1, poni1=0.1, poni2=0.1, detector=detector_factory("Eiger_4M"), wavelength=1e-10) - poni = UtilsTest.getimage("LaB6_5.poni") + poni = UtilsTest.getimage("LaB6_5.poni") fi_2 = load(filename=poni, type_="pyFAI.integrator.fiber.FiberIntegrator") - + cal = get_calibrant("LaB6") data_1 = cal.fake_calibration_image(ai=fi_1) data_file = UtilsTest.getimage("Y6.edf")