Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
#pylint: disable=invalid-name,no-init
from mantid.kernel import *
from mantid.api import *
from mantid.simpleapi import *
import numpy as np
import __builtin__
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_projection_input(self, projection_table):
if isinstance(projection_table, ITableWorkspace):
column_names = set(projection_table.getColumnNames())
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):
logger.warning('You are running algorithm %s that is the beta stage of development' % (self.name()))
to_cut = self.getProperty("InputWorkspace").value
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_property = self.getProperty("P4Bin")
p4_bins = p4_bins_property.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 not p4_bins_property.isDefault:
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[0]
max = p4_bins[2]
step_size = p4_bins[1]
dim_range = max - min
if step_size > dim_range:
step_size = dim_range
nbins = int( dim_range / step_size)
extents.append(min)
extents.append(max)
bins.append(int(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.
n_padding = __builtin__.max(0, ndims-3)
for i in range(0, ndims):
if i <= 2:
label = projection_labels[i]
unit = target_units[i]
vec = list(projection[i]) + ( [0] * n_padding )
# These are always orthogonal to the rest.
else:
orthogonal_dimension = to_cut.getDimension(i)
label = orthogonal_dimension.getName()
unit = orthogonal_dimension.getUnits()
vec = [0] * ndims
vec[i] = 1
value = "%s, %s, %s" % ( label, unit, ",".join(map(str, vec)))
cut_alg.setPropertyValue("BasisVector{0}".format(i) , value)
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)