# pylint: disable=invalid-name
"""
SAS generic computation and sld file readers
"""
import sas.sascalc.calculator.core.sld2i as mod
from sas.sascalc.calculator.BaseComponent import BaseComponent
from periodictable import formula
from periodictable import nsf
import numpy
import os
import copy
import sys
import logging
MFACTOR_AM = 2.853E-12
MFACTOR_MT = 2.3164E-9
METER2ANG = 1.0E+10
#Avogadro constant [1/mol]
NA = 6.02214129e+23
[docs]def mag2sld(mag, v_unit=None):
"""
Convert magnetization to magnatic SLD
sldm = Dm * mag where Dm = gamma * classical elec. radius/(2*Bohr magneton)
Dm ~ 2.853E-12 [A^(-2)] ==> Shouldn't be 2.90636E-12 [A^(-2)]???
"""
if v_unit == "A/m":
factor = MFACTOR_AM
elif v_unit == "mT":
factor = MFACTOR_MT
else:
raise ValueError, "Invalid valueunit"
sld_m = factor * mag
return sld_m
[docs]class GenSAS(BaseComponent):
"""
Generic SAS computation Model based on sld (n & m) arrays
"""
def __init__(self):
"""
Init
:Params sld_data: MagSLD object
"""
# Initialize BaseComponent
BaseComponent.__init__(self)
self.sld_data = None
self.data_pos_unit = None
self.data_x = None
self.data_y = None
self.data_z = None
self.data_sldn = None
self.data_mx = None
self.data_my = None
self.data_mz = None
self.data_vol = None #[A^3]
self.is_avg = False
## Name of the model
self.name = "GenSAS"
## Define parameters
self.params = {}
self.params['scale'] = 1.0
self.params['background'] = 0.0
self.params['solvent_SLD'] = 0.0
self.params['total_volume'] = 1.0
self.params['Up_frac_in'] = 1.0
self.params['Up_frac_out'] = 1.0
self.params['Up_theta'] = 0.0
self.description = 'GenSAS'
## Parameter details [units, min, max]
self.details = {}
self.details['scale'] = ['', 0.0, numpy.inf]
self.details['background'] = ['[1/cm]', 0.0, numpy.inf]
self.details['solvent_SLD'] = ['1/A^(2)', -numpy.inf, numpy.inf]
self.details['total_volume'] = ['A^(3)', 0.0, numpy.inf]
self.details['Up_frac_in'] = ['[u/(u+d)]', 0.0, 1.0]
self.details['Up_frac_out'] = ['[u/(u+d)]', 0.0, 1.0]
self.details['Up_theta'] = ['[deg]', -numpy.inf, numpy.inf]
# fixed parameters
self.fixed = []
[docs] def set_pixel_volumes(self, volume):
"""
Set the volume of a pixel in (A^3) unit
:Param volume: pixel volume [float]
"""
if self.data_vol == None:
raise
self.data_vol = volume
[docs] def set_is_avg(self, is_avg=False):
"""
Sets is_avg: [bool]
"""
self.is_avg = is_avg
def _gen(self, x, y, i):
"""
Evaluate the function
:Param x: array of x-values
:Param y: array of y-values
:Param i: array of initial i-value
:return: function value
"""
pos_x = self.data_x
pos_y = self.data_y
pos_z = self.data_z
len_x = len(pos_x)
if self.is_avg == None:
len_x *= -1
pos_x, pos_y, pos_z = transform_center(pos_x, pos_y, pos_z)
len_q = len(x)
sldn = copy.deepcopy(self.data_sldn)
sldn -= self.params['solvent_SLD']
model = mod.new_GenI(len_x, pos_x, pos_y, pos_z,
sldn, self.data_mx, self.data_my,
self.data_mz, self.data_vol,
self.params['Up_frac_in'],
self.params['Up_frac_out'],
self.params['Up_theta'])
if y == []:
mod.genicom(model, len_q, x, i)
else:
mod.genicomXY(model, len_q, x, y, i)
vol_correction = self.data_total_volume / self.params['total_volume']
return self.params['scale'] * vol_correction * i + \
self.params['background']
[docs] def set_sld_data(self, sld_data=None):
"""
Sets sld_data
"""
self.sld_data = sld_data
self.data_pos_unit = sld_data.pos_unit
self.data_x = sld_data.pos_x
self.data_y = sld_data.pos_y
self.data_z = sld_data.pos_z
self.data_sldn = sld_data.sld_n
self.data_mx = sld_data.sld_mx
self.data_my = sld_data.sld_my
self.data_mz = sld_data.sld_mz
self.data_vol = sld_data.vol_pix
self.data_total_volume = sum(sld_data.vol_pix)
self.params['total_volume'] = sum(sld_data.vol_pix)
[docs] def getProfile(self):
"""
Get SLD profile
: return: sld_data
"""
return self.sld_data
[docs] def run(self, x=0.0):
"""
Evaluate the model
:param x: simple value
:return: (I value)
"""
if x.__class__.__name__ == 'list':
if len(x[1]) > 0:
msg = "Not a 1D."
raise ValueError, msg
i_out = numpy.zeros_like(x[0])
# 1D I is found at y =0 in the 2D pattern
out = self._gen(x[0], [], i_out)
return out
else:
msg = "Q must be given as list of qx's and qy's"
raise ValueError, msg
[docs] def runXY(self, x=0.0):
"""
Evaluate the model
:param x: simple value
:return: I value
:Use this runXY() for the computation
"""
if x.__class__.__name__ == 'list':
i_out = numpy.zeros_like(x[0])
out = self._gen(x[0], x[1], i_out)
return out
else:
msg = "Q must be given as list of qx's and qy's"
raise ValueError, msg
[docs] def evalDistribution(self, qdist):
"""
Evaluate a distribution of q-values.
:param qdist: ndarray of scalar q-values (for 1D) or list [qx,qy]
where qx,qy are 1D ndarrays (for 2D).
"""
if qdist.__class__.__name__ == 'list':
if len(qdist[1]) < 1:
out = self.run(qdist)
else:
out = self.runXY(qdist)
return out
else:
mesg = "evalDistribution is expecting an ndarray of "
mesg += "a list [qx,qy] where qx,qy are arrays."
raise RuntimeError, mesg
[docs]class OMF2SLD(object):
"""
Convert OMFData to MAgData
"""
def __init__(self):
"""
Init
"""
self.pos_x = None
self.pos_y = None
self.pos_z = None
self.mx = None
self.my = None
self.mz = None
self.sld_n = None
self.vol_pix = None
self.output = None
self.omfdata = None
[docs] def set_data(self, omfdata, shape='rectangular'):
"""
Set all data
"""
self.omfdata = omfdata
length = int(omfdata.xnodes * omfdata.ynodes * omfdata.znodes)
pos_x = numpy.arange(omfdata.xmin,
omfdata.xnodes*omfdata.xstepsize + omfdata.xmin,
omfdata.xstepsize)
pos_y = numpy.arange(omfdata.ymin,
omfdata.ynodes*omfdata.ystepsize + omfdata.ymin,
omfdata.ystepsize)
pos_z = numpy.arange(omfdata.zmin,
omfdata.znodes*omfdata.zstepsize + omfdata.zmin,
omfdata.zstepsize)
self.pos_x = numpy.tile(pos_x, int(omfdata.ynodes * omfdata.znodes))
self.pos_y = pos_y.repeat(int(omfdata.xnodes))
self.pos_y = numpy.tile(self.pos_y, int(omfdata.znodes))
self.pos_z = pos_z.repeat(int(omfdata.xnodes * omfdata.ynodes))
self.mx = omfdata.mx
self.my = omfdata.my
self.mz = omfdata.mz
self.sld_n = numpy.zeros(length)
if omfdata.mx == None:
self.mx = numpy.zeros(length)
if omfdata.my == None:
self.my = numpy.zeros(length)
if omfdata.mz == None:
self.mz = numpy.zeros(length)
self._check_data_length(length)
self.remove_null_points(False, False)
mask = numpy.ones(len(self.sld_n), dtype=bool)
if shape.lower() == 'ellipsoid':
try:
# Pixel (step) size included
x_c = max(self.pos_x) + min(self.pos_x)
y_c = max(self.pos_y) + min(self.pos_y)
z_c = max(self.pos_z) + min(self.pos_z)
x_d = max(self.pos_x) - min(self.pos_x)
y_d = max(self.pos_y) - min(self.pos_y)
z_d = max(self.pos_z) - min(self.pos_z)
x_r = (x_d + omfdata.xstepsize) / 2.0
y_r = (y_d + omfdata.ystepsize) / 2.0
z_r = (z_d + omfdata.zstepsize) / 2.0
x_dir2 = ((self.pos_x - x_c / 2.0) / x_r)
x_dir2 *= x_dir2
y_dir2 = ((self.pos_y - y_c / 2.0) / y_r)
y_dir2 *= y_dir2
z_dir2 = ((self.pos_z - z_c / 2.0) / z_r)
z_dir2 *= z_dir2
mask = (x_dir2 + y_dir2 + z_dir2) <= 1.0
except:
logging.error(sys.exc_value)
self.output = MagSLD(self.pos_x[mask], self.pos_y[mask],
self.pos_z[mask], self.sld_n[mask],
self.mx[mask], self.my[mask], self.mz[mask])
self.output.set_pix_type('pixel')
self.output.set_pixel_symbols('pixel')
[docs] def get_omfdata(self):
"""
Return all data
"""
return self.omfdata
[docs] def get_output(self):
"""
Return output
"""
return self.output
def _check_data_length(self, length):
"""
Check if the data lengths are consistent
:Params length: data length
"""
msg = "Error: Inconsistent data length."
if len(self.pos_x) != length:
raise ValueError, msg
if len(self.pos_y) != length:
raise ValueError, msg
if len(self.pos_z) != length:
raise ValueError, msg
if len(self.mx) != length:
raise ValueError, msg
if len(self.my) != length:
raise ValueError, msg
if len(self.mz) != length:
raise ValueError, msg
[docs] def remove_null_points(self, remove=False, recenter=False):
"""
Removes any mx, my, and mz = 0 points
"""
if remove:
is_nonzero = (numpy.fabs(self.mx) + numpy.fabs(self.my) +
numpy.fabs(self.mz)).nonzero()
if len(is_nonzero[0]) > 0:
self.pos_x = self.pos_x[is_nonzero]
self.pos_y = self.pos_y[is_nonzero]
self.pos_z = self.pos_z[is_nonzero]
self.sld_n = self.sld_n[is_nonzero]
self.mx = self.mx[is_nonzero]
self.my = self.my[is_nonzero]
self.mz = self.mz[is_nonzero]
if recenter:
self.pos_x -= (min(self.pos_x) + max(self.pos_x)) / 2.0
self.pos_y -= (min(self.pos_y) + max(self.pos_y)) / 2.0
self.pos_z -= (min(self.pos_z) + max(self.pos_z)) / 2.0
[docs] def get_magsld(self):
"""
return MagSLD
"""
return self.output
[docs]class OMFReader(object):
"""
Class to load omf/ascii files (3 columns w/header).
"""
## File type
type_name = "OMF ASCII"
## Wildcards
type = ["OMF files (*.OMF, *.omf)|*.omf"]
## List of allowed extensions
ext = ['.omf', '.OMF']
[docs] def read(self, path):
"""
Load data file
:param path: file path
:return: x, y, z, sld_n, sld_mx, sld_my, sld_mz
"""
desc = ""
mx = numpy.zeros(0)
my = numpy.zeros(0)
mz = numpy.zeros(0)
try:
input_f = open(path, 'rb')
buff = input_f.read()
lines = buff.split('\n')
input_f.close()
output = OMFData()
valueunit = None
for line in lines:
toks = line.split()
# Read data
try:
_mx = float(toks[0])
_my = float(toks[1])
_mz = float(toks[2])
_mx = mag2sld(_mx, valueunit)
_my = mag2sld(_my, valueunit)
_mz = mag2sld(_mz, valueunit)
mx = numpy.append(mx, _mx)
my = numpy.append(my, _my)
mz = numpy.append(mz, _mz)
except:
# Skip non-data lines
logging.error(sys.exc_value)
#Reading Header; Segment count ignored
s_line = line.split(":", 1)
if s_line[0].lower().count("oommf") > 0:
oommf = s_line[1].lstrip()
if s_line[0].lower().count("title") > 0:
title = s_line[1].lstrip()
if s_line[0].lower().count("desc") > 0:
desc += s_line[1].lstrip()
desc += '\n'
if s_line[0].lower().count("meshtype") > 0:
meshtype = s_line[1].lstrip()
if s_line[0].lower().count("meshunit") > 0:
meshunit = s_line[1].lstrip()
if meshunit.count("m") < 1:
msg = "Error: \n"
msg += "We accept only m as meshunit"
raise ValueError, msg
if s_line[0].lower().count("xbase") > 0:
xbase = s_line[1].lstrip()
if s_line[0].lower().count("ybase") > 0:
ybase = s_line[1].lstrip()
if s_line[0].lower().count("zbase") > 0:
zbase = s_line[1].lstrip()
if s_line[0].lower().count("xstepsize") > 0:
xstepsize = s_line[1].lstrip()
if s_line[0].lower().count("ystepsize") > 0:
ystepsize = s_line[1].lstrip()
if s_line[0].lower().count("zstepsize") > 0:
zstepsize = s_line[1].lstrip()
if s_line[0].lower().count("xnodes") > 0:
xnodes = s_line[1].lstrip()
if s_line[0].lower().count("ynodes") > 0:
ynodes = s_line[1].lstrip()
if s_line[0].lower().count("znodes") > 0:
znodes = s_line[1].lstrip()
if s_line[0].lower().count("xmin") > 0:
xmin = s_line[1].lstrip()
if s_line[0].lower().count("ymin") > 0:
ymin = s_line[1].lstrip()
if s_line[0].lower().count("zmin") > 0:
zmin = s_line[1].lstrip()
if s_line[0].lower().count("xmax") > 0:
xmax = s_line[1].lstrip()
if s_line[0].lower().count("ymax") > 0:
ymax = s_line[1].lstrip()
if s_line[0].lower().count("zmax") > 0:
zmax = s_line[1].lstrip()
if s_line[0].lower().count("valueunit") > 0:
valueunit = s_line[1].lstrip().rstrip()
if s_line[0].lower().count("valuemultiplier") > 0:
valuemultiplier = s_line[1].lstrip()
if s_line[0].lower().count("valuerangeminmag") > 0:
valuerangeminmag = s_line[1].lstrip()
if s_line[0].lower().count("valuerangemaxmag") > 0:
valuerangemaxmag = s_line[1].lstrip()
if s_line[0].lower().count("end") > 0:
output.filename = os.path.basename(path)
output.oommf = oommf
output.title = title
output.desc = desc
output.meshtype = meshtype
output.xbase = float(xbase) * METER2ANG
output.ybase = float(ybase) * METER2ANG
output.zbase = float(zbase) * METER2ANG
output.xstepsize = float(xstepsize) * METER2ANG
output.ystepsize = float(ystepsize) * METER2ANG
output.zstepsize = float(zstepsize) * METER2ANG
output.xnodes = float(xnodes)
output.ynodes = float(ynodes)
output.znodes = float(znodes)
output.xmin = float(xmin) * METER2ANG
output.ymin = float(ymin) * METER2ANG
output.zmin = float(zmin) * METER2ANG
output.xmax = float(xmax) * METER2ANG
output.ymax = float(ymax) * METER2ANG
output.zmax = float(zmax) * METER2ANG
output.valuemultiplier = valuemultiplier
output.valuerangeminmag = mag2sld(float(valuerangeminmag), \
valueunit)
output.valuerangemaxmag = mag2sld(float(valuerangemaxmag), \
valueunit)
output.set_m(mx, my, mz)
return output
except:
msg = "%s is not supported: \n" % path
msg += "We accept only Text format OMF file."
raise RuntimeError, msg
[docs]class PDBReader(object):
"""
PDB reader class: limited for reading the lines starting with 'ATOM'
"""
type_name = "PDB"
## Wildcards
type = ["pdb files (*.PDB, *.pdb)|*.pdb"]
## List of allowed extensions
ext = ['.pdb', '.PDB']
[docs] def read(self, path):
"""
Load data file
:param path: file path
:return: MagSLD
:raise RuntimeError: when the file can't be opened
"""
pos_x = numpy.zeros(0)
pos_y = numpy.zeros(0)
pos_z = numpy.zeros(0)
sld_n = numpy.zeros(0)
sld_mx = numpy.zeros(0)
sld_my = numpy.zeros(0)
sld_mz = numpy.zeros(0)
vol_pix = numpy.zeros(0)
pix_symbol = numpy.zeros(0)
x_line = []
y_line = []
z_line = []
x_lines = []
y_lines = []
z_lines = []
try:
input_f = open(path, 'rb')
buff = input_f.read()
lines = buff.split('\n')
input_f.close()
num = 0
for line in lines:
try:
# check if line starts with "ATOM"
if line[0:6].strip().count('ATM') > 0 or \
line[0:6].strip() == 'ATOM':
# define fields of interest
atom_name = line[12:16].strip()
try:
float(line[12])
atom_name = atom_name[1].upper()
except:
if len(atom_name) == 4:
atom_name = atom_name[0].upper()
elif line[12] != ' ':
atom_name = atom_name[0].upper() + \
atom_name[1].lower()
else:
atom_name = atom_name[0].upper()
_pos_x = float(line[30:38].strip())
_pos_y = float(line[38:46].strip())
_pos_z = float(line[46:54].strip())
pos_x = numpy.append(pos_x, _pos_x)
pos_y = numpy.append(pos_y, _pos_y)
pos_z = numpy.append(pos_z, _pos_z)
try:
val = nsf.neutron_sld(atom_name)[0]
# sld in Ang^-2 unit
val *= 1.0e-6
sld_n = numpy.append(sld_n, val)
atom = formula(atom_name)
# cm to A units
vol = 1.0e+24 * atom.mass / atom.density / NA
vol_pix = numpy.append(vol_pix, vol)
except:
print "Error: set the sld of %s to zero"% atom_name
sld_n = numpy.append(sld_n, 0.0)
sld_mx = numpy.append(sld_mx, 0)
sld_my = numpy.append(sld_my, 0)
sld_mz = numpy.append(sld_mz, 0)
pix_symbol = numpy.append(pix_symbol, atom_name)
elif line[0:6].strip().count('CONECT') > 0:
toks = line.split()
num = int(toks[1]) - 1
val_list = []
for val in toks[2:]:
try:
int_val = int(val)
except:
break
if int_val == 0:
break
val_list.append(int_val)
#need val_list ordered
for val in val_list:
index = val - 1
if (pos_x[index], pos_x[num]) in x_line and \
(pos_y[index], pos_y[num]) in y_line and \
(pos_z[index], pos_z[num]) in z_line:
continue
x_line.append((pos_x[num], pos_x[index]))
y_line.append((pos_y[num], pos_y[index]))
z_line.append((pos_z[num], pos_z[index]))
if len(x_line) > 0:
x_lines.append(x_line)
y_lines.append(y_line)
z_lines.append(z_line)
except:
logging.error(sys.exc_value)
output = MagSLD(pos_x, pos_y, pos_z, sld_n, sld_mx, sld_my, sld_mz)
output.set_conect_lines(x_line, y_line, z_line)
output.filename = os.path.basename(path)
output.set_pix_type('atom')
output.set_pixel_symbols(pix_symbol)
output.set_nodes()
output.set_pixel_volumes(vol_pix)
output.sld_unit = '1/A^(2)'
return output
except:
raise RuntimeError, "%s is not a sld file" % path
[docs] def write(self, path, data):
"""
Write
"""
print "Not implemented... "
[docs]class SLDReader(object):
"""
Class to load ascii files (7 columns).
"""
## File type
type_name = "SLD ASCII"
## Wildcards
type = ["sld files (*.SLD, *.sld)|*.sld",
"txt files (*.TXT, *.txt)|*.txt",
"all files (*.*)|*.*"]
## List of allowed extensions
ext = ['.sld', '.SLD', '.txt', '.TXT', '.*']
[docs] def read(self, path):
"""
Load data file
:param path: file path
:return MagSLD: x, y, z, sld_n, sld_mx, sld_my, sld_mz
:raise RuntimeError: when the file can't be opened
:raise ValueError: when the length of the data vectors are inconsistent
"""
try:
pos_x = numpy.zeros(0)
pos_y = numpy.zeros(0)
pos_z = numpy.zeros(0)
sld_n = numpy.zeros(0)
sld_mx = numpy.zeros(0)
sld_my = numpy.zeros(0)
sld_mz = numpy.zeros(0)
try:
# Use numpy to speed up loading
input_f = numpy.loadtxt(path, dtype='float', skiprows=1,
ndmin=1, unpack=True)
pos_x = numpy.array(input_f[0])
pos_y = numpy.array(input_f[1])
pos_z = numpy.array(input_f[2])
sld_n = numpy.array(input_f[3])
sld_mx = numpy.array(input_f[4])
sld_my = numpy.array(input_f[5])
sld_mz = numpy.array(input_f[6])
ncols = len(input_f)
if ncols == 8:
vol_pix = numpy.array(input_f[7])
elif ncols == 7:
vol_pix = None
except:
# For older version of numpy
input_f = open(path, 'rb')
buff = input_f.read()
lines = buff.split('\n')
input_f.close()
for line in lines:
toks = line.split()
try:
_pos_x = float(toks[0])
_pos_y = float(toks[1])
_pos_z = float(toks[2])
_sld_n = float(toks[3])
_sld_mx = float(toks[4])
_sld_my = float(toks[5])
_sld_mz = float(toks[6])
pos_x = numpy.append(pos_x, _pos_x)
pos_y = numpy.append(pos_y, _pos_y)
pos_z = numpy.append(pos_z, _pos_z)
sld_n = numpy.append(sld_n, _sld_n)
sld_mx = numpy.append(sld_mx, _sld_mx)
sld_my = numpy.append(sld_my, _sld_my)
sld_mz = numpy.append(sld_mz, _sld_mz)
try:
_vol_pix = float(toks[7])
vol_pix = numpy.append(vol_pix, _vol_pix)
except:
vol_pix = None
except:
# Skip non-data lines
logging.error(sys.exc_value)
output = MagSLD(pos_x, pos_y, pos_z, sld_n,
sld_mx, sld_my, sld_mz)
output.filename = os.path.basename(path)
output.set_pix_type('pixel')
output.set_pixel_symbols('pixel')
if vol_pix != None:
output.set_pixel_volumes(vol_pix)
return output
except:
raise RuntimeError, "%s is not a sld file" % path
[docs] def write(self, path, data):
"""
Write sld file
:Param path: file path
:Param data: MagSLD data object
"""
if path == None:
raise ValueError, "Missing the file path."
if data == None:
raise ValueError, "Missing the data to save."
x_val = data.pos_x
y_val = data.pos_y
z_val = data.pos_z
vol_pix = data.vol_pix
length = len(x_val)
sld_n = data.sld_n
if sld_n == None:
sld_n = numpy.zeros(length)
sld_mx = data.sld_mx
if sld_mx == None:
sld_mx = numpy.zeros(length)
sld_my = numpy.zeros(length)
sld_mz = numpy.zeros(length)
else:
sld_my = data.sld_my
sld_mz = data.sld_mz
out = open(path, 'w')
# First Line: Column names
out.write("X Y Z SLDN SLDMx SLDMy SLDMz VOLUMEpix")
for ind in range(length):
out.write("\n%g %g %g %g %g %g %g %g" % \
(x_val[ind], y_val[ind], z_val[ind], sld_n[ind],
sld_mx[ind], sld_my[ind], sld_mz[ind], vol_pix[ind]))
out.close()
[docs]class OMFData(object):
"""
OMF Data.
"""
_meshunit = "A"
_valueunit = "A^(-2)"
def __init__(self):
"""
Init for mag SLD
"""
self.filename = 'default'
self.oommf = ''
self.title = ''
self.desc = ''
self.meshtype = ''
self.meshunit = self._meshunit
self.valueunit = self._valueunit
self.xbase = 0.0
self.ybase = 0.0
self.zbase = 0.0
self.xstepsize = 6.0
self.ystepsize = 6.0
self.zstepsize = 6.0
self.xnodes = 10.0
self.ynodes = 10.0
self.znodes = 10.0
self.xmin = 0.0
self.ymin = 0.0
self.zmin = 0.0
self.xmax = 60.0
self.ymax = 60.0
self.zmax = 60.0
self.mx = None
self.my = None
self.mz = None
self.valuemultiplier = 1.
self.valuerangeminmag = 0
self.valuerangemaxmag = 0
def __str__(self):
"""
doc strings
"""
_str = "Type: %s\n" % self.__class__.__name__
_str += "File: %s\n" % self.filename
_str += "OOMMF: %s\n" % self.oommf
_str += "Title: %s\n" % self.title
_str += "Desc: %s\n" % self.desc
_str += "meshtype: %s\n" % self.meshtype
_str += "meshunit: %s\n" % str(self.meshunit)
_str += "xbase: %s [%s]\n" % (str(self.xbase), self.meshunit)
_str += "ybase: %s [%s]\n" % (str(self.ybase), self.meshunit)
_str += "zbase: %s [%s]\n" % (str(self.zbase), self.meshunit)
_str += "xstepsize: %s [%s]\n" % (str(self.xstepsize),
self.meshunit)
_str += "ystepsize: %s [%s]\n" % (str(self.ystepsize),
self.meshunit)
_str += "zstepsize: %s [%s]\n" % (str(self.zstepsize),
self.meshunit)
_str += "xnodes: %s\n" % str(self.xnodes)
_str += "ynodes: %s\n" % str(self.ynodes)
_str += "znodes: %s\n" % str(self.znodes)
_str += "xmin: %s [%s]\n" % (str(self.xmin), self.meshunit)
_str += "ymin: %s [%s]\n" % (str(self.ymin), self.meshunit)
_str += "zmin: %s [%s]\n" % (str(self.zmin), self.meshunit)
_str += "xmax: %s [%s]\n" % (str(self.xmax), self.meshunit)
_str += "ymax: %s [%s]\n" % (str(self.ymax), self.meshunit)
_str += "zmax: %s [%s]\n" % (str(self.zmax), self.meshunit)
_str += "valueunit: %s\n" % self.valueunit
_str += "valuemultiplier: %s\n" % str(self.valuemultiplier)
_str += "ValueRangeMinMag:%s [%s]\n" % (str(self.valuerangeminmag),
self.valueunit)
_str += "ValueRangeMaxMag:%s [%s]\n" % (str(self.valuerangemaxmag),
self.valueunit)
return _str
[docs] def set_m(self, mx, my, mz):
"""
Set the Mx, My, Mz values
"""
self.mx = mx
self.my = my
self.mz = mz
[docs]class MagSLD(object):
"""
Magnetic SLD.
"""
pos_x = None
pos_y = None
pos_z = None
sld_n = None
sld_mx = None
sld_my = None
sld_mz = None
# Units
_pos_unit = 'A'
_sld_unit = '1/A^(2)'
_pix_type = 'pixel'
def __init__(self, pos_x, pos_y, pos_z, sld_n=None,
sld_mx=None, sld_my=None, sld_mz=None, vol_pix=None):
"""
Init for mag SLD
:params : All should be numpy 1D array
"""
self.is_data = True
self.filename = ''
self.xstepsize = 6.0
self.ystepsize = 6.0
self.zstepsize = 6.0
self.xnodes = 10.0
self.ynodes = 10.0
self.znodes = 10.0
self.has_stepsize = False
self.has_conect = False
self.pos_unit = self._pos_unit
self.sld_unit = self._sld_unit
self.pix_type = 'pixel'
self.pos_x = pos_x
self.pos_y = pos_y
self.pos_z = pos_z
self.sld_n = sld_n
self.line_x = None
self.line_y = None
self.line_z = None
self.sld_mx = sld_mx
self.sld_my = sld_my
self.sld_mz = sld_mz
self.vol_pix = vol_pix
self.sld_m = None
self.sld_phi = None
self.sld_theta = None
self.pix_symbol = None
if sld_mx != None and sld_my != None and sld_mz != None:
self.set_sldms(sld_mx, sld_my, sld_mz)
self.set_nodes()
def __str__(self):
"""
doc strings
"""
_str = "Type: %s\n" % self.__class__.__name__
_str += "File: %s\n" % self.filename
_str += "Axis_unit: %s\n" % self.pos_unit
_str += "SLD_unit: %s\n" % self.sld_unit
return _str
[docs] def set_pix_type(self, pix_type):
"""
Set pixel type
:Param pix_type: string, 'pixel' or 'atom'
"""
self.pix_type = pix_type
[docs] def set_sldn(self, sld_n):
"""
Sets neutron SLD
"""
if sld_n.__class__.__name__ == 'float':
if self.is_data:
# For data, put the value to only the pixels w non-zero M
is_nonzero = (numpy.fabs(self.sld_mx) +
numpy.fabs(self.sld_my) +
numpy.fabs(self.sld_mz)).nonzero()
self.sld_n = numpy.zeros(len(self.pos_x))
if len(self.sld_n[is_nonzero]) > 0:
self.sld_n[is_nonzero] = sld_n
else:
self.sld_n.fill(sld_n)
else:
# For non-data, put the value to all the pixels
self.sld_n = numpy.ones(len(self.pos_x)) * sld_n
else:
self.sld_n = sld_n
[docs] def set_sldms(self, sld_mx, sld_my, sld_mz):
r"""
Sets (\|m\|, m_theta, m_phi)
"""
if sld_mx.__class__.__name__ == 'float':
self.sld_mx = numpy.ones(len(self.pos_x)) * sld_mx
else:
self.sld_mx = sld_mx
if sld_my.__class__.__name__ == 'float':
self.sld_my = numpy.ones(len(self.pos_x)) * sld_my
else:
self.sld_my = sld_my
if sld_mz.__class__.__name__ == 'float':
self.sld_mz = numpy.ones(len(self.pos_x)) * sld_mz
else:
self.sld_mz = sld_mz
sld_m = numpy.sqrt(sld_mx * sld_mx + sld_my * sld_my + \
sld_mz * sld_mz)
self.sld_m = sld_m
[docs] def set_pixel_symbols(self, symbol='pixel'):
"""
Set pixel
:Params pixel: str; pixel or atomic symbol, or array of strings
"""
if self.sld_n == None:
return
if symbol.__class__.__name__ == 'str':
self.pix_symbol = numpy.repeat(symbol, len(self.sld_n))
else:
self.pix_symbol = symbol
[docs] def set_pixel_volumes(self, vol):
"""
Set pixel volumes
:Params pixel: str; pixel or atomic symbol, or array of strings
"""
if self.sld_n == None:
return
if vol.__class__.__name__ == 'ndarray':
self.vol_pix = vol
elif vol.__class__.__name__.count('float') > 0:
self.vol_pix = numpy.repeat(vol, len(self.sld_n))
else:
self.vol_pix = None
[docs] def get_sldn(self):
"""
Returns nuclear sld
"""
return self.sld_n
[docs] def set_nodes(self):
"""
Set xnodes, ynodes, and znodes
"""
self.set_stepsize()
if self.pix_type == 'pixel':
try:
xdist = (max(self.pos_x) - min(self.pos_x)) / self.xstepsize
ydist = (max(self.pos_y) - min(self.pos_y)) / self.ystepsize
zdist = (max(self.pos_z) - min(self.pos_z)) / self.zstepsize
self.xnodes = int(xdist) + 1
self.ynodes = int(ydist) + 1
self.znodes = int(zdist) + 1
except:
self.xnodes = None
self.ynodes = None
self.znodes = None
else:
self.xnodes = None
self.ynodes = None
self.znodes = None
[docs] def set_stepsize(self):
"""
Set xtepsize, ystepsize, and zstepsize
"""
if self.pix_type == 'pixel':
try:
xpos_pre = self.pos_x[0]
ypos_pre = self.pos_y[0]
zpos_pre = self.pos_z[0]
for x_pos in self.pos_x:
if xpos_pre != x_pos:
self.xstepsize = numpy.fabs(x_pos - xpos_pre)
break
for y_pos in self.pos_y:
if ypos_pre != y_pos:
self.ystepsize = numpy.fabs(y_pos - ypos_pre)
break
for z_pos in self.pos_z:
if zpos_pre != z_pos:
self.zstepsize = numpy.fabs(z_pos - zpos_pre)
break
#default pix volume
self.vol_pix = numpy.ones(len(self.pos_x))
vol = self.xstepsize * self.ystepsize * self.zstepsize
self.set_pixel_volumes(vol)
self.has_stepsize = True
except:
self.xstepsize = None
self.ystepsize = None
self.zstepsize = None
self.vol_pix = None
self.has_stepsize = False
else:
self.xstepsize = None
self.ystepsize = None
self.zstepsize = None
self.has_stepsize = True
return self.xstepsize, self.ystepsize, self.zstepsize
[docs] def set_conect_lines(self, line_x, line_y, line_z):
"""
Set bonding line data if taken from pdb
"""
if line_x.__class__.__name__ != 'list' or len(line_x) < 1:
return
if line_y.__class__.__name__ != 'list' or len(line_y) < 1:
return
if line_z.__class__.__name__ != 'list' or len(line_z) < 1:
return
self.has_conect = True
self.line_x = line_x
self.line_y = line_y
self.line_z = line_z
[docs]def test_load():
"""
Test code
"""
from mpl_toolkits.mplot3d import Axes3D
current_dir = os.path.abspath(os.path.curdir)
print current_dir
for i in range(6):
current_dir, _ = os.path.split(current_dir)
tfile = os.path.join(current_dir, "test", "CoreXY_ShellZ.txt")
ofile = os.path.join(current_dir, "test", "A_Raw_Example-1.omf")
if os.path.isfile(tfile):
tfpath = tfile
ofpath = ofile
break
reader = SLDReader()
oreader = OMFReader()
output = reader.read(tfpath)
ooutput = oreader.read(ofpath)
foutput = OMF2SLD()
foutput.set_data(ooutput)
import matplotlib.pyplot as plt
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(output.pos_x, output.pos_y, output.pos_z, '.', c="g",
alpha=0.7, markeredgecolor='gray', rasterized=True)
gap = 7
max_mx = max(output.sld_mx)
max_my = max(output.sld_my)
max_mz = max(output.sld_mz)
max_m = max(max_mx, max_my, max_mz)
x2 = output.pos_x+output.sld_mx/max_m * gap
y2 = output.pos_y+output.sld_my/max_m * gap
z2 = output.pos_z+output.sld_mz/max_m * gap
x_arrow = numpy.column_stack((output.pos_x, x2))
y_arrow = numpy.column_stack((output.pos_y, y2))
z_arrow = numpy.column_stack((output.pos_z, z2))
unit_x2 = output.sld_mx / max_m
unit_y2 = output.sld_my / max_m
unit_z2 = output.sld_mz / max_m
color_x = numpy.fabs(unit_x2 * 0.8)
color_y = numpy.fabs(unit_y2 * 0.8)
color_z = numpy.fabs(unit_z2 * 0.8)
colors = numpy.column_stack((color_x, color_y, color_z))
plt.show()
[docs]def test():
"""
Test code
"""
current_dir = os.path.abspath(os.path.curdir)
for i in range(3):
current_dir, _ = os.path.split(current_dir)
ofile = os.path.join(current_dir, "test", "A_Raw_Example-1.omf")
if os.path.isfile(ofile):
ofpath = ofile
break
oreader = OMFReader()
ooutput = oreader.read(ofpath)
foutput = OMF2SLD()
foutput.set_data(ooutput)
writer = SLDReader()
writer.write(os.path.join(os.path.dirname(ofpath), "out.txt"),
foutput.output)
model = GenSAS()
model.set_sld_data(foutput.output)
x = numpy.arange(1000)/10000. + 1e-5
y = numpy.arange(1000)/10000. + 1e-5
i = numpy.zeros(1000)
model.runXY([x, y, i])
if __name__ == "__main__":
test()
test_load()