Skip to content
Snippets Groups Projects
Commit b66348fe authored by Owen Arnold's avatar Owen Arnold
Browse files

refs #10530. Remove the algorithm parts

I want to get the key parts of this in. The python algorithm can relatively safely be added later.
parent 7dc4cb15
No related branches found
No related tags found
No related merge requests found
from mantid.kernel import *
from mantid.api import *
from mantid.simpleapi import *
import numpy as np
import os.path
import re
class Projection(object):
u = "u"
v = "v"
w = "w"
class ProjectionUnit(object):
r = "r"
a = "a"
class CutMD(DataProcessorAlgorithm):
def category(self):
return 'MDAlgorithms'
def summary(self):
return 'Slices multidimensional workspaces using input projection information and binning limits.'
def PyInit(self):
self.declareProperty(IMDEventWorkspaceProperty('InputWorkspace', '', direction=Direction.Input),
doc='MDWorkspace to slice')
self.declareProperty(ITableWorkspaceProperty('Projection', '', direction=Direction.Input, optional = PropertyMode.Optional), doc='Projection')
self.declareProperty(FloatArrayProperty(name='P1Bin', values=[]),
doc='Projection 1 binning')
self.declareProperty(FloatArrayProperty(name='P2Bin', values=[]),
doc='Projection 2 binning')
self.declareProperty(FloatArrayProperty(name='P3Bin', values=[]),
doc='Projection 3 binning')
self.declareProperty(FloatArrayProperty(name='P4Bin', values=[]),
doc='Projection 4 binning')
self.declareProperty(IMDWorkspaceProperty('OutputWorkspace', '',
direction=Direction.Output),
doc='Output cut workspace')
self.declareProperty(name="NoPix", defaultValue=False, doc="If False creates a full MDEventWorkspaces as output. True to create an MDHistoWorkspace as output. This is DND only in Horace terminology.")
self.declareProperty(name="CheckAxes", defaultValue=True, doc="Check that the axis look to be correct, and abort if not.")
def __calculate_steps(self, extents, horace_binning ):
# Because the step calculations may involve moving the extents, we re-publish the extents.
out_extents = extents
out_n_bins = list()
for i in range(len(horace_binning)):
n_arguments = len(horace_binning[i])
max_extent_index = (i*2) + 1
min_extent_index = (i*2)
dim_range = extents[ max_extent_index ] - extents[ min_extent_index ]
if n_arguments == 0:
raise ValueError("binning parameter cannot be empty")
elif n_arguments == 1:
step_size = horace_binning[i][0]
if step_size > dim_range:
step_size = dim_range
n_bins = int( dim_range / step_size)
# Correct the max extent to fit n * step_size
out_extents[max_extent_index] = extents[min_extent_index] + ( n_bins * step_size )
elif n_arguments == 2:
out_extents[ min_extent_index ] = horace_binning[i][0]
out_extents[ max_extent_index ] = horace_binning[i][1]
n_bins = 1
elif n_arguments == 3:
dim_min = horace_binning[i][0]
dim_max = horace_binning[i][2]
step_size = horace_binning[i][1]
dim_range = dim_max - dim_min
if step_size > dim_range:
step_size = dim_range
n_bins = int( dim_range / step_size)
# Correct the max extent to fit n * step_size
out_extents[ max_extent_index ] = dim_min + ( n_bins * step_size )
out_extents[ min_extent_index ] = dim_min
if n_bins < 1:
raise ValueError("Number of bins calculated to be < 1")
out_n_bins.append( n_bins )
return out_extents, out_n_bins
def __extents_in_current_projection(self, to_cut, dimension_index):
dim = to_cut.getDimension(dimension_index)
dim_min = dim.getMinimum()
dim_max = dim.getMaximum()
return (dim_min, dim_max)
def __calculate_extents(self, v, u, w, limits):
M=np.array([u,v,w])
Minv=np.linalg.inv(M)
# unpack limits
Hrange, Krange, Lrange = limits
# Create a numpy 2D array. Makes finding minimums and maximums for each transformed coordinates over every corner easier.
new_coords = np.empty([8, 3])
counter = 0
for h in Hrange:
for k in Krange:
for l in Lrange:
original_corner=np.array([h,k,l])
new_coords[counter]=np.dot(original_corner, Minv)
counter += 1
# Get the min max extents
extents = list()
for i in range(0,3):
# Vertical slice down each corner for each dimension, then determine the max, min and use as extents
extents.append(np.amin(new_coords[:,i]))
extents.append(np.amax(new_coords[:,i]))
return extents
def __uvw_from_projection_table(self, projection_table):
if not isinstance(projection_table, ITableWorkspace):
I = np.identity(3)
return (I[0, :], I[1, :], I[2, :])
column_names = projection_table.getColumnNames()
u = np.array(projection_table.column(Projection.u))
v = np.array(projection_table.column(Projection.v))
if not Projection.w in column_names:
w = np.cross(v,u)
else:
w = np.array(projection_table.column(Projection.w))
return (u, v, w)
def __units_from_projection_table(self, projection_table):
if not isinstance(projection_table, ITableWorkspace) or not "type" in projection_table.getColumnNames():
units = (ProjectionUnit.r, ProjectionUnit.r, ProjectionUnit.r)
else:
units = tuple(projection_table.column("type"))
return units
def __make_labels(self, projection):
class Mapping:
def __init__(self, replace):
self.__replace = replace
def replace(self, entry):
if np.absolute(entry) == 1:
if entry > 0:
return self.__replace
else:
return "-" + self.__replace
elif entry == 0:
return 0
else:
return "%.2f%s" % ( entry, self.__replace )
return entry
crystallographic_names = ['zeta', 'eta', 'xi' ]
labels = list()
for i in range(len(projection)):
cmapping = Mapping(crystallographic_names[i])
labels.append( [cmapping.replace(x) for x in projection[i] ] )
return labels
def __verify_input_workspace(self, to_cut):
coord_system = to_cut.getSpecialCoordinateSystem()
if not coord_system == SpecialCoordinateSystem.HKL:
raise ValueError("Input Workspace must be in reciprocal lattice dimensions (HKL)")
ndims = to_cut.getNumDims()
if ndims < 3 or ndims > 4:
raise ValueError("Input Workspace should be 3 or 4 dimensional")
# Try to sanity check the order of the dimensions. This is important.
axes_check = self.getProperty("CheckAxes").value
if axes_check:
predicates = ["^(H.*)|(\\[H,0,0\\].*)$","^(K.*)|(\\[0,K,0\\].*)$","^(L.*)|(\\[0,0,L\\].*)$"]
for i in range(ndims):
dimension = to_cut.getDimension(i)
p = re.compile(predicates[i])
if not p.match( dimension.getName() ):
raise ValueError("Dimensions must be in order H, K, L")
def __verify_projection_input(self, projection_table):
if isinstance(projection_table, ITableWorkspace):
column_names = set(projection_table.getColumnNames())
logger.warning(str(column_names))
if not column_names == set([Projection.u, Projection.v, 'type']):
if not column_names == set([Projection.u, Projection.v, 'offsets', 'type']):
if not column_names == set([Projection.u, Projection.v, Projection.w, 'offsets', 'type']):
raise ValueError("Projection table schema is wrong! Column names received: " + str(column_names) )
if projection_table.rowCount() != 3:
raise ValueError("Projection table expects 3 rows")
def __scale_projection(self, (u, v, w), origin_units, target_units, to_cut):
if set(origin_units) == set(target_units):
return (u,v,w) # Nothing to do.
ol = to_cut.getExperimentInfo(0).sample().getOrientedLattice()
projection_scaled = [u, v, w]
to_from_pairs = zip(origin_units, target_units)
for i in range(len(to_from_pairs)) :
proj = projection_scaled[i]
d_star = 2 * np.pi * ol.dstar( float(proj[0]), float(proj[1]), float(proj[2]) )
from_unit, to_unit = to_from_pairs[i]
if from_unit == to_unit:
continue
elif from_unit == ProjectionUnit.a: # From inverse Angstroms to rlu
projection_scaled[i] *= d_star
else: # From rlu to inverse Anstroms
projection_scaled[i] /= d_star
return projection_scaled
def PyExec(self):
to_cut = self.getProperty("InputWorkspace").value
self.__verify_input_workspace(to_cut)
ndims = to_cut.getNumDims()
nopix = self.getProperty("NoPix").value
projection_table = self.getProperty("Projection").value
self.__verify_projection_input(projection_table)
p1_bins = self.getProperty("P1Bin").value
p2_bins = self.getProperty("P2Bin").value
p3_bins = self.getProperty("P3Bin").value
p4_bins = self.getProperty("P4Bin").value
x_extents = self.__extents_in_current_projection(to_cut, 0);
y_extents = self.__extents_in_current_projection(to_cut, 1);
z_extents = self.__extents_in_current_projection(to_cut, 2);
projection = self.__uvw_from_projection_table(projection_table)
target_units = self.__units_from_projection_table(projection_table)
origin_units = (ProjectionUnit.r, ProjectionUnit.r, ProjectionUnit.r) # TODO. This is a hack!
u,v,w = self.__scale_projection(projection, origin_units, target_units, to_cut)
extents = self.__calculate_extents(v, u, w, ( x_extents, y_extents, z_extents ) )
extents, bins = self.__calculate_steps( extents, ( p1_bins, p2_bins, p3_bins ) )
if p4_bins != None:
if (ndims == 4):
n_args = len(p4_bins)
min, max = self.__extents_in_current_projection(to_cut, 3);
d_range = max - min
if n_args == 1:
step_size = p4_bins[0]
nbins = d_range / step_size
elif n_args == 2:
min = p4_bins[0]
max = p4_bins[1]
nbins = 1
elif n_args == 3:
min = p4_bins[i][0]
max = p4_bins[i][2]
step_size = p4_bins[1]
dim_range = max - min
if step_size > dim_range:
step_size = dim_range
n_bins = int( dim_range / step_size)
extents.append(min)
extents.append(max)
bins.append(nbins)
e_units = to_cut.getDimension(3).getUnits()
temp = list(target_units)
temp.append(target_units)
target_units = tuple(temp)
else:
raise ValueError("Cannot specify P4Bins unless the workspace is of sufficient dimensions")
projection_labels = self.__make_labels(projection)
cut_alg_name = "BinMD" if nopix else "SliceMD"
'''
Actually perform the binning operation
'''
cut_alg = self.createChildAlgorithm(name=cut_alg_name, startProgress=0, endProgress=1.0)
cut_alg.initialize()
cut_alg.setProperty("InputWorkspace", to_cut)
cut_alg.setPropertyValue("OutputWorkspace", "sliced")
cut_alg.setProperty("NormalizeBasisVectors", False)
cut_alg.setProperty("AxisAligned", False)
# Now for the basis vectors.
for i in range(0, to_cut.getNumDims()):
if i <= 2:
label = projection_labels[i]
unit = target_units[i]
vec = projection[i]
value = "%s, %s, %s" % ( label, unit, ",".join(map(str, vec)))
cut_alg.setPropertyValue("BasisVector{0}".format(i) , value)
if i > 2:
raise RuntimeError("Not implemented yet for non-crystallographic basis vector generation.")
cut_alg.setProperty("OutputExtents", extents)
cut_alg.setProperty("OutputBins", bins)
cut_alg.execute()
slice = cut_alg.getProperty("OutputWorkspace").value
# Attach the w-matrix (projection matrix)
if slice.getNumExperimentInfo() > 0:
u, v, w = projection
w_matrix = np.array([u, v, w]).flatten().tolist()
info = slice.getExperimentInfo(0)
info.run().addProperty("W_MATRIX", w_matrix, True)
self.setProperty("OutputWorkspace", slice)
AlgorithmFactory.subscribe(CutMD)
......@@ -10,7 +10,6 @@ set ( TEST_PY_FILES
AnalysisDataServiceTest.py
AxisTest.py
CatalogManagerTest.py
CutMDTest.py
DataProcessorAlgorithmTest.py
DeprecatedAlgorithmCheckerTest.py
ExperimentInfoTest.py
......
import unittest
import testhelpers
import numpy as np
from mantid.simpleapi import *
from mantid.api import IMDHistoWorkspace, IMDEventWorkspace
class CutMDTest(unittest.TestCase):
def setUp(self):
# Create a workspace
data_ws = CreateMDWorkspace(Dimensions=3, Extents=[-10,10,-10,10,-10,10], Names="A,B,C", Units="U,U,U")
# Mark the workspace as being in HKL
SetSpecialCoordinates(InputWorkspace=data_ws, SpecialCoordinates='HKL')
# Set the UB
SetUB(Workspace=data_ws, a = 1, b = 1, c = 1, alpha =90, beta=90, gamma = 90)
# Add some data to the workspace
FakeMDEventData(InputWorkspace=data_ws, PeakParams=[10000,0,0,0,1])
self.__in_md = data_ws
def tearDown(self):
DeleteWorkspace(self.__in_md )
def test_exec_throws_if_not_a_hkl_workspace(self):
test_md = CreateMDWorkspace(Dimensions=3, Extents=[-10,10,-10,10,-10,10], Names="A,B,C", Units="U,U,U")
# Explicitly set the coordinate system to lab Q.
SetSpecialCoordinates(InputWorkspace=test_md, SpecialCoordinates='Q (lab frame)')
self.assertRaises(RuntimeError, CutMD, InputWorkspace=test_md, OutputWorkspace="out_ws", P1Bin=[0.1], P2Bin=[0.1], P3Bin=[0.1], CheckAxes=False)
def test_exec_throws_if_set_to_be_a_hkl_workspace_but_with_missaligned_dimension_names(self):
test_md = CreateMDWorkspace(Dimensions=3, Extents=[-10,10,-10,10,-10,10], Names="K,H,L", Units="U,U,U") # K,H,L are the dimension names
SetSpecialCoordinates(InputWorkspace=test_md, SpecialCoordinates='HKL')
self.assertRaises(RuntimeError, CutMD, InputWorkspace=test_md, OutputWorkspace="out_ws", P1Bin=[0.1], P2Bin=[0.1], P3Bin=[0.1], CheckAxes=True)
def test_exec_throws_if_giving_4th_binning_parameter_when_workspace_is_3D(self):
test_md = CreateMDWorkspace(Dimensions=3, Extents=[-10,10,-10,10,-10,10], Names="H,K,L", Units="U,U,U")
# Explicitly set the coordinate system to lab Q.
SetSpecialCoordinates(InputWorkspace=test_md, SpecialCoordinates='HKL')
self.assertRaises(RuntimeError, CutMD, InputWorkspace=test_md, OutputWorkspace="out_ws", P1Bin=[0.1], P2Bin=[0.1], P3Bin=[0.1], P4Bin=[0.1])
def test_slice_to_original(self):
out_md = CutMD(self.__in_md, P1Bin=[0.1], P2Bin=[0.1], P3Bin=[0.1], CheckAxes=False)
self.assertTrue(isinstance(out_md, IMDEventWorkspace), "Should default to producing an IMDEventWorkspace.")
# No rotation. Basis vectors should have been left the same, so no extent changes.
self.assertEquals(self.__in_md.getDimension(0).getMinimum(), out_md.getDimension(0).getMinimum())
self.assertEquals(self.__in_md.getDimension(0).getMaximum(), out_md.getDimension(0).getMaximum())
self.assertEquals(self.__in_md.getDimension(1).getMinimum(), out_md.getDimension(1).getMinimum())
self.assertEquals(self.__in_md.getDimension(1).getMaximum(), out_md.getDimension(1).getMaximum())
self.assertEquals(self.__in_md.getDimension(2).getMinimum(), out_md.getDimension(2).getMinimum())
self.assertEquals(self.__in_md.getDimension(2).getMaximum(), out_md.getDimension(2).getMaximum())
self.assertEquals("['zeta', 0, 0]", out_md.getDimension(0).getName() )
self.assertEquals("[0, 'eta', 0]", out_md.getDimension(1).getName() )
self.assertEquals("[0, 0, 'xi']", out_md.getDimension(2).getName() )
self.assertTrue(isinstance(out_md, IMDEventWorkspace), "nopix defaults to True. Should get an IMDEventWorkspace")
def test_recalculate_extents_with_3_bin_arguments(self):
out_md = CutMD(self.__in_md, P1Bin=[0, 0.3, 0.8], P2Bin=[0.1], P3Bin=[0.1], CheckAxes=False, NoPix=True)
dim = out_md.getDimension(0)
self.assertAlmostEqual(0, dim.getMinimum(), 6, "Wrong minimum")
self.assertEqual(2, dim.getNBins(), "Wrong calculated number of bins")
self.assertAlmostEqual(0.6, dim.getMaximum(), 6, "Wrong calculated maximum")
def test_truncate_extents(self):
out_md = CutMD(self.__in_md, P1Bin=[0, 1.1, 1], P2Bin=[21], P3Bin=[0.1], CheckAxes=False, NoPix=True)
self.assertEqual(1, out_md.getDimension(0).getNBins(), "Step is beyond range. Should just be integrated")
self.assertEqual(1, out_md.getDimension(1).getNBins(), "Step is beyond range. Should just be integrated")
def test_wrong_projection_workspace_format_wrong_column_numbers(self):
projection = CreateEmptyTableWorkspace()
projection.addColumn("double", "u")
# missing other columns
self.assertRaises(RuntimeError, CutMD, InputWorkspace=self.__in_md, Projection=projection, OutputWorkspace="out_ws", P1Bin=[0.1], P2Bin=[0.1], P3Bin=[0.1], CheckAxes=False)
def test_wrong_table_workspace_format_wrong_row_numbers(self):
projection = CreateEmptyTableWorkspace()
# Correct number of columns, and names
projection.addColumn("double", "u")
projection.addColumn("double", "v")
projection.addColumn("double", "w")
projection.addColumn("double", "offset")
projection.addColumn("string", "type")
# Incorrect number of rows i.e. zero in this case as none added.
self.assertRaises(RuntimeError, CutMD, InputWorkspace=self.__in_md, Projection=projection, OutputWorkspace="out_ws", P1Bin=[0.1], P2Bin=[0.1], P3Bin=[0.1], CheckAxes=False)
def test_orthogonal_slice_with_scaling(self):
# We create a fake workspace and check to see that the extents get scaled with the new coordinate system when sliced
to_cut = CreateMDWorkspace(Dimensions=3, Extents=[-1,1,-1,1,-1,1], Names='H,K,L', Units='U,U,U')
# Set the UB
SetUB(Workspace=to_cut, a = 1, b = 1, c = 1, alpha =90, beta=90, gamma = 90)
SetSpecialCoordinates(InputWorkspace=to_cut, SpecialCoordinates='HKL')
scale_x = 2.0
scale_y = 2.0
projection = CreateEmptyTableWorkspace()
# Correct number of columns, and names
projection.addColumn("double", "u")
projection.addColumn("double", "v")
projection.addColumn("str", "type")
projection.addRow([scale_x,0,"r"])
projection.addRow([0,scale_y,"r"])
projection.addRow([0,0,"r"])
out_md = CutMD(to_cut, Projection=projection, P1Bin=[0.1], P2Bin=[0.1], P3Bin=[0.1])
scale_z = np.cross(projection.column(1), projection.column(0))[-1]
'''
Here we check that the corners in HKL end up in the expected positions when transformed into the new scaled basis
provided by the W transform (projection table)
'''
self.assertEquals(-(1/scale_x), out_md.getDimension(0).getMinimum())
self.assertEquals((1/scale_x), out_md.getDimension(0).getMaximum())
self.assertEquals(-(1/scale_y), out_md.getDimension(1).getMinimum())
self.assertEquals((1/scale_y), out_md.getDimension(1).getMaximum())
self.assertEquals((1/scale_z), out_md.getDimension(2).getMinimum())
self.assertEquals(-(1/scale_z), out_md.getDimension(2).getMaximum())
self.assertEquals("['2.00zeta', 0, 0]", out_md.getDimension(0).getName() )
self.assertEquals("[0, '2.00eta', 0]", out_md.getDimension(1).getName() )
self.assertEquals("[0, 0, '-4.00xi']", out_md.getDimension(2).getName() )
def test_non_orthogonal_slice(self):
# We create a fake workspace and check to see that the extents get transformed to the new coordinate system.
to_cut = CreateMDWorkspace(Dimensions=3, Extents=[-1,1,-1,1,-1,1], Names='H,K,L', Units='U,U,U')
# Set the UB
SetUB(Workspace=to_cut, a = 1, b = 1, c = 1, alpha =90, beta=90, gamma = 90)
SetSpecialCoordinates(InputWorkspace=to_cut, SpecialCoordinates='HKL')
projection = CreateEmptyTableWorkspace()
# Correct number of columns, and names
projection.addColumn("double", "u")
projection.addColumn("double", "v")
projection.addColumn("double", "w")
projection.addColumn("double", "offsets")
projection.addColumn("str", "type")
projection.addRow([1,-1, 0, 0, "r"])
projection.addRow([1, 1, 0, 0, "r"])
projection.addRow([0, 0, 1, 0, "r"])
out_md = CutMD(to_cut, Projection=projection, P1Bin=[0.1], P2Bin=[0.1], P3Bin=[0.1], NoPix=True)
'''
Here we check that the corners in HKL end up in the expected positions when transformed into the new scaled basis
provided by the W transform (projection table)
'''
self.assertEquals(-1, out_md.getDimension(0).getMinimum())
self.assertEquals(1, out_md.getDimension(0).getMaximum())
self.assertEquals(-1, out_md.getDimension(1).getMinimum())
self.assertEquals(1, out_md.getDimension(1).getMaximum())
self.assertEquals(-1, out_md.getDimension(2).getMinimum())
self.assertEquals(1, out_md.getDimension(2).getMaximum())
self.assertEquals("['zeta', 'zeta', 0]", out_md.getDimension(0).getName() )
self.assertEquals("['-eta', 'eta', 0]", out_md.getDimension(1).getName() )
self.assertEquals("[0, 0, 'xi']", out_md.getDimension(2).getName() )
self.assertTrue(isinstance(out_md, IMDHistoWorkspace), "Expect that the output was an IMDHistoWorkspace given the NoPix flag.")
run = out_md.getExperimentInfo(0).run()
w_matrix = run.getLogData("W_MATRIX").value
w_matrix = list(w_matrix)
self.assertEquals([1,1,0,-1,1,0,0,0,1], w_matrix, "W-matrix should have been set, but should be an identity matrix")
def test_orthogonal_slice_with_cropping(self):
# We create a fake workspace and check to see that using bin inputs for cropping works
to_cut = CreateMDWorkspace(Dimensions=3, Extents=[-1,1,-1,1,-1,1], Names='H,K,L', Units='U,U,U')
# Set the UB
SetUB(Workspace=to_cut, a = 1, b = 1, c = 1, alpha =90, beta=90, gamma = 90)
SetSpecialCoordinates(InputWorkspace=to_cut, SpecialCoordinates='HKL')
projection = CreateEmptyTableWorkspace()
# Correct number of columns, and names
projection.addColumn("double", "u")
projection.addColumn("double", "v")
projection.addColumn("double", "w")
projection.addColumn("double", "offsets")
projection.addColumn("str", "type")
projection.addRow([1, 0, 0, 0, "r"])
projection.addRow([0, 1, 0, 0, "r"])
projection.addRow([0, 0, 1, 0, "r"])
'''
Specify the cropping boundaries as part of the bin inputs.
'''
out_md = CutMD(to_cut, Projection=projection, P1Bin=[-0.5,0.5], P2Bin=[-0.1,0.1], P3Bin=[-0.3,0.3], NoPix=True)
'''
Here we check that the corners in HKL end up in the expected positions when transformed into the new scaled basis
provided by the W transform (projection table)
'''
self.assertAlmostEqual(-0.5, out_md.getDimension(0).getMinimum(), 6)
self.assertAlmostEqual(0.5, out_md.getDimension(0).getMaximum(), 6)
self.assertAlmostEqual(-0.1, out_md.getDimension(1).getMinimum(), 6)
self.assertAlmostEqual(0.1, out_md.getDimension(1).getMaximum(), 6)
self.assertAlmostEqual(-0.3, out_md.getDimension(2).getMinimum(), 6)
self.assertAlmostEqual(0.3, out_md.getDimension(2).getMaximum(), 6)
self.assertEquals("['zeta', 0, 0]", out_md.getDimension(0).getName() )
self.assertEquals("[0, 'eta', 0]", out_md.getDimension(1).getName() )
self.assertEquals("[0, 0, 'xi']", out_md.getDimension(2).getName() )
self.assertTrue(isinstance(out_md, IMDHistoWorkspace), "Expect that the output was an IMDHistoWorkspace given the NoPix flag.")
if __name__ == '__main__':
unittest.main()
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment