Skip to content

Commit

Permalink
Bugfixes Dec 2023 (#3)
Browse files Browse the repository at this point in the history 
* 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
  • Loading branch information
@mducle
mducle authored Dec 28, 2023
1 parent 051af49 commit b75249d
Show file tree
Hide file tree
Showing 6 changed files with 105 additions and 69 deletions.
80 changes: 39 additions & 41 deletions PyChop/Chop.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,6 @@
"""

import numpy as np
import collections
import warnings

warnings.simplefilter("always", UserWarning)
Expand All @@ -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


Expand All @@ -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


Expand Down Expand Up @@ -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
56 changes: 29 additions & 27 deletions PyChop/Instruments.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down Expand Up @@ -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)
Expand All @@ -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):
"""
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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"])
Expand All @@ -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

Expand Down
2 changes: 1 addition & 1 deletion PyChop/PyChopGui.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down
1 change: 1 addition & 0 deletions PyChop/mari.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
1 change: 1 addition & 0 deletions PyChop/merlin.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
34 changes: 34 additions & 0 deletions tests/PyChopTest.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down

0 comments on commit b75249d

Please sign in to comment.