diff --git a/Code/Mantid/Framework/PythonInterface/plugins/algorithms/LoadNMoldyn3Ascii.py b/Code/Mantid/Framework/PythonInterface/plugins/algorithms/LoadNMoldyn3Ascii.py
index 951eb39e9911ddecb55cded5ec548f190b79c40c..8124fb10de1fa0e48cbb1bc1f5daa5b3eb41b2ca 100644
--- a/Code/Mantid/Framework/PythonInterface/plugins/algorithms/LoadNMoldyn3Ascii.py
+++ b/Code/Mantid/Framework/PythonInterface/plugins/algorithms/LoadNMoldyn3Ascii.py
@@ -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)
diff --git a/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/LoadNMoldyn3AsciiTest.py b/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/LoadNMoldyn3AsciiTest.py
index ba26ed1c9f5acb451c06474f4ccb7caf224e440b..9e65c1569384ecd852eb8205dce7faf667f23588 100644
--- a/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/LoadNMoldyn3AsciiTest.py
+++ b/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/LoadNMoldyn3AsciiTest.py
@@ -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')
diff --git a/Code/Mantid/Testing/Data/UnitTest/CDOS_Croco_total_10K.dat.md5 b/Code/Mantid/Testing/Data/UnitTest/CDOS_Croco_total_10K.dat.md5
new file mode 100644
index 0000000000000000000000000000000000000000..d11a46b00c2a8df75d4c32288b60d06e7c0ebfb4
--- /dev/null
+++ b/Code/Mantid/Testing/Data/UnitTest/CDOS_Croco_total_10K.dat.md5
@@ -0,0 +1 @@
+6f70d4023fb7ddf43f30994cd0c4f6b6
diff --git a/Code/Mantid/Testing/SystemTests/tests/analysis/reference/ISISIndirectSimulation_MolDynDAT.nxs.md5 b/Code/Mantid/Testing/SystemTests/tests/analysis/reference/ISISIndirectSimulation_MolDynDAT.nxs.md5
index a21e268edc40adb5165688c907468f36c0130937..ddc570ba6d30987d200d9cce06962f55263501bc 100644
--- a/Code/Mantid/Testing/SystemTests/tests/analysis/reference/ISISIndirectSimulation_MolDynDAT.nxs.md5
+++ b/Code/Mantid/Testing/SystemTests/tests/analysis/reference/ISISIndirectSimulation_MolDynDAT.nxs.md5
@@ -1 +1 @@
-bc22e1ca81aa608b028c11f7099b9f6c
\ No newline at end of file
+48df154fcd3ae3119874479f929a5666
diff --git a/Code/Mantid/docs/source/algorithms/LoadNMoldyn3Ascii-v1.rst b/Code/Mantid/docs/source/algorithms/LoadNMoldyn3Ascii-v1.rst
index e8e00ea6fb3d8e8ad9fad75c664ae3b6bc8fb0b9..6e0aca39fff33dcc414e7c857e94d67803ef95a0 100644
--- a/Code/Mantid/docs/source/algorithms/LoadNMoldyn3Ascii-v1.rst
+++ b/Code/Mantid/docs/source/algorithms/LoadNMoldyn3Ascii-v1.rst
@@ -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
-----