[14de349] | 1 | #!/usr/bin/env python |
---|
| 2 | # -*- coding: utf-8 -*- |
---|
| 3 | |
---|
| 4 | import sys |
---|
| 5 | |
---|
| 6 | import numpy as np |
---|
| 7 | |
---|
[13d86bc] | 8 | from sasmodels.core import BumpsModel, plot_data, tic, opencl_model, dll_model |
---|
[14de349] | 9 | |
---|
| 10 | def sasview_model(modelname, **pars): |
---|
| 11 | """ |
---|
| 12 | Load a sasview model given the model name. |
---|
| 13 | """ |
---|
| 14 | modelname = modelname+"Model" |
---|
| 15 | sans = __import__('sans.models.'+modelname) |
---|
| 16 | ModelClass = getattr(getattr(sans.models,modelname,None),modelname,None) |
---|
| 17 | if ModelClass is None: |
---|
| 18 | raise ValueError("could not find model %r in sans.models"%modelname) |
---|
| 19 | model = ModelClass() |
---|
| 20 | |
---|
| 21 | for k,v in pars.items(): |
---|
| 22 | if k.endswith("_pd"): |
---|
| 23 | model.dispersion[k[:-3]]['width'] = v |
---|
| 24 | elif k.endswith("_pd_n"): |
---|
| 25 | model.dispersion[k[:-5]]['npts'] = v |
---|
| 26 | elif k.endswith("_pd_nsigma"): |
---|
| 27 | model.dispersion[k[:-10]]['nsigmas'] = v |
---|
[13d86bc] | 28 | elif k.endswith("_pd_type"): |
---|
| 29 | model.dispersion[k[:-8]]['type'] = v |
---|
[14de349] | 30 | else: |
---|
| 31 | model.setParam(k, v) |
---|
| 32 | return model |
---|
| 33 | |
---|
[13d86bc] | 34 | def load_opencl(modelname, dtype='single'): |
---|
| 35 | sasmodels = __import__('sasmodels.models.'+modelname) |
---|
| 36 | module = getattr(sasmodels.models, modelname, None) |
---|
| 37 | kernel = opencl_model(module, dtype=dtype) |
---|
| 38 | return kernel |
---|
| 39 | |
---|
| 40 | def load_dll(modelname, dtype='single'): |
---|
| 41 | sasmodels = __import__('sasmodels.models.'+modelname) |
---|
| 42 | module = getattr(sasmodels.models, modelname, None) |
---|
| 43 | kernel = dll_model(module, dtype=dtype) |
---|
| 44 | return kernel |
---|
| 45 | |
---|
[14de349] | 46 | |
---|
| 47 | def compare(Ncpu, cpuname, cpupars, Ngpu, gpuname, gpupars): |
---|
| 48 | |
---|
| 49 | #from sasmodels.core import load_data |
---|
| 50 | #data = load_data('December/DEC07098.DAT') |
---|
[13d86bc] | 51 | from sasmodels.core import empty_data1D |
---|
| 52 | data = empty_data1D(np.logspace(-4, -1, 128)) |
---|
| 53 | #from sasmodels.core import empty_2D, set_beam_stop |
---|
| 54 | #data = empty_data2D(np.linspace(-0.05, 0.05, 128)) |
---|
| 55 | #set_beam_stop(data, 0.004) |
---|
| 56 | is2D = hasattr(data, 'qx_data') |
---|
[14de349] | 57 | |
---|
| 58 | if Ngpu > 0: |
---|
[13d86bc] | 59 | gpumodel = load_opencl(gpuname, dtype='single') |
---|
[14de349] | 60 | model = BumpsModel(data, gpumodel, **gpupars) |
---|
| 61 | toc = tic() |
---|
| 62 | for i in range(Ngpu): |
---|
| 63 | #pars['scale'] = np.random.rand() |
---|
| 64 | model.update() |
---|
| 65 | gpu = model.theory() |
---|
| 66 | gpu_time = toc()*1000./Ngpu |
---|
[ce27e21] | 67 | print "ocl t=%.1f ms, intensity=%.0f"%(gpu_time, sum(gpu[~np.isnan(gpu)])) |
---|
[14de349] | 68 | #print max(gpu), min(gpu) |
---|
| 69 | |
---|
| 70 | if 0 and Ncpu > 0: # Hack to compare ctypes vs. opencl |
---|
[13d86bc] | 71 | dllmodel = load_dll(gpuname, dtype='double') |
---|
[14de349] | 72 | model = BumpsModel(data, dllmodel, **gpupars) |
---|
| 73 | toc = tic() |
---|
| 74 | for i in range(Ncpu): |
---|
| 75 | model.update() |
---|
| 76 | cpu = model.theory() |
---|
| 77 | cpu_time = toc()*1000./Ncpu |
---|
| 78 | print "dll t=%.1f ms"%cpu_time |
---|
| 79 | |
---|
[ce27e21] | 80 | elif 0: # Hack to check new vs old for GpuCylinder |
---|
[14de349] | 81 | from Models.code_cylinder_f import GpuCylinder as oldgpu |
---|
| 82 | from sasmodel import SasModel |
---|
| 83 | oldmodel = SasModel(data, oldgpu, dtype='single', **cpupars) |
---|
| 84 | toc = tic() |
---|
| 85 | for i in range(Ngpu): |
---|
| 86 | oldmodel.update() |
---|
| 87 | cpu = oldmodel.theory() |
---|
| 88 | cpu_time = toc()*1000./Ngpu |
---|
| 89 | print "old t=%.1f ms"%cpu_time |
---|
| 90 | |
---|
| 91 | elif Ncpu > 0: |
---|
| 92 | cpumodel = sasview_model(cpuname, **cpupars) |
---|
| 93 | toc = tic() |
---|
[13d86bc] | 94 | if is2D: |
---|
| 95 | for i in range(Ncpu): |
---|
| 96 | cpu = cpumodel.evalDistribution([data.qx_data, data.qy_data]) |
---|
| 97 | else: |
---|
| 98 | for i in range(Ncpu): |
---|
| 99 | cpu = cpumodel.evalDistribution(data.x) |
---|
[14de349] | 100 | cpu_time = toc()*1000./Ncpu |
---|
[ce27e21] | 101 | print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[model.index])) |
---|
[14de349] | 102 | |
---|
| 103 | if Ngpu > 0 and Ncpu > 0: |
---|
| 104 | print "gpu/cpu", max(abs(gpu/cpu)), "%.15g"%max(abs(gpu)), "%.15g"%max(abs(cpu)) |
---|
| 105 | #cpu *= max(gpu/cpu) |
---|
| 106 | abserr = (gpu - cpu) |
---|
| 107 | relerr = (gpu - cpu)/cpu |
---|
| 108 | print "max(|ocl-omp|)", max(abs(abserr[model.index])) |
---|
| 109 | print "max(|(ocl-omp)/ocl|)", max(abs(relerr[model.index])) |
---|
| 110 | #return |
---|
| 111 | |
---|
| 112 | import matplotlib.pyplot as plt |
---|
| 113 | if Ncpu > 0: |
---|
| 114 | if Ngpu > 0: plt.subplot(131) |
---|
[13d86bc] | 115 | plot_data(data, cpu, scale='log') |
---|
[14de349] | 116 | plt.title("omp t=%.1f ms"%cpu_time) |
---|
| 117 | if Ngpu > 0: |
---|
| 118 | if Ncpu > 0: plt.subplot(132) |
---|
[13d86bc] | 119 | plot_data(data, gpu, scale='log') |
---|
[14de349] | 120 | plt.title("ocl t=%.1f ms"%gpu_time) |
---|
| 121 | if Ncpu > 0 and Ngpu > 0: |
---|
| 122 | plt.subplot(133) |
---|
[13d86bc] | 123 | plot_data(data, 1e8*relerr, scale='linear') |
---|
[14de349] | 124 | plt.title("max rel err = %.3g"%max(abs(relerr))) |
---|
[13d86bc] | 125 | if is2D: plt.colorbar() |
---|
[14de349] | 126 | plt.show() |
---|
| 127 | |
---|
| 128 | def rename(pars, **names): |
---|
| 129 | newpars = pars.copy() |
---|
| 130 | for new,old in names.items(): |
---|
| 131 | for variant in ("", "_pd", "_pd_n", "_pd_nsigma"): |
---|
| 132 | if old+variant in newpars: |
---|
| 133 | newpars[new+variant] = pars[old+variant] |
---|
| 134 | del newpars[old+variant] |
---|
| 135 | return newpars |
---|
| 136 | |
---|
| 137 | def rescale_sld(pars, names): |
---|
| 138 | newpars = pars.copy() |
---|
| 139 | for p in names: |
---|
| 140 | newpars[p] *= 1e6 |
---|
| 141 | return newpars |
---|
| 142 | |
---|
| 143 | # =========================================================================== |
---|
| 144 | # |
---|
| 145 | MODELS = {} |
---|
| 146 | def model(name): |
---|
| 147 | def gather_function(fn): |
---|
| 148 | MODELS[name] = fn |
---|
| 149 | return fn |
---|
| 150 | return gather_function |
---|
| 151 | |
---|
| 152 | USAGE=""" |
---|
| 153 | usage: compare model [Nopencl] [Nsasview] |
---|
| 154 | |
---|
| 155 | Compare the speed and value for a model between the SasView original and the |
---|
| 156 | OpenCL rewrite. |
---|
| 157 | |
---|
| 158 | * Nopencl is the number of times to run the OpenCL model (default=5) |
---|
| 159 | |
---|
| 160 | * Nsasview is the number of times to run the Sasview model (default=1) |
---|
| 161 | |
---|
| 162 | * model is the name of the model to compare: |
---|
| 163 | |
---|
| 164 | %s |
---|
| 165 | """ |
---|
| 166 | |
---|
| 167 | def main(): |
---|
| 168 | if len(sys.argv) == 1: |
---|
| 169 | models = "\n ".join("%-7s: %s"%(k,v.__name__.replace('_',' ')) |
---|
| 170 | for k,v in sorted(MODELS.items())) |
---|
| 171 | print(USAGE%models) |
---|
| 172 | sys.exit(1) |
---|
| 173 | |
---|
| 174 | cpuname, cpupars, gpuname, gpupars = MODELS[sys.argv[1]]() |
---|
| 175 | Nopencl = int(sys.argv[2]) if len(sys.argv) > 2 else 5 |
---|
| 176 | Nsasview = int(sys.argv[3]) if len(sys.argv) > 3 else 1 |
---|
| 177 | |
---|
| 178 | compare(Nsasview, cpuname, cpupars, Nopencl, gpuname, gpupars) |
---|
| 179 | |
---|
| 180 | @model('cyl') |
---|
| 181 | def cylinder(): |
---|
| 182 | cpupars = dict( |
---|
| 183 | scale=.003, background=.1, |
---|
| 184 | sldCyl=.291e-6, sldSolv=5.77e-6, |
---|
| 185 | radius=264.1, length=66.96, |
---|
| 186 | cyl_theta=85, cyl_phi=0, |
---|
| 187 | radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3, |
---|
| 188 | length_pd=0.1,length_pd_n=1, length_pd_nsigma=3, |
---|
| 189 | cyl_theta_pd=45, cyl_theta_pd_n=50, cyl_theta_pd_nsigma=3, |
---|
| 190 | cyl_phi_pd=0.1, cyl_phi_pd_n=5, cyl_phi_pd_nsigma=3, |
---|
| 191 | ) |
---|
| 192 | cpuname = 'Cylinder' |
---|
| 193 | |
---|
| 194 | gpupars = rename(cpupars, theta='cyl_theta', phi='cyl_phi', sld='sldCyl', solvent_sld='sldSolv') |
---|
| 195 | gpupars = rescale_sld(gpupars, ['sld', 'solvent_sld']) |
---|
| 196 | gpuname = 'cylinder' |
---|
| 197 | return cpuname, cpupars, gpuname, gpupars |
---|
| 198 | |
---|
| 199 | |
---|
| 200 | @model('ell') |
---|
| 201 | def ellipse(): |
---|
| 202 | pars = dict( |
---|
| 203 | scale=.027, background=4.9, |
---|
| 204 | sldEll=.297e-6, sldSolv=5.773e-6, |
---|
| 205 | radius_a=60, radius_b=180, |
---|
| 206 | axis_theta=0, axis_phi=90, |
---|
| 207 | radius_a_pd=0.1, radius_a_pd_n=10, radius_a_pd_nsigma=3, |
---|
| 208 | radius_b_pd=0.1, radius_b_pd_n=10, radius_b_pd_nsigma=3, |
---|
| 209 | axis_theta_pd=0.1, axis_theta_pd_n=6, axis_theta_pd_nsigma=3, |
---|
| 210 | axis_phi_pd=0.1, axis_phi_pd_n=6, axis_phi_pd_nsigma=3, |
---|
| 211 | ) |
---|
| 212 | |
---|
| 213 | from Models.code_ellipse import GpuEllipse as gpumodel |
---|
| 214 | model = sasview_model('Ellipsoid', **pars) |
---|
| 215 | |
---|
| 216 | pars = rename(pars, theta='axis_theta', phi='axis_phi', sld='sldEll', solvent_sld='sldSolv') |
---|
| 217 | pars = rescale_sld(pars, ['sld', 'solvent_sld']) |
---|
| 218 | return model, gpumodel, pars |
---|
| 219 | |
---|
| 220 | |
---|
| 221 | @model('cscyl') |
---|
| 222 | def core_shell_cylinder(N=1): |
---|
| 223 | pars = dict( |
---|
| 224 | scale= 1.77881e-06, background=223.827, |
---|
| 225 | core_sld=1e-6, shell_sld=.291e-6, solvent_sld=7.105e-6, |
---|
| 226 | radius=325, thickness=25, length=34.2709, |
---|
| 227 | axis_theta=90, axis_phi=0, |
---|
| 228 | radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3, |
---|
| 229 | length_pd=0.1, length_pd_n=10, length_pd_nsigma=3, |
---|
| 230 | thickness_pd=0.1, thickness_pd_n=5, thickness_pd_nsigma=3, |
---|
| 231 | axis_theta_pd=15.8, axis_theta_pd_n=20, axis_theta_pd_nsigma=5, |
---|
| 232 | axis_phi_pd=0.0008748, axis_phi_pd_n=5, axis_phi_pd_nsigma=3, |
---|
| 233 | ) |
---|
| 234 | |
---|
| 235 | model = sasview_model('CoreShellCylinder', **pars) |
---|
| 236 | from Models.code_coreshellcyl_f import GpuCoreShellCylinder as gpumodel |
---|
| 237 | |
---|
| 238 | pars = rename(pars, theta='axis_theta', phi='axis_phi') |
---|
| 239 | pars = rescale_sld(pars, ['core_sld', 'shell_sld', 'solvent_sld']) |
---|
| 240 | return model, gpumodel, pars |
---|
| 241 | |
---|
| 242 | |
---|
| 243 | @model('ell3') |
---|
| 244 | def triaxial_ellipse(N=1): |
---|
| 245 | pars = dict( |
---|
| 246 | scale=0.08, background=5, |
---|
| 247 | sldEll=7.105e-6, sldSolv=.291e-6, |
---|
| 248 | axis_theta=0, axis_phi=0, axis_psi=0, |
---|
| 249 | semi_axisA=15, semi_axisB=20, semi_axisC=500, |
---|
| 250 | axis_theta_pd=20, axis_theta_pd_n=10, axis_theta_pd_nsigma=3, |
---|
| 251 | axis_phi_pd=.1, axis_phi_pd_n=10, axis_phi_pd_nsigma=3, |
---|
| 252 | axis_psi_pd=30, axis_psi_pd_n=5, axis_psi_pd_nsigma=3, |
---|
| 253 | semi_axisA_pd=.1, semi_axisA_pd_n=5, semi_axisA_pd_nsigma=3, |
---|
| 254 | semi_axisB_pd=.1, semi_axisB_pd_n=5, semi_axisB_pd_nsigma=3, |
---|
| 255 | semi_axisC_pd=.1, semi_axisC_pd_n=5, semi_axisC_pd_nsigma=3, |
---|
| 256 | ) |
---|
| 257 | |
---|
| 258 | model = sasview_model('TriaxialEllipsoid', **pars) |
---|
| 259 | from Models.code_triaxialellipse import GpuTriEllipse as gpumodel |
---|
| 260 | |
---|
| 261 | pars = rename(pars, |
---|
| 262 | theta='axis_theta', phi='axis_phi', psi='axis_psi', |
---|
| 263 | sld='sldEll', solvent_sld='sldSolv', |
---|
| 264 | radius_a='semi_axisA', radius_b='semi_axisB', |
---|
| 265 | radius_c='semi_axisC', |
---|
| 266 | ) |
---|
| 267 | pars = rescale_sld(pars, ['sld', 'solvent_sld']) |
---|
| 268 | return model, gpumodel, pars |
---|
| 269 | |
---|
| 270 | @model('lam') |
---|
| 271 | def lamellar(N=1): |
---|
| 272 | pars = dict( |
---|
| 273 | scale=0.08, background=0.003, |
---|
| 274 | sld_bi=5.38e-6,sld_sol=7.105e-6, |
---|
| 275 | bi_thick=19.2946, |
---|
| 276 | bi_thick_pd= 0.37765, bi_thick_pd_n=40, bi_thick_pd_nsigma=3, |
---|
| 277 | ) |
---|
| 278 | |
---|
| 279 | model = sasview_model('Lamellar', **pars) |
---|
| 280 | from Models.code_lamellar import GpuLamellar as gpumodel |
---|
| 281 | |
---|
| 282 | pars = rename(pars, sld='sld_bi', solvent_sld='sld_sol', thickness='bi_thick') |
---|
| 283 | pars = rescale_sld(pars, ['sld', 'solvent_sld']) |
---|
| 284 | return model, gpumodel, pars |
---|
| 285 | |
---|
| 286 | @model('capcyl') |
---|
| 287 | def capped_cylinder(N=1): |
---|
| 288 | pars = dict( |
---|
| 289 | scale=.08, background=0, |
---|
| 290 | sld_capcyl=1e-6, sld_solv=6.3e-6, |
---|
| 291 | rad_cyl=20, rad_cap=40, len_cyl=400, |
---|
| 292 | theta=0, phi=0, |
---|
| 293 | rad_cyl_pd=.1, rad_cyl_pd_n=10, rad_cyl_pd_nsigma=3, |
---|
| 294 | rad_cap_pd=.1, rad_cap_pd_n=10, rad_cap_pd_nsigma=3, |
---|
| 295 | len_cyl_pd=.1, len_cyl_pd_n=3, len_cyl_pd_nsigma=3, |
---|
| 296 | theta_pd=.1, theta_pd_n=3, theta_pd_nsigma=3, |
---|
| 297 | phi_pd=.1, phi_pd_n=3, phi_pd_nsigma=3, |
---|
| 298 | ) |
---|
| 299 | |
---|
| 300 | |
---|
| 301 | model = sasview_model('CappedCylinder', **pars) |
---|
| 302 | from Models.code_capcyl import GpuCapCylinder as gpumodel |
---|
| 303 | |
---|
| 304 | pars = rename(pars, |
---|
| 305 | sld='sld_capcyl', solvent_sld='sld_solv', |
---|
| 306 | length='len_cyl', radius='rad_cyl', |
---|
| 307 | cap_radius='rad_cap') |
---|
| 308 | pars = rescale_sld(pars, ['sld', 'solvent_sld']) |
---|
| 309 | return model, gpumodel, pars |
---|
| 310 | |
---|
| 311 | |
---|
| 312 | if __name__ == "__main__": |
---|
| 313 | main() |
---|