Changes in / [43345ab:e3f1fed] in sasview
- Files:
-
- 2 added
- 30 deleted
- 26 edited
Legend:
- Unmodified
- Added
- Removed
-
check_packages.py
r3e1f417 r235f514 72 72 for package_name, test_vals in win_required_package_list.items(): 73 73 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']))) 76 81 except ImportError: 77 82 print('%s NOT INSTALLED'% package_name) -
docs/sphinx-docs/build_sphinx.py
r4abc05d8 r01f1e17 38 38 #/sasview-local-trunk/docs/sphinx-docs/build_sphinx.py 39 39 SASMODELS_SOURCE_PROLOG = os.path.join(CURRENT_SCRIPT_DIR, "..", "..", "..", "sasmodels", "doc") 40 SASMODELS_SOURCE_GPU = os.path.join(CURRENT_SCRIPT_DIR, "..", "..", "..", "sasmodels", "doc", " guide", "gpu")41 SASMODELS_SOURCE_SESANS = os.path.join(CURRENT_SCRIPT_DIR, "..", "..", "..", "sasmodels", "doc", " guide", "sesans")42 SASMODELS_SOURCE_SESANSIMG = os.path.join(CURRENT_SCRIPT_DIR, "..", "..", "..", "sasmodels", "doc", " guide", "sesans", "sesans_img")43 SASMODELS_SOURCE_MAGNETISM = os.path.join(CURRENT_SCRIPT_DIR, "..", "..", "..", "sasmodels", "doc", " guide", "magnetism")44 SASMODELS_SOURCE_MAGIMG = os.path.join(CURRENT_SCRIPT_DIR, "..", "..", "..", "sasmodels", "doc", " guide", "magnetism", "mag_img")45 SASMODELS_SOURCE_REF_MODELS = os.path.join(CURRENT_SCRIPT_DIR, "..", "..", "..", "sasmodels", "doc", " guide", "models")40 SASMODELS_SOURCE_GPU = os.path.join(CURRENT_SCRIPT_DIR, "..", "..", "..", "sasmodels", "doc", "ref", "gpu") 41 SASMODELS_SOURCE_SESANS = os.path.join(CURRENT_SCRIPT_DIR, "..", "..", "..", "sasmodels", "doc", "ref", "sesans") 42 SASMODELS_SOURCE_SESANSIMG = os.path.join(CURRENT_SCRIPT_DIR, "..", "..", "..", "sasmodels", "doc", "ref", "sesans", "sesans_img") 43 SASMODELS_SOURCE_MAGNETISM = os.path.join(CURRENT_SCRIPT_DIR, "..", "..", "..", "sasmodels", "doc", "ref", "magnetism") 44 SASMODELS_SOURCE_MAGIMG = os.path.join(CURRENT_SCRIPT_DIR, "..", "..", "..", "sasmodels", "doc", "ref", "magnetism", "mag_img") 45 SASMODELS_SOURCE_REF_MODELS = os.path.join(CURRENT_SCRIPT_DIR, "..", "..", "..", "sasmodels", "doc", "ref", "models") 46 46 SASMODELS_SOURCE_MODELS = os.path.join(CURRENT_SCRIPT_DIR, "..", "..", "..", "sasmodels", "doc", "model") 47 47 SASMODELS_SOURCE_IMG = os.path.join(CURRENT_SCRIPT_DIR, "..", "..", "..", "sasmodels", "doc", "model", "img") -
run.py
rf94a935 rf36e01f 51 51 def import_package(modname, path): 52 52 """Import a package into a particular point in the python namespace""" 53 #logger.debug("Dynamicly importing: %s", path)54 53 mod = imp.load_source(modname, abspath(joinpath(path, '__init__.py'))) 55 54 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 1 DataFileTitle "Polystyrene of Markus Strobl, Full Sine, ++ only " 2 Sample "Polystyrene 2 um in 53% H2O, 47% D2O " 3 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 " 4 Operator CPD 5 Date do 10 jul 2014 16:37:30 6 ScanType sine one element scan 7 Thickness [cm] 2.00E-01 8 Q_zmax [\A^-1] 0.05 9 Q_ymax [\A^-1] 0.05 10 11 spin 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 12 391.56 0.0041929 0.0036894 19.578 2.11 0.1055 1.0037 0.0032974 13 1564 -0.0046571 0.0038185 78.2 2.11 0.1055 0.99586 0.003386 14 2735.6 -0.017007 0.0038132 136.78 2.11 0.1055 0.98497 0.0033444 15 3907.9 -0.033462 0.0035068 195.39 2.11 0.1055 0.97064 0.0030309 16 5080.2 -0.047483 0.0038208 254.01 2.11 0.1055 0.9586 0.0032613 17 6251.8 -0.070375 0.00376 312.59 2.11 0.1055 0.93926 0.0031446 18 7423.2 -0.092217 0.0037927 371.16 2.11 0.1055 0.92117 0.0031108 19 8595.5 -0.10238 0.004006 429.77 2.11 0.1055 0.91287 0.0032562 20 9767.7 -0.12672 0.0038534 488.39 2.11 0.1055 0.8933 0.0030651 21 10940 -0.1374 0.004243 546.98 2.11 0.1055 0.88484 0.003343 22 12112 -0.16072 0.0045837 605.58 2.11 0.1055 0.86666 0.0035372 23 13284 -0.16623 0.0045613 664.2 2.11 0.1055 0.86242 0.0035027 24 14456 -0.18468 0.0044918 722.79 2.11 0.1055 0.84837 0.0033931 25 15628 -0.19143 0.0048967 781.38 2.11 0.1055 0.84328 0.0036768 26 16800 -0.20029 0.0045421 840.02 2.11 0.1055 0.83666 0.0033837 27 17971 -0.19798 0.0046642 898.56 2.11 0.1055 0.83838 0.0034819 28 19143 -0.21442 0.0047052 957.17 2.11 0.1055 0.82619 0.0034614 29 20316 -0.20885 0.0044931 1015.8 2.11 0.1055 0.8303 0.0033218 30 21488 -0.21393 0.0049186 1074.4 2.11 0.1055 0.82655 0.00362 31 22660 -0.20685 0.004423 1133 2.11 0.1055 0.83179 0.0032758 32 23832 -0.20802 0.0046979 1191.6 2.11 0.1055 0.83092 0.0034758 33 25003 -0.19848 0.0045953 1250.2 2.11 0.1055 0.838 0.0034289 34 26175 -0.21117 0.0044567 1308.8 2.11 0.1055 0.82859 0.0032881 35 27347 -0.21283 0.004137 1367.4 2.11 0.1055 0.82736 0.0030477 36 28520 -0.2042 0.0044587 1426 2.11 0.1055 0.83375 0.0033101 37 29692 -0.2112 0.0042852 1484.6 2.11 0.1055 0.82857 0.0031615 38 30864 -0.20319 0.0043483 1543.2 2.11 0.1055 0.8345 0.003231 39 32036 -0.20752 0.0044297 1601.8 2.11 0.1055 0.83129 0.0032788 40 33207 -0.20654 0.0043188 1660.4 2.11 0.1055 0.83201 0.0031995 41 34380 -0.20126 0.0046375 1719 2.11 0.1055 0.83593 0.0034518 42 35551 -0.20924 0.0042871 1777.6 2.11 0.1055 0.83001 0.0031684 43 36724 -0.21323 0.0045471 1836.2 2.11 0.1055 0.82707 0.0033487 44 37895 -0.21324 0.0045354 1894.7 2.11 0.1055 0.82706 0.00334 45 39067 -0.19905 0.0044141 1953.4 2.11 0.1055 0.83758 0.003292 46 40239 -0.1991 0.0047441 2012 2.11 0.1055 0.83754 0.003538 47 41411 -0.20359 0.0050136 2070.5 2.11 0.1055 0.8342 0.003724 48 42583 -0.21032 0.0049474 2129.1 2.11 0.1055 0.82922 0.0036529 49 43755 -0.20689 0.0048203 2187.8 2.11 0.1055 0.83176 0.00357 50 44927 -0.21075 0.0052337 2246.4 2.11 0.1055 0.8289 0.0038628 51 46099 -0.19956 0.0047827 2304.9 2.11 0.1055 0.8372 0.0035653 -
src/sas/sascalc/data_util/nxsunit.py
r8e9536f rb699768 136 136 sld = { '10^-6 Angstrom^-2': 1e-6, 'Angstrom^-2': 1 } 137 137 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, 139 139 'nm^-1': 0.1, '1/nm': 0.1, 'n_m^-1': 0.1 } 140 140 … … 195 195 _check(1,Converter('n_m^-1')(10,'invA')) # 10 nm^-1 = 1 inv Angstroms 196 196 _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/m198 197 _check(0.003,Converter('microseconds')(3,units='ms')) # 3 us -> 0.003 ms 199 198 _check(45,Converter('nanokelvin')(45)) # 45 nK -> 45 nK -
src/sas/sascalc/data_util/qsmearing.py
r235f514 r235f514 65 65 raise ValueError('one or more of your dx values are negative, please check the data file!') 66 66 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 69 70 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"), 75 72 Rmax = 10000000 76 hankel = SesansTransform(data.x, SElength, 77 data.source.wavelength, 78 zaccept, Rmax) 73 hankel = SesansTransform(data.x, SElength, zaccept, Rmax) 79 74 # Then return the actual transform, as if it were a smearing function 80 75 return PySmear(hankel, model, offset=0) -
src/sas/sascalc/dataloader/manipulations.py
r324e0bf r7432acb 1 from __future__ import division2 1 """ 3 2 Data manipulations for 2D data sets. 4 3 Using the meta data information, various types of averaging 5 4 are performed in Q-space 6 7 To test this module use:8 ```9 cd test10 PYTHONPATH=../src/ python2 -m sasdataloader.test.utest_averaging DataInfoTests.test_sectorphi_quarter11 ```12 5 """ 13 6 ##################################################################### 14 # 15 # 16 # 17 # 18 # 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 19 12 ###################################################################### 20 13 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 23 15 import math 24 import numpy as np 25 import sys 16 import numpy 26 17 27 18 #from data_info import plottable_2D … … 79 70 return phi_out 80 71 72 73 def 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 112 class _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 235 class 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 249 class 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 262 class 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 345 class 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 81 370 def get_pixel_fraction_square(x, xmin, xmax): 82 371 """ … … 101 390 return 1.0 102 391 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 393 class 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 546 class 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 126 648 127 649 def get_pixel_fraction(qmax, q_00, q_01, q_10, q_11): … … 188 710 return frac_max 189 711 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 713 def 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) 255 732 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 738 class _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 338 764 :return: Data1D object 339 765 """ 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 384 844 if frac == 0: 385 845 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: 400 872 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)) 411 878 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 424 931 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 940 class 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 """ 435 948 def __call__(self, data2D): 436 949 """ 437 Compute average I(Qy) for a region of interest950 Perform sector average and return I(phi). 438 951 439 952 :param data2D: Data2D object 440 953 :return: Data1D object 441 954 """ 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 958 class 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 """ 450 967 def __call__(self, data2D): 451 968 """ 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 454 973 :return: Data1D object 455 974 """ 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 978 class 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 1021 class Boxcut(object): 1022 """ 1023 Find a rectangular 2D region of interest. 1024 """ 465 1025 def __init__(self, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0): 466 1026 # Minimum Qx value [A-1] … … 475 1035 def __call__(self, data2D): 476 1036 """ 477 Perform the sum in the region of interest478 479 :param data2D: Data2D object480 :return: number of counts, error on number of counts,481 number of points summed482 """483 y, err_y, y_counts = self._sum(data2D)484 485 # Average the sums486 counts = 0 if y_counts == 0 else y487 error = 0 if y_counts == 0 else math.sqrt(err_y)488 489 # Added y_counts to return, SMK & PDB, 04/03/2013490 return counts, error, y_counts491 492 def _sum(self, data2D):493 """494 Perform the sum in the region of interest495 496 :param data2D: Data2D object497 :return: number of counts,498 error on number of counts, number of entries summed499 """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 data505 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.0511 err_y = 0.0512 y_counts = 0.0513 514 # Average pixelsize in q space515 for npts in range(len(data)):516 # default frac517 frac_x = 0518 frac_y = 0519 520 # get min and max at each points521 qx = qx_data[npts]522 qy = qy_data[npts]523 524 # get the ROI525 if self.x_min <= qx and self.x_max > qx:526 frac_x = 1527 if self.y_min <= qy and self.y_max > qy:528 frac_y = 1529 # Find the fraction along each directions530 frac = frac_x * frac_y531 if frac == 0:532 continue533 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 += frac541 return y, err_y, y_counts542 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 interest556 557 :param data2D: Data2D object558 :return: average counts, error on average counts559 560 """561 y, err_y, y_counts = self._sum(data2D)562 563 # Average the sums564 counts = 0 if y_counts == 0 else y / y_counts565 error = 0 if y_counts == 0 else math.sqrt(err_y) / y_counts566 567 return counts, error568 569 ################################################################################570 571 class CircularAverage(object):572 """573 Perform circular averaging on 2D data574 575 The data returned is the distribution of counts576 as a function of Q577 """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_min582 # Maximum radius included in the average [A-1]583 self.r_max = r_max584 # Bin width (step size) [A-1]585 self.bin_width = bin_width586 587 def __call__(self, data2D, ismask=False):588 """589 Perform circular averaging on the data590 591 :param data2D: Data2D object592 :return: Data1D object593 """594 # Get data W/ finite values595 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 = None601 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_data606 raise RuntimeError(msg)607 608 # Build array of Q intervals609 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 continue621 622 frac = 0623 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 range629 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 = 1634 if frac == 0:635 continue636 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 - 1641 y[i_q] += frac * data_n642 # Take dqs from data to get the q_average643 x[i_q] += frac * q_value644 if err_data is None or err_data[npt] == 0.0:645 if data_n < 0:646 data_n = -data_n647 err_y[i_q] += frac * frac * data_n648 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) because653 # it should not depend on the number of the q points654 # in the qr bins.655 err_x[i_q] += frac * dq_data[npt]656 else:657 err_x = None658 y_counts[i_q] += frac659 660 # Average the sums661 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_counts669 err_y[err_y == 0] = np.average(err_y)670 y = y / y_counts671 x = x / y_counts672 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 = None678 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, and691 the position of the center of the ring.692 693 The data returned is the distribution of counts694 around the ring as a function of phi.695 696 Phi_min and phi_max should be defined between 0 and 2*pi697 in anti-clockwise starting from the x- axis on the left-hand side698 """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 radius703 self.r_min = r_min704 # Maximum radius705 self.r_max = r_max706 # Center of the ring in x707 self.center_x = center_x708 # Center of the ring in y709 self.center_y = center_y710 # Number of angular bins711 self.nbins_phi = nbins712 713 def __call__(self, data2D):714 """715 Apply the ring to the data set.716 Returns the angular distribution for a given q range717 718 :param data2D: Data2D object719 720 :return: Data1D object721 """722 if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:723 raise RuntimeError("Ring averaging only take plottable_2D objects")724 725 Pi = math.pi726 727 # Get data728 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 outputs735 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 order741 # to center first bin at zero742 phi_shift = Pi / self.nbins_phi743 744 for npt in range(len(data)):745 frac = 0746 # 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]) + Pi752 753 if self.r_min <= q_value and q_value <= self.r_max:754 frac = 1755 if frac == 0:756 continue757 # binning758 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 = 0764 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_n769 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] += frac773 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 ring795 where phi_min and phi_max are defined by the right796 and left lines wrt central line797 and phi_max could be less than phi_min.798 799 Phi is defined between 0 and 2*pi in anti-clockwise800 starting from the x- axis on the left-hand side801 """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 number808 '''809 self.r_min = r_min810 self.r_max = r_max811 self.phi_min = phi_min812 self.phi_max = phi_max813 self.nbins = nbins814 self.base = base815 816 def _agv(self, data2D, run='phi'):817 """818 Perform sector averaging.819 820 :param data2D: Data2D object821 :param run: define the varying parameter ('phi' , 'q' , or 'q2')822 823 :return: Data1D object824 """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 & info829 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 = None836 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 outputs840 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 < 2Pi847 phi_min = flip_phi(self.phi_min)848 phi_max = flip_phi(self.phi_max)849 850 # binning object851 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 = False864 865 # phi-value of the pixel (j,i)866 phi_value = math.atan2(qy_data[n], qx_data[n]) + math.pi867 868 # No need to calculate: data outside of the radius869 if self.r_min > q_value or q_value > self.r_max:870 continue871 872 # In case of two ROIs (symmetric major and minor regions)(for 'q2')873 if run.lower() == 'q2':874 # For minor sector wing875 # Calculate the minor wing phis876 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 2pi879 if phi_min_minor > phi_max_minor:880 is_in = (phi_value > phi_min_minor or881 phi_value < phi_max_minor)882 else:883 is_in = (phi_value > phi_min_minor and884 phi_value < phi_max_minor)885 886 # For all cases(i.e.,for 'q', 'q2', and 'phi')887 # Find pixels within ROI888 if phi_min > phi_max:889 is_in = is_in or (phi_value > phi_min or890 phi_value < phi_max)891 else:892 is_in = is_in or (phi_value >= phi_min and893 phi_value < phi_max)894 895 # data oustide of the phi range896 if not is_in:897 continue898 899 # Get the binning index900 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 - 1908 909 # Get the total y910 y[i_bin] += data_n911 x[i_bin] += q_value912 if err_data[n] is None or err_data[n] == 0.0:913 if data_n < 0:914 data_n = -data_n915 y_err[i_bin] += data_n916 else:917 y_err[i_bin] += err_data[n]**2918 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) because922 # it should not depend on the number of the q points923 # in the qr bins.924 x_err[i_bin] += dq_data[n]925 else:926 x_err = None927 y_counts[i_bin] += 1928 929 # Organize the results930 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 q935 # Calculate x[i]should be at the center of the bin936 if run.lower() == 'phi':937 x[i] = (self.phi_max - self.phi_min) / self.nbins * \938 (1.0 * i + 0.5) + self.phi_min939 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.nbins943 # r_inner = self.r_min + delta_r * i944 # r_outer = r_inner + delta_r945 # 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 = None953 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 object976 :return: Data1D object977 """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 object996 997 :return: Data1D object998 """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, and1007 the position of the center of the ring.1008 1009 The data returned is the region inside the ring1010 1011 Phi_min and phi_max should be defined between 0 and 2*pi1012 in anti-clockwise starting from the x- axis on the left-hand side1013 """1014 1015 def __init__(self, r_min=0, r_max=0, center_x=0, center_y=0):1016 # Minimum radius1017 self.r_min = r_min1018 # Maximum radius1019 self.r_max = r_max1020 # Center of the ring in x1021 self.center_x = center_x1022 # Center of the ring in y1023 self.center_y = center_y1024 1025 def __call__(self, data2D):1026 """1027 Apply the ring to the data set.1028 Returns the angular distribution for a given q range1029 1030 :param data2D: Data2D object1031 1032 :return: index array in the range1033 """1034 if data2D.__class__.__name__ not in ["Data2D", "plottable_2D"]:1035 raise RuntimeError("Ring cut only take plottable_2D objects")1036 1037 # Get data1038 qx_data = data2D.qx_data1039 qy_data = data2D.qy_data1040 q_data = np.sqrt(qx_data * qx_data + qy_data * qy_data)1041 1042 # check whether or not the data point is inside ROI1043 out = (self.r_min <= q_data) & (self.r_max >= q_data)1044 return out1045 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_min1056 # Maximum Qx value [A-1]1057 self.x_max = x_max1058 # Minimum Qy value [A-1]1059 self.y_min = y_min1060 # Maximum Qy value [A-1]1061 self.y_max = y_max1062 1063 def __call__(self, data2D):1064 """1065 1037 Find a rectangular 2D region of interest. 1066 1038 … … 1083 1055 """ 1084 1056 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" 1086 1058 # Get qx_ and qy_data 1087 1059 qx_data = data2D.qx_data … … 1094 1066 return outx & outy 1095 1067 1096 ################################################################################1097 1068 1098 1069 class Sectorcut(object): … … 1106 1077 and (phi_max-phi_min) should not be larger than pi 1107 1078 """ 1108 1109 1079 def __init__(self, phi_min=0, phi_max=math.pi): 1110 1080 self.phi_min = phi_min … … 1136 1106 """ 1137 1107 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" 1139 1109 Pi = math.pi 1140 1110 # Get data … … 1143 1113 1144 1114 # get phi from data 1145 phi_data = n p.arctan2(qy_data, qx_data)1115 phi_data = numpy.arctan2(qy_data, qx_data) 1146 1116 1147 1117 # Get the min and max into the region: -pi <= phi < Pi … … 1150 1120 # check for major sector 1151 1121 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) 1154 1123 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) 1157 1125 1158 1126 # minor sector … … 1164 1132 if phi_min_minor > phi_max_minor: 1165 1133 out_minor = (phi_min_minor <= phi_data) + \ 1166 (phi_max_minor >= phi_data)1134 (phi_max_minor >= phi_data) 1167 1135 else: 1168 1136 out_minor = (phi_min_minor <= phi_data) & \ 1169 (phi_max_minor >= phi_data)1137 (phi_max_minor >= phi_data) 1170 1138 out = out_major + out_minor 1171 1139 -
src/sas/sascalc/dataloader/readers/sesans_reader.py
r149b8f6 r9a5097c 1 1 """ 2 2 SESANS reader (based on ASCII reader) 3 3 4 4 Reader for .ses or .sesans file format 5 6 Jurrian Bakker 5 6 Jurrian Bakker 7 7 """ 8 8 import numpy as np … … 18 18 _ZERO = 1e-16 19 19 20 21 20 class Reader: 22 21 """ 23 22 Class to load sesans files (6 columns). 24 23 """ 25 # File type24 ## File type 26 25 type_name = "SESANS" 27 28 # Wildcards26 27 ## Wildcards 29 28 type = ["SESANS files (*.ses)|*.ses", 30 29 "SESANS files (*..sesans)|*.sesans"] 31 # List of allowed extensions30 ## List of allowed extensions 32 31 ext = ['.ses', '.SES', '.sesans', '.SESANS'] 33 34 # Flag to bypass extension check32 33 ## Flag to bypass extension check 35 34 allow_all = True 36 35 37 36 def read(self, path): 37 38 # print "reader triggered" 39 38 40 """ 39 41 Load data file 40 42 41 43 :param path: file path 42 44 43 45 :return: SESANSData1D object, or None 44 46 45 47 :raise RuntimeError: when the file can't be opened 46 48 :raise ValueError: when the length of the data vectors are inconsistent … … 49 51 basename = os.path.basename(path) 50 52 _, 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 54 163 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 65 166 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) 173 171 new_unit = default_unit 174 172 else: 175 173 new_unit = value_unit 176 174 return value, new_unit 177 178 def _unit_fetch(self, unit):179 return self.params[unit+"_unit"] -
src/sas/sascalc/pr/c_extensions/Cinvertor.c
rcb62bd5 r959eb01 294 294 } 295 295 296 const char set_ est_bck_doc[] =296 const char set_has_bck_doc[] = 297 297 "Sets background flag\n"; 298 298 … … 300 300 * Sets the maximum distance 301 301 */ 302 static PyObject * set_ est_bck(Cinvertor *self, PyObject *args) {303 int est_bck;304 305 if (!PyArg_ParseTuple(args, "i", & est_bck)) return NULL;306 self->params. est_bck = est_bck;307 return Py_BuildValue("i", self->params. est_bck);308 } 309 310 const char get_ est_bck_doc[] =302 static PyObject * set_has_bck(Cinvertor *self, PyObject *args) { 303 int has_bck; 304 305 if (!PyArg_ParseTuple(args, "i", &has_bck)) return NULL; 306 self->params.has_bck = has_bck; 307 return Py_BuildValue("i", self->params.has_bck); 308 } 309 310 const char get_has_bck_doc[] = 311 311 "Gets background flag\n"; 312 312 … … 314 314 * Gets the maximum distance 315 315 */ 316 static PyObject * get_ est_bck(Cinvertor *self, PyObject *args) {317 return Py_BuildValue("i", self->params. est_bck);316 static PyObject * get_has_bck(Cinvertor *self, PyObject *args) { 317 return Py_BuildValue("i", self->params.has_bck); 318 318 } 319 319 … … 882 882 sqrt_alpha = sqrt(self->params.alpha); 883 883 pi = acos(-1.0); 884 offset = (self->params. est_bck==1) ? 0 : 1;884 offset = (self->params.has_bck==1) ? 0 : 1; 885 885 886 886 for (j=0; j<nfunc; j++) { … … 892 892 } 893 893 if (accept_q(self, self->params.x[i])){ 894 if (self->params. est_bck==1 && j==0) {894 if (self->params.has_bck==1 && j==0) { 895 895 a[i*nfunc+j] = 1.0/self->params.err[i]; 896 896 } else { … … 906 906 } 907 907 for (i_r=0; i_r<nr; i_r++){ 908 if (self->params. est_bck==1 && j==0) {908 if (self->params.has_bck==1 && j==0) { 909 909 a[(i_r+self->params.npoints)*nfunc+j] = 0.0; 910 910 } else { … … 1029 1029 {"set_slit_height", (PyCFunction)set_slit_height, METH_VARARGS, set_slit_height_doc}, 1030 1030 {"get_slit_height", (PyCFunction)get_slit_height, METH_VARARGS, get_slit_height_doc}, 1031 {"set_ est_bck", (PyCFunction)set_est_bck, METH_VARARGS, set_est_bck_doc},1032 {"get_ est_bck", (PyCFunction)get_est_bck, METH_VARARGS, get_est_bck_doc},1031 {"set_has_bck", (PyCFunction)set_has_bck, METH_VARARGS, set_has_bck_doc}, 1032 {"get_has_bck", (PyCFunction)get_has_bck, METH_VARARGS, get_has_bck_doc}, 1033 1033 {"get_nx", (PyCFunction)get_nx, METH_VARARGS, get_nx_doc}, 1034 1034 {"get_ny", (PyCFunction)get_ny, METH_VARARGS, get_ny_doc}, -
src/sas/sascalc/pr/c_extensions/invertor.c
rcb62bd5 r959eb01 20 20 pars->q_min = -1.0; 21 21 pars->q_max = -1.0; 22 pars-> est_bck = 0;22 pars->has_bck = 0; 23 23 } 24 24 … … 313 313 return sqrt(sum_r2/(2.0*sum)); 314 314 } 315 -
src/sas/sascalc/pr/c_extensions/invertor.h
rcb62bd5 r959eb01 27 27 double q_max; 28 28 /// Flag for whether or not to evalute a constant background while inverting 29 int est_bck;29 int has_bck; 30 30 /// Slit height in units of q [A-1] 31 31 double slit_height; -
src/sas/sascalc/pr/invertor.py
rcb62bd5 r45dffa69 121 121 self.q_min, self.q_max, 122 122 self.x, self.y, 123 self.err, self. est_bck,123 self.err, self.has_bck, 124 124 self.slit_height, self.slit_width) = state 125 125 … … 133 133 self.q_min, self.q_max, 134 134 self.x, self.y, 135 self.err, self. est_bck,135 self.err, self.has_bck, 136 136 self.slit_height, self.slit_width, 137 137 ) … … 175 175 elif name == 'slit_width': 176 176 return self.set_slit_width(value) 177 elif name == ' est_bck':177 elif name == 'has_bck': 178 178 if value == True: 179 return self.set_ est_bck(1)179 return self.set_has_bck(1) 180 180 elif value == False: 181 return self.set_ est_bck(0)181 return self.set_has_bck(0) 182 182 else: 183 raise ValueError, "Invertor: est_bck can only be True or False"183 raise ValueError, "Invertor: has_bck can only be True or False" 184 184 185 185 return Cinvertor.__setattr__(self, name, value) … … 220 220 elif name == 'slit_width': 221 221 return self.get_slit_width() 222 elif name == ' est_bck':223 value = self.get_ est_bck()222 elif name == 'has_bck': 223 value = self.get_has_bck() 224 224 if value == 1: 225 225 return True … … 248 248 invertor.y = self.y 249 249 invertor.err = self.err 250 invertor.est_bck = self.est_bck 251 invertor.background = self.background 250 invertor.has_bck = self.has_bck 252 251 invertor.slit_height = self.slit_height 253 252 invertor.slit_width = self.slit_width … … 291 290 """ 292 291 # Reset the background value before proceeding 293 # self.background = 0.0 294 if not self.est_bck: 295 self.y -= self.background 296 out, cov = self.lstsq(nfunc, nr=nr) 297 if not self.est_bck: 298 self.y += self.background 299 return out, cov 292 self.background = 0.0 293 return self.lstsq(nfunc, nr=nr) 300 294 301 295 def iq(self, out, q): … … 460 454 461 455 # If we need to fit the background, add a term 462 if self. est_bck == True:456 if self.has_bck == True: 463 457 nfunc_0 = nfunc 464 458 nfunc += 1 … … 506 500 507 501 # Keep a copy of the last output 508 if self.est_bck == False: 502 if self.has_bck == False: 503 self.background = 0 509 504 self.out = c 510 505 self.cov = err … … 658 653 file.write("#slit_width=%g\n" % self.slit_width) 659 654 file.write("#background=%g\n" % self.background) 660 if self. est_bck == True:655 if self.has_bck == True: 661 656 file.write("#has_bck=1\n") 662 657 else: … … 739 734 toks = line.split('=') 740 735 if int(toks[1]) == 1: 741 self. est_bck = True736 self.has_bck = True 742 737 else: 743 self. est_bck = False738 self.has_bck = False 744 739 745 740 # Now read in the parameters -
src/sas/sasgui/guiframe/local_perspectives/plotting/Plotter2D.py
r3e5648b r7432acb 361 361 if self.slicer.__class__.__name__ != "BoxSum": 362 362 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') 365 364 wx.EVT_MENU(self, wx_id, self._onEditSlicer) 366 365 slicerpop.AppendSeparator() … … 533 532 534 533 """ 535 # Clear current slicer534 ## Clear current slicer 536 535 if self.slicer is not None: 537 536 self.slicer.clear() 538 # Create a new slicer537 ## Create a new slicer 539 538 self.slicer_z += 1 540 539 self.slicer = slicer(self, self.subplot, zorder=self.slicer_z) 541 540 self.subplot.set_ylim(self.data2D.ymin, self.data2D.ymax) 542 541 self.subplot.set_xlim(self.data2D.xmin, self.data2D.xmax) 543 # Draw slicer542 ## Draw slicer 544 543 self.update() 545 544 self.slicer.update() … … 573 572 npt = math.floor(npt) 574 573 from sas.sascalc.dataloader.manipulations import CircularAverage 575 # compute the maximum radius of data2D574 ## compute the maximum radius of data2D 576 575 self.qmax = max(math.fabs(self.data2D.xmax), 577 576 math.fabs(self.data2D.xmin)) … … 579 578 math.fabs(self.data2D.ymin)) 580 579 self.radius = math.sqrt(math.pow(self.qmax, 2) + math.pow(self.ymax, 2)) 581 # 580 ##Compute beam width 582 581 bin_width = (self.qmax + self.qmax) / npt 583 # Create data1D circular average of data2D582 ## Create data1D circular average of data2D 584 583 Circle = CircularAverage(r_min=0, r_max=self.radius, 585 584 bin_width=bin_width) … … 600 599 new_plot.name = "Circ avg " + self.data2D.name 601 600 new_plot.source = self.data2D.source 602 # 601 #new_plot.info = self.data2D.info 603 602 new_plot.interactive = True 604 603 new_plot.detector = self.data2D.detector 605 604 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... 607 606 new_plot.xaxis("\\rm{Q}", "A^{-1}") 608 607 if hasattr(self.data2D, "scale") and \ … … 616 615 new_plot.id = "Circ avg " + self.data2D.name 617 616 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) 619 619 wx.PostEvent(self.parent, 620 620 NewPlotEvent(plot=new_plot, title=new_plot.name)) … … 630 630 """ 631 631 if self.slicer is not None: 632 from parameters_panel_slicerimport SlicerParameterPanel632 from SlicerParameters import SlicerParameterPanel 633 633 dialog = SlicerParameterPanel(self, -1, "Slicer Parameters") 634 634 dialog.set_slicer(self.slicer.__class__.__name__, … … 668 668 params = self.slicer.get_params() 669 669 ## Create a new panel to display results of summation of Data2D 670 from parameters_panel_boxsumimport SlicerPanel670 from slicerpanel import SlicerPanel 671 671 win = MDIFrame(self.parent, None, 'None', (100, 200)) 672 672 new_panel = SlicerPanel(parent=win, id=-1, … … 758 758 if default_name.count('.') > 0: 759 759 default_name = default_name.split('.')[0] 760 #default_name += "_out" 760 761 if self.parent is not None: 761 762 self.parent.show_data2d(data, default_name) 762 763 763 764 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) 766 766 self.graphApp.setDefaults(self.grid_on, self.legend_on, 767 767 self.xaxis_label, self.yaxis_label, -
src/sas/sasgui/guiframe/local_perspectives/plotting/SectorSlicer.py
r8de66b6 r7432acb 37 37 ## Absolute value of the Angle between the middle line and any side line 38 38 self.phi = math.pi / 12 39 # Binning base for log/lin binning40 self.bin_base = 041 39 ## Middle line 42 40 self.main_line = LineInteractor(self, self.base.subplot, color='blue', … … 153 151 phimin = -self.left_line.phi + self.main_line.theta 154 152 phimax = self.left_line.phi + self.main_line.theta 155 bin_base = self.bin_base156 153 if nbins is None: 157 154 nbins = 20 158 155 sect = SectorQ(r_min=0.0, r_max=radius, 159 156 phi_min=phimin + math.pi, 160 phi_max=phimax + math.pi, nbins=nbins , base=bin_base)157 phi_max=phimax + math.pi, nbins=nbins) 161 158 162 159 sector = sect(self.base.data2D) … … 242 239 params["Delta_Phi [deg]"] = math.fabs(self.left_line.phi * 180 / math.pi) 243 240 params["nbins"] = self.nbins 244 params["binning base"] = self.bin_base245 241 return params 246 242 … … 256 252 phi = math.fabs(params["Delta_Phi [deg]"] * math.pi / 180) 257 253 self.nbins = int(params["nbins"]) 258 self.bin_base = params["binning base"]259 254 self.main_line.theta = main 260 255 ## Reset the slicer parameters -
src/sas/sasgui/perspectives/calculator/model_editor.py
r07ec714 ra1b8fee 643 643 self.name_hsizer = None 644 644 self.name_tcl = None 645 self.overwrite_cb = None646 645 self.desc_sizer = None 647 646 self.desc_tcl = None … … 658 657 self.warning = "" 659 658 #This does not seem to be used anywhere so commenting out for now 660 # -- PDB 2/26/17 659 # -- PDB 2/26/17 661 660 #self._description = "New Plugin Model" 662 661 self.function_tcl = None … … 690 689 #title name [string] 691 690 name_txt = wx.StaticText(self, -1, 'Function Name : ') 692 self.overwrite_cb = wx.CheckBox(self, -1, "Overwrite existing plugin model of this name?", (10, 10))693 self.overwrite_cb.SetValue(False)694 self.overwrite_cb.SetToolTipString("Overwrite it if already exists?")695 wx.EVT_CHECKBOX(self, self.overwrite_cb.GetId(), self.on_over_cb)691 overwrite_cb = wx.CheckBox(self, -1, "Overwrite existing plugin model of this name?", (10, 10)) 692 overwrite_cb.SetValue(False) 693 overwrite_cb.SetToolTipString("Overwrite it if already exists?") 694 wx.EVT_CHECKBOX(self, overwrite_cb.GetId(), self.on_over_cb) 696 695 self.name_tcl = wx.TextCtrl(self, -1, size=(PANEL_WIDTH * 3 / 5, -1)) 697 696 self.name_tcl.Bind(wx.EVT_TEXT_ENTER, self.on_change_name) … … 701 700 self.name_tcl.SetToolTipString(hint_name) 702 701 self.name_hsizer.AddMany([(self.name_tcl, 0, wx.LEFT | wx.TOP, 0), 703 ( self.overwrite_cb, 0, wx.LEFT, 20)])702 (overwrite_cb, 0, wx.LEFT, 20)]) 704 703 self.name_sizer.AddMany([(name_txt, 0, wx.LEFT | wx.TOP, 10), 705 704 (self.name_hsizer, 0, … … 741 740 self.param_sizer.AddMany([(param_txt, 0, wx.LEFT, 10), 742 741 (self.param_tcl, 1, wx.EXPAND | wx.ALL, 10)]) 743 742 744 743 # Parameters with polydispersity 745 744 pd_param_txt = wx.StaticText(self, -1, 'Fit Parameters requiring ' + \ … … 756 755 self.pd_param_tcl.setDisplayLineNumbers(True) 757 756 self.pd_param_tcl.SetToolTipString(pd_param_tip) 758 757 759 758 self.param_sizer.AddMany([(pd_param_txt, 0, wx.LEFT, 10), 760 759 (self.pd_param_tcl, 1, wx.EXPAND | wx.ALL, 10)]) … … 996 995 info = 'Error' 997 996 color = 'red' 998 self.overwrite_cb.SetValue(True)999 self.overwrite_name = True1000 997 else: 1001 998 self._notes = result … … 1033 1030 if has_scipy: 1034 1031 lines.insert(0, 'import scipy') 1035 1036 # Think about 2D later 1032 1033 # Think about 2D later 1037 1034 #self.is_2d = func_str.count("#self.ndim = 2") 1038 1035 #line_2d = '' 1039 1036 #if self.is_2d: 1040 1037 # line_2d = CUSTOM_2D_TEMP.split('\n') 1041 1042 # Also think about test later 1038 1039 # Also think about test later 1043 1040 #line_test = TEST_TEMPLATE.split('\n') 1044 1041 #local_params = '' … … 1046 1043 spaces4 = ' '*4 1047 1044 spaces13 = ' '*13 1048 spaces16 = ' '*16 1045 spaces16 = ' '*16 1049 1046 param_names = [] # to store parameter names 1050 1047 has_scipy = func_str.count("scipy.") … … 1058 1055 out_f.write(line + '\n') 1059 1056 if line.count('#name'): 1060 out_f.write('name = "%s" \n' % name) 1057 out_f.write('name = "%s" \n' % name) 1061 1058 elif line.count('#title'): 1062 out_f.write('title = "User model for %s"\n' % name) 1059 out_f.write('title = "User model for %s"\n' % name) 1063 1060 elif line.count('#description'): 1064 out_f.write('description = "%s"\n' % desc_str) 1061 out_f.write('description = "%s"\n' % desc_str) 1065 1062 elif line.count('#parameters'): 1066 1063 out_f.write('parameters = [ \n') … … 1068 1065 p_line = param_line.lstrip().rstrip() 1069 1066 if p_line: 1070 pname, pvalue , desc= self.get_param_helper(p_line)1067 pname, pvalue = self.get_param_helper(p_line) 1071 1068 param_names.append(pname) 1072 out_f.write("%s['%s', '', %s, [-numpy.inf, numpy.inf], '', ' %s'],\n" % (spaces16, pname, pvalue, desc))1069 out_f.write("%s['%s', '', %s, [-numpy.inf, numpy.inf], '', ''],\n" % (spaces16, pname, pvalue)) 1073 1070 for param_line in pd_param_str.split('\n'): 1074 1071 p_line = param_line.lstrip().rstrip() 1075 1072 if p_line: 1076 pname, pvalue , desc= self.get_param_helper(p_line)1073 pname, pvalue = self.get_param_helper(p_line) 1077 1074 param_names.append(pname) 1078 out_f.write("%s['%s', '', %s, [-numpy.inf, numpy.inf], 'volume', ' %s'],\n" % (spaces16, pname, pvalue, desc))1075 out_f.write("%s['%s', '', %s, [-numpy.inf, numpy.inf], 'volume', ''],\n" % (spaces16, pname, pvalue)) 1079 1076 out_f.write('%s]\n' % spaces13) 1080 1077 1081 1078 # No form_volume or ER available in simple model editor 1082 1079 out_f.write('def form_volume(*arg): \n') … … 1085 1082 out_f.write('def ER(*arg): \n') 1086 1083 out_f.write(' return 1.0 \n') 1087 1084 1088 1085 # function to compute 1089 1086 out_f.write('\n') … … 1094 1091 for func_line in func_str.split('\n'): 1095 1092 out_f.write('%s%s\n' % (spaces4, func_line)) 1096 1093 1097 1094 Iqxy_string = 'return Iq(numpy.sqrt(x**2+y**2) ' 1098 1095 1099 1096 out_f.write('\n') 1100 1097 out_f.write('def Iqxy(x, y ') … … 1116 1113 items = line.split(";") 1117 1114 for item in items: 1118 name = item.split("=")[0].strip() 1119 description = "" 1115 name = item.split("=")[0].lstrip().rstrip() 1120 1116 try: 1121 value = item.split("=")[1].strip() 1122 if value.count("#"): 1123 # If line ends in a comment, remove it before parsing float 1124 index = value.index("#") 1125 description = value[(index + 1):].strip() 1126 value = value[:value.index("#")].strip() 1117 value = item.split("=")[1].lstrip().rstrip() 1127 1118 float(value) 1128 except ValueError:1119 except: 1129 1120 value = 1.0 # default 1130 1121 1131 return name, value , description1122 return name, value 1132 1123 1133 1124 def set_function_helper(self, line): … … 1213 1204 import numpy 1214 1205 1215 #name 1206 #name 1216 1207 1217 1208 #title … … 1219 1210 #description 1220 1211 1221 #parameters 1212 #parameters 1222 1213 1223 1214 """ -
src/sas/sasgui/perspectives/calculator/pyconsole.py
r4627657 r7432acb 37 37 Iqxy = model.evalDistribution([qx, qy]) 38 38 39 # check the model's unit tests run 40 from sasmodels.model_test import run_one 41 result = run_one(path) 39 result = """ 40 Iq(%s) = %s 41 Iqxy(%s, %s) = %s 42 """%(q, Iq, qx, qy, Iqxy) 42 43 43 44 return result … … 88 89 ok = wx.Button(self, wx.ID_OK, "OK") 89 90 90 # Mysterious constraint layouts from 91 # Mysterious constraint layouts from 91 92 # https://www.wxpython.org/docs/api/wx.lib.layoutf.Layoutf-class.html 92 93 lc = layoutf.Layoutf('t=t5#1;b=t5#2;l=l5#1;r=r5#1', (self,ok)) -
src/sas/sasgui/perspectives/fitting/basepage.py
ra1b8fee ra1b8fee 254 254 if not hasattr(self, "model_view"): 255 255 return 256 toggle_mode_on = self.model_view.IsEnabled() or self.data is None256 toggle_mode_on = self.model_view.IsEnabled() 257 257 if toggle_mode_on: 258 258 if self.enable2D and not check_data_validity(self.data): -
src/sas/sasgui/perspectives/fitting/fitting.py
r489f53a ra1b8fee 257 257 toks = os.path.splitext(label) 258 258 path = os.path.join(models.find_plugins_dir(), toks[0]) 259 message = "Are you sure you want to delete the file {}?".format(path)260 dlg = wx.MessageDialog(self.frame, message, '', wx.YES_NO | wx.ICON_QUESTION)261 if not dlg.ShowModal() == wx.ID_YES:262 return263 259 try: 264 260 for ext in ['.py', '.pyc']: 265 261 p_path = path + ext 266 if ext == '.pyc' and not os.path.isfile(path + ext):267 # If model is invalid, .pyc file may not exist as model has268 # never been compiled. Don't try and delete it269 continue270 262 os.remove(p_path) 271 263 self.update_custom_combo() … … 391 383 '(Re)Load all models present in user plugin_models folder') 392 384 wx.EVT_MENU(owner, wx_id, self.load_plugin_models) 393 385 394 386 def set_edit_menu_helper(self, owner=None, menu=None): 395 387 """ … … 1742 1734 @param unsmeared_error: data error, rescaled to unsmeared model 1743 1735 """ 1744 number_finite = np.count_nonzero(np.isfinite(y)) 1745 np.nan_to_num(y) 1746 new_plot = self.create_theory_1D(x, y, page_id, model, data, state, 1747 data_description=model.name, 1748 data_id=str(page_id) + " " + data.name) 1749 if unsmeared_model is not None: 1750 self.create_theory_1D(x, unsmeared_model, page_id, model, data, state, 1751 data_description=model.name + " unsmeared", 1752 data_id=str(page_id) + " " + data.name + " unsmeared") 1753 1754 if unsmeared_data is not None and unsmeared_error is not None: 1755 self.create_theory_1D(x, unsmeared_data, page_id, model, data, state, 1756 data_description="Data unsmeared", 1757 data_id="Data " + data.name + " unsmeared", 1758 dy=unsmeared_error) 1759 # Comment this out until we can get P*S models with correctly populated parameters 1760 if sq_model is not None and pq_model is not None: 1761 self.create_theory_1D(x, sq_model, page_id, model, data, state, 1762 data_description=model.name + " S(q)", 1763 data_id=str(page_id) + " " + data.name + " S(q)") 1764 self.create_theory_1D(x, pq_model, page_id, model, data, state, 1765 data_description=model.name + " P(q)", 1766 data_id=str(page_id) + " " + data.name + " P(q)") 1767 1768 current_pg = self.fit_panel.get_page_by_id(page_id) 1769 title = new_plot.title 1770 batch_on = self.fit_panel.get_page_by_id(page_id).batch_on 1771 if not batch_on: 1772 wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=str(title))) 1773 elif plot_result: 1774 top_data_id = self.fit_panel.get_page_by_id(page_id).data.id 1775 if data.id == top_data_id: 1776 wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=str(title))) 1777 caption = current_pg.window_caption 1778 self.page_finder[page_id].set_fit_tab_caption(caption=caption) 1779 1780 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, 1781 1774 fid=data.id) 1782 if toggle_mode_on:1783 wx.PostEvent(self.parent,1784 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", 1785 1778 action="Hide")) 1786 else:1787 if update_chisqr:1788 wx.PostEvent(current_pg,1789 Chi2UpdateEvent(output=self._cal_chisqr(1779 else: 1780 if update_chisqr: 1781 wx.PostEvent(current_pg, 1782 Chi2UpdateEvent(output=self._cal_chisqr( 1790 1783 data=data, 1791 1784 fid=fid, 1792 1785 weight=weight, 1793 page_id=page_id, 1794 index=index))) 1795 else: 1796 self._plot_residuals(page_id=page_id, data=data, fid=fid, 1797 index=index, weight=weight) 1798 1799 if not number_finite: 1800 logger.error("Using the present parameters the model does not return any finite value. ") 1801 msg = "Computing Error: Model did not return any finite value." 1802 wx.PostEvent(self.parent, StatusEvent(status = msg, info="error")) 1803 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 1804 1792 msg = "Computation completed!" 1805 if number_finite != y.size:1806 msg += ' PROBLEM: For some Q values the model returns non finite intensities!'1807 logger.error("For some Q values the model returns non finite intensities.")1808 1793 wx.PostEvent(self.parent, StatusEvent(status=msg, type="stop")) 1794 except: 1795 raise 1809 1796 1810 1797 def _calc_exception(self, etype, value, tb): … … 1831 1818 that can be plot. 1832 1819 """ 1833 number_finite = np.count_nonzero(np.isfinite(image))1834 1820 np.nan_to_num(image) 1835 1821 new_plot = Data2D(image=image, err_image=data.err_data) … … 1890 1876 self._plot_residuals(page_id=page_id, data=data, fid=fid, 1891 1877 index=index, weight=weight) 1892 1893 if not number_finite: 1894 logger.error("Using the present parameters the model does not return any finite value. ") 1895 msg = "Computing Error: Model did not return any finite value." 1896 wx.PostEvent(self.parent, StatusEvent(status = msg, info="error")) 1897 else: 1898 msg = "Computation completed!" 1899 if number_finite != image.size: 1900 msg += ' PROBLEM: For some Qx,Qy values the model returns non finite intensities!' 1901 logger.error("For some Qx,Qy values the model returns non finite intensities.") 1902 wx.PostEvent(self.parent, StatusEvent(status=msg, type="stop")) 1878 msg = "Computation completed!" 1879 wx.PostEvent(self.parent, StatusEvent(status=msg, type="stop")) 1903 1880 1904 1881 def _draw_model2D(self, model, page_id, qmin, -
src/sas/sasgui/perspectives/fitting/models.py
rb1c2011 ra1b8fee 156 156 try: 157 157 import compileall 158 compileall.compile_dir(dir=dir, ddir=dir, force= 0,158 compileall.compile_dir(dir=dir, ddir=dir, force=1, 159 159 quiet=report_problem) 160 160 except: … … 163 163 164 164 165 def _find _models():165 def _findModels(dir): 166 166 """ 167 167 Find custom models 168 168 """ 169 169 # List of plugin objects 170 dir ectory= find_plugins_dir()170 dir = find_plugins_dir() 171 171 # Go through files in plug-in directory 172 if not os.path.isdir(dir ectory):173 msg = "SasView couldn't locate Model plugin folder %r." % dir ectory172 if not os.path.isdir(dir): 173 msg = "SasView couldn't locate Model plugin folder %r." % dir 174 174 logger.warning(msg) 175 175 return {} 176 176 177 plugin_log("looking for models in: %s" % str(dir ectory))178 # compile_file(directory) #always recompile the folder plugin179 logger.info("plugin model dir: %s" % str(dir ectory))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)) 180 180 181 181 plugins = {} 182 for filename in os.listdir(dir ectory):182 for filename in os.listdir(dir): 183 183 name, ext = os.path.splitext(filename) 184 184 if ext == '.py' and not name == '__init__': 185 path = os.path.abspath(os.path.join(dir ectory, filename))185 path = os.path.abspath(os.path.join(dir, filename)) 186 186 try: 187 187 model = load_custom_model(path) 188 if not model.name.count(PLUGIN_NAME_BASE): 189 model.name = PLUGIN_NAME_BASE + model.name 188 model.name = PLUGIN_NAME_BASE + model.name 190 189 plugins[model.name] = model 191 190 except Exception: … … 194 193 plugin_log(msg) 195 194 logger.warning("Failed to load plugin %r. See %s for details" 196 % (path, PLUGIN_LOG))197 195 % (path, PLUGIN_LOG)) 196 198 197 return plugins 199 198 … … 265 264 temp = {} 266 265 if self.is_changed(): 267 return _find _models()266 return _findModels(dir) 268 267 logger.info("plugin model : %s" % str(temp)) 269 268 return temp … … 298 297 for name, plug in self.stored_plugins.iteritems(): 299 298 self.model_dictionary[name] = plug 300 299 301 300 self._get_multifunc_models() 302 301 … … 340 339 """ 341 340 self.plugins = [] 342 new_plugins = _find _models()341 new_plugins = _findModels(dir) 343 342 for name, plug in new_plugins.iteritems(): 344 343 for stored_name, stored_plug in self.stored_plugins.iteritems(): -
src/sas/sasgui/perspectives/pr/inversion_panel.py
rcb62bd5 r7432acb 70 70 self.rg_ctl = None 71 71 self.iq0_ctl = None 72 self.bck_value = None 73 self.bck_est_ctl = None 74 self.bck_man_ctl = None 75 self.est_bck = True 76 self.bck_input = None 72 self.bck_chk = None 77 73 self.bck_ctl = None 78 74 … … 316 312 # Read the panel's parameters 317 313 flag, alpha, dmax, nfunc, qmin, \ 318 qmax, height, width , bck= self._read_pars()314 qmax, height, width = self._read_pars() 319 315 320 316 state.nfunc = nfunc … … 330 326 331 327 # Background evaluation checkbox 332 state.estimate_bck = self.est_bck 333 state.bck_value = bck 328 state.estimate_bck = self.bck_chk.IsChecked() 334 329 335 330 # Estimates … … 376 371 self.plot_data.SetValue(str(state.file)) 377 372 378 # Background value 379 self.bck_est_ctl.SetValue(state.estimate_bck) 380 self.bck_man_ctl.SetValue(not state.estimate_bck) 381 if not state.estimate_bck: 382 self.bck_input.Enable() 383 self.bck_input.SetValue(str(state.bck_value)) 384 self.est_bck = state.estimate_bck 385 self.bck_value = state.bck_value 373 # Background evaluation checkbox 374 self.bck_chk.SetValue(state.estimate_bck) 386 375 387 376 # Estimates … … 442 431 wx.EXPAND | wx.LEFT | wx.RIGHT | wx.ADJUST_MINSIZE, 15) 443 432 444 radio_sizer = wx.GridBagSizer(5, 5) 445 446 self.bck_est_ctl = wx.RadioButton(self, -1, "Estimate background level", 447 name="estimate_bck", style=wx.RB_GROUP) 448 self.bck_man_ctl = wx.RadioButton(self, -1, "Input manual background level", 449 name="manual_bck") 450 451 self.bck_est_ctl.Bind(wx.EVT_RADIOBUTTON, self._on_bck_changed) 452 self.bck_man_ctl.Bind(wx.EVT_RADIOBUTTON, self._on_bck_changed) 453 454 radio_sizer.Add(self.bck_est_ctl, (0,0), (1,1), wx.LEFT | wx.EXPAND) 455 radio_sizer.Add(self.bck_man_ctl, (0,1), (1,1), wx.RIGHT | wx.EXPAND) 456 433 self.bck_chk = wx.CheckBox(self, -1, "Estimate background level") 434 hint_msg = "Check box to let the fit estimate " 435 hint_msg += "the constant background level." 436 self.bck_chk.SetToolTipString(hint_msg) 437 self.bck_chk.Bind(wx.EVT_CHECKBOX, self._on_pars_changed) 457 438 iy += 1 458 pars_sizer.Add( radio_sizer, (iy, 0), (1, 2),439 pars_sizer.Add(self.bck_chk, (iy, 0), (1, 2), 459 440 wx.LEFT | wx.EXPAND | wx.ADJUST_MINSIZE, 15) 460 461 background_label = wx.StaticText(self, -1, "Background: ")462 self.bck_input = PrTextCtrl(self, -1, style=wx.TE_PROCESS_ENTER,463 size=(60, 20), value="0.0")464 self.bck_input.Disable()465 self.bck_input.Bind(wx.EVT_TEXT, self._read_pars)466 background_units = wx.StaticText(self, -1, "[A^(-1)]", size=(55, 20))467 iy += 1468 469 background_sizer = wx.GridBagSizer(5, 5)470 471 background_sizer.Add(background_label, (0, 0), (1,1),472 wx.LEFT | wx.EXPAND | wx.ADJUST_MINSIZE, 23)473 background_sizer.Add(self.bck_input, (0, 1), (1,1),474 wx.LEFT | wx.ADJUST_MINSIZE, 5)475 background_sizer.Add(background_units, (0, 2), (1,1),476 wx.LEFT | wx.ADJUST_MINSIZE, 5)477 pars_sizer.Add(background_sizer, (iy, 0), (1, 2),478 wx.LEFT | wx.EXPAND | wx.ADJUST_MINSIZE, 15)479 480 441 boxsizer1.Add(pars_sizer, 0, wx.EXPAND) 481 442 vbox.Add(boxsizer1, (iy_vb, 0), (1, 1), … … 803 764 self._on_pars_changed() 804 765 805 def _on_bck_changed(self, evt=None):806 self.est_bck = self.bck_est_ctl.GetValue()807 if self.est_bck:808 self.bck_input.Disable()809 else:810 self.bck_input.Enable()811 812 766 def _on_pars_changed(self, evt=None): 813 767 """ … … 816 770 scenes. 817 771 """ 818 flag, alpha, dmax, nfunc, qmin, qmax, height, width, bck = self._read_pars() 772 flag, alpha, dmax, nfunc, qmin, qmax, height, width = self._read_pars() 773 has_bck = self.bck_chk.IsChecked() 819 774 820 775 # If the pars are valid, estimate alpha … … 828 783 d_max=dmax, 829 784 q_min=qmin, q_max=qmax, 830 est_bck=self.est_bck, 831 bck_val=bck, 785 bck=has_bck, 832 786 height=height, 833 787 width=width) … … 843 797 height = 0 844 798 width = 0 845 background = 0846 799 flag = True 847 800 # Read slit height … … 937 890 self.qmax_ctl.Refresh() 938 891 939 # Read background 940 if not self.est_bck: 941 try: 942 bck_str = self.bck_input.GetValue() 943 if len(bck_str.strip()) == 0: 944 background = 0.0 945 else: 946 background = float(bck_str) 947 self.bck_input.SetBackgroundColour(wx.WHITE) 948 except ValueError: 949 background = 0.0 950 self.bck_input.SetBackgroundColour("pink") 951 self.bck_input.Refresh() 952 953 return flag, alpha, dmax, nfunc, qmin, qmax, height, width, background 892 return flag, alpha, dmax, nfunc, qmin, qmax, height, width 954 893 955 894 def _on_explore(self, evt): … … 976 915 # Push it to the manager 977 916 978 flag, alpha, dmax, nfunc, qmin, qmax, height, width, bck = self._read_pars() 917 flag, alpha, dmax, nfunc, qmin, qmax, height, width = self._read_pars() 918 has_bck = self.bck_chk.IsChecked() 979 919 980 920 if flag: … … 988 928 d_max=dmax, 989 929 q_min=qmin, q_max=qmax, 990 est_bck=self.est_bck, 991 bck_val = bck, 930 bck=has_bck, 992 931 height=height, 993 932 width=width) -
src/sas/sasgui/perspectives/pr/inversion_state.py
ra0e6b1b r7432acb 36 36 ["qmin", "qmin"], 37 37 ["qmax", "qmax"], 38 ["estimate_bck", "estimate_bck"], 39 ["bck_value", "bck_value"]] 38 ["estimate_bck", "estimate_bck"]] 40 39 41 40 ## List of P(r) inversion outputs … … 63 62 self.estimate_bck = False 64 63 self.timestamp = time.time() 65 self.bck_value = 0.066 64 67 65 # Inversion parameters … … 111 109 state += "Timestamp: %s\n" % self.timestamp 112 110 state += "Estimate bck: %s\n" % str(self.estimate_bck) 113 state += "Bck Value: %s\n" % str(self.bck_value)114 111 state += "No. terms: %s\n" % str(self.nfunc) 115 112 state += "D_max: %s\n" % str(self.d_max) … … 299 296 self.coefficients.append(float(c)) 300 297 except: 301 # Bad data, skip. We will count the number of 302 # coefficients at the very end and deal with 298 # Bad data, skip. We will count the number of 299 # coefficients at the very end and deal with 303 300 # inconsistencies then. 304 301 pass … … 332 329 cov_row.append(float(c)) 333 330 except: 334 # Bad data, skip. We will count the number of 335 # coefficients at the very end and deal with 331 # Bad data, skip. We will count the number of 332 # coefficients at the very end and deal with 336 333 # inconsistencies then. 337 334 pass … … 464 461 tree = etree.parse(path, parser=etree.ETCompatXMLParser()) 465 462 # Check the format version number 466 # Specifying the namespace will take care of the file 467 #format version 463 # Specifying the namespace will take care of the file 464 #format version 468 465 root = tree.getroot() 469 466 -
src/sas/sasgui/perspectives/pr/media/pr_help.rst
r1abd19c r1221196 49 49 P(r) inversion requires that the background be perfectly subtracted. This is 50 50 often difficult to do well and thus many data sets will include a background. 51 For those cases, the user should check the "Estimate background level" option 52 and the module will do its best to estimate it. If you know the background value 53 for your data, select the "Input manual background level" option. Note that 54 this value will be treated as having 0 error. 51 For those cases, the user should check the "estimate background" box and the 52 module will do its best to estimate it. 55 53 56 54 The P(r) module is constantly computing in the background what the optimum -
src/sas/sasgui/perspectives/pr/pr.py
rcb62bd5 ra1b8fee 68 68 self.q_min = None 69 69 self.q_max = None 70 self.est_bck = False 71 self.bck_val = 0 70 self.has_bck = False 72 71 self.slit_height = 0 73 72 self.slit_width = 0 … … 829 828 self.control_panel.iq0 = pr.iq0(out) 830 829 self.control_panel.bck = pr.background 831 self.control_panel.bck_input.SetValue("{:.2g}".format(pr.background))832 830 833 831 # Show I(q) fit … … 909 907 910 908 def setup_plot_inversion(self, alpha, nfunc, d_max, q_min=None, q_max=None, 911 est_bck=False, bck_val=0, height=0, width=0):909 bck=False, height=0, width=0): 912 910 """ 913 911 Set up inversion from plotted data … … 918 916 self.q_min = q_min 919 917 self.q_max = q_max 920 self.est_bck = est_bck 921 self.bck_val = bck_val 918 self.has_bck = bck 922 919 self.slit_height = height 923 920 self.slit_width = width … … 933 930 def estimate_plot_inversion(self, alpha, nfunc, d_max, 934 931 q_min=None, q_max=None, 935 est_bck=False, bck_val=0, height=0, width=0):932 bck=False, height=0, width=0): 936 933 """ 937 934 Estimate parameters from plotted data … … 942 939 self.q_min = q_min 943 940 self.q_max = q_max 944 self.est_bck = est_bck 945 self.bck_val = bck_val 941 self.has_bck = bck 946 942 self.slit_height = height 947 943 self.slit_width = width … … 977 973 pr.x = self.current_plottable.x 978 974 pr.y = self.current_plottable.y 979 pr. est_bck = self.est_bck975 pr.has_bck = self.has_bck 980 976 pr.slit_height = self.slit_height 981 977 pr.slit_width = self.slit_width 982 pr.background = self.bck_val983 978 984 979 # Keep track of the plot window title to ensure that … … 1024 1019 self.q_min = q_min 1025 1020 self.q_max = q_max 1026 self. est_bck = bck1021 self.has_bck = bck 1027 1022 self.slit_height = height 1028 1023 self.slit_width = width … … 1047 1042 self.q_min = q_min 1048 1043 self.q_max = q_max 1049 self. est_bck = bck1044 self.has_bck = bck 1050 1045 self.slit_height = height 1051 1046 self.slit_width = width … … 1120 1115 pr.y = y 1121 1116 pr.err = err 1122 pr. est_bck = self.est_bck1117 pr.has_bck = self.has_bck 1123 1118 pr.slit_height = self.slit_height 1124 1119 pr.slit_width = self.slit_width -
test/pr_inversion/test/utest_invertor.py
rcb62bd5 raaf5e49 68 68 self.x_in[i] = 1.0*(i+1) 69 69 70 def test_ est_bck_flag(self):71 """ 72 Tests the est_bck flag operations73 """ 74 self.assertEqual(self.invertor. est_bck, False)75 self.invertor. est_bck=True76 self.assertEqual(self.invertor. est_bck, True)70 def test_has_bck_flag(self): 71 """ 72 Tests the has_bck flag operations 73 """ 74 self.assertEqual(self.invertor.has_bck, False) 75 self.invertor.has_bck=True 76 self.assertEqual(self.invertor.has_bck, True) 77 77 def doit_float(): 78 self.invertor. est_bck = 2.078 self.invertor.has_bck = 2.0 79 79 def doit_str(): 80 self.invertor. est_bck = 'a'80 self.invertor.has_bck = 'a' 81 81 82 82 self.assertRaises(ValueError, doit_float) -
test/sasdataloader/test/utest_averaging.py
re123eb9 r9a5097c 1 1 2 import unittest 2 3 import math 3 import os 4 import unittest 5 from scipy.stats import binned_statistic_2d 4 5 from sas.sascalc.dataloader.loader import Loader 6 from sas.sascalc.dataloader.manipulations import Ring, CircularAverage, SectorPhi, get_q,reader2D_converter 7 6 8 import numpy as np 7 8 9 import sas.sascalc.dataloader.data_info as data_info 9 from sas.sascalc.dataloader.loader import Loader10 from sas.sascalc.dataloader.manipulations import (Boxavg, Boxsum,11 CircularAverage, Ring,12 SectorPhi, SectorQ, SlabX,13 SlabY, get_q,14 reader2D_converter)15 16 10 17 11 class Averaging(unittest.TestCase): … … 19 13 Test averaging manipulations on a flat distribution 20 14 """ 21 22 15 def setUp(self): 23 16 """ … … 25 18 should return the predefined height of the distribution (1.0). 26 19 """ 27 x_0 = np.ones([100,100])28 dx_0 = np.ones([100, 29 20 x_0 = np.ones([100,100]) 21 dx_0 = np.ones([100,100]) 22 30 23 self.data = data_info.Data2D(data=x_0, err_data=dx_0) 31 24 detector = data_info.Detector() 32 detector.distance = 1000.0 # 33 detector.pixel_size.x = 1.0 #mm34 detector.pixel_size.y = 1.0 #mm35 25 detector.distance = 1000.0 #mm 26 detector.pixel_size.x = 1.0 #mm 27 detector.pixel_size.y = 1.0 #mm 28 36 29 # center in pixel position = (len(x_0)-1)/2 37 detector.beam_center.x = (len(x_0) - 1) / 2 #pixel number38 detector.beam_center.y = (len(x_0) - 1) / 2 #pixel number30 detector.beam_center.x = (len(x_0)-1)/2 #pixel number 31 detector.beam_center.y = (len(x_0)-1)/2 #pixel number 39 32 self.data.detector.append(detector) 40 33 41 34 source = data_info.Source() 42 source.wavelength = 10.0 #A35 source.wavelength = 10.0 #A 43 36 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. 47 39 self.qmin = get_q(1.0, 1.0, detector.distance, source.wavelength) 48 40 49 41 self.qmax = get_q(49.5, 49.5, detector.distance, source.wavelength) 50 42 51 43 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 =x61 self.data.y_bins =y44 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 62 54 self.data = reader2D_converter(self.data) 63 55 64 56 def test_ring_flat_distribution(self): 65 57 """ 66 58 Test ring averaging 67 59 """ 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, 70 62 center_y=self.data.detector[0].beam_center.y) 71 63 r.nbins_phi = 20 72 64 73 65 o = r(self.data) 74 66 for i in range(20): 75 67 self.assertEqual(o.y[i], 1.0) 76 68 77 69 def test_sectorphi_full(self): 78 70 """ 79 71 Test sector averaging 80 72 """ 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) 83 75 r.nbins_phi = 20 84 76 o = r(self.data) 85 77 for i in range(7): 86 78 self.assertEqual(o.y[i], 1.0) 87 79 80 88 81 def test_sectorphi_partial(self): 89 82 """ 90 83 """ 91 84 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, 93 86 phi_min=0, phi_max=phi_max) 94 87 self.assertEqual(r.phi_max, phi_max) … … 98 91 for i in range(17): 99 92 self.assertEqual(o.y[i], 1.0) 100 101 102 class DataInfoTests(unittest.TestCase): 103 93 94 95 96 class data_info_tests(unittest.TestCase): 97 104 98 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 109 101 def test_ring(self): 110 102 """ 111 103 Test ring averaging 112 104 """ 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, 115 107 center_y=self.data.detector[0].beam_center.y, 116 nbins =20)108 nbins = 20) 117 109 ##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 124 114 for i in range(r.nbins_phi - 1): 125 115 self.assertAlmostEqual(o.x[i + 1], answer.x[i], 4) 126 116 self.assertAlmostEqual(o.y[i + 1], answer.y[i], 4) 127 117 self.assertAlmostEqual(o.dy[i + 1], answer.dy[i], 4) 128 118 129 119 def test_circularavg(self): 130 120 """ 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') 143 131 for i in range(r.nbins_phi): 144 132 self.assertAlmostEqual(o.x[i], answer.x[i], 4) 145 133 self.assertAlmostEqual(o.y[i], answer.y[i], 4) 146 134 self.assertAlmostEqual(o.dy[i], answer.dy[i], 4) 147 135 148 136 def test_box(self): 149 137 """ … … 151 139 The test data was not generated by IGOR. 152 140 """ 153 141 from sas.sascalc.dataloader.manipulations import Boxsum, Boxavg 142 154 143 r = Boxsum(x_min=.01, x_max=.015, y_min=0.01, y_max=0.015) 155 144 s, ds, npoints = r(self.data) 156 145 self.assertAlmostEqual(s, 34.278990899999997, 4) 157 146 self.assertAlmostEqual(ds, 7.8007981835194293, 4) 158 self.assertAlmostEqual(npoints, 324.0000, 4) 159 147 self.assertAlmostEqual(npoints, 324.0000, 4) 148 160 149 r = Boxavg(x_min=.01, x_max=.015, y_min=0.01, y_max=0.015) 161 150 s, ds = r(self.data) 162 151 self.assertAlmostEqual(s, 0.10579935462962962, 4) 163 152 self.assertAlmostEqual(ds, 0.024076537603455028, 4) 164 153 165 154 def test_slabX(self): 166 155 """ … … 168 157 The test data was not generated by IGOR. 169 158 """ 170 171 r = SlabX(x_min=-.01, x_max=.01, y_min=-0.0002,172 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) 173 162 r.fold = False 174 163 o = r(self.data) 175 164 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 184 171 def test_slabY(self): 185 172 """ … … 187 174 The test data was not generated by IGOR. 188 175 """ 189 190 r = SlabY(x_min=.005, x_max=.01, y_min=-191 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) 192 179 r.fold = False 193 180 o = r(self.data) 194 181 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 203 188 def test_sectorphi_full(self): 204 189 """ … … 208 193 The test data was not generated by IGOR. 209 194 """ 210 195 from sas.sascalc.dataloader.manipulations import SectorPhi 196 import math 197 211 198 nbins = 19 212 199 phi_min = math.pi / (nbins + 1) 213 200 phi_max = math.pi * 2 - phi_min 214 201 215 202 r = SectorPhi(r_min=.005, 216 203 r_max=.01, … … 220 207 o = r(self.data) 221 208 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 230 215 def test_sectorphi_quarter(self): 231 216 """ … … 233 218 The test data was not generated by IGOR. 234 219 """ 235 236 r = SectorPhi(r_min=.005, r_max=.01, phi_min=0, phi_max=math.pi / 2.0)237 r.nbins_phi = 20238 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 248 233 def test_sectorq_full(self): 249 234 """ … … 251 236 The test data was not generated by IGOR. 252 237 """ 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 297 251 298 252 if __name__ == '__main__': -
test/sasguiframe/test/utest_manipulations.py
r959eb01 r959eb01 2 2 Unit tests for data manipulations 3 3 """ 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 6 import unittest 6 7 import math 8 import numpy as np 9 from sas.sascalc.dataloader.loader import Loader 10 from sas.sasgui.guiframe.dataFitting import Data1D, Data2D 11 from sas.sasgui.guiframe.dataFitting import Data1D as Theory1D 12 7 13 import 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 15 class data_info_tests(unittest.TestCase): 16 19 17 def setUp(self): 20 18 data = Loader().load("cansas1d.xml") 21 19 self.data = data[0] 22 20 23 21 def test_clone1D(self): 24 22 """ … … 26 24 """ 27 25 clone = self.data.clone_without_data() 28 26 29 27 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 30 class theory1d_tests(unittest.TestCase): 31 36 32 def setUp(self): 37 33 data = Loader().load("cansas1d.xml") 38 34 self.data = data[0] 39 35 40 36 def test_clone_theory1D(self): 41 37 """ 42 38 Test basic cloning 43 39 """ 44 theory = Data1D(x=[], y=[], dy=None)40 theory = Theory1D(x=[], y=[], dy=None) 45 41 theory.clone_without_data(clone=self.data) 46 42 theory.copy_from_datainfo(data1d=self.data) 47 43 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 51 46 for i in range(len(self.data.x)): 52 47 self.assertEqual(self.data.x[i], theory.x[i]) … … 54 49 self.assertEqual(self.data.dy[i], theory.dy[i]) 55 50 56 57 class ManipTests(unittest.TestCase): 58 51 class manip_tests(unittest.TestCase): 52 59 53 def setUp(self): 60 54 # Create two data sets to play with 61 55 x_0 = np.ones(5) 62 56 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) 67 61 self.data = Data1D(x_0, y_0, dy=dy_0) 68 62 69 63 x = self.data.x 70 64 y = np.ones(5) 71 65 dy = np.ones(5) 72 66 self.data2 = Data1D(x, y, dy=dy) 73 67 68 74 69 def test_load(self): 75 70 """ … … 78 73 # There should be 5 entries in the file 79 74 self.assertEqual(len(self.data.x), 5) 80 75 81 76 for i in range(5): 82 77 # 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 85 80 # 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 88 83 # 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 91 86 def test_add(self): 92 result = self.data2 +self.data87 result = self.data2+self.data 93 88 for i in range(5): 94 89 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 97 92 def test_sub(self): 98 result = self.data2 -self.data93 result = self.data2-self.data 99 94 for i in range(5): 100 95 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 103 98 def test_mul(self): 104 result = self.data2 *self.data99 result = self.data2*self.data 105 100 for i in range(5): 106 101 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 110 104 def test_div(self): 111 result = self.data2 /self.data105 result = self.data2/self.data 112 106 for i in range(5): 113 107 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 117 110 def test_radd(self): 118 result = self.data +3.0111 result = self.data+3.0 119 112 for i in range(5): 120 113 self.assertEqual(result.y[i], 5.0) 121 114 self.assertEqual(result.dy[i], 0.5) 122 123 result = 3.0 +self.data115 116 result = 3.0+self.data 124 117 for i in range(5): 125 118 self.assertEqual(result.y[i], 5.0) 126 119 self.assertEqual(result.dy[i], 0.5) 127 120 128 121 def test_rsub(self): 129 result = self.data -3.0122 result = self.data-3.0 130 123 for i in range(5): 131 124 self.assertEqual(result.y[i], -1.0) 132 125 self.assertEqual(result.dy[i], 0.5) 133 134 result = 3.0 -self.data126 127 result = 3.0-self.data 135 128 for i in range(5): 136 129 self.assertEqual(result.y[i], 1.0) 137 130 self.assertEqual(result.dy[i], 0.5) 138 131 139 132 def test_rmul(self): 140 result = self.data *3.0133 result = self.data*3.0 141 134 for i in range(5): 142 135 self.assertEqual(result.y[i], 6.0) 143 136 self.assertEqual(result.dy[i], 1.5) 144 145 result = 3.0 *self.data137 138 result = 3.0*self.data 146 139 for i in range(5): 147 140 self.assertEqual(result.y[i], 6.0) 148 141 self.assertEqual(result.dy[i], 1.5) 149 142 150 143 def test_rdiv(self): 151 result = self.data /4.0144 result = self.data/4.0 152 145 for i in range(5): 153 146 self.assertEqual(result.y[i], 0.5) 154 147 self.assertEqual(result.dy[i], 0.125) 155 156 result = 6.0 /self.data148 149 result = 6.0/self.data 157 150 for i in range(5): 158 151 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 154 class manip_2D(unittest.TestCase): 155 164 156 def setUp(self): 165 157 # 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) 168 160 qx_0 = np.arange(25) 169 161 qy_0 = np.arange(25) 170 162 mask_0 = np.zeros(25) 171 dqx_0 = np.arange(25) /100172 dqy_0 = np.arange(25) /100163 dqx_0 = np.arange(25)/100 164 dqy_0 = np.arange(25)/100 173 165 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, 176 168 dqx_data=dqx_0, dqy_data=dqy_0) 177 169 178 170 y = np.ones(25) 179 171 dy = np.ones(25) … … 182 174 mask = np.zeros(25) 183 175 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, 185 177 q_data=q, mask=mask) 186 178 179 187 180 def test_load(self): 188 181 """ … … 191 184 # There should be 5 entries in the file 192 185 self.assertEqual(np.size(self.data.data), 25) 193 186 194 187 for i in range(25): 195 188 # 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 198 191 # 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 201 194 def test_add(self): 202 result = self.data2 +self.data195 result = self.data2+self.data 203 196 for i in range(25): 204 197 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 207 200 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 209 231 for i in range(25): 210 232 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 245 236 for i in range(25): 246 237 self.assertEqual(result.data[i], 1.0) 247 238 self.assertEqual(result.err_data[i], 0.5) 248 239 249 240 def test_rmul(self): 250 result = self.data *3.0241 result = self.data*3.0 251 242 for i in range(25): 252 243 self.assertEqual(result.data[i], 6.0) 253 244 self.assertEqual(result.err_data[i], 1.5) 254 255 result = 3.0 *self.data245 246 result = 3.0*self.data 256 247 for i in range(25): 257 248 self.assertEqual(result.data[i], 6.0) 258 249 self.assertEqual(result.err_data[i], 1.5) 259 250 260 251 def test_rdiv(self): 261 result = self.data /4.0252 result = self.data/4.0 262 253 for i in range(25): 263 254 self.assertEqual(result.data[i], 0.5) 264 255 self.assertEqual(result.err_data[i], 0.125) 265 256 266 result = 6.0 /self.data257 result = 6.0/self.data 267 258 for i in range(25): 268 259 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 262 class extra_manip_2D(unittest.TestCase): 263 274 264 def setUp(self): 275 265 # 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) 278 268 qx_0 = np.arange(25) 279 269 qy_0 = np.arange(25) 280 270 mask_0 = np.zeros(25) 281 dqx_0 = np.arange(25) /100282 dqy_0 = np.arange(25) /100271 dqx_0 = np.arange(25)/100 272 dqy_0 = np.arange(25)/100 283 273 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, 286 276 dqx_data=dqx_0, dqy_data=dqy_0) 287 277 288 278 y = np.ones(25) 289 279 dy = np.ones(25) … … 292 282 mask = np.zeros(25) 293 283 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, 295 285 q_data=q, mask=mask) 296 286 287 297 288 def test_load(self): 298 289 """ … … 301 292 # There should be 5 entries in the file 302 293 self.assertEqual(np.size(self.data.data), 25) 303 294 304 295 for i in range(25): 305 296 # 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 308 299 # 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 311 302 def test_add(self): 312 result = self.data2 +self.data303 result = self.data2+self.data 313 304 for i in range(25): 314 305 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 317 308 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 319 339 for i in range(25): 320 340 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 355 344 for i in range(25): 356 345 self.assertEqual(result.data[i], 1.0) 357 346 self.assertEqual(result.err_data[i], 0.5) 358 347 359 348 def test_rmul(self): 360 result = self.data *3.0349 result = self.data*3.0 361 350 for i in range(25): 362 351 self.assertEqual(result.data[i], 6.0) 363 352 self.assertEqual(result.err_data[i], 1.5) 364 365 result = 3.0 *self.data353 354 result = 3.0*self.data 366 355 for i in range(25): 367 356 self.assertEqual(result.data[i], 6.0) 368 357 self.assertEqual(result.err_data[i], 1.5) 369 358 370 359 def test_rdiv(self): 371 result = self.data /4.0360 result = self.data/4.0 372 361 for i in range(25): 373 362 self.assertEqual(result.data[i], 0.5) 374 363 self.assertEqual(result.err_data[i], 0.125) 375 364 376 result = 6.0 /self.data365 result = 6.0/self.data 377 366 for i in range(25): 378 367 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) 381 369 382 370 if __name__ == '__main__': 383 371 unittest.main() 372
Note: See TracChangeset
for help on using the changeset viewer.