Changeset 87985ca in sasmodels for compare.py
- Timestamp:
- Sep 2, 2014 2:35:23 PM (10 years ago)
- Branches:
- master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- 19dcb933
- Parents:
- f4cf580
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
compare.py
- Property mode changed from 100644 to 100755
ra953943 r87985ca 2 2 # -*- coding: utf-8 -*- 3 3 4 from sasmodel import SasModel, load_data, set_beam_stop, plot_data, tic, toc 4 import sys 5 import math 6 5 7 import numpy as np 6 8 9 from sasmodels.bumps_model import BumpsModel, plot_data, tic 10 from sasmodels import gpu, dll 11 from sasmodels.convert import revert_model 12 7 13 8 14 def sasview_model(modelname, **pars): 9 modelname = modelname+"Model" 15 """ 16 Load a sasview model given the model name. 17 """ 18 # convert model parameters from sasmodel form to sasview form 19 #print "old",sorted(pars.items()) 20 modelname, pars = revert_model(modelname, pars) 21 #print "new",sorted(pars.items()) 10 22 sans = __import__('sans.models.'+modelname) 11 23 ModelClass = getattr(getattr(sans.models,modelname,None),modelname,None) … … 21 33 elif k.endswith("_pd_nsigma"): 22 34 model.dispersion[k[:-10]]['nsigmas'] = v 35 elif k.endswith("_pd_type"): 36 model.dispersion[k[:-8]]['type'] = v 23 37 else: 24 38 model.setParam(k, v) 25 39 return model 26 40 27 28 def sasview_eval(model, data): 29 theory = model.evalDistribution([data.qx_data, data.qy_data]) 30 return theory 31 32 def plot(data, cpumodel, gpumodel, N, pars): 33 34 model = SasModel(data, gpumodel, dtype='f', **pars) 35 tic() 36 for i in range(N): 37 #pars['scale'] = np.random.rand() 38 model.update() 39 gpu = model.theory() 40 gpu_time = toc()*1000./N 41 print "ocl t=%.1f ms"%gpu_time 42 tic() 43 for i in range(N): 44 cpu = sasview_eval(cpumodel, data) 45 cpu_time = toc()*1000./N 46 print "omp t=%.1f ms"%cpu_time 47 41 def load_opencl(modelname, dtype='single'): 42 sasmodels = __import__('sasmodels.models.'+modelname) 43 module = getattr(sasmodels.models, modelname, None) 44 kernel = gpu.load_model(module, dtype=dtype) 45 return kernel 46 47 def load_ctypes(modelname, dtype='single'): 48 sasmodels = __import__('sasmodels.models.'+modelname) 49 module = getattr(sasmodels.models, modelname, None) 50 kernel = dll.load_model(module, dtype=dtype) 51 return kernel 52 53 def randomize(p, v): 54 """ 55 Randomizing parameter. 56 57 Guess the parameter type from name. 58 """ 59 if any(p.endswith(s) for s in ('_pd_n','_pd_nsigma','_pd_type')): 60 return v 61 elif any(s in p for s in ('theta','phi','psi')): 62 # orientation in [-180,180], orientation pd in [0,45] 63 if p.endswith('_pd'): 64 return 45*np.random.rand() 65 else: 66 return 360*np.random.rand() - 180 67 elif 'sld' in p: 68 # sld in in [-0.5,10] 69 return 10.5*np.random.rand() - 0.5 70 elif p.endswith('_pd'): 71 # length pd in [0,1] 72 return np.random.rand() 73 else: 74 # length, scale, background in [0,200] 75 return 200*np.random.rand() 76 77 def parlist(pars): 78 return "\n".join("%s: %s"%(p,v) for p,v in sorted(pars.items())) 79 80 def suppress_pd(pars): 81 """ 82 Suppress theta_pd for now until the normalization is resolved. 83 84 May also suppress complete polydispersity of the model to test 85 models more quickly. 86 """ 87 for p in pars: 88 if p.endswith("_pd"): pars[p] = 0 89 90 def compare(name, pars, Ncpu, Ngpu, opts, set_pars): 91 92 # Sort out data 93 qmax = 1.0 if '-highq' in opts else (0.2 if '-midq' in opts else 0.05) 94 if "-1d" in opts: 95 from sasmodels.bumps_model import empty_data1D 96 qmax = math.log10(qmax) 97 data = empty_data1D(np.logspace(qmax-3, qmax, 128)) 98 index = slice(None, None) 99 else: 100 from sasmodels.bumps_model import empty_data2D, set_beam_stop 101 data = empty_data2D(np.linspace(-qmax, qmax, 128)) 102 set_beam_stop(data, 0.004) 103 index = ~data.mask 104 is2D = hasattr(data, 'qx_data') 105 106 # modelling accuracy is determined by dtype and cutoff 107 dtype = 'double' if '-double' in opts else 'single' 108 cutoff_opts = [s for s in opts if s.startswith('-cutoff')] 109 cutoff = float(cutoff_opts[0].split('=')[1]) if cutoff_opts else 1e-5 110 111 # randomize parameters 112 random_opts = [s for s in opts if s.startswith('-random')] 113 if random_opts: 114 if '=' in random_opts[0]: 115 seed = int(random_opts[0].split('=')[1]) 116 else: 117 seed = int(np.random.rand()*10000) 118 np.random.seed(seed) 119 print "Randomize using -random=%i"%seed 120 # Note: the sort guarantees order of calls to random number generator 121 pars = dict((p,randomize(p,v)) for p,v in sorted(pars.items())) 122 # The capped cylinder model has a constraint on its parameters 123 if name == 'capped_cylinder' and pars['cap_radius'] < pars['radius']: 124 pars['radius'],pars['cap_radius'] = pars['cap_radius'],pars['radius'] 125 pars.update(set_pars) 126 127 # parameter selection 128 if '-mono' in opts: 129 suppress_pd(pars) 130 if '-pars' in opts: 131 print "pars",parlist(pars) 132 133 # OpenCl calculation 134 if Ngpu > 0: 135 try: 136 model = load_opencl(name, dtype=dtype) 137 except Exception,exc: 138 print exc 139 print "... trying again with single precision" 140 model = load_opencl(name, dtype='single') 141 problem = BumpsModel(data, model, cutoff=cutoff, **pars) 142 toc = tic() 143 for _ in range(Ngpu): 144 #pars['scale'] = np.random.rand() 145 problem.update() 146 gpu = problem.theory() 147 gpu_time = toc()*1000./Ngpu 148 print "opencl t=%.1f ms, intensity=%.0f"%(gpu_time, sum(gpu[index])) 149 #print max(gpu), min(gpu) 150 151 # ctypes/sasview calculation 152 if Ncpu > 0 and "-ctypes" in opts: 153 model = load_ctypes(name, dtype=dtype) 154 problem = BumpsModel(data, model, cutoff=cutoff, **pars) 155 toc = tic() 156 for _ in range(Ncpu): 157 problem.update() 158 cpu = problem.theory() 159 cpu_time = toc()*1000./Ncpu 160 comp = "ctypes" 161 print "ctypes t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[index])) 162 elif Ncpu > 0: 163 model = sasview_model(name, **pars) 164 toc = tic() 165 for _ in range(Ncpu): 166 if is2D: 167 cpu = model.evalDistribution([data.qx_data, data.qy_data]) 168 else: 169 cpu = model.evalDistribution(data.x) 170 cpu_time = toc()*1000./Ncpu 171 comp = "sasview" 172 print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[index])) 173 174 # Compare, but only if computing both forms 175 if Ngpu > 0 and Ncpu > 0: 176 #print "speedup %.2g"%(cpu_time/gpu_time) 177 #print "max |gpu/cpu|", max(abs(gpu/cpu)), "%.15g"%max(abs(gpu)), "%.15g"%max(abs(cpu)) 178 #cpu *= max(gpu/cpu) 179 resid, relerr = np.zeros_like(gpu), np.zeros_like(gpu) 180 resid[index] = (gpu - cpu)[index] 181 relerr[index] = resid[index]/cpu[index] 182 print "max(|ocl-%s|)"%comp, max(abs(resid[index])) 183 print "max(|(ocl-%s)/ocl|)"%comp, max(abs(relerr[index])) 184 185 # Plot if requested 186 if '-noplot' in opts: return 48 187 import matplotlib.pyplot as plt 49 50 print "gpu/cpu", max(gpu/cpu) 51 cpu *= max(gpu/cpu) 52 relerr = (gpu - cpu)/cpu 53 print "max(|(ocl-omp)/ocl|)", max(abs(relerr)) 54 55 56 plt.subplot(131); plot_data(data, cpu); plt.title("omp t=%.1f ms"%cpu_time) 57 plt.subplot(132); plot_data(data, gpu); plt.title("ocl t=%.1f ms"%gpu_time) 58 plt.subplot(133); plot_data(data, relerr); plt.title("relerr x 10^8"); plt.colorbar() 188 if Ncpu > 0: 189 if Ngpu > 0: plt.subplot(131) 190 plot_data(data, cpu, scale='log') 191 plt.title("%s t=%.1f ms"%(comp,cpu_time)) 192 if Ngpu > 0: 193 if Ncpu > 0: plt.subplot(132) 194 plot_data(data, gpu, scale='log') 195 plt.title("opencl t=%.1f ms"%gpu_time) 196 if Ncpu > 0 and Ngpu > 0: 197 plt.subplot(133) 198 err = resid if '-abs' in opts else relerr 199 errstr = "abs err" if '-abs' in opts else "rel err" 200 #err,errstr = gpu/cpu,"ratio" 201 plot_data(data, err, scale='linear') 202 plt.title("max %s = %.3g"%(errstr, max(abs(err[index])))) 203 if is2D: plt.colorbar() 59 204 plt.show() 60 205 61 def newlen(N): 62 import sys 63 64 if len(sys.argv) > 1: 65 N = int(sys.argv[1]) 66 data = load_data('December/DEC07098.DAT') 67 set_beam_stop(data, 0.004) 68 69 return data, N 70 71 def cyl(N=1): 72 73 dtype = "float" 74 pars = dict( 75 scale=.003, radius=64.1, length=66.96, sldCyl=.291e-6, sldSolv=5.77e-6, background=.1, 76 cyl_theta=90, cyl_phi=0, 77 radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3, 78 length_pd=0.1,length_pd_n=5, length_pd_nsigma=3, 79 cyl_theta_pd=0.1, cyl_theta_pd_n=5, cyl_theta_pd_nsigma=3, 80 cyl_phi_pd=0.1, cyl_phi_pd_n=10, cyl_phi_pd_nsigma=3) 81 82 from Models.code_cylinder_f import GpuCylinder as gpumodel 83 model = sasview_model('Cylinder', **pars) 84 data, N = newlen(N) 85 86 plot(data, model, gpumodel, N, pars) 87 88 89 def ellipse(N=1): 90 91 pars = dict(scale=.027, radius_a=60, radius_b=180, sldEll=.297e-6, sldSolv=5.773e-6, background=4.9, 92 axis_theta=0, axis_phi=90, 93 radius_a_pd=0.1, radius_a_pd_n=10, radius_a_pd_nsigma=3, 94 radius_b_pd=0.1, radius_b_pd_n=10, radius_b_pd_nsigma=3, 95 axis_theta_pd=0.1, axis_theta_pd_n=6, axis_theta_pd_nsigma=3, 96 axis_phi_pd=0.1, axis_phi_pd_n=6, axis_phi_pd_nsigma=3,) 97 98 from Models.code_ellipse import GpuEllipse as gpumodel 99 model = sasview_model('Ellipsoid', **pars) 100 101 data, N = newlen(N) 102 plot(data, model, gpumodel, N, pars) 103 104 def coreshell(N=1): 105 106 data, N = newlen(N) 107 108 pars = dict(scale= 1.77881e-06, radius=325, thickness=25, length=34.2709, 109 core_sld=1e-6, shell_sld=.291e-6, solvent_sld=7.105e-6, 110 background=223.827, axis_theta=90, axis_phi=0, 111 axis_theta_pd=15.8, 112 radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3, 113 length_pd=0.1, length_pd_n=10, length_pd_nsigma=3, 114 thickness_pd=0.1, thickness_pd_n=5, thickness_pd_nsigma=3, 115 axis_theta_pd_n=20, axis_theta_pd_nsigma=5, 116 axis_phi_pd=0.0008748, axis_phi_pd_n=5, axis_phi_pd_nsigma=3,) 117 118 model = sasview_model('CoreShellCylinder', **pars) 119 from Models.code_coreshellcyl_f import GpuCoreShellCylinder as gpumodel 120 121 plot(data, model, gpumodel, N, pars) 122 123 def trellipse(N=1): 124 125 data, N = newlen(N) 126 127 pars = dict(scale=0.08, semi_axisA=15, semi_axisB=20, semi_axisC=500, 128 sldEll=7.105e-6, sldSolv=.291e-6, 129 background=5, axis_theta=0, axis_phi=0, axis_psi=0, 130 axis_theta_pd=20, axis_theta_pd_n=10, axis_theta_pd_nsigma=3, 131 axis_phi_pd=.1, axis_phi_pd_n=10, axis_phi_pd_nsigma=3, 132 axis_psi_pd=30, axis_psi_pd_n=5, axis_psi_pd_nsigma=3, 133 semi_axisA_pd=.1, semi_axisA_pd_n=5, semi_axisA_pd_nsigma=3, 134 semi_axisB_pd=.1, semi_axisB_pd_n=5, semi_axisB_pd_nsigma=3, 135 semi_axisC_pd=.1, semi_axisC_pd_n=5, semi_axisC_pd_nsigma=3,) 136 137 model = sasview_model('TriaxialEllipsoid', **pars) 138 from Models.code_triaxialellipse import GpuTriEllipse as gpumodel 139 140 plot(data, model, gpumodel, N, pars) 141 142 def lamellar(N=1): 143 144 data, N = newlen(N) 145 146 pars = dict(scale=0.08, 147 bi_thick=19.2946, 148 sld_bi=5.38e-6,sld_sol=7.105e-6, 149 background=0.003, 150 bi_thick_pd= 0.37765, bi_thick_pd_n=40, bi_thick_pd_nsigma=3) 151 152 model = sasview_model('Lamellar', **pars) 153 from Models.code_lamellar import GpuLamellar as gpumodel 154 155 plot(data, model, gpumodel, N, pars) 156 157 def cap(N=1): 158 159 data, N = newlen(N) 160 161 pars = dict(scale=.08, rad_cyl=20, rad_cap=40, len_cyl=400, 162 sld_capcyl=1e-6, sld_solv=6.3e-6, 163 background=0, theta=0, phi=0, 164 rad_cyl_pd=.1, rad_cyl_pd_n=10, rad_cyl_pd_nsigma=3, 165 rad_cap_pd=.1, rad_cap_pd_n=10, rad_cap_pd_nsigma=3, 166 len_cyl_pd=.1, len_cyl_pd_n=3, len_cyl_pd_nsigma=3, 167 theta_pd=.1, theta_pd_n=3, theta_pd_nsigma=3, 168 phi_pd=.1, phi_pd_n=3, phi_pd_nsigma=3) 169 170 171 model = sasview_model('CappedCylinder', **pars) 172 from Models.code_capcyl import GpuCapCylinder as gpumodel 173 174 plot(data, model, gpumodel, N, pars) 175 206 # =========================================================================== 207 # 208 USAGE=""" 209 usage: compare.py model [Nopencl] [Nsasview] [options...] [key=val] 210 211 Compare the speed and value for a model between the SasView original and the 212 OpenCL rewrite. 213 214 model is the name of the model to compare (see below). 215 Nopencl is the number of times to run the OpenCL model (default=5) 216 Nsasview is the number of times to run the Sasview model (default=1) 217 218 Options (* for default): 219 220 -plot*/-noplot plots or suppress the plot of the model 221 -single/-double uses double precision for comparison 222 -lowq/-midq/-highq use q values up to 0.05, 0.2 or 1.0 223 -2d/-1d uses 1d or 2d random data 224 -preset/-random[=seed] preset or random parameters 225 -poly/-mono force monodisperse/polydisperse 226 -sasview/-ctypes whether cpu is tested using sasview or ctypes 227 -cutoff=1e-5/value cutoff for including a point in polydispersity 228 -nopars/-pars prints the parameter set or not 229 -rel/-abs plot relative or absolute error 230 231 Key=value pairs allow you to set specific values to any of the model 232 parameters. 233 234 Available models: 235 236 %s 237 """ 238 239 VALID_OPTIONS = [ 240 'plot','noplot', 241 'single','double', 242 'lowq','midq','highq', 243 '2d','1d', 244 'preset','random', 245 'poly','mono', 246 'sasview','ctypes', 247 'nopars','pars', 248 'rel','abs', 249 ] 250 251 def main(): 252 opts = [arg for arg in sys.argv[1:] if arg.startswith('-')] 253 args = [arg for arg in sys.argv[1:] if not arg.startswith('-')] 254 models = "\n ".join("%-7s: %s"%(k,v.__name__.replace('_',' ')) 255 for k,v in sorted(MODELS.items())) 256 if len(args) == 0: 257 print(USAGE%models) 258 sys.exit(1) 259 if args[0] not in MODELS: 260 print "Model %r not available. Use one of:\n %s"%(args[0],models) 261 sys.exit(1) 262 263 valid_opts = set(VALID_OPTIONS) 264 invalid = [o[1:] for o in opts 265 if o[1:] not in valid_opts 266 and not o.startswith('-cutoff=') 267 and not o.startswith('-random=')] 268 if invalid: 269 print "Invalid options: %s"%(", ".join(invalid)) 270 sys.exit(1) 271 272 name, pars = MODELS[args[0]]() 273 Nopencl = int(args[1]) if len(args) > 1 else 5 274 Nsasview = int(args[2]) if len(args) > 3 else 1 275 276 # Fill in default polydispersity parameters 277 pds = set(p.split('_pd')[0] for p in pars if p.endswith('_pd')) 278 for p in pds: 279 if p+"_pd_nsigma" not in pars: pars[p+"_pd_nsigma"] = 3 280 if p+"_pd_type" not in pars: pars[p+"_pd_type"] = "gaussian" 281 282 # Fill in parameters given on the command line 283 set_pars = {} 284 for arg in args[3:]: 285 k,v = arg.split('=') 286 if k not in pars: 287 # extract base name without distribution 288 s = set(p.split('_pd')[0] for p in pars) 289 print "%r invalid; parameters are: %s"%(k,", ".join(sorted(s))) 290 sys.exit(1) 291 set_pars[k] = float(v) if not v.endswith('type') else v 292 293 compare(name, pars, Nsasview, Nopencl, opts, set_pars) 294 295 # =========================================================================== 296 # 297 298 MODELS = {} 299 def model(name): 300 def gather_function(fn): 301 MODELS[name] = fn 302 return fn 303 return gather_function 304 305 306 @model('cyl') 307 def cylinder(): 308 pars = dict( 309 scale=1, background=0, 310 sld=6, solvent_sld=1, 311 #radius=5, length=20, 312 radius=260, length=290, 313 theta=30, phi=0, 314 radius_pd=.2, radius_pd_n=1, 315 length_pd=.2,length_pd_n=1, 316 theta_pd=15, theta_pd_n=45, 317 phi_pd=15, phi_pd_n=1, 318 ) 319 return 'cylinder', pars 320 321 @model('capcyl') 322 def capped_cylinder(): 323 pars = dict( 324 scale=1, background=0, 325 sld=6, solvent_sld=1, 326 radius=260, cap_radius=290, length=290, 327 theta=30, phi=15, 328 radius_pd=.2, radius_pd_n=1, 329 cap_radius_pd=.2, cap_radius_pd_n=1, 330 length_pd=.2, length_pd_n=1, 331 theta_pd=15, theta_pd_n=45, 332 phi_pd=15, phi_pd_n=1, 333 ) 334 return 'capped_cylinder', pars 335 336 337 @model('cscyl') 338 def core_shell_cylinder(): 339 pars = dict( 340 scale=1, background=0, 341 core_sld=6, shell_sld=8, solvent_sld=1, 342 radius=45, thickness=25, length=340, 343 theta=30, phi=15, 344 radius_pd=.2, radius_pd_n=1, 345 length_pd=.2, length_pd_n=1, 346 thickness_pd=.2, thickness_pd_n=1, 347 theta_pd=15, theta_pd_n=45, 348 phi_pd=15, phi_pd_n=1, 349 ) 350 return 'core_shell_cylinder', pars 351 352 353 @model('ell') 354 def ellipsoid(): 355 pars = dict( 356 scale=1, background=0, 357 sld=6, solvent_sld=1, 358 rpolar=50, requatorial=30, 359 theta=30, phi=15, 360 rpolar_pd=.2, rpolar_pd_n=1, 361 requatorial_pd=.2, requatorial_pd_n=1, 362 theta_pd=15, theta_pd_n=45, 363 phi_pd=15, phi_pd_n=1, 364 ) 365 return 'ellipsoid', pars 366 367 368 @model('ell3') 369 def triaxial_ellipsoid(): 370 pars = dict( 371 scale=1, background=0, 372 sld=6, solvent_sld=1, 373 theta=30, phi=15, psi=5, 374 req_minor=25, req_major=36, rpolar=50, 375 req_minor_pd=0, req_minor_pd_n=1, 376 req_major_pd=0, req_major_pd_n=1, 377 rpolar_pd=.2, rpolar_pd_n=1, 378 theta_pd=15, theta_pd_n=45, 379 phi_pd=15, phi_pd_n=1, 380 psi_pd=15, psi_pd_n=1, 381 ) 382 return 'triaxial_ellipsoid', pars 383 384 @model('sph') 385 def sphere(): 386 pars = dict( 387 scale=1, background=0, 388 sld=6, solvent_sld=1, 389 radius=120, 390 radius_pd=.2, radius_pd_n=45, 391 ) 392 return 'sphere', pars 393 394 @model('lam') 395 def lamellar(): 396 pars = dict( 397 scale=1, background=0, 398 sld=6, solvent_sld=1, 399 thickness=40, 400 thickness_pd= 0.2, thickness_pd_n=40, 401 ) 402 return 'lamellar', pars 176 403 177 404 if __name__ == "__main__": 178 coreshell(1) 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 405 main()
Note: See TracChangeset
for help on using the changeset viewer.