sas_gen.py @ a1daf86

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

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

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

SasView

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

Download in other formats: