Changes in / [e61a9b1:fc6651c] in sasview


Ignore:
Files:
2 added
30 deleted
14 edited

Legend:

Unmodified
Added
Removed
  • check_packages.py

    r3e1f417 r235f514  
    7272    for package_name, test_vals in win_required_package_list.items(): 
    7373        try: 
    74             i = __import__(test_vals['import_name'], fromlist=['']) 
    75             print("%s Version Installed: %s"% (package_name, getattr(i, test_vals['test'], "unknown"))) 
     74            if package_name == "pywin": 
     75                import win32api 
     76                fixed_file_info = win32api.GetFileVersionInfo(win32api.__file__, '\\') 
     77                print("%s Version Installed: %s"% (package_name, fixed_file_info['FileVersionLS'] >> 16)) 
     78            else: 
     79                i = __import__(test_vals['import_name'], fromlist=['']) 
     80                print("%s Version Installed: %s"% (package_name, getattr(i, test_vals['test']))) 
    7681        except ImportError: 
    7782            print('%s NOT INSTALLED'% package_name) 
  • run.py

    rf94a935 rf36e01f  
    5151def import_package(modname, path): 
    5252    """Import a package into a particular point in the python namespace""" 
    53     #logger.debug("Dynamicly importing: %s", path) 
    5453    mod = imp.load_source(modname, abspath(joinpath(path, '__init__.py'))) 
    5554    sys.modules[modname] = mod 
  • sasview/test/sesans_data/sphere2micron.ses

    rd3bce8c r2a2b43a  
    1 FileFormatVersion       1.0 
    2 DataFileTitle           Polystyrene of Markus Strobl,  Full Sine, ++ only 
    3 Sample                  Polystyrene 2 um in 53% H2O, 47% D2O 
    4 Settings                D1=D2=20x8 mm,Ds = 16x10 mm (WxH), GF1 =scanning, GF2 = 2.5 A. 2 um polystyrene in 53% H2O, 47% D2O; 8.55% contrast 
    5 Operator                CPD 
    6 Date                    do 10 jul 2014 16:37:30 
    7 ScanType                sine one element scan 
    8 Thickness               2.00E-01 
    9 Thickness_unit          cm 
    10 Theta_zmax              0.0168 
    11 Theta_zmax_unit         radians 
    12 Theta_ymax              0.0168 
    13 Theta_ymax_unit         radians 
    14 Orientation             Z 
    15 SpinEchoLength_unit     A 
    16 Depolarisation_unit     A-2 cm-1 
    17 Wavelength_unit         A 
    18  
    19 BEGIN_DATA 
    20 SpinEchoLength Depolarisation Depolarisation_error SpinEchoLength_error Wavelength Wavelength_error Polarisation  Polarisation_error 
    21 391.56 0.0041929 0.0036894 19.578 2.11 0.1055 1.0037 0.0032974 
    22 1564 -0.0046571 0.0038185 78.2 2.11 0.1055 0.99586 0.003386 
    23 2735.6 -0.017007 0.0038132 136.78 2.11 0.1055 0.98497 0.0033444 
    24 3907.9 -0.033462 0.0035068 195.39 2.11 0.1055 0.97064 0.0030309 
    25 5080.2 -0.047483 0.0038208 254.01 2.11 0.1055 0.9586 0.0032613 
    26 6251.8 -0.070375 0.00376 312.59 2.11 0.1055 0.93926 0.0031446 
    27 7423.2 -0.092217 0.0037927 371.16 2.11 0.1055 0.92117 0.0031108 
    28 8595.5 -0.10238 0.004006 429.77 2.11 0.1055 0.91287 0.0032562 
    29 9767.7 -0.12672 0.0038534 488.39 2.11 0.1055 0.8933 0.0030651 
    30 10940 -0.1374 0.004243 546.98 2.11 0.1055 0.88484 0.003343 
    31 12112 -0.16072 0.0045837 605.58 2.11 0.1055 0.86666 0.0035372 
    32 13284 -0.16623 0.0045613 664.2 2.11 0.1055 0.86242 0.0035027 
    33 14456 -0.18468 0.0044918 722.79 2.11 0.1055 0.84837 0.0033931 
    34 15628 -0.19143 0.0048967 781.38 2.11 0.1055 0.84328 0.0036768 
    35 16800 -0.20029 0.0045421 840.02 2.11 0.1055 0.83666 0.0033837 
    36 17971 -0.19798 0.0046642 898.56 2.11 0.1055 0.83838 0.0034819 
    37 19143 -0.21442 0.0047052 957.17 2.11 0.1055 0.82619 0.0034614 
    38 20316 -0.20885 0.0044931 1015.8 2.11 0.1055 0.8303 0.0033218 
    39 21488 -0.21393 0.0049186 1074.4 2.11 0.1055 0.82655 0.00362 
    40 22660 -0.20685 0.004423 1133 2.11 0.1055 0.83179 0.0032758 
    41 23832 -0.20802 0.0046979 1191.6 2.11 0.1055 0.83092 0.0034758 
    42 25003 -0.19848 0.0045953 1250.2 2.11 0.1055 0.838 0.0034289 
    43 26175 -0.21117 0.0044567 1308.8 2.11 0.1055 0.82859 0.0032881 
    44 27347 -0.21283 0.004137 1367.4 2.11 0.1055 0.82736 0.0030477 
    45 28520 -0.2042 0.0044587 1426 2.11 0.1055 0.83375 0.0033101 
    46 29692 -0.2112 0.0042852 1484.6 2.11 0.1055 0.82857 0.0031615 
    47 30864 -0.20319 0.0043483 1543.2 2.11 0.1055 0.8345 0.003231 
    48 32036 -0.20752 0.0044297 1601.8 2.11 0.1055 0.83129 0.0032788 
    49 33207 -0.20654 0.0043188 1660.4 2.11 0.1055 0.83201 0.0031995 
    50 34380 -0.20126 0.0046375 1719 2.11 0.1055 0.83593 0.0034518 
    51 35551 -0.20924 0.0042871 1777.6 2.11 0.1055 0.83001 0.0031684 
    52 36724 -0.21323 0.0045471 1836.2 2.11 0.1055 0.82707 0.0033487 
    53 37895 -0.21324 0.0045354 1894.7 2.11 0.1055 0.82706 0.00334 
    54 39067 -0.19905 0.0044141 1953.4 2.11 0.1055 0.83758 0.003292 
    55 40239 -0.1991 0.0047441 2012 2.11 0.1055 0.83754 0.003538 
    56 41411 -0.20359 0.0050136 2070.5 2.11 0.1055 0.8342 0.003724 
    57 42583 -0.21032 0.0049474 2129.1 2.11 0.1055 0.82922 0.0036529 
    58 43755 -0.20689 0.0048203 2187.8 2.11 0.1055 0.83176 0.00357 
    59 44927 -0.21075 0.0052337 2246.4 2.11 0.1055 0.8289 0.0038628 
    60 46099 -0.19956 0.0047827 2304.9 2.11 0.1055 0.8372 0.0035653 
     1DataFileTitle   "Polystyrene of Markus Strobl,  Full Sine, ++ only "                                             
     2Sample  "Polystyrene 2 um in 53% H2O, 47% D2O "                                          
     3Settings        "D1=D2=20x8 mm,Ds = 16x10 mm (WxH), GF1 =scanning, GF2 = 2.5 A. 2 um polystyrene in 53% H2O, 47% D2O; 8.55% contrast "                                           
     4Operator        CPD                                              
     5Date    do 10 jul 2014 16:37:30                                                  
     6ScanType        sine one element scan                                            
     7Thickness [cm]  2.00E-01                                                 
     8Q_zmax [\A^-1]  0.05                                             
     9Q_ymax [\A^-1]  0.05                                             
     10                                                         
     11spin echo length [A]     Depolarisation [A-2 cm-1]       error depol [A-2 cm-1]          error SEL [A]   wavelength [A]          error wavelength [A]    polarisation    error pol  
     12391.56  0.0041929       0.0036894       19.578  2.11    0.1055  1.0037  0.0032974 
     131564    -0.0046571      0.0038185       78.2    2.11    0.1055  0.99586 0.003386 
     142735.6  -0.017007       0.0038132       136.78  2.11    0.1055  0.98497 0.0033444 
     153907.9  -0.033462       0.0035068       195.39  2.11    0.1055  0.97064 0.0030309 
     165080.2  -0.047483       0.0038208       254.01  2.11    0.1055  0.9586  0.0032613 
     176251.8  -0.070375       0.00376 312.59  2.11    0.1055  0.93926 0.0031446 
     187423.2  -0.092217       0.0037927       371.16  2.11    0.1055  0.92117 0.0031108 
     198595.5  -0.10238        0.004006        429.77  2.11    0.1055  0.91287 0.0032562 
     209767.7  -0.12672        0.0038534       488.39  2.11    0.1055  0.8933  0.0030651 
     2110940   -0.1374 0.004243        546.98  2.11    0.1055  0.88484 0.003343 
     2212112   -0.16072        0.0045837       605.58  2.11    0.1055  0.86666 0.0035372 
     2313284   -0.16623        0.0045613       664.2   2.11    0.1055  0.86242 0.0035027 
     2414456   -0.18468        0.0044918       722.79  2.11    0.1055  0.84837 0.0033931 
     2515628   -0.19143        0.0048967       781.38  2.11    0.1055  0.84328 0.0036768 
     2616800   -0.20029        0.0045421       840.02  2.11    0.1055  0.83666 0.0033837 
     2717971   -0.19798        0.0046642       898.56  2.11    0.1055  0.83838 0.0034819 
     2819143   -0.21442        0.0047052       957.17  2.11    0.1055  0.82619 0.0034614 
     2920316   -0.20885        0.0044931       1015.8  2.11    0.1055  0.8303  0.0033218 
     3021488   -0.21393        0.0049186       1074.4  2.11    0.1055  0.82655 0.00362 
     3122660   -0.20685        0.004423        1133    2.11    0.1055  0.83179 0.0032758 
     3223832   -0.20802        0.0046979       1191.6  2.11    0.1055  0.83092 0.0034758 
     3325003   -0.19848        0.0045953       1250.2  2.11    0.1055  0.838   0.0034289 
     3426175   -0.21117        0.0044567       1308.8  2.11    0.1055  0.82859 0.0032881 
     3527347   -0.21283        0.004137        1367.4  2.11    0.1055  0.82736 0.0030477 
     3628520   -0.2042 0.0044587       1426    2.11    0.1055  0.83375 0.0033101 
     3729692   -0.2112 0.0042852       1484.6  2.11    0.1055  0.82857 0.0031615 
     3830864   -0.20319        0.0043483       1543.2  2.11    0.1055  0.8345  0.003231 
     3932036   -0.20752        0.0044297       1601.8  2.11    0.1055  0.83129 0.0032788 
     4033207   -0.20654        0.0043188       1660.4  2.11    0.1055  0.83201 0.0031995 
     4134380   -0.20126        0.0046375       1719    2.11    0.1055  0.83593 0.0034518 
     4235551   -0.20924        0.0042871       1777.6  2.11    0.1055  0.83001 0.0031684 
     4336724   -0.21323        0.0045471       1836.2  2.11    0.1055  0.82707 0.0033487 
     4437895   -0.21324        0.0045354       1894.7  2.11    0.1055  0.82706 0.00334 
     4539067   -0.19905        0.0044141       1953.4  2.11    0.1055  0.83758 0.003292 
     4640239   -0.1991 0.0047441       2012    2.11    0.1055  0.83754 0.003538 
     4741411   -0.20359        0.0050136       2070.5  2.11    0.1055  0.8342  0.003724 
     4842583   -0.21032        0.0049474       2129.1  2.11    0.1055  0.82922 0.0036529 
     4943755   -0.20689        0.0048203       2187.8  2.11    0.1055  0.83176 0.00357 
     5044927   -0.21075        0.0052337       2246.4  2.11    0.1055  0.8289  0.0038628 
     5146099   -0.19956        0.0047827       2304.9  2.11    0.1055  0.8372  0.0035653 
  • src/sas/sascalc/data_util/nxsunit.py

    r8e9536f rb699768  
    136136    sld = { '10^-6 Angstrom^-2': 1e-6, 'Angstrom^-2': 1 } 
    137137    Q = { 'invA': 1, 'invAng': 1, 'invAngstroms': 1, '1/A': 1,  
    138           '10^-3 Angstrom^-1': 1e-3, '1/cm': 1e-8, '1/m': 1e-10, 
     138          '10^-3 Angstrom^-1': 1e-3, '1/cm': 1e-8, 
    139139          'nm^-1': 0.1, '1/nm': 0.1, 'n_m^-1': 0.1 } 
    140140 
     
    195195    _check(1,Converter('n_m^-1')(10,'invA')) # 10 nm^-1 = 1 inv Angstroms 
    196196    _check(2,Converter('mm')(2000,'m')) # 2000 mm -> 2 m 
    197     _check(2.011e10,Converter('1/A')(2.011,"1/m")) # 2.011 1/A -> 2.011 * 10^10 1/m 
    198197    _check(0.003,Converter('microseconds')(3,units='ms')) # 3 us -> 0.003 ms 
    199198    _check(45,Converter('nanokelvin')(45))  # 45 nK -> 45 nK 
  • src/sas/sascalc/data_util/qsmearing.py

    r235f514 r235f514  
    6565            raise ValueError('one or more of your dx values are negative, please check the data file!') 
    6666 
    67     if _found_sesans: 
    68         # Pre-compute the Hankel matrix (H) 
     67    if _found_sesans == True: 
     68        #Pre-compute the Hankel matrix (H) 
     69        qmax, qunits = data.sample.zacceptance 
    6970        SElength = Converter(data._xunit)(data.x, "A") 
    70  
    71         theta_max = Converter("radians")(data.sample.zacceptance)[0] 
    72         q_max = 2 * np.pi / np.max(data.source.wavelength) * np.sin(theta_max) 
    73         zaccept = Converter("1/A")(q_max, "1/" + data.source.wavelength_unit), 
    74  
     71        zaccept = Converter(qunits)(qmax, "1/A"), 
    7572        Rmax = 10000000 
    76         hankel = SesansTransform(data.x, SElength, 
    77                                  data.source.wavelength, 
    78                                  zaccept, Rmax) 
     73        hankel = SesansTransform(data.x, SElength, zaccept, Rmax) 
    7974        # Then return the actual transform, as if it were a smearing function 
    8075        return PySmear(hankel, model, offset=0) 
  • src/sas/sascalc/dataloader/manipulations.py

    r324e0bf r7432acb  
    1 from __future__ import division 
    21""" 
    32Data manipulations for 2D data sets. 
    43Using the meta data information, various types of averaging 
    54are performed in Q-space 
    6  
    7 To test this module use: 
    8 ``` 
    9 cd test 
    10 PYTHONPATH=../src/ python2  -m sasdataloader.test.utest_averaging DataInfoTests.test_sectorphi_quarter 
    11 ``` 
    125""" 
    136##################################################################### 
    14 # This software was developed by the University of Tennessee as part of the 
    15 # Distributed Data Analysis of Neutron Scattering Experiments (DANSE) 
    16 # project funded by the US National Science Foundation. 
    17 # See the license text in license.txt 
    18 # copyright 2008, University of Tennessee 
     7#This software was developed by the University of Tennessee as part of the 
     8#Distributed Data Analysis of Neutron Scattering Experiments (DANSE) 
     9#project funded by the US National Science Foundation. 
     10#See the license text in license.txt 
     11#copyright 2008, University of Tennessee 
    1912###################################################################### 
    2013 
    21  
    22 # TODO: copy the meta data from the 2D object to the resulting 1D object 
     14#TODO: copy the meta data from the 2D object to the resulting 1D object 
    2315import math 
    24 import numpy as np 
    25 import sys 
     16import numpy 
    2617 
    2718#from data_info import plottable_2D 
     
    7970    return phi_out 
    8071 
     72 
     73def reader2D_converter(data2d=None): 
     74    """ 
     75    convert old 2d format opened by IhorReader or danse_reader 
     76    to new Data2D format 
     77 
     78    :param data2d: 2d array of Data2D object 
     79    :return: 1d arrays of Data2D object 
     80 
     81    """ 
     82    if data2d.data is None or data2d.x_bins is None or data2d.y_bins is None: 
     83        raise ValueError, "Can't convert this data: data=None..." 
     84    new_x = numpy.tile(data2d.x_bins, (len(data2d.y_bins), 1)) 
     85    new_y = numpy.tile(data2d.y_bins, (len(data2d.x_bins), 1)) 
     86    new_y = new_y.swapaxes(0, 1) 
     87 
     88    new_data = data2d.data.flatten() 
     89    qx_data = new_x.flatten() 
     90    qy_data = new_y.flatten() 
     91    q_data = numpy.sqrt(qx_data * qx_data + qy_data * qy_data) 
     92    if data2d.err_data is None or numpy.any(data2d.err_data <= 0): 
     93        new_err_data = numpy.sqrt(numpy.abs(new_data)) 
     94    else: 
     95        new_err_data = data2d.err_data.flatten() 
     96    mask = numpy.ones(len(new_data), dtype=bool) 
     97 
     98    #TODO: make sense of the following two lines... 
     99    #from sas.sascalc.dataloader.data_info import Data2D 
     100    #output = Data2D() 
     101    output = data2d 
     102    output.data = new_data 
     103    output.err_data = new_err_data 
     104    output.qx_data = qx_data 
     105    output.qy_data = qy_data 
     106    output.q_data = q_data 
     107    output.mask = mask 
     108 
     109    return output 
     110 
     111 
     112class _Slab(object): 
     113    """ 
     114    Compute average I(Q) for a region of interest 
     115    """ 
     116    def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, 
     117                 y_max=0.0, bin_width=0.001): 
     118        # Minimum Qx value [A-1] 
     119        self.x_min = x_min 
     120        # Maximum Qx value [A-1] 
     121        self.x_max = x_max 
     122        # Minimum Qy value [A-1] 
     123        self.y_min = y_min 
     124        # Maximum Qy value [A-1] 
     125        self.y_max = y_max 
     126        # Bin width (step size) [A-1] 
     127        self.bin_width = bin_width 
     128        # If True, I(|Q|) will be return, otherwise, 
     129        # negative q-values are allowed 
     130        self.fold = False 
     131 
     132    def __call__(self, data2D): 
     133        return NotImplemented 
     134 
     135    def _avg(self, data2D, maj): 
     136        """ 
     137        Compute average I(Q_maj) for a region of interest. 
     138        The major axis is defined as the axis of Q_maj. 
     139        The minor axis is the axis that we average over. 
     140 
     141        :param data2D: Data2D object 
     142        :param maj_min: min value on the major axis 
     143        :return: Data1D object 
     144        """ 
     145        if len(data2D.detector) > 1: 
     146            msg = "_Slab._avg: invalid number of " 
     147            msg += " detectors: %g" % len(data2D.detector) 
     148            raise RuntimeError, msg 
     149 
     150        # Get data 
     151        data = data2D.data[numpy.isfinite(data2D.data)] 
     152        err_data = data2D.err_data[numpy.isfinite(data2D.data)] 
     153        qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 
     154        qy_data = data2D.qy_data[numpy.isfinite(data2D.data)] 
     155 
     156        # Build array of Q intervals 
     157        if maj == 'x': 
     158            if self.fold: 
     159                x_min = 0 
     160            else: 
     161                x_min = self.x_min 
     162            nbins = int(math.ceil((self.x_max - x_min) / self.bin_width)) 
     163        elif maj == 'y': 
     164            if self.fold: 
     165                y_min = 0 
     166            else: 
     167                y_min = self.y_min 
     168            nbins = int(math.ceil((self.y_max - y_min) / self.bin_width)) 
     169        else: 
     170            raise RuntimeError, "_Slab._avg: unrecognized axis %s" % str(maj) 
     171 
     172        x = numpy.zeros(nbins) 
     173        y = numpy.zeros(nbins) 
     174        err_y = numpy.zeros(nbins) 
     175        y_counts = numpy.zeros(nbins) 
     176 
     177        # Average pixelsize in q space 
     178        for npts in range(len(data)): 
     179            # default frac 
     180            frac_x = 0 
     181            frac_y = 0 
     182            # get ROI 
     183            if self.x_min <= qx_data[npts] and self.x_max > qx_data[npts]: 
     184                frac_x = 1 
     185            if self.y_min <= qy_data[npts] and self.y_max > qy_data[npts]: 
     186                frac_y = 1 
     187            frac = frac_x * frac_y 
     188 
     189            if frac == 0: 
     190                continue 
     191            # binning: find axis of q 
     192            if maj == 'x': 
     193                q_value = qx_data[npts] 
     194                min_value = x_min 
     195            if maj == 'y': 
     196                q_value = qy_data[npts] 
     197                min_value = y_min 
     198            if self.fold and q_value < 0: 
     199                q_value = -q_value 
     200            # bin 
     201            i_q = int(math.ceil((q_value - min_value) / self.bin_width)) - 1 
     202 
     203            # skip outside of max bins 
     204            if i_q < 0 or i_q >= nbins: 
     205                continue 
     206 
     207            #TODO: find better definition of x[i_q] based on q_data 
     208            # min_value + (i_q + 1) * self.bin_width / 2.0 
     209            x[i_q] += frac * q_value 
     210            y[i_q] += frac * data[npts] 
     211 
     212            if err_data is None or err_data[npts] == 0.0: 
     213                if data[npts] < 0: 
     214                    data[npts] = -data[npts] 
     215                err_y[i_q] += frac * frac * data[npts] 
     216            else: 
     217                err_y[i_q] += frac * frac * err_data[npts] * err_data[npts] 
     218            y_counts[i_q] += frac 
     219 
     220        # Average the sums 
     221        for n in range(nbins): 
     222            err_y[n] = math.sqrt(err_y[n]) 
     223 
     224        err_y = err_y / y_counts 
     225        y = y / y_counts 
     226        x = x / y_counts 
     227        idx = (numpy.isfinite(y) & numpy.isfinite(x)) 
     228 
     229        if not idx.any(): 
     230            msg = "Average Error: No points inside ROI to average..." 
     231            raise ValueError, msg 
     232        return Data1D(x=x[idx], y=y[idx], dy=err_y[idx]) 
     233 
     234 
     235class SlabY(_Slab): 
     236    """ 
     237    Compute average I(Qy) for a region of interest 
     238    """ 
     239    def __call__(self, data2D): 
     240        """ 
     241        Compute average I(Qy) for a region of interest 
     242 
     243        :param data2D: Data2D object 
     244        :return: Data1D object 
     245        """ 
     246        return self._avg(data2D, 'y') 
     247 
     248 
     249class SlabX(_Slab): 
     250    """ 
     251    Compute average I(Qx) for a region of interest 
     252    """ 
     253    def __call__(self, data2D): 
     254        """ 
     255        Compute average I(Qx) for a region of interest 
     256        :param data2D: Data2D object 
     257        :return: Data1D object 
     258        """ 
     259        return self._avg(data2D, 'x') 
     260 
     261 
     262class Boxsum(object): 
     263    """ 
     264    Perform the sum of counts in a 2D region of interest. 
     265    """ 
     266    def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0): 
     267        # Minimum Qx value [A-1] 
     268        self.x_min = x_min 
     269        # Maximum Qx value [A-1] 
     270        self.x_max = x_max 
     271        # Minimum Qy value [A-1] 
     272        self.y_min = y_min 
     273        # Maximum Qy value [A-1] 
     274        self.y_max = y_max 
     275 
     276    def __call__(self, data2D): 
     277        """ 
     278        Perform the sum in the region of interest 
     279 
     280        :param data2D: Data2D object 
     281        :return: number of counts, error on number of counts, 
     282            number of points summed 
     283        """ 
     284        y, err_y, y_counts = self._sum(data2D) 
     285 
     286        # Average the sums 
     287        counts = 0 if y_counts == 0 else y 
     288        error = 0 if y_counts == 0 else math.sqrt(err_y) 
     289 
     290        # Added y_counts to return, SMK & PDB, 04/03/2013 
     291        return counts, error, y_counts 
     292 
     293    def _sum(self, data2D): 
     294        """ 
     295        Perform the sum in the region of interest 
     296 
     297        :param data2D: Data2D object 
     298        :return: number of counts, 
     299            error on number of counts, number of entries summed 
     300        """ 
     301        if len(data2D.detector) > 1: 
     302            msg = "Circular averaging: invalid number " 
     303            msg += "of detectors: %g" % len(data2D.detector) 
     304            raise RuntimeError, msg 
     305        # Get data 
     306        data = data2D.data[numpy.isfinite(data2D.data)] 
     307        err_data = data2D.err_data[numpy.isfinite(data2D.data)] 
     308        qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 
     309        qy_data = data2D.qy_data[numpy.isfinite(data2D.data)] 
     310 
     311        y = 0.0 
     312        err_y = 0.0 
     313        y_counts = 0.0 
     314 
     315        # Average pixelsize in q space 
     316        for npts in range(len(data)): 
     317            # default frac 
     318            frac_x = 0 
     319            frac_y = 0 
     320 
     321            # get min and max at each points 
     322            qx = qx_data[npts] 
     323            qy = qy_data[npts] 
     324 
     325            # get the ROI 
     326            if self.x_min <= qx and self.x_max > qx: 
     327                frac_x = 1 
     328            if self.y_min <= qy and self.y_max > qy: 
     329                frac_y = 1 
     330            #Find the fraction along each directions 
     331            frac = frac_x * frac_y 
     332            if frac == 0: 
     333                continue 
     334            y += frac * data[npts] 
     335            if err_data is None or err_data[npts] == 0.0: 
     336                if data[npts] < 0: 
     337                    data[npts] = -data[npts] 
     338                err_y += frac * frac * data[npts] 
     339            else: 
     340                err_y += frac * frac * err_data[npts] * err_data[npts] 
     341            y_counts += frac 
     342        return y, err_y, y_counts 
     343 
     344 
     345class Boxavg(Boxsum): 
     346    """ 
     347    Perform the average of counts in a 2D region of interest. 
     348    """ 
     349    def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0): 
     350        super(Boxavg, self).__init__(x_min=x_min, x_max=x_max, 
     351                                     y_min=y_min, y_max=y_max) 
     352 
     353    def __call__(self, data2D): 
     354        """ 
     355        Perform the sum in the region of interest 
     356 
     357        :param data2D: Data2D object 
     358        :return: average counts, error on average counts 
     359 
     360        """ 
     361        y, err_y, y_counts = self._sum(data2D) 
     362 
     363        # Average the sums 
     364        counts = 0 if y_counts == 0 else y / y_counts 
     365        error = 0 if y_counts == 0 else math.sqrt(err_y) / y_counts 
     366 
     367        return counts, error 
     368 
     369 
    81370def get_pixel_fraction_square(x, xmin, xmax): 
    82371    """ 
     
    101390        return 1.0 
    102391 
    103 def get_intercept(q, q_0, q_1): 
    104     """ 
    105     Returns the fraction of the side at which the 
    106     q-value intercept the pixel, None otherwise. 
    107     The values returned is the fraction ON THE SIDE 
    108     OF THE LOWEST Q. :: 
    109  
    110             A           B 
    111         +-----------+--------+    <--- pixel size 
    112         0                    1 
    113         Q_0 -------- Q ----- Q_1   <--- equivalent Q range 
    114         if Q_1 > Q_0, A is returned 
    115         if Q_1 < Q_0, B is returned 
    116         if Q is outside the range of [Q_0, Q_1], None is returned 
    117  
    118     """ 
    119     if q_1 > q_0: 
    120         if q > q_0 and q <= q_1: 
    121             return (q - q_0) / (q_1 - q_0) 
    122     else: 
    123         if q > q_1 and q <= q_0: 
    124             return (q - q_1) / (q_0 - q_1) 
    125     return None 
     392 
     393class CircularAverage(object): 
     394    """ 
     395    Perform circular averaging on 2D data 
     396 
     397    The data returned is the distribution of counts 
     398    as a function of Q 
     399    """ 
     400    def __init__(self, r_min=0.0, r_max=0.0, bin_width=0.0005): 
     401        # Minimum radius included in the average [A-1] 
     402        self.r_min = r_min 
     403        # Maximum radius included in the average [A-1] 
     404        self.r_max = r_max 
     405        # Bin width (step size) [A-1] 
     406        self.bin_width = bin_width 
     407 
     408    def __call__(self, data2D, ismask=False): 
     409        """ 
     410        Perform circular averaging on the data 
     411 
     412        :param data2D: Data2D object 
     413        :return: Data1D object 
     414        """ 
     415        # Get data W/ finite values 
     416        data = data2D.data[numpy.isfinite(data2D.data)] 
     417        q_data = data2D.q_data[numpy.isfinite(data2D.data)] 
     418        err_data = data2D.err_data[numpy.isfinite(data2D.data)] 
     419        mask_data = data2D.mask[numpy.isfinite(data2D.data)] 
     420 
     421        dq_data = None 
     422 
     423        # Get the dq for resolution averaging 
     424        if data2D.dqx_data is not None and data2D.dqy_data is not None: 
     425            # The pinholes and det. pix contribution present 
     426            # in both direction of the 2D which must be subtracted when 
     427            # converting to 1D: dq_overlap should calculated ideally at 
     428            # q = 0. Note This method works on only pinhole geometry. 
     429            # Extrapolate dqx(r) and dqy(phi) at q = 0, and take an average. 
     430            z_max = max(data2D.q_data) 
     431            z_min = min(data2D.q_data) 
     432            x_max = data2D.dqx_data[data2D.q_data[z_max]] 
     433            x_min = data2D.dqx_data[data2D.q_data[z_min]] 
     434            y_max = data2D.dqy_data[data2D.q_data[z_max]] 
     435            y_min = data2D.dqy_data[data2D.q_data[z_min]] 
     436            # Find qdx at q = 0 
     437            dq_overlap_x = (x_min * z_max - x_max * z_min) / (z_max - z_min) 
     438            # when extrapolation goes wrong 
     439            if dq_overlap_x > min(data2D.dqx_data): 
     440                dq_overlap_x = min(data2D.dqx_data) 
     441            dq_overlap_x *= dq_overlap_x 
     442            # Find qdx at q = 0 
     443            dq_overlap_y = (y_min * z_max - y_max * z_min) / (z_max - z_min) 
     444            # when extrapolation goes wrong 
     445            if dq_overlap_y > min(data2D.dqy_data): 
     446                dq_overlap_y = min(data2D.dqy_data) 
     447            # get dq at q=0. 
     448            dq_overlap_y *= dq_overlap_y 
     449 
     450            dq_overlap = numpy.sqrt((dq_overlap_x + dq_overlap_y) / 2.0) 
     451            # Final protection of dq 
     452            if dq_overlap < 0: 
     453                dq_overlap = y_min 
     454            dqx_data = data2D.dqx_data[numpy.isfinite(data2D.data)] 
     455            dqy_data = data2D.dqy_data[numpy.isfinite(data2D.data)] - dq_overlap 
     456            # def; dqx_data = dq_r dqy_data = dq_phi 
     457            # Convert dq 2D to 1D here 
     458            dqx = dqx_data * dqx_data 
     459            dqy = dqy_data * dqy_data 
     460            dq_data = numpy.add(dqx, dqy) 
     461            dq_data = numpy.sqrt(dq_data) 
     462 
     463        #q_data_max = numpy.max(q_data) 
     464        if len(data2D.q_data) is None: 
     465            msg = "Circular averaging: invalid q_data: %g" % data2D.q_data 
     466            raise RuntimeError, msg 
     467 
     468        # Build array of Q intervals 
     469        nbins = int(math.ceil((self.r_max - self.r_min) / self.bin_width)) 
     470 
     471        x = numpy.zeros(nbins) 
     472        y = numpy.zeros(nbins) 
     473        err_y = numpy.zeros(nbins) 
     474        err_x = numpy.zeros(nbins) 
     475        y_counts = numpy.zeros(nbins) 
     476 
     477        for npt in range(len(data)): 
     478 
     479            if ismask and not mask_data[npt]: 
     480                continue 
     481 
     482            frac = 0 
     483 
     484            # q-value at the pixel (j,i) 
     485            q_value = q_data[npt] 
     486            data_n = data[npt] 
     487 
     488            ## No need to calculate the frac when all data are within range 
     489            if self.r_min >= self.r_max: 
     490                raise ValueError, "Limit Error: min > max" 
     491 
     492            if self.r_min <= q_value and q_value <= self.r_max: 
     493                frac = 1 
     494            if frac == 0: 
     495                continue 
     496            i_q = int(math.floor((q_value - self.r_min) / self.bin_width)) 
     497 
     498            # Take care of the edge case at phi = 2pi. 
     499            if i_q == nbins: 
     500                i_q = nbins - 1 
     501            y[i_q] += frac * data_n 
     502            # Take dqs from data to get the q_average 
     503            x[i_q] += frac * q_value 
     504            if err_data is None or err_data[npt] == 0.0: 
     505                if data_n < 0: 
     506                    data_n = -data_n 
     507                err_y[i_q] += frac * frac * data_n 
     508            else: 
     509                err_y[i_q] += frac * frac * err_data[npt] * err_data[npt] 
     510            if dq_data is not None: 
     511                # To be consistent with dq calculation in 1d reduction, 
     512                # we need just the averages (not quadratures) because 
     513                # it should not depend on the number of the q points 
     514                # in the qr bins. 
     515                err_x[i_q] += frac * dq_data[npt] 
     516            else: 
     517                err_x = None 
     518            y_counts[i_q] += frac 
     519 
     520        # Average the sums 
     521        for n in range(nbins): 
     522            if err_y[n] < 0: 
     523                err_y[n] = -err_y[n] 
     524            err_y[n] = math.sqrt(err_y[n]) 
     525            #if err_x is not None: 
     526            #    err_x[n] = math.sqrt(err_x[n]) 
     527 
     528        err_y = err_y / y_counts 
     529        err_y[err_y == 0] = numpy.average(err_y) 
     530        y = y / y_counts 
     531        x = x / y_counts 
     532        idx = (numpy.isfinite(y)) & (numpy.isfinite(x)) 
     533 
     534        if err_x is not None: 
     535            d_x = err_x[idx] / y_counts[idx] 
     536        else: 
     537            d_x = None 
     538 
     539        if not idx.any(): 
     540            msg = "Average Error: No points inside ROI to average..." 
     541            raise ValueError, msg 
     542 
     543        return Data1D(x=x[idx], y=y[idx], dy=err_y[idx], dx=d_x) 
     544 
     545 
     546class Ring(object): 
     547    """ 
     548    Defines a ring on a 2D data set. 
     549    The ring is defined by r_min, r_max, and 
     550    the position of the center of the ring. 
     551 
     552    The data returned is the distribution of counts 
     553    around the ring as a function of phi. 
     554 
     555    Phi_min and phi_max should be defined between 0 and 2*pi 
     556    in anti-clockwise starting from the x- axis on the left-hand side 
     557    """ 
     558    #Todo: remove center. 
     559    def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0, nbins=36): 
     560        # Minimum radius 
     561        self.r_min = r_min 
     562        # Maximum radius 
     563        self.r_max = r_max 
     564        # Center of the ring in x 
     565        self.center_x = center_x 
     566        # Center of the ring in y 
     567        self.center_y = center_y 
     568        # Number of angular bins 
     569        self.nbins_phi = nbins 
     570 
     571 
     572    def __call__(self, data2D): 
     573        """ 
     574        Apply the ring to the data set. 
     575        Returns the angular distribution for a given q range 
     576 
     577        :param data2D: Data2D object 
     578 
     579        :return: Data1D object 
     580        """ 
     581        if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: 
     582            raise RuntimeError, "Ring averaging only take plottable_2D objects" 
     583 
     584        Pi = math.pi 
     585 
     586        # Get data 
     587        data = data2D.data[numpy.isfinite(data2D.data)] 
     588        q_data = data2D.q_data[numpy.isfinite(data2D.data)] 
     589        err_data = data2D.err_data[numpy.isfinite(data2D.data)] 
     590        qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 
     591        qy_data = data2D.qy_data[numpy.isfinite(data2D.data)] 
     592 
     593        # Set space for 1d outputs 
     594        phi_bins = numpy.zeros(self.nbins_phi) 
     595        phi_counts = numpy.zeros(self.nbins_phi) 
     596        phi_values = numpy.zeros(self.nbins_phi) 
     597        phi_err = numpy.zeros(self.nbins_phi) 
     598 
     599        # Shift to apply to calculated phi values in order 
     600        # to center first bin at zero 
     601        phi_shift = Pi / self.nbins_phi 
     602 
     603        for npt in range(len(data)): 
     604            frac = 0 
     605            # q-value at the point (npt) 
     606            q_value = q_data[npt] 
     607            data_n = data[npt] 
     608 
     609            # phi-value at the point (npt) 
     610            phi_value = math.atan2(qy_data[npt], qx_data[npt]) + Pi 
     611 
     612            if self.r_min <= q_value and q_value <= self.r_max: 
     613                frac = 1 
     614            if frac == 0: 
     615                continue 
     616            # binning 
     617            i_phi = int(math.floor((self.nbins_phi) * \ 
     618                                   (phi_value + phi_shift) / (2 * Pi))) 
     619 
     620            # Take care of the edge case at phi = 2pi. 
     621            if i_phi >= self.nbins_phi: 
     622                i_phi = 0 
     623            phi_bins[i_phi] += frac * data[npt] 
     624 
     625            if err_data is None or err_data[npt] == 0.0: 
     626                if data_n < 0: 
     627                    data_n = -data_n 
     628                phi_err[i_phi] += frac * frac * math.fabs(data_n) 
     629            else: 
     630                phi_err[i_phi] += frac * frac * err_data[npt] * err_data[npt] 
     631            phi_counts[i_phi] += frac 
     632 
     633        for i in range(self.nbins_phi): 
     634            phi_bins[i] = phi_bins[i] / phi_counts[i] 
     635            phi_err[i] = math.sqrt(phi_err[i]) / phi_counts[i] 
     636            phi_values[i] = 2.0 * math.pi / self.nbins_phi * (1.0 * i) 
     637 
     638        idx = (numpy.isfinite(phi_bins)) 
     639 
     640        if not idx.any(): 
     641            msg = "Average Error: No points inside ROI to average..." 
     642            raise ValueError, msg 
     643        #elif len(phi_bins[idx])!= self.nbins_phi: 
     644        #    print "resulted",self.nbins_phi- len(phi_bins[idx]) 
     645        #,"empty bin(s) due to tight binning..." 
     646        return Data1D(x=phi_values[idx], y=phi_bins[idx], dy=phi_err[idx]) 
     647 
    126648 
    127649def get_pixel_fraction(qmax, q_00, q_01, q_10, q_11): 
     
    188710    return frac_max 
    189711 
    190 def get_dq_data(data2D): 
    191     ''' 
    192     Get the dq for resolution averaging 
    193     The pinholes and det. pix contribution present 
    194     in both direction of the 2D which must be subtracted when 
    195     converting to 1D: dq_overlap should calculated ideally at 
    196     q = 0. Note This method works on only pinhole geometry. 
    197     Extrapolate dqx(r) and dqy(phi) at q = 0, and take an average. 
    198     ''' 
    199     z_max = max(data2D.q_data) 
    200     z_min = min(data2D.q_data) 
    201     dqx_at_z_max = data2D.dqx_data[np.argmax(data2D.q_data)] 
    202     dqx_at_z_min = data2D.dqx_data[np.argmin(data2D.q_data)] 
    203     dqy_at_z_max = data2D.dqy_data[np.argmax(data2D.q_data)] 
    204     dqy_at_z_min = data2D.dqy_data[np.argmin(data2D.q_data)] 
    205     # Find qdx at q = 0 
    206     dq_overlap_x = (dqx_at_z_min * z_max - dqx_at_z_max * z_min) / (z_max - z_min) 
    207     # when extrapolation goes wrong 
    208     if dq_overlap_x > min(data2D.dqx_data): 
    209         dq_overlap_x = min(data2D.dqx_data) 
    210     dq_overlap_x *= dq_overlap_x 
    211     # Find qdx at q = 0 
    212     dq_overlap_y = (dqy_at_z_min * z_max - dqy_at_z_max * z_min) / (z_max - z_min) 
    213     # when extrapolation goes wrong 
    214     if dq_overlap_y > min(data2D.dqy_data): 
    215         dq_overlap_y = min(data2D.dqy_data) 
    216     # get dq at q=0. 
    217     dq_overlap_y *= dq_overlap_y 
    218  
    219     dq_overlap = np.sqrt((dq_overlap_x + dq_overlap_y) / 2.0) 
    220     # Final protection of dq 
    221     if dq_overlap < 0: 
    222         dq_overlap = dqy_at_z_min 
    223     dqx_data = data2D.dqx_data[np.isfinite(data2D.data)] 
    224     dqy_data = data2D.dqy_data[np.isfinite( 
    225         data2D.data)] - dq_overlap 
    226     # def; dqx_data = dq_r dqy_data = dq_phi 
    227     # Convert dq 2D to 1D here 
    228     dq_data = np.sqrt(dqx_data**2 + dqx_data**2) 
    229     return dq_data 
    230  
    231 ################################################################################ 
    232  
    233 def reader2D_converter(data2d=None): 
    234     """ 
    235     convert old 2d format opened by IhorReader or danse_reader 
    236     to new Data2D format 
    237     This is mainly used by the Readers 
    238  
    239     :param data2d: 2d array of Data2D object 
    240     :return: 1d arrays of Data2D object 
    241  
    242     """ 
    243     if data2d.data is None or data2d.x_bins is None or data2d.y_bins is None: 
    244         raise ValueError("Can't convert this data: data=None...") 
    245     new_x = np.tile(data2d.x_bins, (len(data2d.y_bins), 1)) 
    246     new_y = np.tile(data2d.y_bins, (len(data2d.x_bins), 1)) 
    247     new_y = new_y.swapaxes(0, 1) 
    248  
    249     new_data = data2d.data.flatten() 
    250     qx_data = new_x.flatten() 
    251     qy_data = new_y.flatten() 
    252     q_data = np.sqrt(qx_data * qx_data + qy_data * qy_data) 
    253     if data2d.err_data is None or np.any(data2d.err_data <= 0): 
    254         new_err_data = np.sqrt(np.abs(new_data)) 
     712 
     713def get_intercept(q, q_0, q_1): 
     714    """ 
     715    Returns the fraction of the side at which the 
     716    q-value intercept the pixel, None otherwise. 
     717    The values returned is the fraction ON THE SIDE 
     718    OF THE LOWEST Q. :: 
     719 
     720            A           B 
     721        +-----------+--------+    <--- pixel size 
     722        0                    1 
     723        Q_0 -------- Q ----- Q_1   <--- equivalent Q range 
     724        if Q_1 > Q_0, A is returned 
     725        if Q_1 < Q_0, B is returned 
     726        if Q is outside the range of [Q_0, Q_1], None is returned 
     727 
     728    """ 
     729    if q_1 > q_0: 
     730        if q > q_0 and q <= q_1: 
     731            return (q - q_0) / (q_1 - q_0) 
    255732    else: 
    256         new_err_data = data2d.err_data.flatten() 
    257     mask = np.ones(len(new_data), dtype=bool) 
    258  
    259     # TODO: make sense of the following two lines... 
    260     #from sas.sascalc.dataloader.data_info import Data2D 
    261     #output = Data2D() 
    262     output = data2d 
    263     output.data = new_data 
    264     output.err_data = new_err_data 
    265     output.qx_data = qx_data 
    266     output.qy_data = qy_data 
    267     output.q_data = q_data 
    268     output.mask = mask 
    269  
    270     return output 
    271  
    272 ################################################################################ 
    273  
    274 class Binning(object): 
    275     ''' 
    276     This class just creates a binning object 
    277     either linear or log 
    278     ''' 
    279  
    280     def __init__(self, min_value, max_value, n_bins, base=None): 
    281         ''' 
    282         if base is None: Linear binning 
    283         ''' 
    284         self.min = min_value if min_value > 0 else 0.0001 
    285         self.max = max_value 
    286         self.n_bins = n_bins 
    287         self.base = base 
    288  
    289     def get_bin_index(self, value): 
    290         ''' 
    291         The general formula logarithm binning is: 
    292         bin = floor(N * (log(x) - log(min)) / (log(max) - log(min))) 
    293         ''' 
    294         if self.base: 
    295             temp_x = self.n_bins * (math.log(value, self.base) - math.log(self.min, self.base)) 
    296             temp_y = math.log(self.max, self.base) - math.log(self.min, self.base) 
    297         else: 
    298             temp_x = self.n_bins * (value - self.min) 
    299             temp_y = self.max - self.min 
    300         # Bin index calulation 
    301         return int(math.floor(temp_x / temp_y)) 
    302  
    303  
    304 ################################################################################ 
    305  
    306 class _Slab(object): 
    307     """ 
    308     Compute average I(Q) for a region of interest 
    309     """ 
    310  
    311     def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, 
    312                  y_max=0.0, bin_width=0.001): 
    313         # Minimum Qx value [A-1] 
    314         self.x_min = x_min 
    315         # Maximum Qx value [A-1] 
    316         self.x_max = x_max 
    317         # Minimum Qy value [A-1] 
    318         self.y_min = y_min 
    319         # Maximum Qy value [A-1] 
    320         self.y_max = y_max 
    321         # Bin width (step size) [A-1] 
    322         self.bin_width = bin_width 
    323         # If True, I(|Q|) will be return, otherwise, 
    324         # negative q-values are allowed 
    325         self.fold = False 
    326  
    327     def __call__(self, data2D): 
    328         return NotImplemented 
    329  
    330     def _avg(self, data2D, maj): 
    331         """ 
    332         Compute average I(Q_maj) for a region of interest. 
    333         The major axis is defined as the axis of Q_maj. 
    334         The minor axis is the axis that we average over. 
    335  
    336         :param data2D: Data2D object 
    337         :param maj_min: min value on the major axis 
     733        if q > q_1 and q <= q_0: 
     734            return (q - q_1) / (q_0 - q_1) 
     735    return None 
     736 
     737 
     738class _Sector(object): 
     739    """ 
     740    Defines a sector region on a 2D data set. 
     741    The sector is defined by r_min, r_max, phi_min, phi_max, 
     742    and the position of the center of the ring 
     743    where phi_min and phi_max are defined by the right 
     744    and left lines wrt central line 
     745    and phi_max could be less than phi_min. 
     746 
     747    Phi is defined between 0 and 2*pi in anti-clockwise 
     748    starting from the x- axis on the left-hand side 
     749    """ 
     750    def __init__(self, r_min, r_max, phi_min=0, phi_max=2 * math.pi, nbins=20): 
     751        self.r_min = r_min 
     752        self.r_max = r_max 
     753        self.phi_min = phi_min 
     754        self.phi_max = phi_max 
     755        self.nbins = nbins 
     756 
     757    def _agv(self, data2D, run='phi'): 
     758        """ 
     759        Perform sector averaging. 
     760 
     761        :param data2D: Data2D object 
     762        :param run:  define the varying parameter ('phi' , 'q' , or 'q2') 
     763 
    338764        :return: Data1D object 
    339765        """ 
    340         if len(data2D.detector) > 1: 
    341             msg = "_Slab._avg: invalid number of " 
    342             msg += " detectors: %g" % len(data2D.detector) 
    343             raise RuntimeError(msg) 
    344  
    345         # Get data 
    346         data = data2D.data[np.isfinite(data2D.data)] 
    347         err_data = data2D.err_data[np.isfinite(data2D.data)] 
    348         qx_data = data2D.qx_data[np.isfinite(data2D.data)] 
    349         qy_data = data2D.qy_data[np.isfinite(data2D.data)] 
    350  
    351         # Build array of Q intervals 
    352         if maj == 'x': 
    353             if self.fold: 
    354                 x_min = 0 
    355             else: 
    356                 x_min = self.x_min 
    357             nbins = int(math.ceil((self.x_max - x_min) / self.bin_width)) 
    358         elif maj == 'y': 
    359             if self.fold: 
    360                 y_min = 0 
    361             else: 
    362                 y_min = self.y_min 
    363             nbins = int(math.ceil((self.y_max - y_min) / self.bin_width)) 
    364         else: 
    365             raise RuntimeError("_Slab._avg: unrecognized axis %s" % str(maj)) 
    366  
    367         x = np.zeros(nbins) 
    368         y = np.zeros(nbins) 
    369         err_y = np.zeros(nbins) 
    370         y_counts = np.zeros(nbins) 
    371  
    372         # Average pixelsize in q space 
    373         for npts in range(len(data)): 
    374             # default frac 
    375             frac_x = 0 
    376             frac_y = 0 
    377             # get ROI 
    378             if self.x_min <= qx_data[npts] and self.x_max > qx_data[npts]: 
    379                 frac_x = 1 
    380             if self.y_min <= qy_data[npts] and self.y_max > qy_data[npts]: 
    381                 frac_y = 1 
    382             frac = frac_x * frac_y 
    383  
     766        if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: 
     767            raise RuntimeError, "Ring averaging only take plottable_2D objects" 
     768        Pi = math.pi 
     769 
     770        # Get the all data & info 
     771        data = data2D.data[numpy.isfinite(data2D.data)] 
     772        q_data = data2D.q_data[numpy.isfinite(data2D.data)] 
     773        err_data = data2D.err_data[numpy.isfinite(data2D.data)] 
     774        qx_data = data2D.qx_data[numpy.isfinite(data2D.data)] 
     775        qy_data = data2D.qy_data[numpy.isfinite(data2D.data)] 
     776        dq_data = None 
     777 
     778        # Get the dq for resolution averaging 
     779        if data2D.dqx_data is not None and data2D.dqy_data is not None: 
     780            # The pinholes and det. pix contribution present 
     781            # in both direction of the 2D which must be subtracted when 
     782            # converting to 1D: dq_overlap should calculated ideally at 
     783            # q = 0. 
     784            # Extrapolate dqy(perp) at q = 0 
     785            z_max = max(data2D.q_data) 
     786            z_min = min(data2D.q_data) 
     787            x_max = data2D.dqx_data[data2D.q_data[z_max]] 
     788            x_min = data2D.dqx_data[data2D.q_data[z_min]] 
     789            y_max = data2D.dqy_data[data2D.q_data[z_max]] 
     790            y_min = data2D.dqy_data[data2D.q_data[z_min]] 
     791            # Find qdx at q = 0 
     792            dq_overlap_x = (x_min * z_max - x_max * z_min) / (z_max - z_min) 
     793            # when extrapolation goes wrong 
     794            if dq_overlap_x > min(data2D.dqx_data): 
     795                dq_overlap_x = min(data2D.dqx_data) 
     796            dq_overlap_x *= dq_overlap_x 
     797            # Find qdx at q = 0 
     798            dq_overlap_y = (y_min * z_max - y_max * z_min) / (z_max - z_min) 
     799            # when extrapolation goes wrong 
     800            if dq_overlap_y > min(data2D.dqy_data): 
     801                dq_overlap_y = min(data2D.dqy_data) 
     802            # get dq at q=0. 
     803            dq_overlap_y *= dq_overlap_y 
     804 
     805            dq_overlap = numpy.sqrt((dq_overlap_x + dq_overlap_y) / 2.0) 
     806            if dq_overlap < 0: 
     807                dq_overlap = y_min 
     808            dqx_data = data2D.dqx_data[numpy.isfinite(data2D.data)] 
     809            dqy_data = data2D.dqy_data[numpy.isfinite(data2D.data)] - dq_overlap 
     810            # def; dqx_data = dq_r dqy_data = dq_phi 
     811            # Convert dq 2D to 1D here 
     812            dqx = dqx_data * dqx_data 
     813            dqy = dqy_data * dqy_data 
     814            dq_data = numpy.add(dqx, dqy) 
     815            dq_data = numpy.sqrt(dq_data) 
     816 
     817        #set space for 1d outputs 
     818        x = numpy.zeros(self.nbins) 
     819        y = numpy.zeros(self.nbins) 
     820        y_err = numpy.zeros(self.nbins) 
     821        x_err = numpy.zeros(self.nbins) 
     822        y_counts = numpy.zeros(self.nbins) 
     823 
     824        # Get the min and max into the region: 0 <= phi < 2Pi 
     825        phi_min = flip_phi(self.phi_min) 
     826        phi_max = flip_phi(self.phi_max) 
     827 
     828        for n in range(len(data)): 
     829            frac = 0 
     830 
     831            # q-value at the pixel (j,i) 
     832            q_value = q_data[n] 
     833            data_n = data[n] 
     834 
     835            # Is pixel within range? 
     836            is_in = False 
     837 
     838            # phi-value of the pixel (j,i) 
     839            phi_value = math.atan2(qy_data[n], qx_data[n]) + Pi 
     840 
     841            ## No need to calculate the frac when all data are within range 
     842            if self.r_min <= q_value and q_value <= self.r_max: 
     843                frac = 1 
    384844            if frac == 0: 
    385845                continue 
    386             # binning: find axis of q 
    387             if maj == 'x': 
    388                 q_value = qx_data[npts] 
    389                 min_value = x_min 
    390             if maj == 'y': 
    391                 q_value = qy_data[npts] 
    392                 min_value = y_min 
    393             if self.fold and q_value < 0: 
    394                 q_value = -q_value 
    395             # bin 
    396             i_q = int(math.ceil((q_value - min_value) / self.bin_width)) - 1 
    397  
    398             # skip outside of max bins 
    399             if i_q < 0 or i_q >= nbins: 
     846            #In case of two ROIs (symmetric major and minor regions)(for 'q2') 
     847            if run.lower() == 'q2': 
     848                ## For minor sector wing 
     849                # Calculate the minor wing phis 
     850                phi_min_minor = flip_phi(phi_min - Pi) 
     851                phi_max_minor = flip_phi(phi_max - Pi) 
     852                # Check if phis of the minor ring is within 0 to 2pi 
     853                if phi_min_minor > phi_max_minor: 
     854                    is_in = (phi_value > phi_min_minor or \ 
     855                              phi_value < phi_max_minor) 
     856                else: 
     857                    is_in = (phi_value > phi_min_minor and \ 
     858                              phi_value < phi_max_minor) 
     859 
     860            #For all cases(i.e.,for 'q', 'q2', and 'phi') 
     861            #Find pixels within ROI 
     862            if phi_min > phi_max: 
     863                is_in = is_in or (phi_value > phi_min or \ 
     864                                   phi_value < phi_max) 
     865            else: 
     866                is_in = is_in or (phi_value >= phi_min  and \ 
     867                                    phi_value < phi_max) 
     868 
     869            if not is_in: 
     870                frac = 0 
     871            if frac == 0: 
    400872                continue 
    401  
    402             # TODO: find better definition of x[i_q] based on q_data 
    403             # min_value + (i_q + 1) * self.bin_width / 2.0 
    404             x[i_q] += frac * q_value 
    405             y[i_q] += frac * data[npts] 
    406  
    407             if err_data is None or err_data[npts] == 0.0: 
    408                 if data[npts] < 0: 
    409                     data[npts] = -data[npts] 
    410                 err_y[i_q] += frac * frac * data[npts] 
     873            # Check which type of averaging we need 
     874            if run.lower() == 'phi': 
     875                temp_x = (self.nbins) * (phi_value - self.phi_min) 
     876                temp_y = (self.phi_max - self.phi_min) 
     877                i_bin = int(math.floor(temp_x / temp_y)) 
    411878            else: 
    412                 err_y[i_q] += frac * frac * err_data[npts] * err_data[npts] 
    413             y_counts[i_q] += frac 
    414  
    415         # Average the sums 
    416         for n in range(nbins): 
    417             err_y[n] = math.sqrt(err_y[n]) 
    418  
    419         err_y = err_y / y_counts 
    420         y = y / y_counts 
    421         x = x / y_counts 
    422         idx = (np.isfinite(y) & np.isfinite(x)) 
    423  
     879                temp_x = (self.nbins) * (q_value - self.r_min) 
     880                temp_y = (self.r_max - self.r_min) 
     881                i_bin = int(math.floor(temp_x / temp_y)) 
     882 
     883            # Take care of the edge case at phi = 2pi. 
     884            if i_bin == self.nbins: 
     885                i_bin = self.nbins - 1 
     886 
     887            ## Get the total y 
     888            y[i_bin] += frac * data_n 
     889            x[i_bin] += frac * q_value 
     890            if err_data[n] is None or err_data[n] == 0.0: 
     891                if data_n < 0: 
     892                    data_n = -data_n 
     893                y_err[i_bin] += frac * frac * data_n 
     894            else: 
     895                y_err[i_bin] += frac * frac * err_data[n] * err_data[n] 
     896 
     897            if dq_data is not None: 
     898                # To be consistent with dq calculation in 1d reduction, 
     899                # we need just the averages (not quadratures) because 
     900                # it should not depend on the number of the q points 
     901                # in the qr bins. 
     902                x_err[i_bin] += frac * dq_data[n] 
     903            else: 
     904                x_err = None 
     905            y_counts[i_bin] += frac 
     906 
     907        # Organize the results 
     908        for i in range(self.nbins): 
     909            y[i] = y[i] / y_counts[i] 
     910            y_err[i] = math.sqrt(y_err[i]) / y_counts[i] 
     911 
     912            # The type of averaging: phi,q2, or q 
     913            # Calculate x[i]should be at the center of the bin 
     914            if run.lower() == 'phi': 
     915                x[i] = (self.phi_max - self.phi_min) / self.nbins * \ 
     916                    (1.0 * i + 0.5) + self.phi_min 
     917            else: 
     918                # We take the center of ring area, not radius. 
     919                # This is more accurate than taking the radial center of ring. 
     920                #delta_r = (self.r_max - self.r_min) / self.nbins 
     921                #r_inner = self.r_min + delta_r * i 
     922                #r_outer = r_inner + delta_r 
     923                #x[i] = math.sqrt((r_inner * r_inner + r_outer * r_outer) / 2) 
     924                x[i] = x[i] / y_counts[i] 
     925        y_err[y_err == 0] = numpy.average(y_err) 
     926        idx = (numpy.isfinite(y) & numpy.isfinite(y_err)) 
     927        if x_err is not None: 
     928            d_x = x_err[idx] / y_counts[idx] 
     929        else: 
     930            d_x = None 
    424931        if not idx.any(): 
    425             msg = "Average Error: No points inside ROI to average..." 
    426             raise ValueError(msg) 
    427         return Data1D(x=x[idx], y=y[idx], dy=err_y[idx]) 
    428  
    429  
    430 class SlabY(_Slab): 
    431     """ 
    432     Compute average I(Qy) for a region of interest 
    433     """ 
    434  
     932            msg = "Average Error: No points inside sector of ROI to average..." 
     933            raise ValueError, msg 
     934        #elif len(y[idx])!= self.nbins: 
     935        #    print "resulted",self.nbins- len(y[idx]), 
     936        #"empty bin(s) due to tight binning..." 
     937        return Data1D(x=x[idx], y=y[idx], dy=y_err[idx], dx=d_x) 
     938 
     939 
     940class SectorPhi(_Sector): 
     941    """ 
     942    Sector average as a function of phi. 
     943    I(phi) is return and the data is averaged over Q. 
     944 
     945    A sector is defined by r_min, r_max, phi_min, phi_max. 
     946    The number of bin in phi also has to be defined. 
     947    """ 
    435948    def __call__(self, data2D): 
    436949        """ 
    437         Compute average I(Qy) for a region of interest 
     950        Perform sector average and return I(phi). 
    438951 
    439952        :param data2D: Data2D object 
    440953        :return: Data1D object 
    441954        """ 
    442         return self._avg(data2D, 'y') 
    443  
    444  
    445 class SlabX(_Slab): 
    446     """ 
    447     Compute average I(Qx) for a region of interest 
    448     """ 
    449  
     955        return self._agv(data2D, 'phi') 
     956 
     957 
     958class SectorQ(_Sector): 
     959    """ 
     960    Sector average as a function of Q for both symatric wings. 
     961    I(Q) is return and the data is averaged over phi. 
     962 
     963    A sector is defined by r_min, r_max, phi_min, phi_max. 
     964    r_min, r_max, phi_min, phi_max >0. 
     965    The number of bin in Q also has to be defined. 
     966    """ 
    450967    def __call__(self, data2D): 
    451968        """ 
    452         Compute average I(Qx) for a region of interest 
    453         :param data2D: Data2D object 
     969        Perform sector average and return I(Q). 
     970 
     971        :param data2D: Data2D object 
     972 
    454973        :return: Data1D object 
    455974        """ 
    456         return self._avg(data2D, 'x') 
    457  
    458 ################################################################################ 
    459  
    460 class Boxsum(object): 
    461     """ 
    462     Perform the sum of counts in a 2D region of interest. 
    463     """ 
    464  
     975        return self._agv(data2D, 'q2') 
     976 
     977 
     978class Ringcut(object): 
     979    """ 
     980    Defines a ring on a 2D data set. 
     981    The ring is defined by r_min, r_max, and 
     982    the position of the center of the ring. 
     983 
     984    The data returned is the region inside the ring 
     985 
     986    Phi_min and phi_max should be defined between 0 and 2*pi 
     987    in anti-clockwise starting from the x- axis on the left-hand side 
     988    """ 
     989    def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0): 
     990        # Minimum radius 
     991        self.r_min = r_min 
     992        # Maximum radius 
     993        self.r_max = r_max 
     994        # Center of the ring in x 
     995        self.center_x = center_x 
     996        # Center of the ring in y 
     997        self.center_y = center_y 
     998 
     999    def __call__(self, data2D): 
     1000        """ 
     1001        Apply the ring to the data set. 
     1002        Returns the angular distribution for a given q range 
     1003 
     1004        :param data2D: Data2D object 
     1005 
     1006        :return: index array in the range 
     1007        """ 
     1008        if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: 
     1009            raise RuntimeError, "Ring cut only take plottable_2D objects" 
     1010 
     1011        # Get data 
     1012        qx_data = data2D.qx_data 
     1013        qy_data = data2D.qy_data 
     1014        q_data = numpy.sqrt(qx_data * qx_data + qy_data * qy_data) 
     1015 
     1016        # check whether or not the data point is inside ROI 
     1017        out = (self.r_min <= q_data) & (self.r_max >= q_data) 
     1018        return out 
     1019 
     1020 
     1021class Boxcut(object): 
     1022    """ 
     1023    Find a rectangular 2D region of interest. 
     1024    """ 
    4651025    def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0): 
    4661026        # Minimum Qx value [A-1] 
     
    4751035    def __call__(self, data2D): 
    4761036        """ 
    477         Perform the sum in the region of interest 
    478  
    479         :param data2D: Data2D object 
    480         :return: number of counts, error on number of counts, 
    481             number of points summed 
    482         """ 
    483         y, err_y, y_counts = self._sum(data2D) 
    484  
    485         # Average the sums 
    486         counts = 0 if y_counts == 0 else y 
    487         error = 0 if y_counts == 0 else math.sqrt(err_y) 
    488  
    489         # Added y_counts to return, SMK & PDB, 04/03/2013 
    490         return counts, error, y_counts 
    491  
    492     def _sum(self, data2D): 
    493         """ 
    494         Perform the sum in the region of interest 
    495  
    496         :param data2D: Data2D object 
    497         :return: number of counts, 
    498             error on number of counts, number of entries summed 
    499         """ 
    500         if len(data2D.detector) > 1: 
    501             msg = "Circular averaging: invalid number " 
    502             msg += "of detectors: %g" % len(data2D.detector) 
    503             raise RuntimeError(msg) 
    504         # Get data 
    505         data = data2D.data[np.isfinite(data2D.data)] 
    506         err_data = data2D.err_data[np.isfinite(data2D.data)] 
    507         qx_data = data2D.qx_data[np.isfinite(data2D.data)] 
    508         qy_data = data2D.qy_data[np.isfinite(data2D.data)] 
    509  
    510         y = 0.0 
    511         err_y = 0.0 
    512         y_counts = 0.0 
    513  
    514         # Average pixelsize in q space 
    515         for npts in range(len(data)): 
    516             # default frac 
    517             frac_x = 0 
    518             frac_y = 0 
    519  
    520             # get min and max at each points 
    521             qx = qx_data[npts] 
    522             qy = qy_data[npts] 
    523  
    524             # get the ROI 
    525             if self.x_min <= qx and self.x_max > qx: 
    526                 frac_x = 1 
    527             if self.y_min <= qy and self.y_max > qy: 
    528                 frac_y = 1 
    529             # Find the fraction along each directions 
    530             frac = frac_x * frac_y 
    531             if frac == 0: 
    532                 continue 
    533             y += frac * data[npts] 
    534             if err_data is None or err_data[npts] == 0.0: 
    535                 if data[npts] < 0: 
    536                     data[npts] = -data[npts] 
    537                 err_y += frac * frac * data[npts] 
    538             else: 
    539                 err_y += frac * frac * err_data[npts] * err_data[npts] 
    540             y_counts += frac 
    541         return y, err_y, y_counts 
    542  
    543  
    544 class Boxavg(Boxsum): 
    545     """ 
    546     Perform the average of counts in a 2D region of interest. 
    547     """ 
    548  
    549     def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0): 
    550         super(Boxavg, self).__init__(x_min=x_min, x_max=x_max, 
    551                                      y_min=y_min, y_max=y_max) 
    552  
    553     def __call__(self, data2D): 
    554         """ 
    555         Perform the sum in the region of interest 
    556  
    557         :param data2D: Data2D object 
    558         :return: average counts, error on average counts 
    559  
    560         """ 
    561         y, err_y, y_counts = self._sum(data2D) 
    562  
    563         # Average the sums 
    564         counts = 0 if y_counts == 0 else y / y_counts 
    565         error = 0 if y_counts == 0 else math.sqrt(err_y) / y_counts 
    566  
    567         return counts, error 
    568  
    569 ################################################################################ 
    570  
    571 class CircularAverage(object): 
    572     """ 
    573     Perform circular averaging on 2D data 
    574  
    575     The data returned is the distribution of counts 
    576     as a function of Q 
    577     """ 
    578  
    579     def __init__(self, r_min=0.0, r_max=0.0, bin_width=0.0005): 
    580         # Minimum radius included in the average [A-1] 
    581         self.r_min = r_min 
    582         # Maximum radius included in the average [A-1] 
    583         self.r_max = r_max 
    584         # Bin width (step size) [A-1] 
    585         self.bin_width = bin_width 
    586  
    587     def __call__(self, data2D, ismask=False): 
    588         """ 
    589         Perform circular averaging on the data 
    590  
    591         :param data2D: Data2D object 
    592         :return: Data1D object 
    593         """ 
    594         # Get data W/ finite values 
    595         data = data2D.data[np.isfinite(data2D.data)] 
    596         q_data = data2D.q_data[np.isfinite(data2D.data)] 
    597         err_data = data2D.err_data[np.isfinite(data2D.data)] 
    598         mask_data = data2D.mask[np.isfinite(data2D.data)] 
    599  
    600         dq_data = None 
    601         if data2D.dqx_data is not None and data2D.dqy_data is not None: 
    602             dq_data = get_dq_data(data2D) 
    603  
    604         if len(q_data) == 0: 
    605             msg = "Circular averaging: invalid q_data: %g" % data2D.q_data 
    606             raise RuntimeError(msg) 
    607  
    608         # Build array of Q intervals 
    609         nbins = int(math.ceil((self.r_max - self.r_min) / self.bin_width)) 
    610  
    611         x = np.zeros(nbins) 
    612         y = np.zeros(nbins) 
    613         err_y = np.zeros(nbins) 
    614         err_x = np.zeros(nbins) 
    615         y_counts = np.zeros(nbins) 
    616  
    617         for npt in range(len(data)): 
    618  
    619             if ismask and not mask_data[npt]: 
    620                 continue 
    621  
    622             frac = 0 
    623  
    624             # q-value at the pixel (j,i) 
    625             q_value = q_data[npt] 
    626             data_n = data[npt] 
    627  
    628             # No need to calculate the frac when all data are within range 
    629             if self.r_min >= self.r_max: 
    630                 raise ValueError("Limit Error: min > max") 
    631  
    632             if self.r_min <= q_value and q_value <= self.r_max: 
    633                 frac = 1 
    634             if frac == 0: 
    635                 continue 
    636             i_q = int(math.floor((q_value - self.r_min) / self.bin_width)) 
    637  
    638             # Take care of the edge case at phi = 2pi. 
    639             if i_q == nbins: 
    640                 i_q = nbins - 1 
    641             y[i_q] += frac * data_n 
    642             # Take dqs from data to get the q_average 
    643             x[i_q] += frac * q_value 
    644             if err_data is None or err_data[npt] == 0.0: 
    645                 if data_n < 0: 
    646                     data_n = -data_n 
    647                 err_y[i_q] += frac * frac * data_n 
    648             else: 
    649                 err_y[i_q] += frac * frac * err_data[npt] * err_data[npt] 
    650             if dq_data is not None: 
    651                 # To be consistent with dq calculation in 1d reduction, 
    652                 # we need just the averages (not quadratures) because 
    653                 # it should not depend on the number of the q points 
    654                 # in the qr bins. 
    655                 err_x[i_q] += frac * dq_data[npt] 
    656             else: 
    657                 err_x = None 
    658             y_counts[i_q] += frac 
    659  
    660         # Average the sums 
    661         for n in range(nbins): 
    662             if err_y[n] < 0: 
    663                 err_y[n] = -err_y[n] 
    664             err_y[n] = math.sqrt(err_y[n]) 
    665             # if err_x is not None: 
    666             #    err_x[n] = math.sqrt(err_x[n]) 
    667  
    668         err_y = err_y / y_counts 
    669         err_y[err_y == 0] = np.average(err_y) 
    670         y = y / y_counts 
    671         x = x / y_counts 
    672         idx = (np.isfinite(y)) & (np.isfinite(x)) 
    673  
    674         if err_x is not None: 
    675             d_x = err_x[idx] / y_counts[idx] 
    676         else: 
    677             d_x = None 
    678  
    679         if not idx.any(): 
    680             msg = "Average Error: No points inside ROI to average..." 
    681             raise ValueError(msg) 
    682  
    683         return Data1D(x=x[idx], y=y[idx], dy=err_y[idx], dx=d_x) 
    684  
    685 ################################################################################ 
    686  
    687 class Ring(object): 
    688     """ 
    689     Defines a ring on a 2D data set. 
    690     The ring is defined by r_min, r_max, and 
    691     the position of the center of the ring. 
    692  
    693     The data returned is the distribution of counts 
    694     around the ring as a function of phi. 
    695  
    696     Phi_min and phi_max should be defined between 0 and 2*pi 
    697     in anti-clockwise starting from the x- axis on the left-hand side 
    698     """ 
    699     # Todo: remove center. 
    700  
    701     def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0, nbins=36): 
    702         # Minimum radius 
    703         self.r_min = r_min 
    704         # Maximum radius 
    705         self.r_max = r_max 
    706         # Center of the ring in x 
    707         self.center_x = center_x 
    708         # Center of the ring in y 
    709         self.center_y = center_y 
    710         # Number of angular bins 
    711         self.nbins_phi = nbins 
    712  
    713     def __call__(self, data2D): 
    714         """ 
    715         Apply the ring to the data set. 
    716         Returns the angular distribution for a given q range 
    717  
    718         :param data2D: Data2D object 
    719  
    720         :return: Data1D object 
    721         """ 
    722         if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: 
    723             raise RuntimeError("Ring averaging only take plottable_2D objects") 
    724  
    725         Pi = math.pi 
    726  
    727         # Get data 
    728         data = data2D.data[np.isfinite(data2D.data)] 
    729         q_data = data2D.q_data[np.isfinite(data2D.data)] 
    730         err_data = data2D.err_data[np.isfinite(data2D.data)] 
    731         qx_data = data2D.qx_data[np.isfinite(data2D.data)] 
    732         qy_data = data2D.qy_data[np.isfinite(data2D.data)] 
    733  
    734         # Set space for 1d outputs 
    735         phi_bins = np.zeros(self.nbins_phi) 
    736         phi_counts = np.zeros(self.nbins_phi) 
    737         phi_values = np.zeros(self.nbins_phi) 
    738         phi_err = np.zeros(self.nbins_phi) 
    739  
    740         # Shift to apply to calculated phi values in order 
    741         # to center first bin at zero 
    742         phi_shift = Pi / self.nbins_phi 
    743  
    744         for npt in range(len(data)): 
    745             frac = 0 
    746             # q-value at the point (npt) 
    747             q_value = q_data[npt] 
    748             data_n = data[npt] 
    749  
    750             # phi-value at the point (npt) 
    751             phi_value = math.atan2(qy_data[npt], qx_data[npt]) + Pi 
    752  
    753             if self.r_min <= q_value and q_value <= self.r_max: 
    754                 frac = 1 
    755             if frac == 0: 
    756                 continue 
    757             # binning 
    758             i_phi = int(math.floor((self.nbins_phi) * 
    759                                    (phi_value + phi_shift) / (2 * Pi))) 
    760  
    761             # Take care of the edge case at phi = 2pi. 
    762             if i_phi >= self.nbins_phi: 
    763                 i_phi = 0 
    764             phi_bins[i_phi] += frac * data[npt] 
    765  
    766             if err_data is None or err_data[npt] == 0.0: 
    767                 if data_n < 0: 
    768                     data_n = -data_n 
    769                 phi_err[i_phi] += frac * frac * math.fabs(data_n) 
    770             else: 
    771                 phi_err[i_phi] += frac * frac * err_data[npt] * err_data[npt] 
    772             phi_counts[i_phi] += frac 
    773  
    774         for i in range(self.nbins_phi): 
    775             phi_bins[i] = phi_bins[i] / phi_counts[i] 
    776             phi_err[i] = math.sqrt(phi_err[i]) / phi_counts[i] 
    777             phi_values[i] = 2.0 * math.pi / self.nbins_phi * (1.0 * i) 
    778  
    779         idx = (np.isfinite(phi_bins)) 
    780  
    781         if not idx.any(): 
    782             msg = "Average Error: No points inside ROI to average..." 
    783             raise ValueError(msg) 
    784         # elif len(phi_bins[idx])!= self.nbins_phi: 
    785         #    print "resulted",self.nbins_phi- len(phi_bins[idx]) 
    786         #,"empty bin(s) due to tight binning..." 
    787         return Data1D(x=phi_values[idx], y=phi_bins[idx], dy=phi_err[idx]) 
    788  
    789  
    790 class _Sector(object): 
    791     """ 
    792     Defines a sector region on a 2D data set. 
    793     The sector is defined by r_min, r_max, phi_min, phi_max, 
    794     and the position of the center of the ring 
    795     where phi_min and phi_max are defined by the right 
    796     and left lines wrt central line 
    797     and phi_max could be less than phi_min. 
    798  
    799     Phi is defined between 0 and 2*pi in anti-clockwise 
    800     starting from the x- axis on the left-hand side 
    801     """ 
    802  
    803     def __init__(self, r_min, r_max, phi_min=0, phi_max=2 * math.pi, nbins=20, 
    804                  base=None): 
    805         ''' 
    806         :param base: must be a valid base for an algorithm, i.e., 
    807         a positive number 
    808         ''' 
    809         self.r_min = r_min 
    810         self.r_max = r_max 
    811         self.phi_min = phi_min 
    812         self.phi_max = phi_max 
    813         self.nbins = nbins 
    814         self.base = base 
    815  
    816     def _agv(self, data2D, run='phi'): 
    817         """ 
    818         Perform sector averaging. 
    819  
    820         :param data2D: Data2D object 
    821         :param run:  define the varying parameter ('phi' , 'q' , or 'q2') 
    822  
    823         :return: Data1D object 
    824         """ 
    825         if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: 
    826             raise RuntimeError("Ring averaging only take plottable_2D objects") 
    827  
    828         # Get the all data & info 
    829         data = data2D.data[np.isfinite(data2D.data)] 
    830         q_data = data2D.q_data[np.isfinite(data2D.data)] 
    831         err_data = data2D.err_data[np.isfinite(data2D.data)] 
    832         qx_data = data2D.qx_data[np.isfinite(data2D.data)] 
    833         qy_data = data2D.qy_data[np.isfinite(data2D.data)] 
    834  
    835         dq_data = None 
    836         if data2D.dqx_data is not None and data2D.dqy_data is not None: 
    837             dq_data = get_dq_data(data2D) 
    838  
    839         # set space for 1d outputs 
    840         x = np.zeros(self.nbins) 
    841         y = np.zeros(self.nbins) 
    842         y_err = np.zeros(self.nbins) 
    843         x_err = np.zeros(self.nbins) 
    844         y_counts = np.zeros(self.nbins)  # Cycle counts (for the mean) 
    845  
    846         # Get the min and max into the region: 0 <= phi < 2Pi 
    847         phi_min = flip_phi(self.phi_min) 
    848         phi_max = flip_phi(self.phi_max) 
    849  
    850         #  binning object 
    851         if run.lower() == 'phi': 
    852             binning = Binning(self.phi_min, self.phi_max, self.nbins, self.base) 
    853         else: 
    854             binning = Binning(self.r_min, self.r_max, self.nbins, self.base) 
    855  
    856         for n in range(len(data)): 
    857  
    858             # q-value at the pixel (j,i) 
    859             q_value = q_data[n] 
    860             data_n = data[n] 
    861  
    862             # Is pixel within range? 
    863             is_in = False 
    864  
    865             # phi-value of the pixel (j,i) 
    866             phi_value = math.atan2(qy_data[n], qx_data[n]) + math.pi 
    867  
    868             # No need to calculate: data outside of the radius 
    869             if self.r_min > q_value or q_value > self.r_max: 
    870                 continue 
    871  
    872             # In case of two ROIs (symmetric major and minor regions)(for 'q2') 
    873             if run.lower() == 'q2': 
    874                 # For minor sector wing 
    875                 # Calculate the minor wing phis 
    876                 phi_min_minor = flip_phi(phi_min - math.pi) 
    877                 phi_max_minor = flip_phi(phi_max - math.pi) 
    878                 # Check if phis of the minor ring is within 0 to 2pi 
    879                 if phi_min_minor > phi_max_minor: 
    880                     is_in = (phi_value > phi_min_minor or 
    881                              phi_value < phi_max_minor) 
    882                 else: 
    883                     is_in = (phi_value > phi_min_minor and 
    884                              phi_value < phi_max_minor) 
    885  
    886             # For all cases(i.e.,for 'q', 'q2', and 'phi') 
    887             # Find pixels within ROI 
    888             if phi_min > phi_max: 
    889                 is_in = is_in or (phi_value > phi_min or 
    890                                   phi_value < phi_max) 
    891             else: 
    892                 is_in = is_in or (phi_value >= phi_min and 
    893                                   phi_value < phi_max) 
    894  
    895             # data oustide of the phi range 
    896             if not is_in: 
    897                 continue 
    898  
    899             # Get the binning index 
    900             if run.lower() == 'phi': 
    901                 i_bin = binning.get_bin_index(phi_value) 
    902             else: 
    903                 i_bin = binning.get_bin_index(q_value) 
    904  
    905             # Take care of the edge case at phi = 2pi. 
    906             if i_bin == self.nbins: 
    907                 i_bin = self.nbins - 1 
    908  
    909             # Get the total y 
    910             y[i_bin] += data_n 
    911             x[i_bin] += q_value 
    912             if err_data[n] is None or err_data[n] == 0.0: 
    913                 if data_n < 0: 
    914                     data_n = -data_n 
    915                 y_err[i_bin] += data_n 
    916             else: 
    917                 y_err[i_bin] += err_data[n]**2 
    918  
    919             if dq_data is not None: 
    920                 # To be consistent with dq calculation in 1d reduction, 
    921                 # we need just the averages (not quadratures) because 
    922                 # it should not depend on the number of the q points 
    923                 # in the qr bins. 
    924                 x_err[i_bin] += dq_data[n] 
    925             else: 
    926                 x_err = None 
    927             y_counts[i_bin] += 1 
    928  
    929         # Organize the results 
    930         for i in range(self.nbins): 
    931             y[i] = y[i] / y_counts[i] 
    932             y_err[i] = math.sqrt(y_err[i]) / y_counts[i] 
    933  
    934             # The type of averaging: phi,q2, or q 
    935             # Calculate x[i]should be at the center of the bin 
    936             if run.lower() == 'phi': 
    937                 x[i] = (self.phi_max - self.phi_min) / self.nbins * \ 
    938                     (1.0 * i + 0.5) + self.phi_min 
    939             else: 
    940                 # We take the center of ring area, not radius. 
    941                 # This is more accurate than taking the radial center of ring. 
    942                 # delta_r = (self.r_max - self.r_min) / self.nbins 
    943                 # r_inner = self.r_min + delta_r * i 
    944                 # r_outer = r_inner + delta_r 
    945                 # x[i] = math.sqrt((r_inner * r_inner + r_outer * r_outer) / 2) 
    946                 x[i] = x[i] / y_counts[i] 
    947         y_err[y_err == 0] = np.average(y_err) 
    948         idx = (np.isfinite(y) & np.isfinite(y_err)) 
    949         if x_err is not None: 
    950             d_x = x_err[idx] / y_counts[idx] 
    951         else: 
    952             d_x = None 
    953         if not idx.any(): 
    954             msg = "Average Error: No points inside sector of ROI to average..." 
    955             raise ValueError(msg) 
    956         # elif len(y[idx])!= self.nbins: 
    957         #    print "resulted",self.nbins- len(y[idx]), 
    958         # "empty bin(s) due to tight binning..." 
    959         return Data1D(x=x[idx], y=y[idx], dy=y_err[idx], dx=d_x) 
    960  
    961  
    962 class SectorPhi(_Sector): 
    963     """ 
    964     Sector average as a function of phi. 
    965     I(phi) is return and the data is averaged over Q. 
    966  
    967     A sector is defined by r_min, r_max, phi_min, phi_max. 
    968     The number of bin in phi also has to be defined. 
    969     """ 
    970  
    971     def __call__(self, data2D): 
    972         """ 
    973         Perform sector average and return I(phi). 
    974  
    975         :param data2D: Data2D object 
    976         :return: Data1D object 
    977         """ 
    978         return self._agv(data2D, 'phi') 
    979  
    980  
    981 class SectorQ(_Sector): 
    982     """ 
    983     Sector average as a function of Q for both symatric wings. 
    984     I(Q) is return and the data is averaged over phi. 
    985  
    986     A sector is defined by r_min, r_max, phi_min, phi_max. 
    987     r_min, r_max, phi_min, phi_max >0. 
    988     The number of bin in Q also has to be defined. 
    989     """ 
    990  
    991     def __call__(self, data2D): 
    992         """ 
    993         Perform sector average and return I(Q). 
    994  
    995         :param data2D: Data2D object 
    996  
    997         :return: Data1D object 
    998         """ 
    999         return self._agv(data2D, 'q2') 
    1000  
    1001 ################################################################################ 
    1002  
    1003 class Ringcut(object): 
    1004     """ 
    1005     Defines a ring on a 2D data set. 
    1006     The ring is defined by r_min, r_max, and 
    1007     the position of the center of the ring. 
    1008  
    1009     The data returned is the region inside the ring 
    1010  
    1011     Phi_min and phi_max should be defined between 0 and 2*pi 
    1012     in anti-clockwise starting from the x- axis on the left-hand side 
    1013     """ 
    1014  
    1015     def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0): 
    1016         # Minimum radius 
    1017         self.r_min = r_min 
    1018         # Maximum radius 
    1019         self.r_max = r_max 
    1020         # Center of the ring in x 
    1021         self.center_x = center_x 
    1022         # Center of the ring in y 
    1023         self.center_y = center_y 
    1024  
    1025     def __call__(self, data2D): 
    1026         """ 
    1027         Apply the ring to the data set. 
    1028         Returns the angular distribution for a given q range 
    1029  
    1030         :param data2D: Data2D object 
    1031  
    1032         :return: index array in the range 
    1033         """ 
    1034         if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: 
    1035             raise RuntimeError("Ring cut only take plottable_2D objects") 
    1036  
    1037         # Get data 
    1038         qx_data = data2D.qx_data 
    1039         qy_data = data2D.qy_data 
    1040         q_data = np.sqrt(qx_data * qx_data + qy_data * qy_data) 
    1041  
    1042         # check whether or not the data point is inside ROI 
    1043         out = (self.r_min <= q_data) & (self.r_max >= q_data) 
    1044         return out 
    1045  
    1046 ################################################################################ 
    1047  
    1048 class Boxcut(object): 
    1049     """ 
    1050     Find a rectangular 2D region of interest. 
    1051     """ 
    1052  
    1053     def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0): 
    1054         # Minimum Qx value [A-1] 
    1055         self.x_min = x_min 
    1056         # Maximum Qx value [A-1] 
    1057         self.x_max = x_max 
    1058         # Minimum Qy value [A-1] 
    1059         self.y_min = y_min 
    1060         # Maximum Qy value [A-1] 
    1061         self.y_max = y_max 
    1062  
    1063     def __call__(self, data2D): 
    1064         """ 
    10651037       Find a rectangular 2D region of interest. 
    10661038 
     
    10831055        """ 
    10841056        if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: 
    1085             raise RuntimeError("Boxcut take only plottable_2D objects") 
     1057            raise RuntimeError, "Boxcut take only plottable_2D objects" 
    10861058        # Get qx_ and qy_data 
    10871059        qx_data = data2D.qx_data 
     
    10941066        return outx & outy 
    10951067 
    1096 ################################################################################ 
    10971068 
    10981069class Sectorcut(object): 
     
    11061077    and (phi_max-phi_min) should not be larger than pi 
    11071078    """ 
    1108  
    11091079    def __init__(self, phi_min=0, phi_max=math.pi): 
    11101080        self.phi_min = phi_min 
     
    11361106        """ 
    11371107        if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]: 
    1138             raise RuntimeError("Sectorcut take only plottable_2D objects") 
     1108            raise RuntimeError, "Sectorcut take only plottable_2D objects" 
    11391109        Pi = math.pi 
    11401110        # Get data 
     
    11431113 
    11441114        # get phi from data 
    1145         phi_data = np.arctan2(qy_data, qx_data) 
     1115        phi_data = numpy.arctan2(qy_data, qx_data) 
    11461116 
    11471117        # Get the min and max into the region: -pi <= phi < Pi 
     
    11501120        # check for major sector 
    11511121        if phi_min_major > phi_max_major: 
    1152             out_major = (phi_min_major <= phi_data) + \ 
    1153                 (phi_max_major > phi_data) 
     1122            out_major = (phi_min_major <= phi_data) + (phi_max_major > phi_data) 
    11541123        else: 
    1155             out_major = (phi_min_major <= phi_data) & ( 
    1156                 phi_max_major > phi_data) 
     1124            out_major = (phi_min_major <= phi_data) & (phi_max_major > phi_data) 
    11571125 
    11581126        # minor sector 
     
    11641132        if phi_min_minor > phi_max_minor: 
    11651133            out_minor = (phi_min_minor <= phi_data) + \ 
    1166                 (phi_max_minor >= phi_data) 
     1134                            (phi_max_minor >= phi_data) 
    11671135        else: 
    11681136            out_minor = (phi_min_minor <= phi_data) & \ 
    1169                 (phi_max_minor >= phi_data) 
     1137                            (phi_max_minor >= phi_data) 
    11701138        out = out_major + out_minor 
    11711139 
  • src/sas/sascalc/dataloader/readers/sesans_reader.py

    r149b8f6 r9a5097c  
    11""" 
    22    SESANS reader (based on ASCII reader) 
    3  
     3     
    44    Reader for .ses or .sesans file format 
    5  
    6     Jurrian Bakker 
     5     
     6    Jurrian Bakker  
    77""" 
    88import numpy as np 
     
    1818_ZERO = 1e-16 
    1919 
    20  
    2120class Reader: 
    2221    """ 
    2322    Class to load sesans files (6 columns). 
    2423    """ 
    25     # File type 
     24    ## File type 
    2625    type_name = "SESANS" 
    27  
    28     # Wildcards 
     26     
     27    ## Wildcards 
    2928    type = ["SESANS files (*.ses)|*.ses", 
    3029            "SESANS files (*..sesans)|*.sesans"] 
    31     # List of allowed extensions 
     30    ## List of allowed extensions 
    3231    ext = ['.ses', '.SES', '.sesans', '.SESANS'] 
    33  
    34     # Flag to bypass extension check 
     32     
     33    ## Flag to bypass extension check 
    3534    allow_all = True 
    36  
     35     
    3736    def read(self, path): 
     37         
     38#        print "reader triggered" 
     39         
    3840        """ 
    3941        Load data file 
    40  
     42         
    4143        :param path: file path 
    42  
     44         
    4345        :return: SESANSData1D object, or None 
    44  
     46         
    4547        :raise RuntimeError: when the file can't be opened 
    4648        :raise ValueError: when the length of the data vectors are inconsistent 
     
    4951            basename = os.path.basename(path) 
    5052            _, extension = os.path.splitext(basename) 
    51             if not (self.allow_all or extension.lower() in self.ext): 
    52                 raise RuntimeError( 
    53                     "{} has an unrecognized file extension".format(path)) 
     53            if self.allow_all or extension.lower() in self.ext: 
     54                try: 
     55                    # Read in binary mode since GRASP frequently has no-ascii 
     56                    # characters that brakes the open operation 
     57                    input_f = open(path,'rb') 
     58                except: 
     59                    raise  RuntimeError, "sesans_reader: cannot open %s" % path 
     60                buff = input_f.read() 
     61                lines = buff.splitlines() 
     62                x  = np.zeros(0) 
     63                y  = np.zeros(0) 
     64                dy = np.zeros(0) 
     65                lam  = np.zeros(0) 
     66                dlam = np.zeros(0) 
     67                dx = np.zeros(0) 
     68                 
     69               #temp. space to sort data 
     70                tx  = np.zeros(0) 
     71                ty  = np.zeros(0) 
     72                tdy = np.zeros(0) 
     73                tlam  = np.zeros(0) 
     74                tdlam = np.zeros(0) 
     75                tdx = np.zeros(0) 
     76                output = Data1D(x=x, y=y, lam=lam, dy=dy, dx=dx, dlam=dlam, isSesans=True) 
     77                self.filename = output.filename = basename 
     78 
     79                paramnames=[] 
     80                paramvals=[] 
     81                zvals=[] 
     82                dzvals=[] 
     83                lamvals=[] 
     84                dlamvals=[] 
     85                Pvals=[] 
     86                dPvals=[] 
     87 
     88                for line in lines: 
     89                    # Initial try for CSV (split on ,) 
     90                    line=line.strip() 
     91                    toks = line.split('\t') 
     92                    if len(toks)==2: 
     93                        paramnames.append(toks[0]) 
     94                        paramvals.append(toks[1]) 
     95                    if len(toks)>5: 
     96                        zvals.append(toks[0]) 
     97                        dzvals.append(toks[3]) 
     98                        lamvals.append(toks[4]) 
     99                        dlamvals.append(toks[5]) 
     100                        Pvals.append(toks[1]) 
     101                        dPvals.append(toks[2]) 
     102                    else: 
     103                        continue 
     104 
     105                x=[] 
     106                y=[] 
     107                lam=[] 
     108                dx=[] 
     109                dy=[] 
     110                dlam=[] 
     111                lam_header = lamvals[0].split() 
     112                data_conv_z = None 
     113                default_z_unit = "A" 
     114                data_conv_P = None 
     115                default_p_unit = " " # Adjust unit for axis (L^-3) 
     116                lam_unit = lam_header[1].replace("[","").replace("]","") 
     117                if lam_unit == 'AA': 
     118                    lam_unit = 'A' 
     119                varheader=[zvals[0],dzvals[0],lamvals[0],dlamvals[0],Pvals[0],dPvals[0]] 
     120                valrange=range(1, len(zvals)) 
     121                for i in valrange: 
     122                    x.append(float(zvals[i])) 
     123                    y.append(float(Pvals[i])) 
     124                    lam.append(float(lamvals[i])) 
     125                    dy.append(float(dPvals[i])) 
     126                    dx.append(float(dzvals[i])) 
     127                    dlam.append(float(dlamvals[i])) 
     128 
     129                x,y,lam,dy,dx,dlam = [ 
     130                    np.asarray(v, 'double') 
     131                   for v in (x,y,lam,dy,dx,dlam) 
     132                ] 
     133 
     134                input_f.close() 
     135 
     136                output.x, output.x_unit = self._unit_conversion(x, lam_unit, default_z_unit) 
     137                output.y = y 
     138                output.y_unit = r'\AA^{-2} cm^{-1}'  # output y_unit added 
     139                output.dx, output.dx_unit = self._unit_conversion(dx, lam_unit, default_z_unit) 
     140                output.dy = dy 
     141                output.lam, output.lam_unit = self._unit_conversion(lam, lam_unit, default_z_unit) 
     142                output.dlam, output.dlam_unit = self._unit_conversion(dlam, lam_unit, default_z_unit) 
     143                 
     144                output.xaxis(r"\rm{z}", output.x_unit) 
     145                output.yaxis(r"\rm{ln(P)/(t \lambda^2)}", output.y_unit)  # Adjust label to ln P/(lam^2 t), remove lam column refs 
     146 
     147                # Store loading process information 
     148                output.meta_data['loader'] = self.type_name 
     149                #output.sample.thickness = float(paramvals[6]) 
     150                output.sample.name = paramvals[1] 
     151                output.sample.ID = paramvals[0] 
     152                zaccept_unit_split = paramnames[7].split("[") 
     153                zaccept_unit = zaccept_unit_split[1].replace("]","") 
     154                if zaccept_unit.strip() == r'\AA^-1' or zaccept_unit.strip() == r'\A^-1': 
     155                    zaccept_unit = "1/A" 
     156                output.sample.zacceptance=(float(paramvals[7]),zaccept_unit) 
     157                output.vars = varheader 
     158 
     159                if len(output.x) < 1: 
     160                    raise RuntimeError, "%s is empty" % path 
     161                return output 
     162 
    54163        else: 
    55             raise RuntimeError("{} is not a file".format(path)) 
    56         with open(path, 'r') as input_f: 
    57             line = input_f.readline() 
    58             params = {} 
    59             while not line.startswith("BEGIN_DATA"): 
    60                 terms = line.split() 
    61                 if len(terms) >= 2: 
    62                     params[terms[0]] = " ".join(terms[1:]) 
    63                 line = input_f.readline() 
    64             self.params = params 
     164            raise RuntimeError, "%s is not a file" % path 
     165        return None 
    65166 
    66             if "FileFormatVersion" not in self.params: 
    67                 raise RuntimeError("SES file missing FileFormatVersion") 
    68             if float(self.params["FileFormatVersion"]) >= 2.0: 
    69                 raise RuntimeError("SASView only supports SES version 1") 
    70  
    71             if "SpinEchoLength_unit" not in self.params: 
    72                 raise RuntimeError("SpinEchoLength has no units") 
    73             if "Wavelength_unit" not in self.params: 
    74                 raise RuntimeError("Wavelength has no units") 
    75             if params["SpinEchoLength_unit"] != params["Wavelength_unit"]: 
    76                 raise RuntimeError("The spin echo data has rudely used " 
    77                                    "different units for the spin echo length " 
    78                                    "and the wavelength.  While sasview could " 
    79                                    "handle this instance, it is a violation " 
    80                                    "of the file format and will not be " 
    81                                    "handled by other software.") 
    82  
    83             headers = input_f.readline().split() 
    84  
    85             self._insist_header(headers, "SpinEchoLength") 
    86             self._insist_header(headers, "Depolarisation") 
    87             self._insist_header(headers, "Depolarisation_error") 
    88             self._insist_header(headers, "Wavelength") 
    89  
    90             data = np.loadtxt(input_f) 
    91  
    92             if data.shape[1] != len(headers): 
    93                 raise RuntimeError( 
    94                     "File has {} headers, but {} columns".format( 
    95                         len(headers), 
    96                         data.shape[1])) 
    97  
    98             if not data.size: 
    99                 raise RuntimeError("{} is empty".format(path)) 
    100             x = data[:, headers.index("SpinEchoLength")] 
    101             if "SpinEchoLength_error" in headers: 
    102                 dx = data[:, headers.index("SpinEchoLength_error")] 
    103             else: 
    104                 dx = x * 0.05 
    105             lam = data[:, headers.index("Wavelength")] 
    106             if "Wavelength_error" in headers: 
    107                 dlam = data[:, headers.index("Wavelength_error")] 
    108             else: 
    109                 dlam = lam * 0.05 
    110             y = data[:, headers.index("Depolarisation")] 
    111             dy = data[:, headers.index("Depolarisation_error")] 
    112  
    113             lam_unit = self._unit_fetch("Wavelength") 
    114             x, x_unit = self._unit_conversion(x, "A", 
    115                                               self._unit_fetch( 
    116                                                   "SpinEchoLength")) 
    117             dx, dx_unit = self._unit_conversion( 
    118                 dx, lam_unit, 
    119                 self._unit_fetch("SpinEchoLength")) 
    120             dlam, dlam_unit = self._unit_conversion( 
    121                 dlam, lam_unit, 
    122                 self._unit_fetch("Wavelength")) 
    123             y_unit = self._unit_fetch("Depolarisation") 
    124  
    125             output = Data1D(x=x, y=y, lam=lam, dy=dy, dx=dx, dlam=dlam, 
    126                             isSesans=True) 
    127  
    128             output.y_unit = y_unit 
    129             output.x_unit = x_unit 
    130             output.source.wavelength_unit = lam_unit 
    131             output.source.wavelength = lam 
    132             self.filename = output.filename = basename 
    133             output.xaxis(r"\rm{z}", x_unit) 
    134             # Adjust label to ln P/(lam^2 t), remove lam column refs 
    135             output.yaxis(r"\rm{ln(P)/(t \lambda^2)}", y_unit) 
    136             # Store loading process information 
    137             output.meta_data['loader'] = self.type_name 
    138             output.sample.name = params["Sample"] 
    139             output.sample.ID = params["DataFileTitle"] 
    140             output.sample.thickness = self._unit_conversion( 
    141                 float(params["Thickness"]), "cm", 
    142                 self._unit_fetch("Thickness"))[0] 
    143  
    144             output.sample.zacceptance = ( 
    145                 float(params["Theta_zmax"]), 
    146                 self._unit_fetch("Theta_zmax")) 
    147  
    148             output.sample.yacceptance = ( 
    149                 float(params["Theta_ymax"]), 
    150                 self._unit_fetch("Theta_ymax")) 
    151             return output 
    152  
    153     @staticmethod 
    154     def _insist_header(headers, name): 
    155         if name not in headers: 
    156             raise RuntimeError( 
    157                 "Missing {} column in spin echo data".format(name)) 
    158  
    159     @staticmethod 
    160     def _unit_conversion(value, value_unit, default_unit): 
    161         """ 
    162         Performs unit conversion on a measurement. 
    163  
    164         :param value: The magnitude of the measurement 
    165         :param value_unit: a string containing the final desired unit 
    166         :param default_unit: string with the units of the original measurement 
    167         :return: The magnitude of the measurement in the new units 
    168         """ 
    169         # (float, string, string) -> float 
    170         if has_converter and value_unit != default_unit: 
    171             data_conv_q = Converter(default_unit) 
    172             value = data_conv_q(value, units=value_unit) 
     167    def _unit_conversion(self, value, value_unit, default_unit): 
     168        if has_converter == True and value_unit != default_unit: 
     169            data_conv_q = Converter(value_unit) 
     170            value = data_conv_q(value, units=default_unit) 
    173171            new_unit = default_unit 
    174172        else: 
    175173            new_unit = value_unit 
    176174        return value, new_unit 
    177  
    178     def _unit_fetch(self, unit): 
    179         return self.params[unit+"_unit"] 
  • src/sas/sasgui/guiframe/local_perspectives/plotting/Plotter2D.py

    r3e5648b r7432acb  
    361361                if self.slicer.__class__.__name__ != "BoxSum": 
    362362                    wx_id = ids.next() 
    363                     name = '&Edit Slicer Parameters and Batch Slicing' 
    364                     slicerpop.Append(wx_id, name) 
     363                    slicerpop.Append(wx_id, '&Edit Slicer Parameters') 
    365364                    wx.EVT_MENU(self, wx_id, self._onEditSlicer) 
    366365            slicerpop.AppendSeparator() 
     
    533532 
    534533        """ 
    535         # Clear current slicer 
     534        ## Clear current slicer 
    536535        if self.slicer is not None: 
    537536            self.slicer.clear() 
    538         # Create a new slicer 
     537        ## Create a new slicer 
    539538        self.slicer_z += 1 
    540539        self.slicer = slicer(self, self.subplot, zorder=self.slicer_z) 
    541540        self.subplot.set_ylim(self.data2D.ymin, self.data2D.ymax) 
    542541        self.subplot.set_xlim(self.data2D.xmin, self.data2D.xmax) 
    543         # Draw slicer 
     542        ## Draw slicer 
    544543        self.update() 
    545544        self.slicer.update() 
     
    573572        npt = math.floor(npt) 
    574573        from sas.sascalc.dataloader.manipulations import CircularAverage 
    575         # compute the maximum radius of data2D 
     574        ## compute the maximum radius of data2D 
    576575        self.qmax = max(math.fabs(self.data2D.xmax), 
    577576                        math.fabs(self.data2D.xmin)) 
     
    579578                        math.fabs(self.data2D.ymin)) 
    580579        self.radius = math.sqrt(math.pow(self.qmax, 2) + math.pow(self.ymax, 2)) 
    581         # Compute beam width 
     580        ##Compute beam width 
    582581        bin_width = (self.qmax + self.qmax) / npt 
    583         # Create data1D circular average of data2D 
     582        ## Create data1D circular average of data2D 
    584583        Circle = CircularAverage(r_min=0, r_max=self.radius, 
    585584                                 bin_width=bin_width) 
     
    600599        new_plot.name = "Circ avg " + self.data2D.name 
    601600        new_plot.source = self.data2D.source 
    602         # new_plot.info = self.data2D.info 
     601        #new_plot.info = self.data2D.info 
    603602        new_plot.interactive = True 
    604603        new_plot.detector = self.data2D.detector 
    605604 
    606         # If the data file does not tell us what the axes are, just assume... 
     605        ## If the data file does not tell us what the axes are, just assume... 
    607606        new_plot.xaxis("\\rm{Q}", "A^{-1}") 
    608607        if hasattr(self.data2D, "scale") and \ 
     
    616615        new_plot.id = "Circ avg " + self.data2D.name 
    617616        new_plot.is_data = True 
    618         self.parent.update_theory(data_id=self.data2D.id, theory=new_plot) 
     617        self.parent.update_theory(data_id=self.data2D.id, \ 
     618                                       theory=new_plot) 
    619619        wx.PostEvent(self.parent, 
    620620                     NewPlotEvent(plot=new_plot, title=new_plot.name)) 
     
    630630        """ 
    631631        if self.slicer is not None: 
    632             from parameters_panel_slicer import SlicerParameterPanel 
     632            from SlicerParameters import SlicerParameterPanel 
    633633            dialog = SlicerParameterPanel(self, -1, "Slicer Parameters") 
    634634            dialog.set_slicer(self.slicer.__class__.__name__, 
     
    668668        params = self.slicer.get_params() 
    669669        ## Create a new panel to display results of summation of Data2D 
    670         from parameters_panel_boxsum import SlicerPanel 
     670        from slicerpanel import SlicerPanel 
    671671        win = MDIFrame(self.parent, None, 'None', (100, 200)) 
    672672        new_panel = SlicerPanel(parent=win, id=-1, 
     
    758758        if default_name.count('.') > 0: 
    759759            default_name = default_name.split('.')[0] 
     760        #default_name += "_out" 
    760761        if self.parent is not None: 
    761762            self.parent.show_data2d(data, default_name) 
    762763 
    763764    def modifyGraphAppearance(self, e): 
    764         self.graphApp = graphAppearance(self, 'Modify graph appearance', 
    765                                         legend=False) 
     765        self.graphApp = graphAppearance(self, 'Modify graph appearance', legend=False) 
    766766        self.graphApp.setDefaults(self.grid_on, self.legend_on, 
    767767                                  self.xaxis_label, self.yaxis_label, 
  • src/sas/sasgui/guiframe/local_perspectives/plotting/SectorSlicer.py

    r8de66b6 r7432acb  
    3737        ## Absolute value of the Angle between the middle line and any side line 
    3838        self.phi = math.pi / 12 
    39         # Binning base for log/lin binning 
    40         self.bin_base = 0 
    4139        ## Middle line 
    4240        self.main_line = LineInteractor(self, self.base.subplot, color='blue', 
     
    153151        phimin = -self.left_line.phi + self.main_line.theta 
    154152        phimax = self.left_line.phi + self.main_line.theta 
    155         bin_base = self.bin_base 
    156153        if nbins is None: 
    157154            nbins = 20 
    158155        sect = SectorQ(r_min=0.0, r_max=radius, 
    159156                       phi_min=phimin + math.pi, 
    160                        phi_max=phimax + math.pi, nbins=nbins, base=bin_base) 
     157                       phi_max=phimax + math.pi, nbins=nbins) 
    161158 
    162159        sector = sect(self.base.data2D) 
     
    242239        params["Delta_Phi [deg]"] = math.fabs(self.left_line.phi * 180 / math.pi) 
    243240        params["nbins"] = self.nbins 
    244         params["binning base"] = self.bin_base 
    245241        return params 
    246242 
     
    256252        phi = math.fabs(params["Delta_Phi [deg]"] * math.pi / 180) 
    257253        self.nbins = int(params["nbins"]) 
    258         self.bin_base = params["binning base"] 
    259254        self.main_line.theta = main 
    260255        ## Reset the slicer parameters 
  • src/sas/sasgui/perspectives/fitting/basepage.py

    ra1b8fee ra1b8fee  
    254254        if not hasattr(self, "model_view"): 
    255255            return 
    256         toggle_mode_on = self.model_view.IsEnabled() or self.data is None 
     256        toggle_mode_on = self.model_view.IsEnabled() 
    257257        if toggle_mode_on: 
    258258            if self.enable2D and not check_data_validity(self.data): 
  • src/sas/sasgui/perspectives/fitting/fitting.py

    ra534432 ra1b8fee  
    383383          '(Re)Load all models present in user plugin_models folder') 
    384384        wx.EVT_MENU(owner, wx_id, self.load_plugin_models) 
    385  
     385    
    386386    def set_edit_menu_helper(self, owner=None, menu=None): 
    387387        """ 
     
    17341734            @param unsmeared_error: data error, rescaled to unsmeared model 
    17351735        """ 
    1736         number_finite = np.count_nonzero(np.isfinite(y)) 
    1737         np.nan_to_num(y) 
    1738         new_plot = self.create_theory_1D(x, y, page_id, model, data, state, 
    1739                                          data_description=model.name, 
    1740                                          data_id=str(page_id) + " " + data.name) 
    1741         if unsmeared_model is not None: 
    1742             self.create_theory_1D(x, unsmeared_model, page_id, model, data, state, 
    1743                                   data_description=model.name + " unsmeared", 
    1744                                   data_id=str(page_id) + " " + data.name + " unsmeared") 
    1745  
    1746             if unsmeared_data is not None and unsmeared_error is not None: 
    1747                 self.create_theory_1D(x, unsmeared_data, page_id, model, data, state, 
    1748                                       data_description="Data unsmeared", 
    1749                                       data_id="Data  " + data.name + " unsmeared", 
    1750                                       dy=unsmeared_error) 
    1751         # Comment this out until we can get P*S models with correctly populated parameters 
    1752         if sq_model is not None and pq_model is not None: 
    1753             self.create_theory_1D(x, sq_model, page_id, model, data, state, 
    1754                                   data_description=model.name + " S(q)", 
    1755                                   data_id=str(page_id) + " " + data.name + " S(q)") 
    1756             self.create_theory_1D(x, pq_model, page_id, model, data, state, 
    1757                                   data_description=model.name + " P(q)", 
    1758                                   data_id=str(page_id) + " " + data.name + " P(q)") 
    1759  
    1760         current_pg = self.fit_panel.get_page_by_id(page_id) 
    1761         title = new_plot.title 
    1762         batch_on = self.fit_panel.get_page_by_id(page_id).batch_on 
    1763         if not batch_on: 
    1764             wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=str(title))) 
    1765         elif plot_result: 
    1766             top_data_id = self.fit_panel.get_page_by_id(page_id).data.id 
    1767             if data.id == top_data_id: 
    1768                 wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=str(title))) 
    1769         caption = current_pg.window_caption 
    1770         self.page_finder[page_id].set_fit_tab_caption(caption=caption) 
    1771  
    1772         self.page_finder[page_id].set_theory_data(data=new_plot, 
     1736        try: 
     1737            np.nan_to_num(y) 
     1738            new_plot = self.create_theory_1D(x, y, page_id, model, data, state, 
     1739                                             data_description=model.name, 
     1740                                             data_id=str(page_id) + " " + data.name) 
     1741            if unsmeared_model is not None: 
     1742                self.create_theory_1D(x, unsmeared_model, page_id, model, data, state, 
     1743                                      data_description=model.name + " unsmeared", 
     1744                                      data_id=str(page_id) + " " + data.name + " unsmeared") 
     1745 
     1746                if unsmeared_data is not None and unsmeared_error is not None: 
     1747                    self.create_theory_1D(x, unsmeared_data, page_id, model, data, state, 
     1748                                          data_description="Data unsmeared", 
     1749                                          data_id="Data  " + data.name + " unsmeared", 
     1750                                          dy=unsmeared_error) 
     1751            if sq_model is not None and pq_model is not None: 
     1752                self.create_theory_1D(x, sq_model, page_id, model, data, state, 
     1753                                      data_description=model.name + " S(q)", 
     1754                                      data_id=str(page_id) + " " + data.name + " S(q)") 
     1755                self.create_theory_1D(x, pq_model, page_id, model, data, state, 
     1756                                      data_description=model.name + " P(q)", 
     1757                                      data_id=str(page_id) + " " + data.name + " P(q)") 
     1758 
     1759            current_pg = self.fit_panel.get_page_by_id(page_id) 
     1760            title = new_plot.title 
     1761            batch_on = self.fit_panel.get_page_by_id(page_id).batch_on 
     1762            if not batch_on: 
     1763                wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, 
     1764                                            title=str(title))) 
     1765            elif plot_result: 
     1766                top_data_id = self.fit_panel.get_page_by_id(page_id).data.id 
     1767                if data.id == top_data_id: 
     1768                    wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, 
     1769                                            title=str(title))) 
     1770            caption = current_pg.window_caption 
     1771            self.page_finder[page_id].set_fit_tab_caption(caption=caption) 
     1772 
     1773            self.page_finder[page_id].set_theory_data(data=new_plot, 
    17731774                                                      fid=data.id) 
    1774         if toggle_mode_on: 
    1775             wx.PostEvent(self.parent, 
    1776                          NewPlotEvent(group_id=str(page_id) + " Model2D", 
     1775            if toggle_mode_on: 
     1776                wx.PostEvent(self.parent, 
     1777                             NewPlotEvent(group_id=str(page_id) + " Model2D", 
    17771778                                          action="Hide")) 
    1778         else: 
    1779             if update_chisqr: 
    1780                 wx.PostEvent(current_pg, 
    1781                              Chi2UpdateEvent(output=self._cal_chisqr( 
     1779            else: 
     1780                if update_chisqr: 
     1781                    wx.PostEvent(current_pg, 
     1782                                 Chi2UpdateEvent(output=self._cal_chisqr( 
    17821783                                                                data=data, 
    17831784                                                                fid=fid, 
    17841785                                                                weight=weight, 
    1785                                                                 page_id=page_id, 
    1786                                                                 index=index))) 
    1787             else: 
    1788                 self._plot_residuals(page_id=page_id, data=data, fid=fid, 
    1789                                      index=index, weight=weight) 
    1790  
    1791         if not number_finite: 
    1792             logger.error("Using the present parameters the model does not return any finite value. ") 
    1793             msg = "Computing Error: Model did not return any finite value." 
    1794             wx.PostEvent(self.parent, StatusEvent(status = msg, info="error")) 
    1795         else: 
     1786                                                            page_id=page_id, 
     1787                                                            index=index))) 
     1788                else: 
     1789                    self._plot_residuals(page_id=page_id, data=data, fid=fid, 
     1790                                         index=index, weight=weight) 
     1791 
    17961792            msg = "Computation  completed!" 
    1797             if number_finite != y.size: 
    1798                 msg += ' PROBLEM: For some Q values the model returns non finite intensities!' 
    1799                 logger.error("For some Q values the model returns non finite intensities.") 
    18001793            wx.PostEvent(self.parent, StatusEvent(status=msg, type="stop")) 
     1794        except: 
     1795            raise 
    18011796 
    18021797    def _calc_exception(self, etype, value, tb): 
     
    18231818        that can be plot. 
    18241819        """ 
    1825         number_finite = np.count_nonzero(np.isfinite(image)) 
    18261820        np.nan_to_num(image) 
    18271821        new_plot = Data2D(image=image, err_image=data.err_data) 
     
    18821876                self._plot_residuals(page_id=page_id, data=data, fid=fid, 
    18831877                                      index=index, weight=weight) 
    1884  
    1885         if not number_finite: 
    1886             logger.error("Using the present parameters the model does not return any finite value. ") 
    1887             msg = "Computing Error: Model did not return any finite value." 
    1888             wx.PostEvent(self.parent, StatusEvent(status = msg, info="error")) 
    1889         else: 
    1890             msg = "Computation  completed!" 
    1891             if number_finite != image.size: 
    1892                 msg += ' PROBLEM: For some Qx,Qy values the model returns non finite intensities!' 
    1893                 logger.error("For some Qx,Qy values the model returns non finite intensities.") 
    1894             wx.PostEvent(self.parent, StatusEvent(status=msg, type="stop")) 
     1878        msg = "Computation  completed!" 
     1879        wx.PostEvent(self.parent, StatusEvent(status=msg, type="stop")) 
    18951880 
    18961881    def _draw_model2D(self, model, page_id, qmin, 
  • src/sas/sasgui/perspectives/fitting/models.py

    r56a282c ra1b8fee  
    156156    try: 
    157157        import compileall 
    158         compileall.compile_dir(dir=dir, ddir=dir, force=0, 
     158        compileall.compile_dir(dir=dir, ddir=dir, force=1, 
    159159                               quiet=report_problem) 
    160160    except: 
     
    163163 
    164164 
    165 def _find_models(): 
     165def _findModels(dir): 
    166166    """ 
    167167    Find custom models 
    168168    """ 
    169169    # List of plugin objects 
    170     directory = find_plugins_dir() 
     170    dir = find_plugins_dir() 
    171171    # Go through files in plug-in directory 
    172     if not os.path.isdir(directory): 
    173         msg = "SasView couldn't locate Model plugin folder %r." % directory 
     172    if not os.path.isdir(dir): 
     173        msg = "SasView couldn't locate Model plugin folder %r." % dir 
    174174        logger.warning(msg) 
    175175        return {} 
    176176 
    177     plugin_log("looking for models in: %s" % str(directory)) 
    178     # compile_file(directory)  #always recompile the folder plugin 
    179     logger.info("plugin model dir: %s" % str(directory)) 
     177    plugin_log("looking for models in: %s" % str(dir)) 
     178    #compile_file(dir)  #always recompile the folder plugin 
     179    logger.info("plugin model dir: %s" % str(dir)) 
    180180 
    181181    plugins = {} 
    182     for filename in os.listdir(directory): 
     182    for filename in os.listdir(dir): 
    183183        name, ext = os.path.splitext(filename) 
    184184        if ext == '.py' and not name == '__init__': 
    185             path = os.path.abspath(os.path.join(directory, filename)) 
     185            path = os.path.abspath(os.path.join(dir, filename)) 
    186186            try: 
    187187                model = load_custom_model(path) 
     
    193193                plugin_log(msg) 
    194194                logger.warning("Failed to load plugin %r. See %s for details" 
    195                                % (path, PLUGIN_LOG)) 
    196  
     195                                % (path, PLUGIN_LOG)) 
     196             
    197197    return plugins 
    198198 
     
    264264        temp = {} 
    265265        if self.is_changed(): 
    266             return  _find_models() 
     266            return  _findModels(dir) 
    267267        logger.info("plugin model : %s" % str(temp)) 
    268268        return temp 
     
    339339        """ 
    340340        self.plugins = [] 
    341         new_plugins = _find_models() 
     341        new_plugins = _findModels(dir) 
    342342        for name, plug in  new_plugins.iteritems(): 
    343343            for stored_name, stored_plug in self.stored_plugins.iteritems(): 
  • test/sasdataloader/test/utest_averaging.py

    re123eb9 r9a5097c  
    11 
     2import unittest 
    23import math 
    3 import os 
    4 import unittest 
    5 from scipy.stats import binned_statistic_2d 
     4 
     5from sas.sascalc.dataloader.loader import  Loader 
     6from sas.sascalc.dataloader.manipulations import Ring, CircularAverage, SectorPhi, get_q,reader2D_converter 
     7 
    68import numpy as np 
    7  
    89import sas.sascalc.dataloader.data_info as data_info 
    9 from sas.sascalc.dataloader.loader import Loader 
    10 from sas.sascalc.dataloader.manipulations import (Boxavg, Boxsum, 
    11                                                   CircularAverage, Ring, 
    12                                                   SectorPhi, SectorQ, SlabX, 
    13                                                   SlabY, get_q, 
    14                                                   reader2D_converter) 
    15  
    1610 
    1711class Averaging(unittest.TestCase): 
     
    1913        Test averaging manipulations on a flat distribution 
    2014    """ 
    21  
    2215    def setUp(self): 
    2316        """ 
     
    2518            should return the predefined height of the distribution (1.0). 
    2619        """ 
    27         x_0 = np.ones([100, 100]) 
    28         dx_0 = np.ones([100, 100]) 
    29  
     20        x_0  = np.ones([100,100]) 
     21        dx_0 = np.ones([100,100]) 
     22         
    3023        self.data = data_info.Data2D(data=x_0, err_data=dx_0) 
    3124        detector = data_info.Detector() 
    32         detector.distance = 1000.0  # mm 
    33         detector.pixel_size.x = 1.0  # mm 
    34         detector.pixel_size.y = 1.0  # mm 
    35  
     25        detector.distance = 1000.0  #mm 
     26        detector.pixel_size.x = 1.0 #mm 
     27        detector.pixel_size.y = 1.0 #mm 
     28         
    3629        # center in pixel position = (len(x_0)-1)/2 
    37         detector.beam_center.x = (len(x_0) - 1) / 2  # pixel number 
    38         detector.beam_center.y = (len(x_0) - 1) / 2  # pixel number 
     30        detector.beam_center.x = (len(x_0)-1)/2 #pixel number 
     31        detector.beam_center.y = (len(x_0)-1)/2 #pixel number 
    3932        self.data.detector.append(detector) 
    40  
     33         
    4134        source = data_info.Source() 
    42         source.wavelength = 10.0  # A 
     35        source.wavelength = 10.0 #A 
    4336        self.data.source = source 
    44  
    45         # get_q(dx, dy, det_dist, wavelength) where units are mm,mm,mm,and A 
    46         # respectively. 
     37         
     38        # get_q(dx, dy, det_dist, wavelength) where units are mm,mm,mm,and A respectively. 
    4739        self.qmin = get_q(1.0, 1.0, detector.distance, source.wavelength) 
    4840 
    4941        self.qmax = get_q(49.5, 49.5, detector.distance, source.wavelength) 
    50  
     42         
    5143        self.qstep = len(x_0) 
    52         x = np.linspace(start=-1 * self.qmax, 
    53                         stop=self.qmax, 
    54                         num=self.qstep, 
    55                         endpoint=True) 
    56         y = np.linspace(start=-1 * self.qmax, 
    57                         stop=self.qmax, 
    58                         num=self.qstep, 
    59                         endpoint=True) 
    60         self.data.x_bins = x 
    61         self.data.y_bins = y 
     44        x=  np.linspace(start= -1*self.qmax, 
     45                               stop= self.qmax, 
     46                               num= self.qstep, 
     47                               endpoint=True )   
     48        y = np.linspace(start= -1*self.qmax, 
     49                               stop= self.qmax, 
     50                               num= self.qstep, 
     51                               endpoint=True ) 
     52        self.data.x_bins=x 
     53        self.data.y_bins=y 
    6254        self.data = reader2D_converter(self.data) 
    63  
     55             
    6456    def test_ring_flat_distribution(self): 
    6557        """ 
    6658            Test ring averaging 
    6759        """ 
    68         r = Ring(r_min=2 * self.qmin, r_max=5 * self.qmin, 
    69                  center_x=self.data.detector[0].beam_center.x, 
     60        r = Ring(r_min=2*self.qmin, r_max=5*self.qmin,  
     61                 center_x=self.data.detector[0].beam_center.x,  
    7062                 center_y=self.data.detector[0].beam_center.y) 
    7163        r.nbins_phi = 20 
    72  
     64         
    7365        o = r(self.data) 
    7466        for i in range(20): 
    7567            self.assertEqual(o.y[i], 1.0) 
    76  
     68             
    7769    def test_sectorphi_full(self): 
    7870        """ 
    7971            Test sector averaging 
    8072        """ 
    81         r = SectorPhi(r_min=self.qmin, r_max=3 * self.qmin, 
    82                       phi_min=0, phi_max=math.pi * 2.0) 
     73        r = SectorPhi(r_min=self.qmin, r_max=3*self.qmin,  
     74                      phi_min=0, phi_max=math.pi*2.0) 
    8375        r.nbins_phi = 20 
    8476        o = r(self.data) 
    8577        for i in range(7): 
    8678            self.assertEqual(o.y[i], 1.0) 
    87  
     79             
     80             
    8881    def test_sectorphi_partial(self): 
    8982        """ 
    9083        """ 
    9184        phi_max = math.pi * 1.5 
    92         r = SectorPhi(r_min=self.qmin, r_max=3 * self.qmin, 
     85        r = SectorPhi(r_min=self.qmin, r_max=3*self.qmin,  
    9386                      phi_min=0, phi_max=phi_max) 
    9487        self.assertEqual(r.phi_max, phi_max) 
     
    9891        for i in range(17): 
    9992            self.assertEqual(o.y[i], 1.0) 
    100  
    101  
    102 class DataInfoTests(unittest.TestCase): 
    103  
     93             
     94             
     95 
     96class data_info_tests(unittest.TestCase): 
     97     
    10498    def setUp(self): 
    105         filepath = os.path.join(os.path.dirname( 
    106             os.path.realpath(__file__)), 'MAR07232_rest.ASC') 
    107         self.data = Loader().load(filepath) 
    108  
     99        self.data = Loader().load('MAR07232_rest.ASC') 
     100         
    109101    def test_ring(self): 
    110102        """ 
    111103            Test ring averaging 
    112104        """ 
    113         r = Ring(r_min=.005, r_max=.01, 
    114                  center_x=self.data.detector[0].beam_center.x, 
     105        r = Ring(r_min=.005, r_max=.01,  
     106                 center_x=self.data.detector[0].beam_center.x,  
    115107                 center_y=self.data.detector[0].beam_center.y, 
    116                  nbins=20) 
     108                 nbins = 20) 
    117109        ##r.nbins_phi = 20 
    118  
    119         o = r(self.data) 
    120         filepath = os.path.join(os.path.dirname( 
    121             os.path.realpath(__file__)), 'ring_testdata.txt') 
    122         answer = Loader().load(filepath) 
    123  
     110         
     111        o = r(self.data) 
     112        answer = Loader().load('ring_testdata.txt') 
     113         
    124114        for i in range(r.nbins_phi - 1): 
    125115            self.assertAlmostEqual(o.x[i + 1], answer.x[i], 4) 
    126116            self.assertAlmostEqual(o.y[i + 1], answer.y[i], 4) 
    127117            self.assertAlmostEqual(o.dy[i + 1], answer.dy[i], 4) 
    128  
     118             
    129119    def test_circularavg(self): 
    130120        """ 
    131         Test circular averaging 
    132         The test data was not generated by IGOR. 
    133         """ 
    134         r = CircularAverage(r_min=.00, r_max=.025, 
    135                             bin_width=0.0003) 
    136         r.nbins_phi = 20 
    137  
    138         o = r(self.data) 
    139  
    140         filepath = os.path.join(os.path.dirname( 
    141             os.path.realpath(__file__)), 'avg_testdata.txt') 
    142         answer = Loader().load(filepath) 
     121            Test circular averaging 
     122            The test data was not generated by IGOR. 
     123        """ 
     124        r = CircularAverage(r_min=.00, r_max=.025,  
     125                 bin_width=0.0003) 
     126        r.nbins_phi = 20 
     127         
     128        o = r(self.data) 
     129 
     130        answer = Loader().load('avg_testdata.txt') 
    143131        for i in range(r.nbins_phi): 
    144132            self.assertAlmostEqual(o.x[i], answer.x[i], 4) 
    145133            self.assertAlmostEqual(o.y[i], answer.y[i], 4) 
    146134            self.assertAlmostEqual(o.dy[i], answer.dy[i], 4) 
    147  
     135             
    148136    def test_box(self): 
    149137        """ 
     
    151139            The test data was not generated by IGOR. 
    152140        """ 
    153  
     141        from sas.sascalc.dataloader.manipulations import Boxsum, Boxavg 
     142         
    154143        r = Boxsum(x_min=.01, x_max=.015, y_min=0.01, y_max=0.015) 
    155144        s, ds, npoints = r(self.data) 
    156145        self.assertAlmostEqual(s, 34.278990899999997, 4) 
    157146        self.assertAlmostEqual(ds, 7.8007981835194293, 4) 
    158         self.assertAlmostEqual(npoints, 324.0000, 4) 
    159  
     147        self.assertAlmostEqual(npoints, 324.0000, 4)         
     148     
    160149        r = Boxavg(x_min=.01, x_max=.015, y_min=0.01, y_max=0.015) 
    161150        s, ds = r(self.data) 
    162151        self.assertAlmostEqual(s, 0.10579935462962962, 4) 
    163152        self.assertAlmostEqual(ds, 0.024076537603455028, 4) 
    164  
     153             
    165154    def test_slabX(self): 
    166155        """ 
     
    168157            The test data was not generated by IGOR. 
    169158        """ 
    170  
    171         r = SlabX(x_min=-.01, x_max=.01, y_min=-0.0002, 
    172                   y_max=0.0002, bin_width=0.0004) 
     159        from sas.sascalc.dataloader.manipulations import SlabX 
     160         
     161        r = SlabX(x_min=-.01, x_max=.01, y_min=-0.0002, y_max=0.0002, bin_width=0.0004) 
    173162        r.fold = False 
    174163        o = r(self.data) 
    175164 
    176         filepath = os.path.join(os.path.dirname( 
    177             os.path.realpath(__file__)), 'slabx_testdata.txt') 
    178         answer = Loader().load(filepath) 
    179         for i in range(len(o.x)): 
    180             self.assertAlmostEqual(o.x[i], answer.x[i], 4) 
    181             self.assertAlmostEqual(o.y[i], answer.y[i], 4) 
    182             self.assertAlmostEqual(o.dy[i], answer.dy[i], 4) 
    183  
     165        answer = Loader().load('slabx_testdata.txt') 
     166        for i in range(len(o.x)): 
     167            self.assertAlmostEqual(o.x[i], answer.x[i], 4) 
     168            self.assertAlmostEqual(o.y[i], answer.y[i], 4) 
     169            self.assertAlmostEqual(o.dy[i], answer.dy[i], 4) 
     170             
    184171    def test_slabY(self): 
    185172        """ 
     
    187174            The test data was not generated by IGOR. 
    188175        """ 
    189  
    190         r = SlabY(x_min=.005, x_max=.01, y_min=- 
    191                   0.01, y_max=0.01, bin_width=0.0004) 
     176        from sas.sascalc.dataloader.manipulations import SlabY 
     177         
     178        r = SlabY(x_min=.005, x_max=.01, y_min=-0.01, y_max=0.01, bin_width=0.0004) 
    192179        r.fold = False 
    193180        o = r(self.data) 
    194181 
    195         filepath = os.path.join(os.path.dirname( 
    196             os.path.realpath(__file__)), 'slaby_testdata.txt') 
    197         answer = Loader().load(filepath) 
    198         for i in range(len(o.x)): 
    199             self.assertAlmostEqual(o.x[i], answer.x[i], 4) 
    200             self.assertAlmostEqual(o.y[i], answer.y[i], 4) 
    201             self.assertAlmostEqual(o.dy[i], answer.dy[i], 4) 
    202  
     182        answer = Loader().load('slaby_testdata.txt') 
     183        for i in range(len(o.x)): 
     184            self.assertAlmostEqual(o.x[i], answer.x[i], 4) 
     185            self.assertAlmostEqual(o.y[i], answer.y[i], 4) 
     186            self.assertAlmostEqual(o.dy[i], answer.dy[i], 4) 
     187             
    203188    def test_sectorphi_full(self): 
    204189        """ 
     
    208193            The test data was not generated by IGOR. 
    209194        """ 
    210  
     195        from sas.sascalc.dataloader.manipulations import SectorPhi 
     196        import math 
     197         
    211198        nbins = 19 
    212199        phi_min = math.pi / (nbins + 1) 
    213200        phi_max = math.pi * 2 - phi_min 
    214  
     201         
    215202        r = SectorPhi(r_min=.005, 
    216203                      r_max=.01, 
     
    220207        o = r(self.data) 
    221208 
    222         filepath = os.path.join(os.path.dirname( 
    223             os.path.realpath(__file__)), 'ring_testdata.txt') 
    224         answer = Loader().load(filepath) 
    225         for i in range(len(o.x)): 
    226             self.assertAlmostEqual(o.x[i], answer.x[i], 4) 
    227             self.assertAlmostEqual(o.y[i], answer.y[i], 4) 
    228             self.assertAlmostEqual(o.dy[i], answer.dy[i], 4) 
    229  
     209        answer = Loader().load('ring_testdata.txt') 
     210        for i in range(len(o.x)): 
     211            self.assertAlmostEqual(o.x[i], answer.x[i], 4) 
     212            self.assertAlmostEqual(o.y[i], answer.y[i], 4) 
     213            self.assertAlmostEqual(o.dy[i], answer.dy[i], 4) 
     214             
    230215    def test_sectorphi_quarter(self): 
    231216        """ 
     
    233218            The test data was not generated by IGOR. 
    234219        """ 
    235  
    236         r = SectorPhi(r_min=.005, r_max=.01, phi_min=0, phi_max=math.pi / 2.0) 
    237         r.nbins_phi = 20 
    238         o = r(self.data) 
    239  
    240         filepath = os.path.join(os.path.dirname( 
    241             os.path.realpath(__file__)), 'sectorphi_testdata.txt') 
    242         answer = Loader().load(filepath) 
    243         for i in range(len(o.x)): 
    244             self.assertAlmostEqual(o.x[i], answer.x[i], 4) 
    245             self.assertAlmostEqual(o.y[i], answer.y[i], 4) 
    246             self.assertAlmostEqual(o.dy[i], answer.dy[i], 4) 
    247  
     220        from sas.sascalc.dataloader.manipulations import SectorPhi 
     221        import math 
     222         
     223        r = SectorPhi(r_min=.005, r_max=.01, phi_min=0, phi_max=math.pi/2.0) 
     224        r.nbins_phi = 20 
     225        o = r(self.data) 
     226 
     227        answer = Loader().load('sectorphi_testdata.txt') 
     228        for i in range(len(o.x)): 
     229            self.assertAlmostEqual(o.x[i], answer.x[i], 4) 
     230            self.assertAlmostEqual(o.y[i], answer.y[i], 4) 
     231            self.assertAlmostEqual(o.dy[i], answer.dy[i], 4) 
     232             
    248233    def test_sectorq_full(self): 
    249234        """ 
     
    251236            The test data was not generated by IGOR. 
    252237        """ 
    253  
    254         r = SectorQ(r_min=.005, r_max=.01, phi_min=0, phi_max=math.pi / 2.0) 
    255         r.nbins_phi = 20 
    256         o = r(self.data) 
    257  
    258         filepath = os.path.join(os.path.dirname( 
    259             os.path.realpath(__file__)), 'sectorq_testdata.txt') 
    260         answer = Loader().load(filepath) 
    261         for i in range(len(o.x)): 
    262             self.assertAlmostEqual(o.x[i], answer.x[i], 4) 
    263             self.assertAlmostEqual(o.y[i], answer.y[i], 4) 
    264             self.assertAlmostEqual(o.dy[i], answer.dy[i], 4) 
    265  
    266     def test_sectorq_log(self): 
    267         """ 
    268             Test sector averaging I(q) 
    269             The test data was not generated by IGOR. 
    270         """ 
    271  
    272         r = SectorQ(r_min=.005, r_max=.01, phi_min=0, phi_max=math.pi / 2.0, base=10) 
    273         r.nbins_phi = 20 
    274         o = r(self.data) 
    275  
    276         expected_binning = np.logspace(np.log10(0.005), np.log10(0.01), 20, base=10) 
    277         for i in range(len(o.x)): 
    278             self.assertAlmostEqual(o.x[i], expected_binning[i], 3) 
    279          
    280         # TODO: Test for Y values (o.y) 
    281         # print len(self.data.x_bins) 
    282         # print len(self.data.y_bins) 
    283         # print self.data.q_data.shape 
    284         # data_to_bin = np.array(self.data.q_data) 
    285         # data_to_bin = data_to_bin.reshape(len(self.data.x_bins), len(self.data.y_bins)) 
    286         # H, xedges, yedges, binnumber = binned_statistic_2d(self.data.x_bins, self.data.y_bins, data_to_bin, 
    287         #     bins=[[0, math.pi / 2.0], expected_binning], statistic='mean') 
    288         # xedges_width = (xedges[1] - xedges[0]) 
    289         # xedges_center = xedges[1:] - xedges_width / 2 
    290          
    291         # yedges_width = (yedges[1] - yedges[0]) 
    292         # yedges_center = yedges[1:] - yedges_width / 2 
    293          
    294         # print H.flatten().shape 
    295         # print o.y.shape 
    296          
     238        from sas.sascalc.dataloader.manipulations import SectorQ 
     239        import math 
     240         
     241        r = SectorQ(r_min=.005, r_max=.01, phi_min=0, phi_max=math.pi/2.0) 
     242        r.nbins_phi = 20 
     243        o = r(self.data) 
     244 
     245        answer = Loader().load('sectorq_testdata.txt') 
     246        for i in range(len(o.x)): 
     247            self.assertAlmostEqual(o.x[i], answer.x[i], 4) 
     248            self.assertAlmostEqual(o.y[i], answer.y[i], 4) 
     249            self.assertAlmostEqual(o.dy[i], answer.dy[i], 4) 
     250             
    297251 
    298252if __name__ == '__main__': 
  • test/sasguiframe/test/utest_manipulations.py

    r959eb01 r959eb01  
    22    Unit tests for data manipulations 
    33""" 
    4 # TODO: what happens if you add a Data1D to a Data2D? 
    5  
     4#TODO: what happens if you add a Data1D to a Data2D? 
     5 
     6import unittest 
    67import math 
     8import numpy as np 
     9from sas.sascalc.dataloader.loader import  Loader 
     10from sas.sasgui.guiframe.dataFitting import Data1D, Data2D 
     11from sas.sasgui.guiframe.dataFitting import Data1D as Theory1D 
     12  
    713import os.path 
    8 import unittest 
    9  
    10 import numpy as np 
    11  
    12 from sas.sascalc.dataloader.loader import Loader 
    13 from sas.sasgui.guiframe.dataFitting import Data1D 
    14 from sas.sasgui.guiframe.dataFitting import Data2D 
    15  
    16  
    17 class DataInfoTests(unittest.TestCase): 
    18  
     14 
     15class data_info_tests(unittest.TestCase): 
     16     
    1917    def setUp(self): 
    2018        data = Loader().load("cansas1d.xml") 
    2119        self.data = data[0] 
    22  
     20         
    2321    def test_clone1D(self): 
    2422        """ 
     
    2624        """ 
    2725        clone = self.data.clone_without_data() 
    28  
     26         
    2927        for i in range(len(self.data.detector)): 
    30             self.assertEqual( 
    31                 self.data.detector[i].distance, clone.detector[i].distance) 
    32  
    33  
    34 class Theory1DTests(unittest.TestCase): 
    35  
     28            self.assertEqual(self.data.detector[i].distance, clone.detector[i].distance) 
     29             
     30class theory1d_tests(unittest.TestCase): 
     31     
    3632    def setUp(self): 
    3733        data = Loader().load("cansas1d.xml") 
    3834        self.data = data[0] 
    39  
     35         
    4036    def test_clone_theory1D(self): 
    4137        """ 
    4238            Test basic cloning 
    4339        """ 
    44         theory = Data1D(x=[], y=[], dy=None) 
     40        theory = Theory1D(x=[], y=[], dy=None) 
    4541        theory.clone_without_data(clone=self.data) 
    4642        theory.copy_from_datainfo(data1d=self.data) 
    4743        for i in range(len(self.data.detector)): 
    48             self.assertEqual( 
    49                 self.data.detector[i].distance, theory.detector[i].distance) 
    50  
     44            self.assertEqual(self.data.detector[i].distance, theory.detector[i].distance) 
     45             
    5146        for i in range(len(self.data.x)): 
    5247            self.assertEqual(self.data.x[i], theory.x[i]) 
     
    5449            self.assertEqual(self.data.dy[i], theory.dy[i]) 
    5550 
    56  
    57 class ManipTests(unittest.TestCase): 
    58  
     51class manip_tests(unittest.TestCase): 
     52     
    5953    def setUp(self): 
    6054        # Create two data sets to play with 
    6155        x_0 = np.ones(5) 
    6256        for i in range(5): 
    63             x_0[i] = x_0[i] * (i + 1.0) 
    64  
    65         y_0 = 2.0 * np.ones(5) 
    66         dy_0 = 0.5 * np.ones(5) 
     57            x_0[i] = x_0[i]*(i+1.0) 
     58             
     59        y_0 = 2.0*np.ones(5) 
     60        dy_0 = 0.5*np.ones(5) 
    6761        self.data = Data1D(x_0, y_0, dy=dy_0) 
    68  
     62         
    6963        x = self.data.x 
    7064        y = np.ones(5) 
    7165        dy = np.ones(5) 
    7266        self.data2 = Data1D(x, y, dy=dy) 
    73  
     67         
     68         
    7469    def test_load(self): 
    7570        """ 
     
    7873        # There should be 5 entries in the file 
    7974        self.assertEqual(len(self.data.x), 5) 
    80  
     75         
    8176        for i in range(5): 
    8277            # The x values should be from 1 to 5 
    83             self.assertEqual(self.data.x[i], float(i + 1)) 
    84  
     78            self.assertEqual(self.data.x[i], float(i+1)) 
     79         
    8580            # All y-error values should be 0.5 
    86             self.assertEqual(self.data.dy[i], 0.5) 
    87  
     81            self.assertEqual(self.data.dy[i], 0.5)     
     82             
    8883            # All y values should be 2.0 
    89             self.assertEqual(self.data.y[i], 2.0) 
    90  
     84            self.assertEqual(self.data.y[i], 2.0)     
     85         
    9186    def test_add(self): 
    92         result = self.data2 + self.data 
     87        result = self.data2+self.data 
    9388        for i in range(5): 
    9489            self.assertEqual(result.y[i], 3.0) 
    95             self.assertEqual(result.dy[i], math.sqrt(0.5**2 + 1.0)) 
    96  
     90            self.assertEqual(result.dy[i], math.sqrt(0.5**2+1.0)) 
     91         
    9792    def test_sub(self): 
    98         result = self.data2 - self.data 
     93        result = self.data2-self.data 
    9994        for i in range(5): 
    10095            self.assertEqual(result.y[i], -1.0) 
    101             self.assertEqual(result.dy[i], math.sqrt(0.5**2 + 1.0)) 
    102  
     96            self.assertEqual(result.dy[i], math.sqrt(0.5**2+1.0)) 
     97         
    10398    def test_mul(self): 
    104         result = self.data2 * self.data 
     99        result = self.data2*self.data 
    105100        for i in range(5): 
    106101            self.assertEqual(result.y[i], 2.0) 
    107             self.assertEqual(result.dy[i], math.sqrt( 
    108                 (0.5 * 1.0)**2 + (1.0 * 2.0)**2)) 
    109  
     102            self.assertEqual(result.dy[i], math.sqrt((0.5*1.0)**2+(1.0*2.0)**2)) 
     103         
    110104    def test_div(self): 
    111         result = self.data2 / self.data 
     105        result = self.data2/self.data 
    112106        for i in range(5): 
    113107            self.assertEqual(result.y[i], 0.5) 
    114             self.assertEqual(result.dy[i], math.sqrt( 
    115                 (1.0 / 2.0)**2 + (0.5 * 1.0 / 4.0)**2)) 
    116  
     108            self.assertEqual(result.dy[i], math.sqrt((1.0/2.0)**2+(0.5*1.0/4.0)**2)) 
     109         
    117110    def test_radd(self): 
    118         result = self.data + 3.0 
     111        result = self.data+3.0 
    119112        for i in range(5): 
    120113            self.assertEqual(result.y[i], 5.0) 
    121114            self.assertEqual(result.dy[i], 0.5) 
    122  
    123         result = 3.0 + self.data 
     115             
     116        result = 3.0+self.data 
    124117        for i in range(5): 
    125118            self.assertEqual(result.y[i], 5.0) 
    126119            self.assertEqual(result.dy[i], 0.5) 
    127  
     120             
    128121    def test_rsub(self): 
    129         result = self.data - 3.0 
     122        result = self.data-3.0 
    130123        for i in range(5): 
    131124            self.assertEqual(result.y[i], -1.0) 
    132125            self.assertEqual(result.dy[i], 0.5) 
    133  
    134         result = 3.0 - self.data 
     126             
     127        result = 3.0-self.data 
    135128        for i in range(5): 
    136129            self.assertEqual(result.y[i], 1.0) 
    137130            self.assertEqual(result.dy[i], 0.5) 
    138  
     131             
    139132    def test_rmul(self): 
    140         result = self.data * 3.0 
     133        result = self.data*3.0 
    141134        for i in range(5): 
    142135            self.assertEqual(result.y[i], 6.0) 
    143136            self.assertEqual(result.dy[i], 1.5) 
    144  
    145         result = 3.0 * self.data 
     137             
     138        result = 3.0*self.data 
    146139        for i in range(5): 
    147140            self.assertEqual(result.y[i], 6.0) 
    148141            self.assertEqual(result.dy[i], 1.5) 
    149  
     142             
    150143    def test_rdiv(self): 
    151         result = self.data / 4.0 
     144        result = self.data/4.0 
    152145        for i in range(5): 
    153146            self.assertEqual(result.y[i], 0.5) 
    154147            self.assertEqual(result.dy[i], 0.125) 
    155  
    156         result = 6.0 / self.data 
     148             
     149        result = 6.0/self.data 
    157150        for i in range(5): 
    158151            self.assertEqual(result.y[i], 3.0) 
    159             self.assertEqual(result.dy[i], 6.0 * 0.5 / 4.0) 
    160  
    161  
    162 class Manin2DTests(unittest.TestCase): 
    163  
     152            self.assertEqual(result.dy[i], 6.0*0.5/4.0) 
     153             
     154class manip_2D(unittest.TestCase): 
     155     
    164156    def setUp(self): 
    165157        # Create two data sets to play with 
    166         x_0 = 2.0 * np.ones(25) 
    167         dx_0 = 0.5 * np.ones(25) 
     158        x_0 = 2.0*np.ones(25) 
     159        dx_0 = 0.5*np.ones(25) 
    168160        qx_0 = np.arange(25) 
    169161        qy_0 = np.arange(25) 
    170162        mask_0 = np.zeros(25) 
    171         dqx_0 = np.arange(25) / 100 
    172         dqy_0 = np.arange(25) / 100 
     163        dqx_0 = np.arange(25)/100 
     164        dqy_0 = np.arange(25)/100 
    173165        q_0 = np.sqrt(qx_0 * qx_0 + qy_0 * qy_0) 
    174         self.data = Data2D(image=x_0, err_image=dx_0, qx_data=qx_0, 
    175                            qy_data=qy_0, q_data=q_0, mask=mask_0, 
     166        self.data = Data2D(image=x_0, err_image=dx_0, qx_data=qx_0,  
     167                           qy_data=qy_0, q_data=q_0, mask=mask_0,  
    176168                           dqx_data=dqx_0, dqy_data=dqy_0) 
    177  
     169         
    178170        y = np.ones(25) 
    179171        dy = np.ones(25) 
     
    182174        mask = np.zeros(25) 
    183175        q = np.sqrt(qx * qx + qy * qy) 
    184         self.data2 = Data2D(image=y, err_image=dy, qx_data=qx, qy_data=qy, 
     176        self.data2 = Data2D(image=y, err_image=dy, qx_data=qx, qy_data=qy,  
    185177                            q_data=q, mask=mask) 
    186  
     178         
     179         
    187180    def test_load(self): 
    188181        """ 
     
    191184        # There should be 5 entries in the file 
    192185        self.assertEqual(np.size(self.data.data), 25) 
    193  
     186         
    194187        for i in range(25): 
    195188            # All y-error values should be 0.5 
    196             self.assertEqual(self.data.err_data[i], 0.5) 
    197  
     189            self.assertEqual(self.data.err_data[i], 0.5)     
     190             
    198191            # All y values should be 2.0 
    199             self.assertEqual(self.data.data[i], 2.0) 
    200  
     192            self.assertEqual(self.data.data[i], 2.0)     
     193         
    201194    def test_add(self): 
    202         result = self.data2 + self.data 
     195        result = self.data2+self.data 
    203196        for i in range(25): 
    204197            self.assertEqual(result.data[i], 3.0) 
    205             self.assertEqual(result.err_data[i], math.sqrt(0.5**2 + 1.0)) 
    206  
     198            self.assertEqual(result.err_data[i], math.sqrt(0.5**2+1.0)) 
     199         
    207200    def test_sub(self): 
    208         result = self.data2 - self.data 
     201        result = self.data2-self.data 
     202        for i in range(25): 
     203                self.assertEqual(result.data[i], -1.0) 
     204                self.assertEqual(result.err_data[i], math.sqrt(0.5**2+1.0)) 
     205         
     206    def test_mul(self): 
     207        result = self.data2*self.data 
     208        for i in range(25): 
     209            self.assertEqual(result.data[i], 2.0) 
     210            self.assertEqual(result.err_data[i], math.sqrt((0.5*1.0)**2+(1.0*2.0)**2)) 
     211         
     212    def test_div(self): 
     213        result = self.data2/self.data 
     214        for i in range(25): 
     215            self.assertEqual(result.data[i], 0.5) 
     216            self.assertEqual(result.err_data[i], math.sqrt((1.0/2.0)**2+(0.5*1.0/4.0)**2)) 
     217         
     218    def test_radd(self): 
     219        result = self.data+3.0 
     220        for i in range(25): 
     221            self.assertEqual(result.data[i], 5.0) 
     222            self.assertEqual(result.err_data[i], 0.5) 
     223    
     224        result = 3.0+self.data 
     225        for i in range(25): 
     226            self.assertEqual(result.data[i], 5.0) 
     227            self.assertEqual(result.err_data[i], 0.5) 
     228             
     229    def test_rsub(self): 
     230        result = self.data-3.0 
    209231        for i in range(25): 
    210232            self.assertEqual(result.data[i], -1.0) 
    211             self.assertEqual(result.err_data[i], math.sqrt(0.5**2 + 1.0)) 
    212  
    213     def test_mul(self): 
    214         result = self.data2 * self.data 
    215         for i in range(25): 
    216             self.assertEqual(result.data[i], 2.0) 
    217             self.assertEqual(result.err_data[i], math.sqrt( 
    218                 (0.5 * 1.0)**2 + (1.0 * 2.0)**2)) 
    219  
    220     def test_div(self): 
    221         result = self.data2 / self.data 
    222         for i in range(25): 
    223             self.assertEqual(result.data[i], 0.5) 
    224             self.assertEqual(result.err_data[i], math.sqrt( 
    225                 (1.0 / 2.0)**2 + (0.5 * 1.0 / 4.0)**2)) 
    226  
    227     def test_radd(self): 
    228         result = self.data + 3.0 
    229         for i in range(25): 
    230             self.assertEqual(result.data[i], 5.0) 
    231             self.assertEqual(result.err_data[i], 0.5) 
    232  
    233         result = 3.0 + self.data 
    234         for i in range(25): 
    235             self.assertEqual(result.data[i], 5.0) 
    236             self.assertEqual(result.err_data[i], 0.5) 
    237  
    238     def test_rsub(self): 
    239         result = self.data - 3.0 
    240         for i in range(25): 
    241             self.assertEqual(result.data[i], -1.0) 
    242             self.assertEqual(result.err_data[i], 0.5) 
    243  
    244         result = 3.0 - self.data 
     233            self.assertEqual(result.err_data[i], 0.5) 
     234     
     235        result = 3.0-self.data 
    245236        for i in range(25): 
    246237            self.assertEqual(result.data[i], 1.0) 
    247238            self.assertEqual(result.err_data[i], 0.5) 
    248  
     239             
    249240    def test_rmul(self): 
    250         result = self.data * 3.0 
     241        result = self.data*3.0 
    251242        for i in range(25): 
    252243            self.assertEqual(result.data[i], 6.0) 
    253244            self.assertEqual(result.err_data[i], 1.5) 
    254  
    255         result = 3.0 * self.data 
     245  
     246        result = 3.0*self.data 
    256247        for i in range(25): 
    257248            self.assertEqual(result.data[i], 6.0) 
    258249            self.assertEqual(result.err_data[i], 1.5) 
    259  
     250             
    260251    def test_rdiv(self): 
    261         result = self.data / 4.0 
     252        result = self.data/4.0 
    262253        for i in range(25): 
    263254            self.assertEqual(result.data[i], 0.5) 
    264255            self.assertEqual(result.err_data[i], 0.125) 
    265256 
    266         result = 6.0 / self.data 
     257        result = 6.0/self.data 
    267258        for i in range(25): 
    268259            self.assertEqual(result.data[i], 3.0) 
    269             self.assertEqual(result.err_data[i], 6.0 * 0.5 / 4.0) 
    270  
    271  
    272 class ExtraManip2DTests(unittest.TestCase): 
    273  
     260            self.assertEqual(result.err_data[i], 6.0*0.5/4.0) 
     261     
     262class extra_manip_2D(unittest.TestCase): 
     263     
    274264    def setUp(self): 
    275265        # Create two data sets to play with 
    276         x_0 = 2.0 * np.ones(25) 
    277         dx_0 = 0.5 * np.ones(25) 
     266        x_0 = 2.0*np.ones(25) 
     267        dx_0 = 0.5*np.ones(25) 
    278268        qx_0 = np.arange(25) 
    279269        qy_0 = np.arange(25) 
    280270        mask_0 = np.zeros(25) 
    281         dqx_0 = np.arange(25) / 100 
    282         dqy_0 = np.arange(25) / 100 
     271        dqx_0 = np.arange(25)/100 
     272        dqy_0 = np.arange(25)/100 
    283273        q_0 = np.sqrt(qx_0 * qx_0 + qy_0 * qy_0) 
    284         self.data = Data2D(image=x_0, err_image=dx_0, qx_data=qx_0, 
    285                            qy_data=qy_0, q_data=q_0, mask=mask_0, 
     274        self.data = Data2D(image=x_0, err_image=dx_0, qx_data=qx_0,  
     275                           qy_data=qy_0, q_data=q_0, mask=mask_0,  
    286276                           dqx_data=dqx_0, dqy_data=dqy_0) 
    287  
     277         
    288278        y = np.ones(25) 
    289279        dy = np.ones(25) 
     
    292282        mask = np.zeros(25) 
    293283        q = np.sqrt(qx * qx + qy * qy) 
    294         self.data2 = Data2D(image=y, err_image=dy, qx_data=qx, qy_data=qy, 
     284        self.data2 = Data2D(image=y, err_image=dy, qx_data=qx, qy_data=qy,  
    295285                            q_data=q, mask=mask) 
    296  
     286         
     287         
    297288    def test_load(self): 
    298289        """ 
     
    301292        # There should be 5 entries in the file 
    302293        self.assertEqual(np.size(self.data.data), 25) 
    303  
     294         
    304295        for i in range(25): 
    305296            # All y-error values should be 0.5 
    306             self.assertEqual(self.data.err_data[i], 0.5) 
    307  
     297            self.assertEqual(self.data.err_data[i], 0.5)     
     298             
    308299            # All y values should be 2.0 
    309             self.assertEqual(self.data.data[i], 2.0) 
    310  
     300            self.assertEqual(self.data.data[i], 2.0)     
     301         
    311302    def test_add(self): 
    312         result = self.data2 + self.data 
     303        result = self.data2+self.data 
    313304        for i in range(25): 
    314305            self.assertEqual(result.data[i], 3.0) 
    315             self.assertEqual(result.err_data[i], math.sqrt(0.5**2 + 1.0)) 
    316  
     306            self.assertEqual(result.err_data[i], math.sqrt(0.5**2+1.0)) 
     307         
    317308    def test_sub(self): 
    318         result = self.data2 - self.data 
     309        result = self.data2-self.data 
     310        for i in range(25): 
     311                self.assertEqual(result.data[i], -1.0) 
     312                self.assertEqual(result.err_data[i], math.sqrt(0.5**2+1.0)) 
     313         
     314    def test_mul(self): 
     315        result = self.data2*self.data 
     316        for i in range(25): 
     317            self.assertEqual(result.data[i], 2.0) 
     318            self.assertEqual(result.err_data[i], math.sqrt((0.5*1.0)**2+(1.0*2.0)**2)) 
     319         
     320    def test_div(self): 
     321        result = self.data2/self.data 
     322        for i in range(25): 
     323            self.assertEqual(result.data[i], 0.5) 
     324            self.assertEqual(result.err_data[i], math.sqrt((1.0/2.0)**2+(0.5*1.0/4.0)**2)) 
     325         
     326    def test_radd(self): 
     327        result = self.data+3.0 
     328        for i in range(25): 
     329            self.assertEqual(result.data[i], 5.0) 
     330            self.assertEqual(result.err_data[i], 0.5) 
     331    
     332        result = 3.0+self.data 
     333        for i in range(25): 
     334            self.assertEqual(result.data[i], 5.0) 
     335            self.assertEqual(result.err_data[i], 0.5) 
     336             
     337    def test_rsub(self): 
     338        result = self.data-3.0 
    319339        for i in range(25): 
    320340            self.assertEqual(result.data[i], -1.0) 
    321             self.assertEqual(result.err_data[i], math.sqrt(0.5**2 + 1.0)) 
    322  
    323     def test_mul(self): 
    324         result = self.data2 * self.data 
    325         for i in range(25): 
    326             self.assertEqual(result.data[i], 2.0) 
    327             self.assertEqual(result.err_data[i], math.sqrt( 
    328                 (0.5 * 1.0)**2 + (1.0 * 2.0)**2)) 
    329  
    330     def test_div(self): 
    331         result = self.data2 / self.data 
    332         for i in range(25): 
    333             self.assertEqual(result.data[i], 0.5) 
    334             self.assertEqual(result.err_data[i], math.sqrt( 
    335                 (1.0 / 2.0)**2 + (0.5 * 1.0 / 4.0)**2)) 
    336  
    337     def test_radd(self): 
    338         result = self.data + 3.0 
    339         for i in range(25): 
    340             self.assertEqual(result.data[i], 5.0) 
    341             self.assertEqual(result.err_data[i], 0.5) 
    342  
    343         result = 3.0 + self.data 
    344         for i in range(25): 
    345             self.assertEqual(result.data[i], 5.0) 
    346             self.assertEqual(result.err_data[i], 0.5) 
    347  
    348     def test_rsub(self): 
    349         result = self.data - 3.0 
    350         for i in range(25): 
    351             self.assertEqual(result.data[i], -1.0) 
    352             self.assertEqual(result.err_data[i], 0.5) 
    353  
    354         result = 3.0 - self.data 
     341            self.assertEqual(result.err_data[i], 0.5) 
     342     
     343        result = 3.0-self.data 
    355344        for i in range(25): 
    356345            self.assertEqual(result.data[i], 1.0) 
    357346            self.assertEqual(result.err_data[i], 0.5) 
    358  
     347             
    359348    def test_rmul(self): 
    360         result = self.data * 3.0 
     349        result = self.data*3.0 
    361350        for i in range(25): 
    362351            self.assertEqual(result.data[i], 6.0) 
    363352            self.assertEqual(result.err_data[i], 1.5) 
    364  
    365         result = 3.0 * self.data 
     353  
     354        result = 3.0*self.data 
    366355        for i in range(25): 
    367356            self.assertEqual(result.data[i], 6.0) 
    368357            self.assertEqual(result.err_data[i], 1.5) 
    369  
     358             
    370359    def test_rdiv(self): 
    371         result = self.data / 4.0 
     360        result = self.data/4.0 
    372361        for i in range(25): 
    373362            self.assertEqual(result.data[i], 0.5) 
    374363            self.assertEqual(result.err_data[i], 0.125) 
    375364 
    376         result = 6.0 / self.data 
     365        result = 6.0/self.data 
    377366        for i in range(25): 
    378367            self.assertEqual(result.data[i], 3.0) 
    379             self.assertEqual(result.err_data[i], 6.0 * 0.5 / 4.0) 
    380  
     368            self.assertEqual(result.err_data[i], 6.0*0.5/4.0) 
    381369 
    382370if __name__ == '__main__': 
    383371    unittest.main() 
     372    
Note: See TracChangeset for help on using the changeset viewer.