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

magnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since aab23e3 was a1daf86, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

hack around broken isfinite/isnan in tinycc

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