sas_gen.py @ 7af652d

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

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

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

SasView

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

Download in other formats: