diff --git a/Code/Mantid/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/CutMD.py b/Code/Mantid/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/CutMD.py index 9585d3595bc5..feaeb730d7dc 100644 --- a/Code/Mantid/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/CutMD.py +++ b/Code/Mantid/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/CutMD.py @@ -1,6 +1,7 @@ from mantid.kernel import * from mantid.api import * from mantid.simpleapi import * +import numpy as np import os.path @@ -26,11 +27,110 @@ def PyInit(self): 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') + + def __to_mantid_slicing_binning(self, horace_binning, to_cut, dimension_index): + dim = to_cut.getDimension(dimension_index) + min = dim.getMinimum() + max = dim.getMaximum() + n_arguments = len(horace_binning) + if n_arguments == 0: + raise ValueError("binning parameter cannot be empty") + elif n_arguments == 1: + step_size = horace_binning[0] + n_bins = (max - min) / step_size + elif n_arguments == 2: + min = horace_binning[0] + max = horace_binning[1] + n_bins = 1 + elif n_arguments == 3: + step_size = horace_binning[1] + n_bins = (max - min) / step_size + min = horace_binning[0] + max = horace_binning[2] + pass + else: + raise ValueError("Too many arguments given to the binning parameter") + if min >= max: + raise ValueError("Dimension Min >= Max value. Min %.2f Max %.2f", min, max) + if n_bins < 1: + raise ValueError("Number of bins calculated to be < 1") + return (min, max, n_bins) + + + def __calculate_extents(self, v, u, w, ws): + M=np.array([u,v,w]) + Minv=np.linalg.inv(M) + + # We are assuming that the workspace has dimensions H, K, L in that order. The workspace MUST have 3 dimensions at least for the following to work. + Hrange=[ws.getDimension(0).getMinimum(), ws.getDimension(0).getMaximum()] + Krange=[ws.getDimension(1).getMinimum(), ws.getDimension(1).getMaximum()] + Lrange=[ws.getDimension(2).getMinimum(), ws.getDimension(2).getMaximum()] + + # 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]=original_corner.dot(Minv) + counter += 1 + + # Get the min max extents + extents = list() + print new_coords + 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])) + + # Copy extents for non crystallographic dimensions + non_crystallographic_dimensions = ws.getNumDims() - 3 + if non_crystallographic_dimensions > 0: + for i in range(0, non_crystallographic_dimensions): + extents.append(ws.getDimension(i + 3).getMinimum()) + extents.append(ws.getDimension(i + 3).getMaximum()) + + 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, :]) + u = np.array(projection_table.column('u')) + v = np.array(projection_table.column('v')) + w = np.cross(v,u) + return (u, v, w) + + def __make_labels(self, projection): + + class Mapping: + + def __init__(self, replace): + self.__replace = replace + + def replace(self, entry): + if entry > 0: + return 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 PyExec(self): to_cut = self.getProperty("InputWorkspace").value @@ -39,25 +139,67 @@ def PyExec(self): raise ValueError("Input Workspace must be in reciprocal lattice dimensions (HKL)") - projection = self.getProperty("Projection").value - if isinstance(projection, ITableWorkspace): - column_names = set(projection.getColumnNames()) + projection_table = self.getProperty("Projection").value + if isinstance(projection_table, ITableWorkspace): + column_names = set(projection_table.getColumnNames()) logger.warning(str(column_names)) - if not column_names == set(['u', 'v']): - if not column_names == set(['u', 'v', 'offsets']): - if not column_names == set(['u', 'v', 'w', 'offsets']): + if not column_names == set(['u', 'v', 'type']): + if not column_names == set(['u', 'v', 'offsets', 'type']): + if not column_names == set(['u', 'v', 'w', 'offsets', 'type']): raise ValueError("Projection table schema is wrong") - if projection.rowCount() != 3: + if projection_table.rowCount() != 3: raise ValueError("Projection table expects 3 rows") + + p1_bins = self.getProperty("P1Bin").value + p2_bins = self.getProperty("P2Bin").value + p3_bins = self.getProperty("P3Bin").value + p4_bins = self.getProperty("P4Bin").value # TODO handle 3D only slicing. + + # TODO. THESE ARE WRONG. Need to consider the acutal transformed extents as part of this. + xbins = self.__to_mantid_slicing_binning(p1_bins, to_cut, 0); + ybins = self.__to_mantid_slicing_binning(p1_bins, to_cut, 1); + zbins = self.__to_mantid_slicing_binning(p1_bins, to_cut, 2); + #ebins = self.__to_mantid_slicing_binning(p1_bins, to_cut, 3); # TODO. cannot guarantee this one is here + + # TODO. check that workspace is x=H, y=K, z=L before using the projections and slicing. + projection = self.__uvw_from_projection_table(projection_table) + u,v,w = projection + + extents = self.__calculate_extents(v, u, w, to_cut) + nopix = True # TODO from input parameter. + + cut_alg_name = "SliceMD" + if nopix: + cut_alg_name= "BinMD" + + projection_labels = self.__make_labels(projection) - clone_alg = AlgorithmManager.Instance().create("CloneMDWorkspace") - clone_alg.setChild(True) - clone_alg.initialize() - clone_alg.setProperty("InputWorkspace", to_cut) - clone_alg.setPropertyValue("OutputWorkspace", "cloned") - clone_alg.execute() - cloned = clone_alg.getProperty("OutputWorkspace").value - self.setProperty("OutputWorkspace", cloned) + cut_alg = AlgorithmManager.Instance().create(cut_alg_name) + cut_alg.setChild(True) + 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()): + logger.warning("WARNING....") + if i <= 2: + label = projection_labels[i] + unit = "TODO" # Haven't figured out how to do this yet. + vec = projection[i] + value = "%s, %s, %s" % ( label, unit, ",".join(map(str, vec))) + cut_alg.setPropertyValue("BasisVector{0}".format(i) , value) + logger.warning("BasisVector{0}".format(i)) + if i > 2: + raise RuntimeError("Not implmented yet for non-crystallographic basis vector generation.") + cut_alg.setProperty("OutputExtents", extents) + bins = [ int(xbins[2]), int(ybins[2]), int(zbins[2]) ] # Again this is a hack for 3 dimensional data only. + cut_alg.setProperty("OutputBins", bins) + + cut_alg.execute() + slice = cut_alg.getProperty("OutputWorkspace").value + self.setProperty("OutputWorkspace", slice) AlgorithmFactory.subscribe(CutMD) diff --git a/Code/Mantid/Framework/PythonInterface/test/python/mantid/api/CutMDTest.py b/Code/Mantid/Framework/PythonInterface/test/python/mantid/api/CutMDTest.py index e4c167c08b93..97c6e62d43d2 100644 --- a/Code/Mantid/Framework/PythonInterface/test/python/mantid/api/CutMDTest.py +++ b/Code/Mantid/Framework/PythonInterface/test/python/mantid/api/CutMDTest.py @@ -2,6 +2,7 @@ import testhelpers from mantid.simpleapi import * +from mantid.api import IMDHistoWorkspace class CutMDTest(unittest.TestCase): @@ -23,18 +24,27 @@ 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") + self.assertRaises(RuntimeError, CutMD, InputWorkspace=test_md, OutputWorkspace="out_ws", P1Bin=[0.1], P2Bin=[0.1], P3Bin=[0.1]) def test_slice_to_original(self): - out_md = CutMD(self.__in_md) - equals, result = CompareMDWorkspaces(Workspace1=self.__in_md, Workspace2=out_md) - self.assertTrue(equals, "Input and output workspaces should be identical" ) + out_md = CutMD(self.__in_md, P1Bin=[0.1], P2Bin=[0.1], P3Bin=[0.1]) + self.assertTrue(isinstance(out_md, IMDHistoWorkspace)) + # 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, 0.0]", out_md.getDimension(0).getName() ) + self.assertEquals("[0.0, 'eta', 0.0]", out_md.getDimension(1).getName() ) + self.assertEquals("[0.0, 0.0, 'xi']", out_md.getDimension(2).getName() ) 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") + self.assertRaises(RuntimeError, CutMD, InputWorkspace=self.__in_md, Projection=projection, OutputWorkspace="out_ws", P1Bin=[0.1], P2Bin=[0.1], P3Bin=[0.1]) def test_wrong_table_workspace_format_wrong_row_numbers(self): projection = CreateEmptyTableWorkspace() @@ -43,8 +53,9 @@ def test_wrong_table_workspace_format_wrong_row_numbers(self): 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") + self.assertRaises(RuntimeError, CutMD, InputWorkspace=self.__in_md, Projection=projection, OutputWorkspace="out_ws", P1Bin=[0.1], P2Bin=[0.1], P3Bin=[0.1]) def test_run_it(self): pass