Changeset d5b5b71 in sasmodels
- Timestamp:
- Mar 1, 2017 7:29:46 AM (8 years ago)
- Branches:
- master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- b216880
- Parents:
- e09d1e0 (diff), e1aa129 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent. - Files:
-
- 1 added
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/sesans.py
rb397165 r94d13f1 14 14 import numpy as np # type: ignore 15 15 from numpy import pi, exp # type: ignore 16 from scipy.special import jv as besselj 17 #import direct_model.DataMixin as model 18 19 def make_q(q_max, Rmax): 20 r""" 21 Return a $q$ vector suitable for SESANS covering from $2\pi/ (10 R_{\max})$ 22 to $q_max$. This is the integration range of the Hankel transform; bigger range and 23 more points makes a better numerical integration. 24 Smaller q_min will increase reliable spin echo length range. 25 Rmax is the "radius" of the largest expected object and can be set elsewhere. 26 q_max is determined by the acceptance angle of the SESANS instrument. 16 from scipy.special import j0 17 18 class SesansTransform(object): 27 19 """ 28 from sas.sascalc.data_util.nxsunit import Converter 20 Spin-Echo SANS transform calculator. Similar to a resolution function, 21 the SesansTransform object takes I(q) for the set of *q_calc* values and 22 produces a transformed dataset 29 23 30 q_min = dq = 0.1 * 2*pi / Rmax31 return np.arange(q_min, 32 Converter(q_max[1])(q_max[0],33 units="1/A"),34 dq) 35 36 def make_all_q(data): 24 *SElength* (A) is the set of spin-echo lengths in the measured data. 25 26 *zaccept* (1/A) is the maximum acceptance of scattering vector in the spin 27 echo encoding dimension (for ToF: Q of min(R) and max(lam)). 28 29 *Rmax* (A) is the maximum size sensitivity; larger radius requires more 30 computation time. 37 31 """ 38 Return a $q$ vector suitable for calculating the total scattering cross section for 39 calculating the effect of finite acceptance angles on Time of Flight SESANS instruments. 40 If no acceptance is given, or unwanted (set "unwanted" flag in paramfile), no all_q vector is needed. 41 If the instrument has a rectangular acceptance, 2 all_q vectors are needed. 42 If the instrument has a circular acceptance, 1 all_q vector is needed 43 44 """ 45 if not data.has_no_finite_acceptance: 46 return [] 47 elif data.has_yz_acceptance(data): 48 # compute qx, qy 49 Qx, Qy = np.meshgrid(qx, qy) 50 return [Qx, Qy] 51 else: 52 # else only need q 53 # data.has_z_acceptance 54 return [q] 32 #: SElength from the data in the original data units; not used by transform 33 #: but the GUI uses it, so make sure that it is present. 34 q = None # type: np.ndarray 55 35 56 def transform(data, q_calc, Iq_calc, qmono, Iq_mono): 57 """ 58 Decides which transform type is to be used, based on the experiment data file contents (header) 59 (2016-03-19: currently controlled from parameters script) 60 nqmono is the number of q vectors to be used for the detector integration 61 """ 62 nqmono = len(qmono) 63 if nqmono == 0: 64 result = call_hankel(data, q_calc, Iq_calc) 65 elif nqmono == 1: 66 q = qmono[0] 67 result = call_HankelAccept(data, q_calc, Iq_calc, q, Iq_mono) 68 else: 69 Qx, Qy = [qmono[0], qmono[1]] 70 Qx = np.reshape(Qx, nqx, nqy) 71 Qy = np.reshape(Qy, nqx, nqy) 72 Iq_mono = np.reshape(Iq_mono, nqx, nqy) 73 qx = Qx[0, :] 74 qy = Qy[:, 0] 75 result = call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono) 36 #: q values to calculate when computing transform 37 q_calc = None # type: np.ndarray 76 38 77 return result 39 # transform arrays 40 _H = None # type: np.ndarray 41 _H0 = None # type: np.ndarray 78 42 79 def call_hankel(data, q_calc, Iq_calc): 80 return hankel((data.x, data.x_unit), 81 (data.lam, data.lam_unit), 82 (data.sample.thickness, 83 data.sample.thickness_unit), 84 q_calc, Iq_calc) 85 86 def call_HankelAccept(data, q_calc, Iq_calc, q_mono, Iq_mono): 87 return hankel(data.x, data.lam * 1e-9, 88 data.sample.thickness / 10, 89 q_calc, Iq_calc) 90 91 def call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono): 92 return hankel(data.x, data.y, data.lam * 1e-9, 93 data.sample.thickness / 10, 94 q_calc, Iq_calc) 95 96 def TotalScatter(model, parameters): #Work in progress!! 97 # Calls a model with existing model parameters already in place, then integrate the product of q and I(q) from 0 to (4*pi/lambda) 98 allq = np.linspace(0,4*pi/wavelength,1000) 99 allIq = 1 100 integral = allq*allIq 101 43 def __init__(self, z, SElength, zaccept, Rmax): 44 # type: (np.ndarray, float, float) -> None 45 #import logging; logging.info("creating SESANS transform") 46 self.q = z 47 self._set_hankel(SElength, zaccept, Rmax) 102 48 49 def apply(self, Iq): 50 # tye: (np.ndarray) -> np.ndarray 51 G0 = np.dot(self._H0, Iq) 52 G = np.dot(self._H.T, Iq) 53 P = G - G0 54 return P 103 55 104 def Cosine2D(wavelength, magfield, thickness, qy, qz, Iqy, Iqz, modelname): #Work in progress!! Needs to call model still 105 #============================================================================== 106 # 2D Cosine Transform if "wavelength" is a vector 107 #============================================================================== 108 #allq is the q-space needed to create the total scattering cross-section 56 def _set_hankel(self, SElength, zaccept, Rmax): 57 # type: (np.ndarray, float, float) -> None 58 # Force float32 arrays, otherwise run into memory problems on some machines 59 SElength = np.asarray(SElength, dtype='float32') 109 60 110 Gprime = np.zeros_like(wavelength, 'd') 111 s = np.zeros_like(wavelength, 'd') 112 sd = np.zeros_like(wavelength, 'd') 113 Gprime = np.zeros_like(wavelength, 'd') 114 f = np.zeros_like(wavelength, 'd') 115 for i, wavelength_i in enumerate(wavelength): 116 z = magfield*wavelength_i 117 allq=np.linspace() #for calculating the Q-range of the scattering power integral 118 allIq=np.linspace() # This is the model applied to the allq q-space. Needs to refference the model somehow 119 alldq = (allq[1]-allq[0])*1e10 120 sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) 121 s[i]=1-exp(-sigma) 122 for j, Iqy_j, qy_j in enumerate(qy): 123 for k, Iqz_k, qz_k in enumerate(qz): 124 Iq = np.sqrt(Iqy_j^2+Iqz_k^2) 125 q = np.sqrt(qy_j^2 + qz_k^2) 126 Gintegral = Iq*cos(z*Qz_k) 127 Gprime[i] += Gintegral 128 # sigma = wavelength^2*thickness/2/pi* allq[i]*allIq[i] 129 # s[i] += 1-exp(Totalscatter(modelname)*thickness) 130 # For now, work with standard 2-phase scatter 61 #Rmax = #value in text box somewhere in FitPage? 62 q_max = 2*pi / (SElength[1] - SElength[0]) 63 q_min = 0.1 * 2*pi / (np.size(SElength) * SElength[-1]) 64 q = np.arange(q_min, q_max, q_min, dtype='float32') 65 dq = q_min 131 66 67 H0 = np.float32(dq/(2*pi)) * q 132 68 133 sd[i] += Iq134 f[i] = 1-s[i]+sd[i]135 P[i] = (1-sd[i]/f[i])+1/f[i]*Gprime[i]69 repq = np.tile(q, (SElength.size, 1)).T 70 repSE = np.tile(SElength, (q.size, 1)) 71 H = np.float32(dq/(2*pi)) * j0(repSE*repq) * repq 136 72 137 138 139 140 def HankelAccept(wavelength, magfield, thickness, q, Iq, theta, modelname): 141 #============================================================================== 142 # HankelTransform with fixed circular acceptance angle (circular aperture) for Time of Flight SESANS 143 #============================================================================== 144 #acceptq is the q-space needed to create limited acceptance effect 145 SElength= wavelength*magfield 146 G = np.zeros_like(SElength, 'd') 147 threshold=2*pi*theta/wavelength 148 for i, SElength_i in enumerate(SElength): 149 allq=np.linspace() #for calculating the Q-range of the scattering power integral 150 allIq=np.linspace() # This is the model applied to the allq q-space. Needs to refference the model somehow 151 alldq = (allq[1]-allq[0])*1e10 152 sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) 153 s[i]=1-exp(-sigma) 154 155 dq = (q[1]-q[0])*1e10 156 a = (x<threshold) 157 acceptq = a*q 158 acceptIq = a*Iq 159 160 G[i] = np.sum(besselj(0, acceptq*SElength_i)*acceptIq*acceptq*dq) 161 162 # G[i]=np.sum(integral) 163 164 G *= dq*1e10*2*pi 165 166 P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0])) 167 168 def hankel(SElength, wavelength, thickness, q, Iq): 169 r""" 170 Compute the expected SESANS polarization for a given SANS pattern. 171 172 Uses the hankel transform followed by the exponential. The values for *zz* 173 (or spin echo length, or delta), wavelength and sample thickness should 174 come from the dataset. $q$ should be chosen such that the oscillations 175 in $I(q)$ are well sampled (e.g., $5 \cdot 2 \pi/d_{\max}$). 176 177 *SElength* [A] is the set of $z$ points at which to compute the 178 Hankel transform 179 180 *wavelength* [m] is the wavelength of each individual point *zz* 181 182 *thickness* [cm] is the sample thickness. 183 184 *q* [A$^{-1}$] is the set of $q$ points at which the model has been 185 computed. These should be equally spaced. 186 187 *I* [cm$^{-1}$] is the value of the SANS model at *q* 188 """ 189 190 from sas.sascalc.data_util.nxsunit import Converter 191 wavelength = Converter(wavelength[1])(wavelength[0],"A") 192 thickness = Converter(thickness[1])(thickness[0],"A") 193 Iq = Converter("1/cm")(Iq,"1/A") # All models default to inverse centimeters 194 SElength = Converter(SElength[1])(SElength[0],"A") 195 196 G = np.zeros_like(SElength, 'd') 197 #============================================================================== 198 # Hankel Transform method if "wavelength" is a scalar; mono-chromatic SESANS 199 #============================================================================== 200 for i, SElength_i in enumerate(SElength): 201 integral = besselj(0, q*SElength_i)*Iq*q 202 G[i] = np.sum(integral) 203 G0 = np.sum(Iq*q) 204 205 # [m^-1] step size in q, needed for integration 206 dq = (q[1]-q[0]) 207 208 # integration step, convert q into [m**-1] and 2 pi circle integration 209 G *= dq*2*pi 210 G0 = np.sum(Iq*q)*dq*2*np.pi 211 212 P = exp(thickness*wavelength**2/(4*pi**2)*(G-G0)) 213 214 return P 73 self.q_calc = q 74 self._H, self._H0 = H, H0 -
sasmodels/compare_many.py
r424fe00 r5124c969 106 106 header = ('\n"Model","%s","Count","%d","Dimension","%s"' 107 107 % (name, N, "2D" if is_2d else "1D")) 108 if not mono: header += ',"Cutoff",%g'%(cutoff,) 108 if not mono: 109 header += ',"Cutoff",%g'%(cutoff,) 109 110 print(header) 110 111 … … 161 162 max_diff = [0] 162 163 for k in range(N): 163 print(" %s %d"%(name, k), file=sys.stderr)164 print("Model %s %d"%(name, k+1), file=sys.stderr) 164 165 seed = np.random.randint(1e6) 165 166 pars_i = randomize_pars(model_info, pars, seed) 166 167 constrain_pars(model_info, pars_i) 167 constrain_new_to_old(model_info, pars_i) 168 if 'sasview' in (base, comp): 169 constrain_new_to_old(model_info, pars_i) 168 170 if mono: 169 171 pars_i = suppress_pd(pars_i) … … 204 206 print("""\ 205 207 206 MODEL is the model name of the model or "all" for all the models 207 in alphabetical order. 208 MODEL is the model name of the model or one of the model types listed in 209 sasmodels.core.list_models (all, py, c, double, single, opencl, 1d, 2d, 210 nonmagnetic, magnetic). Model types can be combined, such as 2d+single. 208 211 209 212 COUNT is the number of randomly generated parameter sets to try. A value … … 237 240 return 238 241 239 model = argv[0] 240 if not (model in MODELS) and (model != "all"): 241 print('Bad model %s. Use "all" or one of:'%model) 242 target = argv[0] 243 try: 244 model_list = [target] if target in MODELS else core.list_models(target) 245 except ValueError: 246 print('Bad model %s. Use model type or one of:'%model) 242 247 print_models() 248 print('model types: all, py, c, double, single, opencl, 1d, 2d, nonmagnetic, magnetic') 243 249 return 244 250 try: … … 258 264 data, index = make_data({'qmax':1.0, 'is2d':is2D, 'nq':Nq, 'res':0., 259 265 'accuracy': 'Low', 'view':'log', 'zero': False}) 260 model_list = [model] if model != "all" else MODELS261 266 for model in model_list: 262 267 compare_instance(model, data, index, N=count, mono=mono, -
sasmodels/core.py
r52e9a45 r5124c969 69 69 * magnetic: models with an sld 70 70 * nommagnetic: models without an sld 71 """ 72 if kind and kind not in KINDS: 71 72 For multiple conditions, combine with plus. For example, *c+single+2d* 73 would return all oriented models implemented in C which can be computed 74 accurately with single precision arithmetic. 75 """ 76 if kind and any(k not in KINDS for k in kind.split('+')): 73 77 raise ValueError("kind not in " + ", ".join(KINDS)) 74 78 files = sorted(glob(joinpath(generate.MODEL_PATH, "[a-zA-Z]*.py"))) 75 79 available_models = [basename(f)[:-3] for f in files] 76 selected = [name for name in available_models if _matches(name, kind)] 80 if kind and '+' in kind: 81 all_kinds = kind.split('+') 82 condition = lambda name: all(_matches(name, k) for k in all_kinds) 83 else: 84 condition = lambda name: _matches(name, kind) 85 selected = [name for name in available_models if condition(name)] 77 86 78 87 return selected -
sasmodels/model_test.py
r479d0f3 re09d1e0 97 97 is_py = callable(model_info.Iq) 98 98 99 # Some OpenCL drivers seem to be flaky, and are not producing the 100 # expected result. Since we don't have known test values yet for 101 # all of our models, we are instead going to compare the results 102 # for the 'smoke test' (that is, evaluation at q=0.1 for the default 103 # parameters just to see that the model runs to completion) between 104 # the OpenCL and the DLL. To do this, we define a 'stash' which is 105 # shared between OpenCL and DLL tests. This is just a list. If the 106 # list is empty (which it will be when DLL runs, if the DLL runs 107 # first), then the results are appended to the list. If the list 108 # is not empty (which it will be when OpenCL runs second), the results 109 # are compared to the results stored in the first element of the list. 110 # This is a horrible stateful hack which only makes sense because the 111 # test suite is thrown away after being run once. 112 stash = [] 113 99 114 if is_py: # kernel implemented in python 100 115 test_name = "Model: %s, Kernel: python"%model_name … … 103 118 test_method_name, 104 119 platform="dll", # so that 105 dtype="double") 120 dtype="double", 121 stash=stash) 106 122 suite.addTest(test) 107 123 else: # kernel implemented in C 124 125 # test using dll if desired 126 if 'dll' in loaders or not core.HAVE_OPENCL: 127 test_name = "Model: %s, Kernel: dll"%model_name 128 test_method_name = "test_%s_dll" % model_info.id 129 test = ModelTestCase(test_name, model_info, 130 test_method_name, 131 platform="dll", 132 dtype="double", 133 stash=stash) 134 suite.addTest(test) 135 108 136 # test using opencl if desired and available 109 137 if 'opencl' in loaders and core.HAVE_OPENCL: … … 116 144 test = ModelTestCase(test_name, model_info, 117 145 test_method_name, 118 platform="ocl", dtype=None) 146 platform="ocl", dtype=None, 147 stash=stash) 119 148 #print("defining", test_name) 120 suite.addTest(test)121 122 # test using dll if desired123 if 'dll' in loaders or not core.HAVE_OPENCL:124 test_name = "Model: %s, Kernel: dll"%model_name125 test_method_name = "test_%s_dll" % model_info.id126 test = ModelTestCase(test_name, model_info,127 test_method_name,128 platform="dll",129 dtype="double")130 149 suite.addTest(test) 131 150 … … 144 163 """ 145 164 def __init__(self, test_name, model_info, test_method_name, 146 platform, dtype ):147 # type: (str, ModelInfo, str, str, DType ) -> None165 platform, dtype, stash): 166 # type: (str, ModelInfo, str, str, DType, List[Any]) -> None 148 167 self.test_name = test_name 149 168 self.info = model_info 150 169 self.platform = platform 151 170 self.dtype = dtype 171 self.stash = stash # container for the results of the first run 152 172 153 173 setattr(self, test_method_name, self.run_all) … … 174 194 ] 175 195 176 tests = s elf.info.tests196 tests = smoke_tests + self.info.tests 177 197 try: 178 198 model = build_model(self.info, dtype=self.dtype, 179 199 platform=self.platform) 180 for test in smoke_tests + tests: 181 self.run_one(model, test) 182 183 if not tests and self.platform == "dll": 184 ## Uncomment the following to make forgetting the test 185 ## values an error. Only do so for the "dll" tests 186 ## to reduce noise from both opencl and dll, and because 187 ## python kernels use platform="dll". 188 #raise Exception("No test cases provided") 189 pass 200 results = [self.run_one(model, test) for test in tests] 201 if self.stash: 202 for test, target, actual in zip(tests, self.stash[0], results): 203 assert np.all(abs(target-actual)<2e-5*abs(actual)),\ 204 "expected %s but got %s for %s"%(target, actual, test[0]) 205 else: 206 self.stash.append(results) 207 208 # Check for missing tests. Only do so for the "dll" tests 209 # to reduce noise from both opencl and dll, and because 210 # python kernels use platform="dll". 211 if self.platform == "dll": 212 missing = [] 213 ## Uncomment the following to require test cases 214 #missing = self._find_missing_tests() 215 if missing: 216 raise ValueError("Missing tests for "+", ".join(missing)) 190 217 191 218 except: 192 219 annotate_exception(self.test_name) 193 220 raise 221 222 def _find_missing_tests(self): 223 # type: () -> None 224 """make sure there are 1D, 2D, ER and VR tests as appropriate""" 225 model_has_VR = callable(self.info.VR) 226 model_has_ER = callable(self.info.ER) 227 model_has_1D = True 228 model_has_2D = any(p.type == 'orientation' 229 for p in self.info.parameters.kernel_parameters) 230 231 # Lists of tests that have a result that is not None 232 single = [test for test in self.info.tests 233 if not isinstance(test[2], list) and test[2] is not None] 234 tests_has_VR = any(test[1] == 'VR' for test in single) 235 tests_has_ER = any(test[1] == 'ER' for test in single) 236 tests_has_1D_single = any(isinstance(test[1], float) for test in single) 237 tests_has_2D_single = any(isinstance(test[1], tuple) for test in single) 238 239 multiple = [test for test in self.info.tests 240 if isinstance(test[2], list) 241 and not all(result is None for result in test[2])] 242 tests_has_1D_multiple = any(isinstance(test[1][0], float) 243 for test in multiple) 244 tests_has_2D_multiple = any(isinstance(test[1][0], tuple) 245 for test in multiple) 246 247 missing = [] 248 if model_has_VR and not tests_has_VR: 249 missing.append("VR") 250 if model_has_ER and not tests_has_ER: 251 missing.append("ER") 252 if model_has_1D and not (tests_has_1D_single or tests_has_1D_multiple): 253 missing.append("1D") 254 if model_has_2D and not (tests_has_2D_single or tests_has_2D_multiple): 255 missing.append("2D") 256 257 return missing 194 258 195 259 def run_one(self, model, test): … … 207 271 208 272 if x[0] == 'ER': 209 actual = [call_ER(model.info, pars)]273 actual = np.array([call_ER(model.info, pars)]) 210 274 elif x[0] == 'VR': 211 actual = [call_VR(model.info, pars)]275 actual = np.array([call_VR(model.info, pars)]) 212 276 elif isinstance(x[0], tuple): 213 277 qx, qy = zip(*x) … … 238 302 'f(%s); expected:%s; actual:%s' 239 303 % (xi, yi, actual_yi)) 304 return actual 240 305 241 306 return ModelTestCase -
sasmodels/modelinfo.py
r85fe7f8 r5124c969 730 730 info.docs = kernel_module.__doc__ 731 731 info.category = getattr(kernel_module, 'category', None) 732 info.single = getattr(kernel_module, 'single', True)733 info.opencl = getattr(kernel_module, 'opencl', True)734 732 info.structure_factor = getattr(kernel_module, 'structure_factor', False) 735 733 info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) … … 745 743 info.profile = getattr(kernel_module, 'profile', None) # type: ignore 746 744 info.sesans = getattr(kernel_module, 'sesans', None) # type: ignore 745 # Default single and opencl to True for C models. Python models have callable Iq. 746 info.opencl = getattr(kernel_module, 'opencl', not callable(info.Iq)) 747 info.single = getattr(kernel_module, 'single', not callable(info.Iq)) 747 748 748 749 # multiplicity info
Note: See TracChangeset
for help on using the changeset viewer.