sans_gen.py @ cd89c6f

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload

Last change on this file since cd89c6f was cd89c6f, checked in by Jae Cho <jhjcho@…>, 11 years ago
cleaned up a little
Property mode set to `100644`
File size: 39.9 KB

Line
1	"""
2	SANS generic computation and sld file readers
3	"""
4	from sans.models.BaseComponent import BaseComponent
5	import sans.models.sans_extension.sld2i as mod
6	from periodictable import formula
7	from periodictable import nsf
8	import numpy
9	import os
10	import copy
11
12	MFactor_AM = 2.853E-12
13	MFactor_mT = 2.3164E-9
14	METER2ANG = 1.0E+10
15	#Avogadro constant [1/mol]
16	NA = 6.02214129e+23
17
18	def mag2sld(mag, v_unit=None):
19	"""
20	Convert magnetization to magnatic SLD
21	sldm = Dm * mag where Dm = gamma * classical elec. radius/(2*Bohr magneton)
22	Dm ~ 2.853E-12 [A^(-2)] ==> Shouldn't be 2.90636E-12 [A^(-2)]???
23	"""
24	if v_unit == "A/m":
25	factor = MFactor_AM
26	elif v_unit == "mT":
27	factor = MFactor_mT
28	else:
29	raise ValueError, "Invalid valueunit"
30	sld_m = factor * mag
31	return sld_m
32
33	class GenSAS(BaseComponent):
34	"""
35	Generic SAS computation Model based on sld (n & m) arrays
36	"""
37	def __init__(self):
38	"""
39	Init
40
41	:Params sld_data: MagSLD object
42	"""
43	# Initialize BaseComponent
44	BaseComponent.__init__(self)
45	self.sld_data = None
46	self.data_pos_unit = None
47	self.data_x = None
48	self.data_y = None
49	self.data_z = None
50	self.data_sldn = None
51	self.data_mx = None
52	self.data_my = None
53	self.data_mz = None
54	self.data_vol = None #[A^3]
55	self.is_avg = False
56	## Name of the model
57	self.name = "GenSAS"
58	## Define parameters
59	self.params = {}
60	self.params['scale'] = 1.0
61	self.params['background'] = 0.0
62	self.params['solvent_SLD'] = 0.0
63	self.params['Up_frac_in'] = 1.0
64	self.params['Up_frac_out'] = 1.0
65	self.params['Up_theta'] = 0.0
66	self.description = 'GenSAS'
67	## Parameter details [units, min, max]
68	self.details = {}
69	self.details['scale'] = ['', None, None]
70	self.details['background'] = ['[1/cm]', None, None]
71	self.details['solvent_SLD'] = ['1/A^(2)', None, None]
72	self.details['Up_frac_in'] = ['[u/(u+d)]', None, None]
73	self.details['Up_frac_out'] = ['[u/(u+d)]', None, None]
74	self.details['Up_theta'] = ['[deg]', None, None]
75	# fixed parameters
76	self.fixed = []
77
78	def set_pixel_volumes(self, volume):
79	"""
80	Set the volume of a pixel in (A^3) unit
81	:Param volume: pixel volume [float]
82	"""
83	if self.data_vol == None:
84	raise
85	self.data_vol = volume
86
87	def set_is_avg(self, is_avg=False):
88	"""
89	Sets is_avg: [bool]
90	"""
91	self.is_avg = is_avg
92
93	def _gen(self, x, y, i):
94	"""
95	Evaluate the function
96	:Param x: array of x-values
97	:Param y: array of y-values
98	:Param i: array of initial i-value
99	:return: function value
100	"""
101	len_x = len(self.data_x)
102	if self.is_avg:
103	len_x *= -1
104	len_q = len(x)
105	sldn = copy.deepcopy(self.data_sldn)
106	sldn -= self.params['solvent_SLD']
107	model = mod.new_GenI(len_x, self.data_x, self.data_y, self.data_z,
108	sldn, self.data_mx, self.data_my,
109	self.data_mz, self.data_vol,
110	self.params['Up_frac_in'],
111	self.params['Up_frac_out'],
112	self.params['Up_theta'])
113
114	mod.genicom(model, len_q, x, y, i)
115	return self.params['scale'] * i + self.params['background']
116
117	def set_sld_data(self, sld_data=None):
118	"""
119	Sets sld_data
120	"""
121	self.sld_data = sld_data
122	self.data_pos_unit = sld_data.pos_unit
123	self.data_x = sld_data.pos_x
124	self.data_y = sld_data.pos_y
125	self.data_z = sld_data.pos_z
126	self.data_sldn = sld_data.sld_n
127	self.data_mx = sld_data.sld_mx
128	self.data_my = sld_data.sld_my
129	self.data_mz = sld_data.sld_mz
130	self.data_vol = sld_data.vol_pix
131
132	def getProfile(self):
133	"""
134	Get SLD profile
135
136	: return: sld_data
137	"""
138	return self.sld_data
139
140	def run(self, x = 0.0):
141	"""
142	Evaluate the model
143	:param x: simple value
144	:return: (I value)
145	"""
146	if x.__class__.__name__ == 'list':
147	i_out = numpy.zeros_like(x[0])
148	y_in = numpy.zeros_like(x[0])
149	# 1D I is found at y =0 in the 2D pattern
150	return self._gen(x[0], y_in, i_out )
151	else:
152	msg = "Q must be given as list of qx's and qy's"
153	raise ValueError, msg
154
155	def runXY(self, x = 0.0):
156	"""
157	Evaluate the model
158	:param x: simple value
159	:return: I value
160	:Use this runXY() for the computation
161	"""
162	if x.__class__.__name__ == 'list':
163	i_out = numpy.zeros_like(x[0])
164	import time
165	s = time.time()
166	out = self._gen(x[0], x[1], i_out)
167	print "i_out_time", time.time() - s
168	return out
169	else:
170	msg = "Q must be given as list of qx's and qy's"
171	raise ValueError, msg
172
173	def evalDistribution(self, qdist):
174	"""
175	Evaluate a distribution of q-values.
176	* For 1D, a numpy array is expected as input:
177	evalDistribution(q)
178	where q is a numpy array.
179	* For 2D, a list of numpy arrays are expected: [qx_prime,qy_prime],
180	where 1D arrays,
181	:param qdist: ndarray of scalar q-values or list [qx,qy]
182	where qx,qy are 1D ndarrays
183	"""
184	if qdist.__class__.__name__ == 'list':
185	out = self.runXY(qdist)
186	return out
187	else:
188	mesg = "evalDistribution is expecting an ndarray of "
189	mesg += "a list [qx,qy] where qx,qy are arrays."
190	raise RuntimeError, mesg
191
192	class OMF2SLD:
193	"""
194	Convert OMFData to MAgData
195	"""
196	def __init__(self):
197	"""
198	Init
199	"""
200	self.pos_x = None
201	self.pos_y = None
202	self.pos_z = None
203	self.mx = None
204	self.my = None
205	self.mz = None
206	self.sld_n = None
207	self.vol_pix = None
208	self.output = None
209	self.omfdata = None
210
211	def set_data(self, omfdata, shape='rectangular'):
212	"""
213	Set all data
214	"""
215	self.omfdata = omfdata
216	length = int(omfdata.xnodes * omfdata.ynodes * omfdata.znodes)
217	pos_x = numpy.arange(omfdata.xmin,
218	omfdata.xnodes*omfdata.xstepsize + omfdata.xmin,
219	omfdata.xstepsize)
220	pos_y = numpy.arange(omfdata.ymin,
221	omfdata.ynodes*omfdata.ystepsize + omfdata.ymin,
222	omfdata.ystepsize)
223	pos_z = numpy.arange(omfdata.zmin,
224	omfdata.znodes*omfdata.zstepsize + omfdata.zmin,
225	omfdata.zstepsize)
226	self.pos_x = numpy.tile(pos_x, int(omfdata.ynodes * omfdata.znodes))
227	self.pos_y = pos_y.repeat(int(omfdata.xnodes))
228	self.pos_y = numpy.tile(self.pos_y, int(omfdata.znodes))
229	self.pos_z = pos_z.repeat(int(omfdata.xnodes * omfdata.ynodes))
230	self.mx = omfdata.mx
231	self.my = omfdata.my
232	self.mz = omfdata.mz
233	self.sld_n = numpy.zeros(length)
234
235	if omfdata.mx == None:
236	self.mx = numpy.zeros(length)
237	if omfdata.my == None:
238	self.my = numpy.zeros(length)
239	if omfdata.mz == None:
240	self.mz = numpy.zeros(length)
241
242	self._check_data_length(length)
243	self.remove_null_points(False, False)
244	mask = numpy.ones(len(self.sld_n), dtype=bool)
245	if shape.lower() == 'ellipsoid':
246	try:
247	# Pixel (step) size included
248	x_c = max(self.pos_x) + min(self.pos_x)
249	y_c = max(self.pos_y) + min(self.pos_y)
250	z_c = max(self.pos_z) + min(self.pos_z)
251	x_d = max(self.pos_x) - min(self.pos_x)
252	y_d = max(self.pos_y) - min(self.pos_y)
253	z_d = max(self.pos_z) - min(self.pos_z)
254	x_r = (x_d + omfdata.xstepsize) / 2.0
255	y_r = (y_d + omfdata.ystepsize) / 2.0
256	z_r = (z_d + omfdata.zstepsize) / 2.0
257	x_dir2 = ((self.pos_x - x_c / 2.0) / x_r)
258	x_dir2 *= x_dir2
259	y_dir2 = ((self.pos_y - y_c / 2.0) / y_r)
260	y_dir2 *= y_dir2
261	z_dir2 = ((self.pos_z - z_c / 2.0) / z_r)
262	z_dir2 *= z_dir2
263	mask = (x_dir2 + y_dir2 + z_dir2) <= 1.0
264	except:
265	pass
266	self.output = MagSLD(self.pos_x[mask], self.pos_y[mask],
267	self.pos_z[mask], self.sld_n[mask],
268	self.mx[mask], self.my[mask], self.mz[mask])
269	self.output.set_pix_type('pixel')
270	self.output.set_pixel_symbols('pixel')
271
272	def get_omfdata(self):
273	"""
274	Return all data
275	"""
276	return self.omfdata
277
278	def get_output(self):
279	"""
280	Return output
281	"""
282	return self.output
283
284	def _check_data_length(self, length):
285	"""
286	Check if the data lengths are consistent
287	:Params length: data length
288	"""
289	msg = "Error: Inconsistent data length."
290	if len(self.pos_x) != length:
291	raise ValueError, msg
292	if len(self.pos_y) != length:
293	raise ValueError, msg
294	if len(self.pos_z) != length:
295	raise ValueError, msg
296	if len(self.mx) != length:
297	raise ValueError, msg
298	if len(self.my) != length:
299	raise ValueError, msg
300	if len(self.mz) != length:
301	raise ValueError, msg
302
303	def remove_null_points(self, remove=False, recenter=False):
304	"""
305	Removes any mx, my, and mz = 0 points
306	"""
307	if remove:
308	is_nonzero = (numpy.fabs(self.mx) + numpy.fabs(self.my) +
309	numpy.fabs(self.mz)).nonzero()
310	if len(is_nonzero[0]) > 0:
311	self.pos_x = self.pos_x[is_nonzero]
312	self.pos_y = self.pos_y[is_nonzero]
313	self.pos_z = self.pos_z[is_nonzero]
314	self.sld_n = self.sld_n[is_nonzero]
315	self.mx = self.mx[is_nonzero]
316	self.my = self.my[is_nonzero]
317	self.mz = self.mz[is_nonzero]
318	if recenter:
319	self.pos_x -= (min(self.pos_x) + max(self.pos_x)) / 2.0
320	self.pos_y -= (min(self.pos_y) + max(self.pos_y)) / 2.0
321	self.pos_z -= (min(self.pos_z) + max(self.pos_z)) / 2.0
322
323	def get_magsld(self):
324	"""
325	return MagSLD
326	"""
327	return self.output
328
329
330	class OMFReader:
331	"""
332	Class to load omf/ascii files (3 columns w/header).
333	"""
334	## File type
335	type_name = "OMF ASCII"
336
337	## Wildcards
338	type = ["OMF files (.OMF, .omf)\|*.omf"]
339	## List of allowed extensions
340	ext = ['.omf', '.OMF']
341
342	def read(self, path):
343	"""
344	Load data file
345	:param path: file path
346	:return: x, y, z, sld_n, sld_mx, sld_my, sld_mz
347	"""
348	#data_conv_m = None
349	desc = ""
350	mx = numpy.zeros(0)
351	my = numpy.zeros(0)
352	mz = numpy.zeros(0)
353	try:
354	data_conv_val = None
355	data_conv_mesh = None
356	input_f = open(path, 'rb')
357	buff = input_f.read()
358	lines = buff.split('\n')
359	input_f.close()
360	output = OMFData()
361	valueunit = None
362	for line in lines:
363	toks = line.split()
364	# Read data
365	try:
366	_mx = float(toks[0])
367	_my = float(toks[1])
368	_mz = float(toks[2])
369	_mx = mag2sld(_mx, valueunit)
370	_my = mag2sld(_my, valueunit)
371	_mz = mag2sld(_mz, valueunit)
372	mx = numpy.append(mx, _mx)
373	my = numpy.append(my, _my)
374	mz = numpy.append(mz, _mz)
375	except:
376	# Skip non-data lines
377	pass
378	#Reading Header; Segment count ignored
379	s_line = line.split(":", 1)
380	if s_line[0].lower().count("oommf") > 0:
381	oommf = s_line[1].lstrip()
382	if s_line[0].lower().count("title") > 0:
383	title = s_line[1].lstrip()
384	if s_line[0].lower().count("desc") > 0:
385	desc += s_line[1].lstrip()
386	desc += '\n'
387	if s_line[0].lower().count("meshtype") > 0:
388	meshtype = s_line[1].lstrip()
389	if s_line[0].lower().count("meshunit") > 0:
390	meshunit = s_line[1].lstrip()
391	if meshunit.count("m") < 1:
392	msg = "Error: \n"
393	msg += "We accept only m as meshunit"
394	raise ValueError, msg
395	if s_line[0].lower().count("xbase") > 0:
396	xbase = s_line[1].lstrip()
397	if s_line[0].lower().count("ybase") > 0:
398	ybase = s_line[1].lstrip()
399	if s_line[0].lower().count("zbase") > 0:
400	zbase = s_line[1].lstrip()
401	if s_line[0].lower().count("xstepsize") > 0:
402	xstepsize = s_line[1].lstrip()
403	if s_line[0].lower().count("ystepsize") > 0:
404	ystepsize = s_line[1].lstrip()
405	if s_line[0].lower().count("zstepsize") > 0:
406	zstepsize = s_line[1].lstrip()
407	if s_line[0].lower().count("xnodes") > 0:
408	xnodes = s_line[1].lstrip()
409	if s_line[0].lower().count("ynodes") > 0:
410	ynodes = s_line[1].lstrip()
411	if s_line[0].lower().count("znodes") > 0:
412	znodes = s_line[1].lstrip()
413	if s_line[0].lower().count("xmin") > 0:
414	xmin = s_line[1].lstrip()
415	if s_line[0].lower().count("ymin") > 0:
416	ymin = s_line[1].lstrip()
417	if s_line[0].lower().count("zmin") > 0:
418	zmin = s_line[1].lstrip()
419	if s_line[0].lower().count("xmax") > 0:
420	xmax = s_line[1].lstrip()
421	if s_line[0].lower().count("ymax") > 0:
422	ymax = s_line[1].lstrip()
423	if s_line[0].lower().count("zmax") > 0:
424	zmax = s_line[1].lstrip()
425	if s_line[0].lower().count("valueunit") > 0:
426	valueunit = s_line[1].lstrip().rstrip()
427	if s_line[0].lower().count("valuemultiplier") > 0:
428	valuemultiplier = s_line[1].lstrip()
429	if s_line[0].lower().count("valuerangeminmag") > 0:
430	valuerangeminmag = s_line[1].lstrip()
431	if s_line[0].lower().count("valuerangemaxmag") > 0:
432	valuerangemaxmag = s_line[1].lstrip()
433	if s_line[0].lower().count("end") > 0:
434	#output.set_sldms(mx, my, mz)
435	output.filename = os.path.basename(path)
436	output.oommf = oommf
437	output.title = title
438	output.desc = desc
439	output.meshtype = meshtype
440	output.xbase = float(xbase) * METER2ANG
441	output.ybase = float(ybase) * METER2ANG
442	output.zbase = float(zbase) * METER2ANG
443	output.xstepsize = float(xstepsize) * METER2ANG
444	output.ystepsize = float(ystepsize) * METER2ANG
445	output.zstepsize = float(zstepsize) * METER2ANG
446	output.xnodes = float(xnodes)
447	output.ynodes = float(ynodes)
448	output.znodes = float(znodes)
449	output.xmin = float(xmin) * METER2ANG
450	output.ymin = float(ymin) * METER2ANG
451	output.zmin = float(zmin) * METER2ANG
452	output.xmax = float(xmax) * METER2ANG
453	output.ymax = float(ymax) * METER2ANG
454	output.zmax = float(zmax) * METER2ANG
455	output.valuemultiplier = valuemultiplier
456	output.valuerangeminmag = mag2sld(float(valuerangeminmag), \
457	valueunit)
458	output.valuerangemaxmag = mag2sld(float(valuerangemaxmag), \
459	valueunit)
460	output.set_m(mx, my, mz)
461	return output
462	except:
463	msg = "%s is not supported: \n" % path
464	msg += "We accept only Text format OMF file."
465	raise RuntimeError, msg
466
467	class PDBReader:
468	"""
469	PDB reader class: limited for reading the lines starting with 'ATOM'
470	"""
471	type_name = "PDB"
472	## Wildcards
473	type = ["pdb files (.PDB, .pdb)\|*.pdb"]
474	## List of allowed extensions
475	ext = ['.pdb', '.PDB']
476
477	def read(self, path):
478	"""
479	Load data file
480
481	:param path: file path
482
483	:return: MagSLD
484
485	:raise RuntimeError: when the file can't be opened
486	"""
487	pos_x = numpy.zeros(0)
488	pos_y = numpy.zeros(0)
489	pos_z = numpy.zeros(0)
490	sld_n = numpy.zeros(0)
491	sld_mx = numpy.zeros(0)
492	sld_my = numpy.zeros(0)
493	sld_mz = numpy.zeros(0)
494	vol_pix = numpy.zeros(0)
495	pix_symbol = numpy.zeros(0)
496	try:
497	input_f = open(path, 'rb')
498	buff = input_f.read()
499	lines = buff.split('\n')
500	input_f.close()
501	for line in lines:
502	try:
503	#toks = line.split()
504	# check if line starts with "ATOM"
505	if line[0:6].strip().count('ATM') > 0 or \
506	line[0:6].strip() == 'ATOM':
507	# define fields of interest
508	atom_id = line[6:11].strip()
509	atom_name = line[12:16].strip()
510	try:
511	float(line[12])
512	atom_name = atom_name[1].upper()
513	except:
514	if len(atom_name) == 4:
515	atom_name = atom_name[0].upper()
516	elif line[12] != ' ':
517	atom_name = atom_name[0].upper() + \
518	atom_name[1].lower()
519	else:
520	atom_name = atom_name[0].upper()
521	#res_name = line[17:20].strip()
522	#chain_name = line[21:22].strip()
523	_pos_x = float(line[30:38].strip())
524	_pos_y = float(line[38:46].strip())
525	_pos_z = float(line[46:54].strip())
526
527	pos_x = numpy.append(pos_x, _pos_x)
528	pos_y = numpy.append(pos_y, _pos_y)
529	pos_z = numpy.append(pos_z, _pos_z)
530	try:
531	# sld in Ang unit (from fm)
532	#val = formula(atom_name).atoms.keys()[0].neutron.b_c
533	#val *= 1.0e-5
534	val = nsf.neutron_sld(atom_name)[0]
535	# sld in Ang^-2 unit
536	val *= 1.0e-6
537	sld_n = numpy.append(sld_n, val)
538	atom = formula(atom_name)
539	# cm to A units
540	vol = 1.0e+24 * atom.mass / atom.density / NA
541	vol_pix = numpy.append(vol_pix, vol)
542	except:
543	print "Error: set the sld of %s to zero"% atom_name
544	sld_n = numpy.append(sld_n, 0.0)
545	sld_mx = numpy.append(sld_mx, 0)
546	sld_my = numpy.append(sld_my, 0)
547	sld_mz = numpy.append(sld_mz, 0)
548	pix_symbol = numpy.append(pix_symbol, atom_name)
549	except:
550	pass
551	#re-centering
552	#pos_x -= (min(pos_x) + max(pos_x)) / 2.0
553	#pos_y -= (min(pos_y) + max(pos_y)) / 2.0
554	#pos_z -= (min(pos_z) + max(pos_z)) / 2.0
555
556	output = MagSLD(pos_x, pos_y, pos_z, sld_n, sld_mx, sld_my, sld_mz)
557	output.filename = os.path.basename(path)
558	output.set_pix_type('atom')
559	output.set_pixel_symbols(pix_symbol)
560	output.set_nodes()
561	output.set_pixel_volumes(vol_pix)
562	output.sld_unit = '1/A^(2)'
563	return output
564	except:
565	RuntimeError, "%s is not a sld file" % path
566
567	def write(self, path, data):
568	"""
569	Write
570	"""
571	#Not implemented
572	print "Not implemented... "
573
574	class SLDReader:
575	"""
576	Class to load ascii files (7 columns).
577	"""
578	## File type
579	type_name = "SLD ASCII"
580
581	## Wildcards
582	type = ["sld files (.SLD, .sld)\|*.sld",
583	"txt files (.TXT, .txt)\|*.txt"]
584	## List of allowed extensions
585	ext = ['.txt', '.TXT', '.sld', 'SLD', '.*']
586
587	def read(self, path):
588	"""
589	Load data file
590
591	:param path: file path
592
593	:return MagSLD: x, y, z, sld_n, sld_mx, sld_my, sld_mz
594
595	:raise RuntimeError: when the file can't be opened
596	:raise ValueError: when the length of the data vectors are inconsistent
597	"""
598	try:
599	pos_x = numpy.zeros(0)
600	pos_y = numpy.zeros(0)
601	pos_z = numpy.zeros(0)
602	sld_n = numpy.zeros(0)
603	sld_mx = numpy.zeros(0)
604	sld_my = numpy.zeros(0)
605	sld_mz = numpy.zeros(0)
606	try:
607	# Use numpy to speed up loading
608	input_f = numpy.loadtxt(path, dtype='float', skiprows=1,
609	ndmin=1, unpack=True)
610	pos_x = numpy.array(input_f[0])
611	pos_y = numpy.array(input_f[1])
612	pos_z = numpy.array(input_f[2])
613	sld_n = numpy.array(input_f[3])
614	sld_mx = numpy.array(input_f[4])
615	sld_my = numpy.array(input_f[5])
616	sld_mz = numpy.array(input_f[6])
617	ncols = len(input_f)
618	if ncols == 8:
619	vol_pix = numpy.array(input_f[7])
620	elif ncols == 7:
621	vol_pix = None
622	except:
623	# For older version of numpy
624	input_f = open(path, 'rb')
625	buff = input_f.read()
626	lines = buff.split('\n')
627	input_f.close()
628	for line in lines:
629	toks = line.split()
630	try:
631	_pos_x = float(toks[0])
632	_pos_y = float(toks[1])
633	_pos_z = float(toks[2])
634	_sld_n = float(toks[3])
635	_sld_mx = float(toks[4])
636	_sld_my = float(toks[5])
637	_sld_mz = float(toks[6])
638
639	pos_x = numpy.append(pos_x, _pos_x)
640	pos_y = numpy.append(pos_y, _pos_y)
641	pos_z = numpy.append(pos_z, _pos_z)
642	sld_n = numpy.append(sld_n, _sld_n)
643	sld_mx = numpy.append(sld_mx, _sld_mx)
644	sld_my = numpy.append(sld_my, _sld_my)
645	sld_mz = numpy.append(sld_mz, _sld_mz)
646	try:
647	_vol_pix = float(toks[7])
648	vol_pix = numpy.append(vol_pix, _vol_pix)
649	except:
650	vol_pix = None
651	except:
652	# Skip non-data lines
653	pass
654
655	# re-center
656	#pos_x -= (min(pos_x) + max(pos_x)) / 2.0
657	#pos_y -= (min(pos_y) + max(pos_y)) / 2.0
658	#pos_z -= (min(pos_z) + max(pos_z)) / 2.0
659
660	output = MagSLD(pos_x, pos_y, pos_z, sld_n,
661	sld_mx, sld_my, sld_mz)
662
663	output.filename = os.path.basename(path)
664	output.set_pix_type('pixel')
665	output.set_pixel_symbols('pixel')
666	if vol_pix != None:
667	output.set_pixel_volumes(vol_pix)
668	return output
669	except:
670	RuntimeError, "%s is not a sld file" % path
671
672	def write(self, path, data):
673	"""
674	Write sld file
675
676	:Param path: file path
677	:Param data: MagSLD data object
678	"""
679	if path == None:
680	raise ValueError, "Missing the file path."
681	if data == None:
682	raise ValueError, "Missing the data to save."
683
684	x_val = data.pos_x
685	y_val = data.pos_y
686	z_val = data.pos_z
687	vol_pix = data.vol_pix
688
689	length = len(x_val)
690
691	sld_n = data.sld_n
692	if sld_n == None:
693	sld_n = numpy.zerros(length)
694
695	sld_mx = data.sld_mx
696	if sld_mx == None:
697	sld_mx = numpy.zerros(length)
698	sld_my = numpy.zerros(length)
699	sld_mz = numpy.zerros(length)
700	else:
701	sld_my = data.sld_my
702	sld_mz = data.sld_mz
703
704	out = open(path, 'w')
705	# First Line: Column names
706	out.write("X Y Z SLDN SLDMx SLDMy SLDMz VOLUMEpix")
707	for ind in range(length):
708	out.write("\n%g %g %g %g %g %g %g %g" % (x_val[ind], y_val[ind],
709	z_val[ind], sld_n[ind],
710	sld_mx[ind], sld_my[ind],
711	sld_mz[ind], vol_pix[ind]))
712	out.close()
713
714
715	class OMFData:
716	"""
717	OMF Data.
718	"""
719	_meshunit = "A"
720	_valueunit = "A^(-2)"
721	def __init__(self):
722	"""
723	Init for mag SLD
724	"""
725	self.filename = 'default'
726	self.oommf = ''
727	self.title = ''
728	self.desc = ''
729	self.meshtype = ''
730	self.meshunit = self._meshunit
731	self.valueunit = self._valueunit
732	self.xbase = 0.0
733	self.ybase = 0.0
734	self.zbase = 0.0
735	self.xstepsize = 6.0
736	self.ystepsize = 6.0
737	self.zstepsize = 6.0
738	self.xnodes = 10.0
739	self.ynodes = 10.0
740	self.znodes = 10.0
741	self.xmin = 0.0
742	self.ymin = 0.0
743	self.zmin = 0.0
744	self.xmax = 60.0
745	self.ymax = 60.0
746	self.zmax = 60.0
747	self.mx = None
748	self.my = None
749	self.mz = None
750	self.valuemultiplier = 1.
751	self.valuerangeminmag = 0
752	self.valuerangemaxmag = 0
753
754	def __str__(self):
755	"""
756	doc strings
757	"""
758	_str = "Type: %s\n" % self.__class__.__name__
759	_str += "File: %s\n" % self.filename
760	_str += "OOMMF: %s\n" % self.oommf
761	_str += "Title: %s\n" % self.title
762	_str += "Desc: %s\n" % self.desc
763	_str += "meshtype: %s\n" % self.meshtype
764	_str += "meshunit: %s\n" % str(self.meshunit)
765	_str += "xbase: %s [%s]\n" % (str(self.xbase), self.meshunit)
766	_str += "ybase: %s [%s]\n" % (str(self.ybase), self.meshunit)
767	_str += "zbase: %s [%s]\n" % (str(self.zbase), self.meshunit)
768	_str += "xstepsize: %s [%s]\n" % (str(self.xstepsize),
769	self.meshunit)
770	_str += "ystepsize: %s [%s]\n" % (str(self.ystepsize),
771	self.meshunit)
772	_str += "zstepsize: %s [%s]\n" % (str(self.zstepsize),
773	self.meshunit)
774	_str += "xnodes: %s\n" % str(self.xnodes)
775	_str += "ynodes: %s\n" % str(self.ynodes)
776	_str += "znodes: %s\n" % str(self.znodes)
777	_str += "xmin: %s [%s]\n" % (str(self.xmin), self.meshunit)
778	_str += "ymin: %s [%s]\n" % (str(self.ymin), self.meshunit)
779	_str += "zmin: %s [%s]\n" % (str(self.zmin), self.meshunit)
780	_str += "xmax: %s [%s]\n" % (str(self.xmax), self.meshunit)
781	_str += "ymax: %s [%s]\n" % (str(self.ymax), self.meshunit)
782	_str += "zmax: %s [%s]\n" % (str(self.zmax), self.meshunit)
783	_str += "valueunit: %s\n" % self.valueunit
784	_str += "valuemultiplier: %s\n" % str(self.valuemultiplier)
785	_str += "ValueRangeMinMag:%s [%s]\n" % (str(self.valuerangeminmag),
786	self.valueunit)
787	_str += "ValueRangeMaxMag:%s [%s]\n" % (str(self.valuerangemaxmag),
788	self.valueunit)
789	return _str
790
791	def set_m(self, mx, my, mz):
792	"""
793	Set the Mx, My, Mz values
794	"""
795	self.mx = mx
796	self.my = my
797	self.mz = mz
798
799	class MagSLD:
800	"""
801	Magnetic SLD.
802	"""
803	pos_x = None
804	pos_y = None
805	pos_z = None
806	sld_n = None
807	sld_mx = None
808	sld_my = None
809	sld_mz = None
810	# Units
811	_pos_unit = 'A'
812	_sld_unit = '1/A^(2)'
813	_pix_type = 'pixel'
814
815	def __init__(self, pos_x, pos_y, pos_z, sld_n=None,
816	sld_mx=None, sld_my=None, sld_mz=None, vol_pix=None):
817	"""
818	Init for mag SLD
819	:params : All should be numpy 1D array
820	"""
821	self.is_data = True
822	self.filename = ''
823	self.xstepsize = 6.0
824	self.ystepsize = 6.0
825	self.zstepsize = 6.0
826	self.xnodes = 10.0
827	self.ynodes = 10.0
828	self.znodes = 10.0
829	self.has_stepsize = False
830	self.pos_unit = self._pos_unit
831	self.sld_unit = self._sld_unit
832	self.pix_type = 'pixel'
833	self.pos_x = pos_x
834	self.pos_y = pos_y
835	self.pos_z = pos_z
836	self.sld_n = sld_n
837	self.sld_mx = sld_mx
838	self.sld_my = sld_my
839	self.sld_mz = sld_mz
840	self.vol_pix = vol_pix
841	self.sld_m = None
842	self.sld_phi = None
843	self.sld_theta = None
844	self.pix_symbol = None
845	if sld_mx != None and sld_my != None and sld_mz != None:
846	self.set_sldms(sld_mx, sld_my, sld_mz)
847	self.set_nodes()
848
849	def __str__(self):
850	"""
851	doc strings
852	"""
853	_str = "Type: %s\n" % self.__class__.__name__
854	_str += "File: %s\n" % self.filename
855	_str += "Axis_unit: %s\n" % self.pos_unit
856	_str += "SLD_unit: %s\n" % self.sld_unit
857	return _str
858
859	def set_pix_type(self, pix_type):
860	"""
861	Set pixel type
862	:Param pix_type: string, 'pixel' or 'atom'
863	"""
864	self.pix_type = pix_type
865
866	def set_sldn(self, sld_n):
867	"""
868	Sets neutron SLD
869	"""
870	if sld_n.__class__.__name__ == 'float':
871	if self.is_data:
872	# For data, put the value to only the pixels w non-zero M
873	is_nonzero = (numpy.fabs(self.sld_mx) +
874	numpy.fabs(self.sld_my) +
875	numpy.fabs(self.sld_mz)).nonzero()
876	self.sld_n = numpy.zeros(len(self.pos_x))
877	if len(self.sld_n[is_nonzero]) > 0:
878	self.sld_n[is_nonzero] = sld_n
879	else:
880	self.sld_n.fill(sld_n)
881	else:
882	# For non-data, put the value to all the pixels
883	self.sld_n = numpy.ones(len(self.pos_x)) * sld_n
884	else:
885	self.sld_n = sld_n
886
887	def set_sldms(self, sld_mx, sld_my, sld_mz):
888	"""
889	Sets (\|m\|, m_theta, m_phi)
890	"""
891	if sld_mx.__class__.__name__ == 'float':
892	self.sld_mx = numpy.ones(len(self.pos_x)) * sld_mx
893	else:
894	self.sld_mx = sld_mx
895	if sld_my.__class__.__name__ == 'float':
896	self.sld_my = numpy.ones(len(self.pos_x)) * sld_my
897	else:
898	self.sld_my = sld_my
899	if sld_mz.__class__.__name__ == 'float':
900	self.sld_mz = numpy.ones(len(self.pos_x)) * sld_mz
901	else:
902	self.sld_mz = sld_mz
903
904	sld_m = numpy.sqrt(sld_mx * sld_mx + sld_my * sld_my + \
905	sld_mz * sld_mz)
906	self.sld_m = sld_m
907
908	def set_pixel_symbols(self, symbol='pixel'):
909	"""
910	Set pixel
911	:Params pixel: str; pixel or atomic symbol, or array of strings
912	"""
913	if self.sld_n == None:
914	return
915	if symbol.__class__.__name__ == 'str':
916	self.pix_symbol = numpy.repeat(symbol, len(self.sld_n))
917	else:
918	self.pix_symbol = symbol
919
920	def set_pixel_volumes(self, vol):
921	"""
922	Set pixel volumes
923	:Params pixel: str; pixel or atomic symbol, or array of strings
924	"""
925	if self.sld_n == None:
926	return
927	if vol.__class__.__name__ == 'ndarray':
928	self.vol_pix = vol
929	elif vol.__class__.__name__.count('float') > 0:
930	self.vol_pix = numpy.repeat(vol, len(self.sld_n))
931	else:
932	self.vol_pix = None
933
934	def get_sldn(self):
935	"""
936	Returns nuclear sld
937	"""
938	return self.sld_n
939
940	def get_sld_mxyz(self):
941	"""
942	Returns (sld_mx, sldmy, sld_mz)
943	"""
944	return (self.sld_mx, self.sldmy, self.sld_mz)
945
946	def get_sld_m(self):
947	"""
948	Returns (sldm, sld_theta, sld_phi)
949	"""
950	return (self.sldm, self.sld_theta, self.sld_phi)
951
952	def set_nodes(self):
953	"""
954	Set xnodes, ynodes, and znodes
955	"""
956	self.set_stepsize()
957	if self.pix_type == 'pixel':
958	try:
959	xdist = (max(self.pos_x) - min(self.pos_x)) / self.xstepsize
960	ydist = (max(self.pos_y) - min(self.pos_y)) / self.ystepsize
961	zdist = (max(self.pos_z) - min(self.pos_z)) / self.zstepsize
962	self.xnodes = int(xdist) + 1
963	self.ynodes = int(ydist) + 1
964	self.znodes = int(zdist) + 1
965	except:
966	self.xnodes = None
967	self.ynodes = None
968	self.znodes = None
969	else:
970	self.xnodes = None
971	self.ynodes = None
972	self.znodes = None
973
974	def set_stepsize(self):
975	"""
976	Set xtepsize, ystepsize, and zstepsize
977	"""
978	if self.pix_type == 'pixel':
979	try:
980	xpos_pre = self.pos_x[0]
981	ypos_pre = self.pos_y[0]
982	zpos_pre = self.pos_z[0]
983	for x_pos in self.pos_x:
984	if xpos_pre != x_pos:
985	self.xstepsize = numpy.fabs(x_pos - xpos_pre)
986	break
987	for y_pos in self.pos_y:
988	if ypos_pre != y_pos:
989	self.ystepsize = numpy.fabs(y_pos - ypos_pre)
990	break
991	for z_pos in self.pos_z:
992	if zpos_pre != z_pos:
993	self.zstepsize = numpy.fabs(z_pos - zpos_pre)
994	break
995	#default pix volume
996	self.vol_pix = numpy.ones(len(self.pos_x))
997	vol = self.xstepsize * self.ystepsize * self.zstepsize
998	self.set_pixel_volumes(vol)
999	self.has_stepsize = True
1000	except:
1001	self.xstepsize = None
1002	self.ystepsize = None
1003	self.zstepsize = None
1004	self.vol_pix = None
1005	self.has_stepsize = False
1006	else:
1007	self.xstepsize = None
1008	self.ystepsize = None
1009	self.zstepsize = None
1010	self.has_stepsize = True
1011	return self.xstepsize, self.ystepsize, self.zstepsize
1012
1013
1014	def test_load():
1015	from danse.common.plottools.arrow3d import Arrow3D
1016	import os
1017	dir = os.path.abspath(os.path.curdir)
1018	print dir
1019	for i in range(6):
1020	dir, _ = os.path.split(dir)
1021	tfile = os.path.join(dir, "test", "CoreXY_ShellZ.txt")
1022	ofile = os.path.join(dir, "test", "A_Raw_Example-1.omf")
1023	if os.path.isfile(tfile):
1024	tfpath = tfile
1025	ofpath = ofile
1026	break
1027	reader = SLDReader()
1028	oreader = OMFReader()
1029	output = reader.read(tfpath)
1030	ooutput = oreader.read(ofpath)
1031	foutput = OMF2SLD()
1032	foutput.set_data(ooutput)
1033
1034	import matplotlib.pyplot as plt
1035	from mpl_toolkits.mplot3d import Axes3D
1036	fig = plt.figure()
1037	ax = fig.gca(projection='3d')
1038	ax.plot(output.pos_x, output.pos_y, output.pos_z, '.', c="g",
1039	alpha = 0.7, markeredgecolor='gray',rasterized=True)
1040	gap = 7
1041	max_mx = max(output.sld_mx)
1042	max_my = max(output.sld_my)
1043	max_mz = max(output.sld_mz)
1044	max_m = max(max_mx,max_my,max_mz)
1045	x2 = output.pos_x+output.sld_mx/max_m * gap
1046	y2 = output.pos_y+output.sld_my/max_m * gap
1047	z2 = output.pos_z+output.sld_mz/max_m * gap
1048	x_arrow = numpy.column_stack((output.pos_x,x2))
1049	y_arrow = numpy.column_stack((output.pos_y,y2))
1050	z_arrow = numpy.column_stack((output.pos_z,z2))
1051	unit_x2 = output.sld_mx / max_m
1052	unit_y2 = output.sld_my / max_m
1053	unit_z2 = output.sld_mz / max_m
1054	color_x = numpy.fabs(unit_x2 * 0.8)
1055	color_y = numpy.fabs(unit_y2 * 0.8)
1056	color_z = numpy.fabs(unit_z2 * 0.8)
1057	colors = numpy.column_stack((color_x, color_y, color_z))
1058	a = Arrow3D(None, x_arrow,y_arrow,z_arrow, colors, mutation_scale=10, lw=1,
1059	arrowstyle="->", color="y", alpha = 0.5)
1060
1061	ax.add_artist(a)
1062	plt.show()
1063
1064	def test():
1065	import os
1066	dir = os.path.abspath(os.path.curdir)
1067	for i in range(3):
1068	dir, _ = os.path.split(dir)
1069	#tfile = os.path.join(dir, "test", "C_Example_Converted.txt")
1070	ofile = os.path.join(dir, "test", "A_Raw_Example-1.omf")
1071	if os.path.isfile(ofile):
1072	#tfpath = tfile
1073	ofpath = ofile
1074	break
1075	#reader = SLDReader()
1076	oreader = OMFReader()
1077	#output = reader.read(tfpath)
1078	ooutput = oreader.read(ofpath)
1079	foutput = OMF2SLD()
1080	foutput.set_data(ooutput)
1081	writer = SLDReader()
1082	writer.write(os.path.join(os.path.dirname(ofpath), "out.txt"),
1083	foutput.output)
1084	model = GenSAS()
1085	model.set_sld_data(foutput.output)
1086	x = numpy.arange(1000)/10000. + 1e-5
1087	y = numpy.arange(1000)/10000. + 1e-5
1088	#z = numpy.arange(1000)/10000. + 1e-5
1089	i = numpy.zeros(1000)
1090	out = model.runXY([x,y,i])
1091
1092	if __name__ == "__main__":
1093	test()
1094	test_load()

Note: See TracBrowser for help on using the repository browser.

SasView

source: sasview/sanscalculator/src/sans/calculator/sans_gen.py @ cd89c6f

Download in other formats: