source: sasview/src/sas/calculator/sans_gen.py @ 79492222

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 79492222 was 79492222, checked in by krzywon, 9 years ago

Changed the file and folder names to remove all SANS references.

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