Skip to content
Snippets Groups Projects
Commit 022facbc authored by Martyn Gigg's avatar Martyn Gigg
Browse files

Merge pull request #13496 from mantidproject/13434_loadnmoldyn3ascii_load_any_function 

Allow LoadNMoldyn3Ascii to load 2D functions exported from the plot
parents 4749f34b ca2d5dac
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 239 additions and 195 deletions
Code/Mantid/Framework/PythonInterface/plugins/algorithms/LoadNMoldyn3Ascii.py
+ 218
− 190
View file @ 022facbc
......@@ -4,18 +4,13 @@ from mantid.simpleapi import *
from mantid.kernel import *
from mantid.api import *
import ast
import re
import math
import os
import numpy as np
import math
def _split_line(a):
elements = a.split() # split line on character
extracted = []
for n in elements:
extracted.append(float(n))
return extracted # values as list
#==============================================================================
def _find_starts(data, c, l1):
for l in range(l1, len(data)):
......@@ -25,6 +20,7 @@ def _find_starts(data, c, l1):
break
return line
#==============================================================================
def _find_tab_starts(data, c, l1):
for l in range(l1, len(data)):
......@@ -34,6 +30,7 @@ def _find_tab_starts(data, c, l1):
break
return line
#==============================================================================
def _find_ends(data, c, l1):
for l in range(l1, len(data)):
......@@ -43,15 +40,7 @@ def _find_ends(data, c, l1):
break
return line
def _find_char(data, c, l1):
for l in range(l1, len(data)):
char = data[l]
if char.find(c):
line = l
break
return line
#==============================================================================
def _make_list(a, l1, l2):
data = ''
......@@ -60,21 +49,52 @@ def _make_list(a, l1, l2):
alist = data.split(',')
return alist
#==============================================================================
def _cdl_find_dimensions(data):
"""
Gets the number of Q, time and frequency values in given raw data.
@param data Raw data to search
"""
num_q_values = _find_tab_starts(data, 'NQVALUES', 0)
num_time_values = _find_tab_starts(data, 'NTIMES', 0)
num_freq_values = _find_tab_starts(data, 'NFREQUENCIES', 0)
q_el = data[num_q_values].split()
num_q = int(q_el[2])
t_el = data[num_time_values].split()
num_t = int(t_el[2])
f_el = data[num_freq_values].split()
num_f = int(f_el[2])
logger.debug(data[2][1:-1])
logger.debug(data[3][1:-1])
logger.debug(data[6][1:-1])
return num_q, num_t, num_f
#==============================================================================
class LoadNMoldyn3Ascii(PythonAlgorithm):
_file_name = None
_file_type = None
_functions = None
_out_ws = None
#-------------------------------------------------------------------------------
def category(self):
return 'PythonAlgorithms;Inelastic;Simulation'
#-------------------------------------------------------------------------------
def summary(self):
return 'Imports functions from CDL and ASCII files output by nMOLDYN 3.'
#-------------------------------------------------------------------------------
def PyInit(self):
self.declareProperty(FileProperty('Filename', '',
......@@ -89,6 +109,7 @@ class LoadNMoldyn3Ascii(PythonAlgorithm):
direction=Direction.Output),
doc='Output workspace name')
#-------------------------------------------------------------------------------
def validateInputs(self):
issues = dict()
......@@ -105,124 +126,84 @@ class LoadNMoldyn3Ascii(PythonAlgorithm):
return issues
#-------------------------------------------------------------------------------
def PyExec(self):
# Do setup
self._setup()
# Load data file
data, name, ext = self._load_file()
loaded_ws = None
# Run nMOLDYN import
if ext == 'cdl':
self._cdl_import(data, name)
elif ext == 'dat':
self._ascii_import(data, name)
if self._file_type == 'cdl':
loaded_ws = self._cdl_import()
elif self._file_type == 'dat':
loaders = [self._ascii_3d_import, self._ascii_2d_import]
for loader in loaders:
try:
logger.information('Attempting to load with loader {0}'.format(loader))
loaded_ws = loader()
break
except (ValueError, SyntaxError) as err:
logger.information('Loader {0} failed to load data: {1}'.format(loader, err))
else:
raise RuntimeError('Unrecognised file format: %s' % ext)
raise RuntimeError('Unrecognised file extension: %s' % self._file_type)
if loaded_ws is None:
raise RuntimeError('Failed to load any data, check file format')
# Set the output workspace
self.setProperty('OutputWorkspace', self._out_ws)
self.setProperty('OutputWorkspace', loaded_ws)
#-------------------------------------------------------------------------------
def _setup(self):
"""
Gets algorithm properties.
"""
self._file_name = self.getPropertyValue('Filename')
self._file_type = os.path.splitext(self._file_name)[1][1:]
self._out_ws = self.getPropertyValue('OutputWorkspace')
raw_functions = self.getProperty('Functions').value
self._functions = [x.strip() for x in raw_functions]
#-------------------------------------------------------------------------------
def _load_file(self):
def _cdl_import(self):
"""
Attempts to load the sample file.
@returns A tuple with the ASCII data, sample file name and file extension
Import data from CDL file.
"""
logger.notice('Loading CDL file')
# Get some data about the file
path = self._file_name
base = os.path.basename(path)
name = os.path.splitext(base)[0]
ext = os.path.splitext(path)[1]
# Remove dot from extension
if len(ext) > 1:
ext = ext[1:]
logger.debug('Base filename for %s is %s' % (self._file_name, name))
logger.debug('File type of %s is %s' % (self._file_name, ext))
if not os.path.isfile(path):
path = FileFinder.getFullPath(path)
logger.information('Got file path for %s: %s' % (self._file_name, path))
# Get file base name
base_filename = os.path.basename(self._file_name)
base_name = os.path.splitext(base_filename)[0]
# Open file and get data
handle = open(path, 'r')
data = []
for line in handle:
line = line.rstrip()
data.append(line)
handle.close()
return data, name, ext
def _find_dimensions(self, data):
"""
Gets the number of Q, time and frequency values in given raw data.
@param data Raw data to search
"""
num_q_values = _find_tab_starts(data, 'NQVALUES', 0)
num_time_values = _find_tab_starts(data, 'NTIMES', 0)
num_freq_values = _find_tab_starts(data, 'NFREQUENCIES', 0)
q_el = data[num_q_values].split()
num_q = int(q_el[2])
t_el = data[num_time_values].split()
num_t = int(t_el[2])
f_el = data[num_freq_values].split()
num_f = int(f_el[2])
logger.debug(data[2][1:-1])
logger.debug(data[3][1:-1])
logger.debug(data[6][1:-1])
return num_q, num_t, num_f
def _cdl_import(self, data, name):
"""
Import data from CDL file.
@param data Raw data
@param name Name of data file
"""
logger.notice('Loading CDL file: %s' % name)
with open(self._file_name, 'r') as handle:
for line in handle:
line = line.rstrip()
data.append(line)
len_data = len(data)
# raw head
nQ, nT, nF = self._find_dimensions(data)
num_spec, nT, nF = _cdl_find_dimensions(data)
ldata = _find_starts(data, 'data:', 0)
lq1 = _find_starts(data, ' q =', ldata) # start Q values
lq2 = _find_starts(data, ' q =', lq1 - 1)
Qlist = _make_list(data, lq1, lq2)
if nQ != len(Qlist):
if num_spec != len(Qlist):
raise RUntimeError('Error reading Q values')
Qf = Qlist[0].split()
Q = [float(Qf[2]) / 10.0]
for m in range(1, nQ - 1):
for m in range(1, num_spec - 1):
Q.append(float(Qlist[m]) / 10.0)
Q.append(float(Qlist[nQ - 1][:-1]) / 10.0)
Q.append(float(Qlist[num_spec - 1][:-1]) / 10.0)
logger.information('Q values = ' + str(Q))
lt1 = _find_starts(data, ' time =', lq2) # start T values
......@@ -263,17 +244,17 @@ class LoadNMoldyn3Ascii(PythonAlgorithm):
if func[:3] == 'Fqt':
nP = nT
xEn = np.array(T)
eZero = np.zeros(nT)
xUnit = 'TOF'
zero_error = np.zeros(nT)
x_unit = 'TOF'
elif func[:3] == 'Sqw':
nP = nF
xEn = np.array(F)
eZero = np.zeros(nF)
xUnit = 'Energy'
zero_error = np.zeros(nF)
x_unit = 'Energy'
else:
raise RuntimeError('Failed to parse function string ' + func)
for n in range(0, nQ):
for n in range(0, num_spec):
for m in range(lstart, len_data):
char = data[m]
if char.startswith(' // ' + func):
......@@ -287,12 +268,12 @@ class LoadNMoldyn3Ascii(PythonAlgorithm):
raise RuntimeError('Failed to parse function string ' + func)
Qaxis = ''
for n in range(0, nQ):
for n in range(0, num_spec):
logger.information(str(start))
logger.information('Reading : ' + data[start[n]])
Slist = _make_list(data, start[n] + 1, start[n + 1] - 1)
if n == nQ - 1:
if n == num_spec - 1:
Slist[nP - 1] = Slist[nP - 1][:-1]
S = []
for m in range(0, nP):
......@@ -303,111 +284,158 @@ class LoadNMoldyn3Ascii(PythonAlgorithm):
logger.information('S values = ' + str(S[:2]) + ' to ' + str(S[-2:]))
if n == 0:
Qaxis += str(Q[n])
xDat = xEn
yDat = np.array(S)
eDat = eZero
data_x = xEn
data_y = np.array(S)
data_e = zero_error
else:
Qaxis += ',' + str(Q[n])
xDat = np.append(xDat, xEn)
yDat = np.append(yDat, np.array(S))
eDat = np.append(eDat, eZero)
outWS = name + '_' + func
CreateWorkspace(OutputWorkspace=outWS,
DataX=xDat,
DataY=yDat,
DataE=eDat,
Nspec=nQ,
UnitX=xUnit,
data_x = np.append(data_x, xEn)
data_y = np.append(data_y, np.array(S))
data_e = np.append(data_e, zero_error)
function_ws_name = base_name + '_' + func
CreateWorkspace(OutputWorkspace=function_ws_name,
DataX=data_x,
DataY=data_y,
DataE=data_e,
Nspec=num_spec,
UnitX=x_unit,
VerticalAxisUnit='MomentumTransfer',
VerticalAxisValues=Qaxis)
output_ws_list.append(outWS)
output_ws_list.append(function_ws_name)
GroupWorkspaces(InputWorkspaces=output_ws_list,
OutputWorkspace=self._out_ws)
out_ws = GroupWorkspaces(InputWorkspaces=output_ws_list,
OutputWorkspace=self._out_ws)
return out_ws
def _ascii_import(self, data, name):
"""
Import ASCII data.
#-------------------------------------------------------------------------------
@param data Raw ASCII data
@param name Name of data file
def _ascii_3d_import(self):
"""
Import 3D ASCII data (e.g. I(Q, t)).
"""
logger.notice('Loading ASCII data')
from IndirectCommon import getEfixed
from IndirectNeutron import ChangeAngles, InstrParas
logger.notice('Loading ASCII data: %s' % name)
val = _split_line(data[3])
Q = []
for n in range(1, len(val)):
Q.append(val[n])
nQ = len(Q)
x = []
y = []
for n in range(4, len(data)):
val = _split_line(data[n])
x.append(val[0])
yval = val[1:]
y.append(yval)
nX = len(x)
logger.information('nQ = ' + str(nQ))
logger.information('nT = ' + str(nX))
xT = np.array(x)
eZero = np.zeros(nX)
#Qaxis = ''
for m in range(0, nQ):
logger.information('Q[' + str(m + 1) + '] : ' + str(Q[m]))
S = []
for n in range(0, nX):
S.append(y[n][m])
if m == 0:
#Qaxis += str(Q[m])
xDat = xT
yDat = np.array(S)
eDat = eZero
else:
#Qaxis += ',' + str(Q[m])
xDat = np.append(xDat, xT)
yDat = np.append(yDat, np.array(S))
eDat = np.append(eDat, eZero)
CreateWorkspace(OutputWorkspace=self._out_ws,
DataX=xDat,
DataY=yDat,
DataE=eDat,
Nspec=nQ,
UnitX='TOF')
Qmax = Q[nQ - 1]
instr = 'MolDyn'
ana = 'simul'
if Qmax <= 2.0:
refl = '2'
else:
refl = '4'
InstrParas(self._out_ws, instr, ana, refl)
efixed = getEfixed(self._out_ws)
logger.information('Qmax = ' + str(Qmax) + ' ; efixed = ' + str(efixed))
pi4 = 4.0 * math.pi
# Read file
data = []
x_axis = ('time', 'ns')
v_axis = ('q', 'ang^-1')
with open(self._file_name, 'r') as handle:
for line in handle:
line = line.strip()
# Ignore empty lines
if line == '':
continue
# Data line (if not comment)
elif line.strip()[0] != '#':
line_values = np.array([ast.literal_eval(t.strip()) if 'nan' not in t.lower() else np.nan for t in line.split()])
data.append(line_values)
if x_axis is None or v_axis is None:
raise ValueError('Data is not in expected format for 3D data')
logger.debug('X axis: {0}'.format(x_axis))
logger.debug('V axis: {0}'.format(v_axis))
# Get axis and Y values
data = np.swapaxes(np.array(data), 0, 1)
x_axis_values = data[0,1:]
v_axis_values = data[1:,0]
y_values = np.ravel(data[1:,1:])
# Create the workspace
wks = CreateWorkspace(OutputWorkspace=self._out_ws,
DataX=x_axis_values,
DataY=y_values,
NSpec=v_axis_values.size,
UnitX=x_axis[1],
EnableLogging=False)
# Load the MolDyn instrument
q_max = v_axis_values[-1]
instrument = 'MolDyn'
reflection = '2' if q_max <= 2.0 else '4'
InstrParas(wks.name(), instrument, 'simul', reflection)
# Process angles
efixed = getEfixed(wks.name())
logger.information('Qmax={0}, Efixed={1}'.format(q_max, efixed))
wave = 1.8 * math.sqrt(25.2429 / efixed)
theta = []
for n in range(0, nQ):
qw = wave * Q[n] / pi4
ang = 2.0 * math.degrees(math.asin(qw))
theta.append(ang)
qw = wave * v_axis_values / (4.0 * math.pi)
theta = 2.0 * np.degrees(np.arcsin(qw))
ChangeAngles(wks.name(), instrument, theta)
return wks
ChangeAngles(self._out_ws, instr, theta)
#-------------------------------------------------------------------------------
def _ascii_2d_import(self):
"""
Import 2D ASCII data (e.g. DoS).
"""
logger.notice('Loading ASCII data')
# Regex
x_axis_regex = re.compile(r"\s*columns-1\s*=\s*([A-z0-9\-\s]+)\s*\(([A-z0-9-\s*]+)\)")
y_axis_regex = re.compile(r"\s*columns-2\s*=\s*([A-z\-\s]+)")
# Read file
data = []
x_axis = None
y_axis = None
with open(self._file_name, 'r') as handle:
for line in handle:
line = line.strip()
# Ignore empty lines
if line == "":
continue
x_axis_match = x_axis_regex.match(line)
y_axis_match = y_axis_regex.match(line)
# Line (X) header
if x_axis_match:
x_axis = (x_axis_match.group(1).strip(), x_axis_match.group(2))
# Data (Y) header
elif y_axis_match:
y_axis = y_axis_match.group(1)
# Data line (if not comment)
elif line.strip()[0] != '#':
line_values = np.array([ast.literal_eval(t.strip()) if 'nan' not in t.lower() else np.nan for t in line.split()])
data.append(line_values)
if x_axis is None or y_axis is None:
raise ValueError('Data is not in expected format for 2D data')
logger.debug('X axis: {0}'.format(x_axis))
logger.debug('Y axis: {0}'.format(y_axis))
# Get axis and Y values
data = np.array(data)
x_axis_values = data[:,0]
y_values = data[:,1]
# Create the workspace
wks = CreateWorkspace(OutputWorkspace=self._out_ws,
DataX=x_axis_values,
DataY=y_values,
NSpec=1,
UnitX=x_axis[1],
WorkspaceTitle=y_axis,
YUnitLabel=y_axis,
EnableLogging=False)
return wks
#==============================================================================
# Register algorithm with Mantid
AlgorithmFactory.subscribe(LoadNMoldyn3Ascii)
Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/LoadNMoldyn3AsciiTest.py
+ 15
− 4
View file @ 022facbc
......@@ -14,7 +14,6 @@ class LoadNMoldyn3AsciiTest(unittest.TestCase):
"""
Load an F(Q, t) function from an nMOLDYN 3 .cdl file
"""
moldyn_group = LoadNMoldyn3Ascii(Filename=self._cdl_filename,
Functions=['Fqt-total'],
OutputWorkspace='__LoadNMoldyn3Ascii_test')
......@@ -34,9 +33,8 @@ class LoadNMoldyn3AsciiTest(unittest.TestCase):
def test_load_sqw_from_cdl(self):
"""
Load an S(Q, w) function from a nMOLDYN 3 .cdl file
Load an S(Q, w) function from an nMOLDYN 3 .cdl file
"""
moldyn_group = LoadNMoldyn3Ascii(Filename=self._cdl_filename,
Functions=['Sqw-total'],
OutputWorkspace='__LoadNMoldyn3Ascii_test')
......@@ -53,11 +51,24 @@ class LoadNMoldyn3AsciiTest(unittest.TestCase):
self.assertEqual(units, 'Energy')
def test_load_dos_from_dat(self):
"""
Load a density of states from an nMOLDYN 3 .dat file.
"""
moldyn_ws = LoadNMoldyn3Ascii(Filename='CDOS_Croco_total_10K.dat',
OutputWorkspace='__LoadNMoldyn3Ascii_test')
self.assertTrue(isinstance(moldyn_ws, MatrixWorkspace))
self.assertTrue(moldyn_ws.getNumberHistograms(), 1)
self.assertEqual(moldyn_ws.getAxis(0).getUnit().label(), 'THz')
self.assertEqual(moldyn_ws.YUnitLabel(), 'dos-total vs frequency')
def test_load_from_dat(self):
"""
Load a function from an nMOLDYN 3 .dat file
"""
moldyn_ws = LoadNMoldyn3Ascii(Filename=self._dat_filename,
OutputWorkspace='__LoadNMoldyn3Ascii_test')
......
Code/Mantid/Testing/Data/UnitTest/CDOS_Croco_total_10K.dat.md5 0 → 100644
+ 1
− 0
View file @ 022facbc
6f70d4023fb7ddf43f30994cd0c4f6b6
Code/Mantid/Testing/SystemTests/tests/analysis/reference/ISISIndirectSimulation_MolDynDAT.nxs.md5
+ 1
− 1
View file @ 022facbc
bc22e1ca81aa608b028c11f7099b9f6c
\ No newline at end of file
48df154fcd3ae3119874479f929a5666
Code/Mantid/docs/source/algorithms/LoadNMoldyn3Ascii-v1.rst
+ 4
− 0
View file @ 022facbc
......@@ -13,6 +13,10 @@ This algorithm is a loader from simulation data saved by version 3 of the
nMOLDYN package, simulations can be loaded form either a plain ASCII (``.dat``)
file or ``.cdl`` file.
Currently this supports loading :math:`S(Q, \omega)` and :math:`F(Q, t)`
functions from ``.cdl`` files, :math:`F(Q, t)` from ``.dat`` files and any 2D
function (such as density of states) from ``.dat`` files.
Usage
-----
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment