Changes in / [81bb668:002adb6] in sasmodels
- Files:
-
- 2 added
- 26 edited
Legend:
- Unmodified
- Added
- Removed
-
doc/guide/index.rst
r19dcb933 rbb6f0f3 1 ********************** 2 SAS Model Organization 3 ********************** 1 ********** 2 SAS Models 3 ********** 4 5 Small angle X-ray and Neutron (SAXS and SANS) scattering examines the 6 scattering patterns produced by a beam travelling through the sample 7 and scattering at low angles. The scattering is computed as a function 8 of $q_x$ and $q_y$, which for a given beam wavelength corresponds to 9 particular scattering angles. Each pixel on the detector corresponds to 10 a different scattering angle. If the sample is unoriented, the scattering 11 pattern will appear as rings on the detector. In this case, a circular 12 average can be taken with 1-dimension data at $q = \surd (q_x^2 + q_y^2)$ 13 compared to the orientationally averaged SAS scattering pattern. 4 14 5 15 Models have certain features in common. -
sasmodels/convert.py
rd15a908 r5054e80 13 13 'gauss_lorentz_gel', 14 14 'be_polyelectrolyte', 15 'correlation_length', 15 16 ] 16 17 -
sasmodels/core.py
reafc9fa rd18582e 73 73 return True 74 74 75 def load_model(model_definition, dtype= "single", platform="ocl"):75 def load_model(model_definition, dtype=None, platform="ocl"): 76 76 """ 77 77 Prepare the model for the default execution platform. … … 87 87 for the calculation. Any valid numpy single or double precision identifier 88 88 is valid, such as 'single', 'f', 'f32', or np.float32 for single, or 89 'double', 'd', 'f64' and np.float64 for double. 89 'double', 'd', 'f64' and np.float64 for double. If *None*, then use 90 'single' unless the model defines single=False. 90 91 91 92 *platform* should be "dll" to force the dll to be used for C models, … … 94 95 if isstr(model_definition): 95 96 model_definition = load_model_definition(model_definition) 97 if dtype is None: 98 dtype = 'single' if getattr(model_definition, 'single', True) else 'double' 96 99 source, info = generate.make(model_definition) 97 100 if callable(info.get('Iq', None)): -
sasmodels/data.py
r5c962df rd18582e 242 242 243 243 244 def empty_data1D(q, resolution=0.0 5):244 def empty_data1D(q, resolution=0.0): 245 245 """ 246 246 Create empty 1D data using the given *q* as the x value. … … 252 252 #dIq = np.sqrt(Iq) 253 253 Iq, dIq = None, None 254 q = np.asarray(q) 254 255 data = Data1D(q, Iq, dx=resolution * q, dy=dIq) 255 256 data.filename = "fake data" … … 257 258 258 259 259 def empty_data2D(qx, qy=None, resolution=0.0 5):260 def empty_data2D(qx, qy=None, resolution=0.0): 260 261 """ 261 262 Create empty 2D data using the given mesh. … … 267 268 if qy is None: 268 269 qy = qx 270 qx, qy = np.asarray(qx), np.asarray(qy) 269 271 # 5% dQ/Q resolution 270 272 Qx, Qy = np.meshgrid(qx, qy) -
sasmodels/direct_model.py
reafc9fa rd18582e 234 234 235 235 model_definition = load_model_definition(model_name) 236 model = load_model(model_definition , dtype='single')236 model = load_model(model_definition) 237 237 calculator = DirectModel(data, model) 238 238 pars = dict((k, float(v)) -
sasmodels/kernel_template.c
r9c79c32 rcaf768d 16 16 using namespace std; 17 17 #if defined(_MSC_VER) 18 # define kernel extern "C" __declspec( dllexport ) 18 #include <float.h> 19 #define kernel extern "C" __declspec( dllexport ) 19 20 inline double trunc(double x) { return x>=0?floor(x):-floor(-x); } 20 inline double fmin(double x, double y) { return x>y ? y : x; } 21 inline double fmax(double x, double y) { return x<y ? y : x; } 21 inline double fmin(double x, double y) { return x>y ? y : x; } 22 inline double fmax(double x, double y) { return x<y ? y : x; } 23 inline double isnan(double x) { return _isnan(x); } 22 24 #else 23 #define kernel extern "C"25 #define kernel extern "C" 24 26 #endif 25 27 inline void SINCOS(double angle, double &svar, double &cvar) { svar=sin(angle); cvar=cos(angle); } -
sasmodels/kernelcl.py
reafc9fa re6a5556 172 172 #self.data_boundary = max(d.min_data_type_align_size 173 173 # for d in self.context.devices) 174 self.queues = [cl.CommandQueue( self.context, d)175 for d in self.context.devices]174 self.queues = [cl.CommandQueue(context, context.devices[0]) 175 for context in self.context] 176 176 self.compiled = {} 177 177 … … 181 181 """ 182 182 dtype = generate.F32 if dtype == 'fast' else np.dtype(dtype) 183 return all(has_type(d, dtype) for d in self.context.devices) 183 return any(has_type(d, dtype) 184 for context in self.context 185 for d in context.devices) 186 187 def get_queue(self, dtype): 188 """ 189 Return a command queue for the kernels of type dtype. 190 """ 191 for context, queue in zip(self.context, self.queues): 192 if all(has_type(d, dtype) for d in context.devices): 193 return queue 194 195 def get_context(self, dtype): 196 """ 197 Return a OpenCL context for the kernels of type dtype. 198 """ 199 for context, queue in zip(self.context, self.queues): 200 if all(has_type(d, dtype) for d in context.devices): 201 return context 184 202 185 203 def _create_some_context(self): … … 190 208 """ 191 209 try: 192 self.context = cl.create_some_context(interactive=False)210 self.context = [cl.create_some_context(interactive=False)] 193 211 except Exception as exc: 194 212 warnings.warn(str(exc)) … … 204 222 #print("compiling",name) 205 223 dtype = np.dtype(dtype) 206 program = compile_model(self. context, source, dtype, fast)224 program = compile_model(self.get_context(dtype), source, dtype, fast) 207 225 self.compiled[key] = program 208 226 return self.compiled[key] … … 218 236 def _get_default_context(): 219 237 """ 220 Get an OpenCL context, preferring GPU over CPU. 221 """ 222 default = None 238 Get an OpenCL context, preferring GPU over CPU, and preferring Intel 239 drivers over AMD drivers. 240 """ 241 # Note: on mobile devices there is automatic clock scaling if either the 242 # CPU or the GPU is underutilized; probably doesn't affect us, but we if 243 # it did, it would mean that putting a busy loop on the CPU while the GPU 244 # is running may increase throughput. 245 # 246 # Macbook pro, base install: 247 # {'Apple': [Intel CPU, NVIDIA GPU]} 248 # Macbook pro, base install: 249 # {'Apple': [Intel CPU, Intel GPU]} 250 # 2 x nvidia 295 with Intel and NVIDIA opencl drivers installed 251 # {'Intel': [CPU], 'NVIDIA': [GPU, GPU, GPU, GPU]} 252 gpu, cpu = None, None 223 253 for platform in cl.get_platforms(): 254 # AMD provides a much weaker CPU driver than Intel/Apple, so avoid it. 255 # If someone has bothered to install the AMD/NVIDIA drivers, prefer them over the integrated 256 # graphics driver that may have been supplied with the CPU chipset. 257 preferred_cpu = platform.vendor.startswith('Intel') or platform.vendor.startswith('Apple') 258 preferred_gpu = platform.vendor.startswith('Advanced') or platform.vendor.startswith('NVIDIA') 224 259 for device in platform.get_devices(): 225 260 if device.type == cl.device_type.GPU: 226 return cl.Context([device]) 227 if default is None: 228 default = device 229 230 if not default: 231 raise RuntimeError("OpenCL device not found") 232 233 return cl.Context([default]) 261 # If the existing type is not GPU then it will be CUSTOM or ACCELERATOR, 262 # so don't override it. 263 if gpu is None or (preferred_gpu and gpu.type == cl.device_type.GPU): 264 gpu = device 265 elif device.type == cl.device_type.CPU: 266 if cpu is None or preferred_cpu: 267 cpu = device 268 else: 269 # System has cl.device_type.ACCELERATOR or cl.device_type.CUSTOM 270 # Intel Phi for example registers as an accelerator 271 # Since the user installed a custom device on their system and went through the 272 # pain of sorting out OpenCL drivers for it, lets assume they really do want to 273 # use it as their primary compute device. 274 gpu = device 275 276 # order the devices by gpu then by cpu; when searching for an available device by data type they 277 # will be checked in this order, which means that if the gpu supports double then the cpu will never 278 # be used (though we may make it possible to explicitly request the cpu at some point). 279 devices = [] 280 if gpu is not None: 281 devices.append(gpu) 282 if cpu is not None: 283 devices.append(cpu) 284 return [cl.Context([d]) for d in devices] 234 285 235 286 … … 314 365 # architectures tested so far. 315 366 self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors] 367 context = env.get_context(self.dtype) 316 368 self.q_buffers = [ 317 cl.Buffer( env.context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q)369 cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q) 318 370 for q in self.q_vectors 319 371 ] … … 363 415 # Note: res may be shorter than res_b if global_size != nq 364 416 env = environment() 365 self.loops_b = [cl.Buffer(env.context, mf.READ_WRITE, 366 2 * MAX_LOOPS * q_input.dtype.itemsize) 367 for _ in env.queues] 368 self.res_b = [cl.Buffer(env.context, mf.READ_WRITE, 369 q_input.global_size[0] * q_input.dtype.itemsize) 370 for _ in env.queues] 417 self.queue = env.get_queue(dtype) 418 self.loops_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 419 2 * MAX_LOOPS * q_input.dtype.itemsize) 420 self.res_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 421 q_input.global_size[0] * q_input.dtype.itemsize) 371 422 self.q_input = q_input 423 424 self._need_release = [self.loops_b, self.res_b, self.q_input] 372 425 373 426 def __call__(self, fixed_pars, pd_pars, cutoff=1e-5): … … 377 430 else np.float32) # will never get here, so use np.float32 378 431 379 device_num = 0 380 queuei = environment().queues[device_num] 381 res_bi = self.res_b[device_num] 432 res_bi = self.res_b 382 433 nq = np.uint32(self.q_input.nq) 383 434 if pd_pars: … … 394 445 raise ValueError("too many polydispersity points") 395 446 396 loops_bi = self.loops_b [device_num]397 cl.enqueue_copy( queuei, loops_bi, loops)447 loops_bi = self.loops_b 448 cl.enqueue_copy(self.queue, loops_bi, loops) 398 449 loops_l = cl.LocalMemory(len(loops.data)) 399 450 #ctx = environment().context … … 404 455 fixed = [real(p) for p in fixed_pars] 405 456 args = self.q_input.q_buffers + [res_bi, nq] + dispersed + fixed 406 self.kernel( queuei, self.q_input.global_size, None, *args)407 cl.enqueue_copy( queuei, self.res, res_bi)457 self.kernel(self.queue, self.q_input.global_size, None, *args) 458 cl.enqueue_copy(self.queue, self.res, res_bi) 408 459 409 460 return self.res … … 413 464 Release resources associated with the kernel. 414 465 """ 415 for b in self.loops_b: 416 b.release() 417 self.loops_b = [] 418 for b in self.res_b: 419 b.release() 420 self.res_b = [] 421 self.q_input.release() 466 for v in self._need_release: 467 v.release() 468 self._need_release = [] 422 469 423 470 def __del__(self): -
sasmodels/model_test.py
r5c962df r13ed84c 100 100 test_name = "Model: %s, Kernel: OpenCL"%model_name 101 101 test_method_name = "test_%s_opencl" % model_name 102 # Using dtype=None so that the models that are only 103 # correct for double precision are not tested using 104 # single precision. The choice is determined by the 105 # presence of *single=False* in the model file. 102 106 test = ModelTestCase(test_name, model_definition, 103 107 test_method_name, 104 platform="ocl", dtype= 'single')108 platform="ocl", dtype=None) 105 109 #print("defining", test_name) 106 110 suite.addTest(test) … … 157 161 ## values an error. Only do so for the "dll" tests 158 162 ## to reduce noise from both opencl and dll, and because 159 ## python kernels us 163 ## python kernels use platform="dll". 160 164 #raise Exception("No test cases provided") 161 165 pass … … 198 202 'invalid f(%s): %s' % (xi, actual_yi)) 199 203 else: 200 err = abs(yi - actual_yi) 201 nrm = abs(yi) 202 self.assertLess(err * 10**5, nrm, 204 self.assertTrue(is_near(yi, actual_yi, 5), 203 205 'f(%s); expected:%s; actual:%s' 204 206 % (xi, yi, actual_yi)) … … 206 208 return ModelTestCase 207 209 208 210 def is_near(target, actual, digits=5): 211 """ 212 Returns true if *actual* is within *digits* significant digits of *target*. 213 """ 214 import math 215 shift = 10**math.ceil(math.log10(abs(target))) 216 return abs(target-actual)/shift < 1.5*10**-digits 209 217 210 218 def main(): … … 217 225 218 226 models = sys.argv[1:] 227 if models and models[0] == '-v': 228 verbosity = 2 229 models = models[1:] 230 else: 231 verbosity = 1 219 232 if models and models[0] == 'opencl': 220 233 if not HAVE_OPENCL: … … 235 248 print("""\ 236 249 usage: 237 python -m sasmodels.model_test [opencl|dll|opencl_and_dll] model1 model2 ... 250 python -m sasmodels.model_test [-v] [opencl|dll] model1 model2 ... 251 252 If -v is included on the 253 If neither opencl nor dll is specified, then models will be tested with 254 both opencl and dll; the compute target is ignored for pure python models. 238 255 239 256 If model1 is 'all', then all except the remaining models will be tested. 240 If no compute target is specified, then models will be tested with both opencl 241 and dll; the compute target is ignored for pure python models.""")257 258 """) 242 259 243 260 return 1 244 261 245 262 #runner = unittest.TextTestRunner() 246 runner = xmlrunner.XMLTestRunner(output='logs' )263 runner = xmlrunner.XMLTestRunner(output='logs', verbosity=verbosity) 247 264 result = runner.run(make_suite(loaders, models)) 248 265 return 1 if result.failures or result.errors else 0 -
sasmodels/models/HayterMSAsq.py
r7f47777 r13ed84c 56 56 parameters used in P(Q). 57 57 """ 58 single = False # double precision only for now 58 59 # [ "name", "units", default, [lower, upper], "type", "description" ], 59 60 parameters = [["effect_radius", "Ang", 20.75, [0, inf], "volume", -
sasmodels/models/bcc.py
rdcdf29d r13ed84c 116 116 """ 117 117 category = "shape:paracrystal" 118 119 single = False 120 118 121 # pylint: disable=bad-whitespace, line-too-long 119 122 # ["name", "units", default, [lower, upper], "type","description" ], -
sasmodels/models/core_shell_ellipsoid.py
r81dd619 r177c1a1 96 96 category = "shape:ellipsoid" 97 97 98 single = False # TODO: maybe using sph_j1c inside gfn would help? 98 99 # pylint: disable=bad-whitespace, line-too-long 99 100 # ["name", "units", default, [lower, upper], "type", "description"], … … 111 112 # pylint: enable=bad-whitespace, line-too-long 112 113 113 source = ["lib/ gfn.c", "lib/gauss76.c", "core_shell_ellipsoid.c"]114 source = ["lib/sph_j1c.c", "lib/gfn.c", "lib/gauss76.c", "core_shell_ellipsoid.c"] 114 115 115 116 demo = dict(scale=1, background=0.001, -
sasmodels/models/fcc.c
r82d239a reeb8bac 101 101 102 102 double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z, 103 double q_z;103 // double q_z; 104 104 double cos_val_b3, cos_val_b2, cos_val_b1; 105 105 double a1_dot_q, a2_dot_q,a3_dot_q; … … 124 124 const double latticescale = 2.0*(4.0/3.0)*M_PI*(radius*radius*radius)/(s1*s1*s1); 125 125 // q vector 126 q_z = 0.0; // for SANS; assuming qz is negligible126 // q_z = 0.0; // for SANS; assuming qz is negligible 127 127 /// Angles here are respect to detector coordinate 128 128 /// instead of against q coordinate(PRB 36(46), 3(6), 1754(3854)) -
sasmodels/models/fcc.py
reb69cce r13ed84c 112 112 category = "shape:paracrystal" 113 113 114 single = False 115 114 116 # ["name", "units", default, [lower, upper], "type","description"], 115 117 parameters = [["dnn", "Ang", 220, [-inf, inf], "", "Nearest neighbour distance"], -
sasmodels/models/flexible_cylinder.py
r168052c r13ed84c 78 78 79 79 category = "shape:cylinder" 80 single = False 80 81 81 82 # pylint: disable=bad-whitespace, line-too-long -
sasmodels/models/flexible_cylinder_ex.py
r504abee r13ed84c 98 98 during model fitting. 99 99 """ 100 single = False 100 101 101 102 category = "shape:cylinder" -
sasmodels/models/gaussian_peak.py
reb69cce r13ed84c 42 42 category = "shape-independent" 43 43 44 single = False 44 45 # ["name", "units", default, [lower, upper], "type","description"], 45 46 parameters = [["q0", "1/Ang", 0.05, [-inf, inf], "", "Peak position"], -
sasmodels/models/hardsphere.py
r7f47777 r13ed84c 55 55 "volume fraction of hard spheres"], 56 56 ] 57 single = False 57 58 58 59 # No volume normalization despite having a volume parameter -
sasmodels/models/lamellarCaille.py
r7f47777 r13ed84c 87 87 category = "shape:lamellae" 88 88 89 single = False 90 89 91 # ["name", "units", default, [lower, upper], "type","description"], 90 92 parameters = [["thickness", "Ang", 30.0, [0, inf], "volume", "sheet thickness"], -
sasmodels/models/lamellarCailleHG.py
r7f47777 r13ed84c 91 91 category = "shape:lamellae" 92 92 93 single = False 93 94 parameters = [ 94 95 # [ "name", "units", default, [lower, upper], "type", -
sasmodels/models/lamellarPC.py
r7f47777 r13ed84c 111 111 category = "shape:lamellae" 112 112 113 single = False 114 113 115 # ["name", "units", default, [lower, upper], "type","description"], 114 116 parameters = [["thickness", "Ang", 33.0, [0, inf], "volume", -
sasmodels/models/lib/gfn.c
r81dd619 r177c1a1 7 7 // function gfn4 for oblate ellipsoids 8 8 double 9 gfn4(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq); 10 double 9 11 gfn4(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq) 10 12 { 11 // local variables 12 double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,tgfn,gfn4,pi43,Pi; 13 // local variables 14 const double pi43=4.0/3.0*M_PI; 15 const double aa = crmaj; 16 const double bb = crmin; 17 const double u2 = (bb*bb*xx*xx + aa*aa*(1.0-xx*xx)); 18 const double uq = sqrt(u2)*qq; 19 // changing to more accurate sph_j1c since the following inexplicably fails on Radeon Nano. 20 //const double siq = (uq == 0.0 ? 1.0 : 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq); 21 const double siq = sph_j1c(uq); 22 const double vc = pi43*aa*aa*bb; 23 const double gfnc = siq*vc*delpc; 13 24 14 Pi = 4.0*atan(1.0); 15 pi43=4.0/3.0*Pi; 16 aa = crmaj; 17 bb = crmin; 18 u2 = (bb*bb*xx*xx + aa*aa*(1.0-xx*xx)); 19 ut2 = (trmin*trmin*xx*xx + trmaj*trmaj*(1.0-xx*xx)); 20 uq = sqrt(u2)*qq; 21 ut= sqrt(ut2)*qq; 22 vc = pi43*aa*aa*bb; 23 vt = pi43*trmaj*trmaj*trmin; 24 if (uq == 0.0){ 25 siq = 1.0/3.0; 26 }else{ 27 siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq; 28 } 29 if (ut == 0.0){ 30 sit = 1.0/3.0; 31 }else{ 32 sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut; 33 } 34 gfnc = 3.0*siq*vc*delpc; 35 gfnt = 3.0*sit*vt*delps; 36 tgfn = gfnc+gfnt; 37 gfn4 = tgfn*tgfn; 25 const double ut2 = (trmin*trmin*xx*xx + trmaj*trmaj*(1.0-xx*xx)); 26 const double ut= sqrt(ut2)*qq; 27 const double vt = pi43*trmaj*trmaj*trmin; 28 //const double sit = (ut == 0.0 ? 1.0 : 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut); 29 const double sit = sph_j1c(ut); 30 const double gfnt = sit*vt*delps; 38 31 39 return (gfn4); 32 const double tgfn = gfnc + gfnt; 33 const double result = tgfn*tgfn; 34 35 return (result); 40 36 } -
sasmodels/models/lib/wrc_cyl.c
r504abee r13ed84c 379 379 } 380 380 381 double Sk_WR(double q, double L, double b); 381 382 double Sk_WR(double q, double L, double b) 382 383 { -
sasmodels/models/pearl_necklace.py
rf12357f rd18582e 96 96 97 97 source = ["lib/Si.c", "pearl_necklace.c"] 98 # new flag to let the compiler know to never use single precision 99 single = False 98 single = False # use double precision unless told otherwise 100 99 101 100 def volume(radius, edge_separation, string_thickness, number_of_pearls): -
sasmodels/models/rpa.c
r82c299f r13ed84c 205 205 const double Kbb = 0.0; 206 206 const double Kcc = 0.0; 207 const double Kdd = 0.0;207 //const double Kdd = 0.0; 208 208 const double Zaa = Kaa - Kad - Kad; 209 209 const double Zab = Kab - Kad - Kbd; … … 278 278 const double Q12 = (-Mab*Mcc + Mac*Mcb)/DenQ; 279 279 const double Q13 = ( Mab*Mbc - Mac*Mbb)/DenQ; 280 const double Q21 = (-Mba*Mcc + Mbc*Mca)/DenQ;280 //const double Q21 = (-Mba*Mcc + Mbc*Mca)/DenQ; 281 281 const double Q22 = ( Maa*Mcc - Mac*Mca)/DenQ; 282 282 const double Q23 = (-Maa*Mbc + Mac*Mba)/DenQ; 283 const double Q31 = ( Mba*Mcb - Mbb*Mca)/DenQ;284 const double Q32 = (-Maa*Mcb + Mab*Mca)/DenQ;283 //const double Q31 = ( Mba*Mcb - Mbb*Mca)/DenQ; 284 //const double Q32 = (-Maa*Mcb + Mab*Mca)/DenQ; 285 285 const double Q33 = ( Maa*Mbb - Mab*Mba)/DenQ; 286 286 -
sasmodels/models/star_polymer.py
r168052c r13ed84c 55 55 """ 56 56 category = "shape-independent" 57 single = False 57 58 # pylint: disable=bad-whitespace, line-too-long 58 59 # ["name", "units", default, [lower, upper], "type","description"], -
sasmodels/models/stickyhardsphere.py
r7f47777 r13ed84c 85 85 category = "structure-factor" 86 86 87 single = False 87 88 # ["name", "units", default, [lower, upper], "type","description"], 88 89 parameters = [
Note: See TracChangeset
for help on using the changeset viewer.