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

magnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1249
Last change on this file since 7af652d was 952ea1f, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

move c extensions from core/module.so to _module.so

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