From b75249daadc4c996078c73ae60bab7d17ea6828d Mon Sep 17 00:00:00 2001 From: Duc Le Date: Thu, 28 Dec 2023 15:06:18 +0000 Subject: [PATCH] Bugfixes Dec 2023 (#3) * Fix mantid bug 34213 - show data ascii changes internal state * Fix mantid feat req 36462 - add checks on max/min Ei for Gd choppers * Comment out unused path in Chop.py * Refactor Chop.py to remove unnecessary import * Add new numerics test * Remove unused resolution polynomial approx in Instrument.calculate() * Remove scipy imports for web interface --- PyChop/Chop.py | 80 +++++++++++++++++++++---------------------- PyChop/Instruments.py | 56 +++++++++++++++--------------- PyChop/PyChopGui.py | 2 +- PyChop/mari.yaml | 1 + PyChop/merlin.yaml | 1 + tests/PyChopTest.py | 34 ++++++++++++++++++ 6 files changed, 105 insertions(+), 69 deletions(-) diff --git a/PyChop/Chop.py b/PyChop/Chop.py index 775ebb7..bbee2a9 100644 --- a/PyChop/Chop.py +++ b/PyChop/Chop.py @@ -24,7 +24,6 @@ """ import numpy as np -import collections import warnings warnings.simplefilter("always", UserWarning) @@ -45,24 +44,24 @@ def tchop(freq, Ei, pslit, radius, rho): veloc = 437.3920 * np.sqrt(Ei) gamm = (2.00 * (R**2) / p) * abs(1.00 / rho - 2.00 * w / veloc) # Find regime and calculate variance: - if hasattr(gamm, "__len__"): - tausqr = np.zeros(len(gamm)) - pre = (p / (2.00 * R * w)) ** 2 / 6.00 - idx = np.where((gamm <= 1.0)) - tausqr[idx] = pre * (1.00 - (gamm[idx] ** 2) ** 2 / 10.00) / (1.00 - (gamm[idx] ** 2) / 6.00) - idx = np.where((gamm > 1.0) * (gamm < 4.0)) - groot = np.sqrt(gamm[idx]) - # area[idx] = pre * 0.60 * gamm[idx] * ((groot-2.00)**2) * (groot+8.00) / (groot+4.00) + #if hasattr(gamm, "__len__"): + # tausqr = np.zeros(len(gamm)) + # pre = (p / (2.00 * R * w)) ** 2 / 6.00 + # idx = np.where((gamm <= 1.0)) + # tausqr[idx] = pre * (1.00 - (gamm[idx] ** 2) ** 2 / 10.00) / (1.00 - (gamm[idx] ** 2) / 6.00) + # idx = np.where((gamm > 1.0) * (gamm < 4.0)) + # groot = np.sqrt(gamm[idx]) + # tausqr[idx] = pre * 0.60 * gamm[idx] * ((groot-2.00)**2) * (groot+8.00) / (groot+4.00) + #else: + if gamm >= 4.00: + warnings.warn("PyChop: tchop(): No transmission at %5.3f meV at %3d Hz" % (Ei, freq)) + return np.nan + if gamm <= 1.00: + gsqr = (1.00 - (gamm**2) ** 2 / 10.00) / (1.00 - (gamm**2) / 6.00) else: - if gamm >= 4.00: - warnings.warn("PyChop: tchop(): No transmission at %5.3f meV at %3d Hz" % (Ei, freq)) - return np.nan - if gamm <= 1.00: - gsqr = (1.00 - (gamm**2) ** 2 / 10.00) / (1.00 - (gamm**2) / 6.00) - else: - groot = np.sqrt(gamm) - gsqr = 0.60 * gamm * ((groot - 2.00) ** 2) * (groot + 8.00) / (groot + 4.00) - tausqr = ((p / (2.00 * R * w)) ** 2 / 6.00) * gsqr + groot = np.sqrt(gamm) + gsqr = 0.60 * gamm * ((groot - 2.00) ** 2) * (groot + 8.00) / (groot + 4.00) + tausqr = ((p / (2.00 * R * w)) ** 2 / 6.00) * gsqr return tausqr @@ -86,25 +85,25 @@ def achop(Ei, freq, dslat, pslit, radius, rho): vela = 437.3920 * np.sqrt(Ei) gamm = (2.00 * (R1**2) / p1) * abs(1.00 / rho1 - 2.00 * w1 / vela) # Find regime and calculate variance: - if hasattr(gamm, "__len__"): - area = np.zeros(len(gamm)) - pre = (p1**2) / (2.00 * R1 * w1) - idx = np.where(gamm <= 1.0) - area[idx] = pre * (1.0 - (gamm[idx] ** 2) / 6.0) - idx = np.where((gamm > 1.0) * (gamm < 4.0)) - groot = np.sqrt(gamm[idx]) - area[idx] = pre * groot * ((groot - 2.0) ** 2) * (groot + 4.0) / 6.0 + #if hasattr(gamm, "__len__"): + # area = np.zeros(len(gamm)) + # pre = (p1**2) / (2.00 * R1 * w1) + # idx = np.where(gamm <= 1.0) + # area[idx] = pre * (1.0 - (gamm[idx] ** 2) / 6.0) + # idx = np.where((gamm > 1.0) * (gamm < 4.0)) + # groot = np.sqrt(gamm[idx]) + # area[idx] = pre * groot * ((groot - 2.0) ** 2) * (groot + 4.0) / 6.0 + #else: + if gamm >= 4.00: + warnings.warn("PyChop: achop(): No transmission at %5.3f meV at %3d Hz" % (Ei, freq), UserWarning) + return np.nan else: - if gamm >= 4.00: - warnings.warn("PyChop: achop(): No transmission at %5.3f meV at %3d Hz" % (Ei, freq), UserWarning) - return np.nan + if gamm <= 1.00: + f1 = 1.00 - (gamm**2) / 6.00 else: - if gamm <= 1.00: - f1 = 1.00 - (gamm**2) / 6.00 - else: - groot = np.sqrt(gamm) - f1 = groot * ((groot - 2.00) ** 2) * (groot + 4.00) / 6.00 - area = ((p1**2) / (2.00 * R1 * w1)) * f1 + groot = np.sqrt(gamm) + f1 = groot * ((groot - 2.00) ** 2) * (groot + 4.00) / 6.00 + area = ((p1**2) / (2.00 * R1 * w1)) * f1 return area @@ -462,13 +461,12 @@ def sam0(sx, sy, sz, isam): # RAE code (assumes plate sample, so correct for MAPS vanadium) # A more sophisticated version would do different things depending on sample type but usually # this contribution is small, and in any case this will be close enough for most geometries + def sample_shape_scaling_factors(typ): + # Sample type: 0==flat plate, 1==ellipse, 2==annulus, 3==sphere, 4==solid cylinder + # Factor is only correct for plate (1/12) and annulus (1/8) + return 1. / 8 if typ == 2 else 1. / 12 varx = 0 # vary = ((sx)**2 + (sy)**2) # WRONG - vary = (sy**2) * sample_shape_scaling_factors[isam] + vary = (sy**2) * sample_shape_scaling_factors(isam) varz = 0 return varx, vary, varz - - -# Sample type: 0==flat plate, 1==ellipse, 2==annulus, 3==sphere, 4==solid cylinder -sample_shape_scaling_factors = collections.defaultdict(lambda: 1.0 / 12) -sample_shape_scaling_factors[2] = 1.0 / 8 diff --git a/PyChop/Instruments.py b/PyChop/Instruments.py index 7032190..126d114 100644 --- a/PyChop/Instruments.py +++ b/PyChop/Instruments.py @@ -14,17 +14,26 @@ import warnings import copy from . import Chop, MulpyRep -from scipy.interpolate import interp1d -from scipy.special import erf -from scipy import constants -from scipy.optimize import curve_fit +from math import erf # Some global constants SIGMA2FWHM = 2 * np.sqrt(2 * np.log(2)) SIGMA2FWHMSQ = SIGMA2FWHM**2 -E2V = np.sqrt((constants.e / 1000) * 2 / constants.neutron_mass) # v = E2V * sqrt(E) veloc in m/s, E in meV -E2L = 1.0e23 * constants.h**2 / (2 * constants.m_n * constants.e) # lam = sqrt(E2L / E) lam in Angst, E in meV -E2K = constants.e * 2 * constants.m_n / constants.hbar**2 / 1e23 # k = sqrt(E2K * E) k in 1/Angst, E in meV +#from scipy import constants +#E2V = np.sqrt((constants.e / 1000) * 2 / constants.neutron_mass) # v = E2V * sqrt(E) veloc in m/s, E in meV +#E2L = 1.0e23 * constants.h**2 / (2 * constants.m_n * constants.e) # lam = sqrt(E2L / E) lam in Angst, E in meV +#E2K = constants.e * 2 * constants.m_n / constants.hbar**2 / 1e23 # k = sqrt(E2K * E) k in 1/Angst, E in meV +E2V = 437.393362604208619 +E2L = 81.8042103582802156 +E2K = 0.48259640220781652 + +try: + from scipy.interpolate import interp1d +except ModuleNotFoundError: + def interp1d(xp, yp, **kwargs): + # Emulates scipy's interp1d using numpy's linear interpolation (ignore 'kind' kwargs) + # Note that this is only accurate to about 1% compared with scipy's spline interpolation + return lambda x: np.interp(x, xp, yp) def wrap_attributes(obj, inval, allowed_var_names): @@ -81,7 +90,7 @@ class FermiChopper(object): Class which represents a Fermi chopper package """ - __allowed_var_names = ["name", "pslit", "pslat", "radius", "rho", "tjit", "fluxcorr", "isPi"] + __allowed_var_names = ["name", "pslit", "pslat", "radius", "rho", "tjit", "fluxcorr", "isPi", "ei_limits"] def __init__(self, inval=None): wrap_attributes(self, inval, self.__allowed_var_names) @@ -101,6 +110,14 @@ def getTransmission(self, Ei, freq): dslat = (self.pslit + self.pslat) / 1000 return Chop.achop(Ei, freq, dslat, self.pslit / 1000.0, self.radius / 1000.0, self.rho / 1000.0) / dslat + @property + def emin(self): + return self.ei_limits[0] if hasattr(self, "ei_limits") else 0.1 + + @property + def emax(self): + return self.ei_limits[1] if hasattr(self, "ei_limits") else 1000. + class ChopperSystem(object): """ @@ -263,6 +280,10 @@ def getFrequency(self): def setEi(self, Ei): """Sets the (focussed) incident energy""" + emin = max(self.emin, self.packages[self.package].emin if self.isFermi else 0) + emax = min(self.emax, self.packages[self.package].emax if self.isFermi else np.inf) + if Ei < emin or Ei > emax: + raise ValueError(f'Ei={Ei} is outside limits [{emin}, {emax}]') self.ei = Ei def getEi(self): @@ -941,13 +962,6 @@ def calculate(cls, *args, **kwargs): ! Instrument.calculate('merlin', 'g', 450., 60., range(55)) ! Instrument.calculate('maps', 'a', 450., 600., etrans=np.linspace(0,550,55)) ! - ! For fast calculations, one can return a polynomial approximation (cubic) of the - ! resolution function. By passing etrans='polynomial', the calculator estimates the - ! resolution for etrans=np.arange(-Ei, Ei, Ei*0.01) then fits it to a cubic polynomial. - ! The resolution is then an array with coefficients, from the lowest power. - ! - ! res, flux = Instrument.calculate(inst='cncs', variant='High flux', freq=240, ei=1.5, etrans='polynomial') - ! ! The results are returned as tuple: (resolution, flux) """ argdict = argparser(args, kwargs, ["inst", "package", "frequency", "ei", "etrans", "variant"]) @@ -959,26 +973,14 @@ def calculate(cls, *args, **kwargs): if argdict["variant"]: obj.variant = argdict["variant"] etrans = argdict["etrans"] - return_polynomial = False if etrans is None: etrans = 0.0 else: - if etrans == "polynomial": - return_polynomial = True - etrans = np.arange(-obj.ei, obj.ei, obj.ei * 0.01) try: etrans = float(etrans) except TypeError: etrans = np.asfarray(etrans) res = obj.getResolution(etrans) - - if return_polynomial: - - def cubic(x, x_0, x_1, x_2, x_3): - return x_0 + x_1 * x + x_2 * x**2 + x_3 * x**3 - - popt, pcov = curve_fit(cubic, etrans, res) - res = popt flux = obj.getFlux() return res, flux diff --git a/PyChop/PyChopGui.py b/PyChop/PyChopGui.py index 2b98078..4d7131e 100644 --- a/PyChop/PyChopGui.py +++ b/PyChop/PyChopGui.py @@ -609,7 +609,7 @@ def plot_frame(self): self.repcanvas.draw() def _gen_text_ei(self, ei, obj_in): - obj = Instrument(obj_in) + obj = copy.deepcopy(obj_in) obj.setEi(ei) en = np.linspace(0, 0.95 * ei, 10) # ValueErrors here will be caught in showText() or writeText() diff --git a/PyChop/mari.yaml b/PyChop/mari.yaml index f479877..d6c58c4 100644 --- a/PyChop/mari.yaml +++ b/PyChop/mari.yaml @@ -65,6 +65,7 @@ chopper_system: tjit: 0.0 # Jitter time (us) fluxcorr: 3.2 # (Empirical/Fudge) factor to scale calculated flux by isPi: True # Should the PI pulse (at 180 deg rotation) be transmitted by this package? + ei_limits: [0, 181] # Limits on ei for this chopper (setting Ei outside this will give error) R: name: MARI R (500meV) pslit: 1.143 # Neutron transparent slit width (mm) diff --git a/PyChop/merlin.yaml b/PyChop/merlin.yaml index 1f8f320..87fc17b 100644 --- a/PyChop/merlin.yaml +++ b/PyChop/merlin.yaml @@ -37,6 +37,7 @@ chopper_system: tjit: 0.0 # Jitter time (us) fluxcorr: 3. # (Empirical/Fudge) factor to scale calculated flux by isPi: True # Should the PI pulse (at 180 deg rotation) be transmitted by this package? + ei_limits: [0, 181] # Limits on ei for this chopper (setting Ei outside this will give error) S: name: MERLIN S (Sloppy) pslit: 2.280 # Neutron transparent slit width (mm) diff --git a/tests/PyChopTest.py b/tests/PyChopTest.py index 8cdbf64..f294303 100644 --- a/tests/PyChopTest.py +++ b/tests/PyChopTest.py @@ -87,6 +87,21 @@ def test_pychop_invalid_ei(self): assert "Cannot calculate for energy transfer greater than Ei" in str(w[0].message) assert np.isnan(res[0]) + def test_pychop_numerics(self): + # Tests all instruments resolution and flux values against reference + instruments = ["ARCS", "CNCS", "HYSPEC", "LET", "MAPS", "MARI", "MERLIN", "SEQUOIA"] + choppers = ["ARCS-100-1.5-AST", "High Flux", "OnlyOne", "High Flux", "S", "S", "G", "SEQ-100-2.0-AST"] + freqs = [ [300], [300, 60], [180], [240, 120], [400, 50], [400], [400], [300] ] + eis = [120, 3.7, 45, 3.7, 120, 80, 120, 120] + ref_res = [10.278744237772832, 0.13188102618129077, 3.6751279831313703, 0.08079912729715726, + 4.9611687063450995, 2.6049587487601764, 6.8755979524827255, 5.396705255853653] + ref_flux = [2055.562054927915, 128986.24972543867, 0.014779264739956933, 45438.33797146135, + 24196.496233770937, 5747.118187298609, 22287.647098883135, 4063.3113893387676] + for inst, ch, frq, ei, res0, flux0 in zip(instruments, choppers, freqs, eis, ref_res, ref_flux): + res, flux = PyChop2.calculate(inst, ch, frq, ei, 0) + np.testing.assert_allclose(res[0], res0, rtol=1e-7, atol=0) + np.testing.assert_allclose(flux[0], flux0, rtol=1e-7, atol=0) + class MockedModule(mock.MagicMock): # A class which is meant to act like a module @@ -232,6 +247,22 @@ def __getattr__(self, attribute): # noqa: E306 self.__dict__[attribute] = mock.MagicMock() return self.__dict__[attribute] + class fake_QLineEdit(mock.MagicMock): # noqa: E306 + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + self.value = '1.0' + + def setText(self, value): # noqa: E306 + self.value = value + + def text(self): # noqa: E306 + return self.value + + def __getattr__(self, attribute): # noqa: E306 + if attribute not in self.__dict__: + self.__dict__[attribute] = mock.MagicMock() + return self.__dict__[attribute] + class fake_Line(mock.MagicMock): # noqa: E306 def __init__(self, parent, *args, **kwargs): super().__init__(*args, **kwargs) @@ -267,6 +298,7 @@ def __init__(self, parent, label, valmin, valmax, **kwargs): cls.mock_modules = MockImports(include=['qtpy', 'matplotlib', 'mantidqt', 'mantid.plots'], replace={'QMainWindow':fake_QMainWindow, 'QComboBox':MockedModule(mock_class=fake_QCombo), + 'QLineEdit':MockedModule(mock_class=fake_QLineEdit), 'Figure':MockedModule(mock_class=fake_Figure), 'Slider':MockedModule(mock_class=fake_Slider)}) # Mess around with import mechanism _inside_ PyChopGui so GUI libs not really imported @@ -282,6 +314,8 @@ def test_hyspec(self): self.window.calc_callback() setS2.assert_not_called() self.window.setInstrument('HYSPEC') + self.window.widgets['EiEdit']['Edit'].setText('50') + self.window.setEi() self.window.calc_callback() setS2.assert_called() # Test that the value of S2 is set correctly