sas_gen.py @ 574adc7

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalcmagnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload

Last change on this file since 574adc7 was 574adc7, checked in by Paul Kienzle <pkienzle@…>, 7 years ago
convert sascalc to python 2/3 syntax
Property mode set to `100644`
File size: 39.2 KB

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

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

SasView

source: sasview/src/sas/sascalc/calculator/sas_gen.py @ 574adc7

Download in other formats: