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

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 9c0f3c17 was 9c0f3c17, checked in by Ricardo Ferraz Leal <ricleal@…>, 7 years ago

After merge conflict

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