sans_gen.py @ b11e127

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

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

SasView

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

Download in other formats: