sans_gen.py @ 9624cda

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

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

SasView

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

Download in other formats: