sas_gen.py @ 6bd4235

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

Last change on this file since 6bd4235 was 144e032a, checked in by Paul Kienzle <pkienzle@…>, 7 years ago
tinycc doesn't return structures, so must pass return structure as pointer
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 .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	# ** 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 = mod.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	mod.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	mod.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:
307	logger.error(sys.exc_value)
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:
603	logger.error(sys.exc_value)
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:
694	vol_pix = None
695	except Exception:
696	# Skip non-data lines
697	logger.error(sys.exc_value)
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 @ 6bd4235

Download in other formats: