sans_gen.py @ b11e127

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload

Last change on this file since b11e127 was 11363ee, checked in by Jae Cho <jhjcho@…>, 12 years ago
fixes of connecting lines
Property mode set to `100644`
File size: 42.4 KB

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

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

SasView

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

Download in other formats: