Changeset ec7e360 in sasmodels
- Timestamp:
- Dec 23, 2015 12:17:49 PM (8 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:
- e21cc31
- Parents:
- ce166d3
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
compare.sh
r5753e4e rec7e360 12 12 #SAS_OPENMP=1; export SAS_OPENMP 13 13 14 python-m sasmodels.compare $*14 ${PYTHON:-python} -m sasmodels.compare $* -
sasmodels/bumps_model.py
r9404dd3 rec7e360 12 12 """ 13 13 14 import datetime15 14 import warnings 16 15 17 16 import numpy as np 18 17 19 from . import sesans20 from . import weights21 18 from .data import plot_theory 22 19 from .direct_model import DataMixin … … 31 28 return experiment 32 29 30 def create_parameters(model_info, **kwargs): 31 # lazy import; this allows the doc builder and nosetests to run even 32 # when bumps is not on the path. 33 from bumps.names import Parameter 34 35 partype = model_info['partype'] 36 37 pars = {} 38 for p in model_info['parameters']: 39 name, default, limits = p[0], p[2], p[3] 40 value = kwargs.pop(name, default) 41 pars[name] = Parameter.default(value, name=name, limits=limits) 42 for name in partype['pd-2d']: 43 for xpart, xdefault, xlimits in [ 44 ('_pd', 0., limits), 45 ('_pd_n', 35., (0, 1000)), 46 ('_pd_nsigma', 3., (0, 10)), 47 ]: 48 xname = name + xpart 49 xvalue = kwargs.pop(xname, xdefault) 50 pars[xname] = Parameter.default(xvalue, name=xname, limits=xlimits) 51 52 pd_types = {} 53 for name in partype['pd-2d']: 54 xname = name + '_pd_type' 55 xvalue = kwargs.pop(xname, 'gaussian') 56 pd_types[xname] = xvalue 57 58 if kwargs: # args not corresponding to parameters 59 raise TypeError("unexpected parameters: %s" 60 % (", ".join(sorted(kwargs.keys())))) 61 62 return pars, pd_types 33 63 34 64 class Model(object): 35 def __init__(self, model, **kw): 36 # lazy import; this allows the doc builder and nosetests to run even 37 # when bumps is not on the path. 38 from bumps.names import Parameter 39 65 def __init__(self, model, **kwargs): 40 66 self._sasmodel = model 41 partype = model.info['partype'] 42 43 pars = [] 44 for p in model.info['parameters']: 45 name, default, limits = p[0], p[2], p[3] 46 value = kw.pop(name, default) 47 setattr(self, name, Parameter.default(value, name=name, limits=limits)) 48 pars.append(name) 49 for name in partype['pd-2d']: 50 for xpart, xdefault, xlimits in [ 51 ('_pd', 0, limits), 52 ('_pd_n', 35, (0, 1000)), 53 ('_pd_nsigma', 3, (0, 10)), 54 ('_pd_type', 'gaussian', None), 55 ]: 56 xname = name + xpart 57 xvalue = kw.pop(xname, xdefault) 58 if xlimits is not None: 59 xvalue = Parameter.default(xvalue, name=xname, limits=xlimits) 60 pars.append(xname) 61 setattr(self, xname, xvalue) 62 self._parameter_names = pars 63 if kw: 64 raise TypeError("unexpected parameters: %s" 65 % (", ".join(sorted(kw.keys())))) 67 pars, pd_types = create_parameters(model.info, **kwargs) 68 for k,v in pars.items(): 69 setattr(self, k, v) 70 for k,v in pd_types.items(): 71 setattr(self, k, v) 72 self._parameter_names = list(pars.keys()) 73 self._pd_type_names = list(pd_types.keys()) 66 74 67 75 def parameters(self): … … 71 79 return dict((k, getattr(self, k)) for k in self._parameter_names) 72 80 81 def state(self): 82 pars = dict((k, getattr(self, k).value) for k in self._parameter_names) 83 pars.update((k, getattr(self, k)) for k in self._pd_type_names) 84 return pars 73 85 74 86 class Experiment(DataMixin): … … 113 125 def theory(self): 114 126 if 'theory' not in self._cache: 115 pars = dict((k, v.value) for k,v in self.model.parameters().items())127 pars = self.model.state() 116 128 self._cache['theory'] = self._calc_theory(pars, cutoff=self.cutoff) 117 """118 if self._fn is None:119 q_input = self.model.kernel.make_input(self._kernel_inputs)120 self._fn = self.model.kernel(q_input)121 122 fixed_pars = [getattr(self.model, p).value for p in self._fn.fixed_pars]123 pd_pars = [self._get_weights(p) for p in self._fn.pd_pars]124 #print(fixed_pars,pd_pars)125 Iq_calc = self._fn(fixed_pars, pd_pars, self.cutoff)126 #self._theory[:] = self._fn.eval(pars, pd_pars)127 if self.data_type == 'sesans':128 result = sesans.hankel(self.data.x, self.data.lam * 1e-9,129 self.data.sample.thickness / 10,130 self._kernel_inputs[0], Iq_calc)131 self._cache['theory'] = result132 else:133 Iq = self.resolution.apply(Iq_calc)134 self._cache['theory'] = Iq135 """136 129 return self._cache['theory'] 137 130 … … 162 155 pass 163 156 164 def remove_get_weights(self, name):165 """166 Get parameter dispersion weights167 """168 info = self.model.kernel.info169 relative = name in info['partype']['pd-rel']170 limits = info['limits'][name]171 disperser, value, npts, width, nsigma = [172 getattr(self.model, name + ext)173 for ext in ('_pd_type', '', '_pd_n', '_pd', '_pd_nsigma')]174 value, weight = weights.get_weights(175 disperser, int(npts.value), width.value, nsigma.value,176 value.value, limits, relative)177 return value, weight / np.sum(weight)178 179 157 def __getstate__(self): 180 158 # Can't pickle gpu functions, so instead make them lazy -
sasmodels/compare.py
rf4f3919 rec7e360 66 66 67 67 68 def sasview_model(model_definition, **pars): 69 """ 70 Load a sasview model given the model name. 71 """ 68 def parameter_range(p, v): 69 """ 70 Choose a parameter range based on parameter name and initial value. 71 """ 72 if p.endswith('_pd_n'): 73 return [0, 100] 74 elif p.endswith('_pd_nsigma'): 75 return [0, 5] 76 elif p.endswith('_pd_type'): 77 return v 78 elif any(s in p for s in ('theta','phi','psi')): 79 # orientation in [-180,180], orientation pd in [0,45] 80 if p.endswith('_pd'): 81 return [0,45] 82 else: 83 return [-180, 180] 84 elif 'sld' in p: 85 return [-0.5, 10] 86 elif p.endswith('_pd'): 87 return [0, 1] 88 elif p in ['background', 'scale']: 89 return [0, 1e3] 90 else: 91 return [0, (2*v if v>0 else 1)] 92 93 def _randomize_one(p, v): 94 """ 95 Randomizing parameter. 96 """ 97 if any(p.endswith(s) for s in ('_pd_n','_pd_nsigma','_pd_type')): 98 return v 99 else: 100 return np.random.uniform(*parameter_range(p, v)) 101 102 def randomize_pars(pars, seed=None): 103 np.random.seed(seed) 104 # Note: the sort guarantees order `of calls to random number generator 105 pars = dict((p,_randomize_one(p,v)) 106 for p,v in sorted(pars.items())) 107 return pars 108 109 def constrain_pars(model_definition, pars): 110 """ 111 Restrict parameters to valid values. 112 """ 113 name = model_definition.name 114 if name == 'capped_cylinder' and pars['cap_radius'] < pars['radius']: 115 pars['radius'],pars['cap_radius'] = pars['cap_radius'],pars['radius'] 116 if name == 'barbell' and pars['bell_radius'] < pars['radius']: 117 pars['radius'],pars['bell_radius'] = pars['bell_radius'],pars['radius'] 118 119 # Limit guinier to an Rg such that Iq > 1e-30 (single precision cutoff) 120 if name == 'guinier': 121 #q_max = 0.2 # mid q maximum 122 q_max = 1.0 # high q maximum 123 rg_max = np.sqrt(90*np.log(10) + 3*np.log(pars['scale']))/q_max 124 pars['rg'] = min(pars['rg'],rg_max) 125 126 def parlist(pars): 127 return "\n".join("%s: %s"%(p,v) for p,v in sorted(pars.items())) 128 129 def suppress_pd(pars): 130 """ 131 Suppress theta_pd for now until the normalization is resolved. 132 133 May also suppress complete polydispersity of the model to test 134 models more quickly. 135 """ 136 pars = pars.copy() 137 for p in pars: 138 if p.endswith("_pd"): pars[p] = 0 139 return pars 140 141 def eval_sasview(model_definition, data): 142 # importing sas here so that the error message will be that sas failed to 143 # import rather than the more obscure smear_selection not imported error 144 import sas 145 from sas.models.qsmearing import smear_selection 146 72 147 # convert model parameters from sasmodel form to sasview form 73 148 #print("old",sorted(pars.items())) 74 modelname, pars = revert_model(model_definition, pars)149 modelname, pars = revert_model(model_definition, {}) 75 150 #print("new",sorted(pars.items())) 76 151 sas = __import__('sas.models.'+modelname) … … 79 154 raise ValueError("could not find model %r in sas.models"%modelname) 80 155 model = ModelClass() 81 82 for k,v in pars.items(): 83 parts = k.split('.') # polydispersity components 84 if len(parts) == 2: 85 model.dispersion[parts[0]][parts[1]] = v 156 smearer = smear_selection(data, model=model) 157 158 if hasattr(data, 'qx_data'): 159 q = np.sqrt(data.qx_data**2 + data.qy_data**2) 160 index = ((~data.mask) & (~np.isnan(data.data)) 161 & (q >= data.qmin) & (q <= data.qmax)) 162 if smearer is not None: 163 smearer.model = model # because smear_selection has a bug 164 smearer.accuracy = data.accuracy 165 smearer.set_index(index) 166 theory = lambda: smearer.get_value() 86 167 else: 87 model.setParam(k, v) 88 return model 89 90 def randomize(p, v): 91 """ 92 Randomizing parameter. 93 94 Guess the parameter type from name. 95 """ 96 if any(p.endswith(s) for s in ('_pd_n','_pd_nsigma','_pd_type')): 97 return v 98 elif any(s in p for s in ('theta','phi','psi')): 99 # orientation in [-180,180], orientation pd in [0,45] 100 if p.endswith('_pd'): 101 return 45*np.random.rand() 102 else: 103 return 360*np.random.rand() - 180 104 elif 'sld' in p: 105 # sld in in [-0.5,10] 106 return 10.5*np.random.rand() - 0.5 107 elif p.endswith('_pd'): 108 # length pd in [0,1] 109 return np.random.rand() 110 else: 111 # values from 0 to 2*x for all other parameters 112 return 2*np.random.rand()*(v if v != 0 else 1) 113 114 def randomize_model(pars, seed=None): 115 if seed is None: 116 seed = np.random.randint(1e9) 117 np.random.seed(seed) 118 # Note: the sort guarantees order of calls to random number generator 119 pars = dict((p,randomize(p,v)) for p,v in sorted(pars.items())) 120 121 return pars, seed 122 123 def constrain_pars(model_definition, pars): 124 """ 125 Restrict parameters to valid values. 126 """ 127 name = model_definition.name 128 if name == 'capped_cylinder' and pars['cap_radius'] < pars['radius']: 129 pars['radius'],pars['cap_radius'] = pars['cap_radius'],pars['radius'] 130 if name == 'barbell' and pars['bell_radius'] < pars['radius']: 131 pars['radius'],pars['bell_radius'] = pars['bell_radius'],pars['radius'] 132 133 # Limit guinier to an Rg such that Iq > 1e-30 (single precision cutoff) 134 if name == 'guinier': 135 #q_max = 0.2 # mid q maximum 136 q_max = 1.0 # high q maximum 137 rg_max = np.sqrt(90*np.log(10) + 3*np.log(pars['scale']))/q_max 138 pars['rg'] = min(pars['rg'],rg_max) 139 140 def parlist(pars): 141 return "\n".join("%s: %s"%(p,v) for p,v in sorted(pars.items())) 142 143 def suppress_pd(pars): 144 """ 145 Suppress theta_pd for now until the normalization is resolved. 146 147 May also suppress complete polydispersity of the model to test 148 models more quickly. 149 """ 150 pars = pars.copy() 151 for p in pars: 152 if p.endswith("_pd"): pars[p] = 0 153 return pars 154 155 def eval_sasview(model_definition, pars, data, Nevals=1): 156 # importing sas here so that the error message will be that sas failed to 157 # import rather than the more obscure smear_selection not imported error 158 import sas 159 from sas.models.qsmearing import smear_selection 160 model = sasview_model(model_definition, **pars) 161 smearer = smear_selection(data, model=model) 162 value = None # silence the linter 163 toc = tic() 164 for _ in range(max(Nevals, 1)): # make sure there is at least one eval 165 if hasattr(data, 'qx_data'): 166 q = np.sqrt(data.qx_data**2 + data.qy_data**2) 167 index = ((~data.mask) & (~np.isnan(data.data)) 168 & (q >= data.qmin) & (q <= data.qmax)) 169 if smearer is not None: 170 smearer.model = model # because smear_selection has a bug 171 smearer.accuracy = data.accuracy 172 smearer.set_index(index) 173 value = smearer.get_value() 168 theory = lambda: model.evalDistribution([data.qx_data[index], data.qy_data[index]]) 169 elif smearer is not None: 170 theory = lambda: smearer(model.evalDistribution(data.x)) 171 else: 172 theory = lambda: model.evalDistribution(data.x) 173 174 def calculator(**pars): 175 # paying for parameter conversion each time to keep life simple, if not fast 176 _, pars = revert_model(model_definition, pars) 177 for k,v in pars.items(): 178 parts = k.split('.') # polydispersity components 179 if len(parts) == 2: 180 model.dispersion[parts[0]][parts[1]] = v 174 181 else: 175 value = model.evalDistribution([data.qx_data[index], data.qy_data[index]]) 176 else: 177 value = model.evalDistribution(data.x) 178 if smearer is not None: 179 value = smearer(value) 180 average_time = toc()*1000./Nevals 181 return value, average_time 182 183 def eval_opencl(model_definition, pars, data, dtype='single', Nevals=1, 184 cutoff=0., fast=False): 182 model.setParam(k, v) 183 return theory() 184 185 calculator.engine = "sasview" 186 return calculator 187 188 DTYPE_MAP = { 189 'half': '16', 190 'fast': 'fast', 191 'single': '32', 192 'double': '64', 193 'quad': '128', 194 'f16': '16', 195 'f32': '32', 196 'f64': '64', 197 'longdouble': '128', 198 } 199 def eval_opencl(model_definition, data, dtype='single', cutoff=0.): 185 200 try: 186 model = core.load_model(model_definition, dtype=dtype, 187 platform="ocl", fast=fast) 201 model = core.load_model(model_definition, dtype=dtype, platform="ocl") 188 202 except Exception as exc: 189 203 print(exc) 190 204 print("... trying again with single precision") 191 model = core.load_model(model_definition, dtype='single',192 platform="ocl", fast=fast)205 dtype = 'single' 206 model = core.load_model(model_definition, dtype=dtype, platform="ocl") 193 207 calculator = DirectModel(data, model, cutoff=cutoff) 194 # Run the calculator once before starting the timing loop 208 calculator.engine = "OCL%s"%DTYPE_MAP[dtype] 209 return calculator 210 211 def eval_ctypes(model_definition, data, dtype='double', cutoff=0.): 212 if dtype=='quad': 213 dtype = 'longdouble' 214 model = core.load_model(model_definition, dtype=dtype, platform="dll") 215 calculator = DirectModel(data, model, cutoff=cutoff) 216 calculator.engine = "OMP%s"%DTYPE_MAP[dtype] 217 return calculator 218 219 def time_calculation(calculator, pars, Nevals=1): 220 # initialize the code so time is more accurate 195 221 value = calculator(**suppress_pd(pars)) 196 # Now run the timing loop197 222 toc = tic() 198 for _ in range(max(Nevals, 1)): # forceat least one eval223 for _ in range(max(Nevals, 1)): # make sure there is at least one eval 199 224 value = calculator(**pars) 200 225 average_time = toc()*1000./Nevals 201 226 return value, average_time 202 227 203 204 def eval_ctypes(model_definition, pars, data, dtype='double', Nevals=1, cutoff=0.): 205 model = core.load_model(model_definition, dtype=dtype, platform="dll") 206 calculator = DirectModel(data, model, cutoff=cutoff) 207 # Run the calculator once before starting the timing loop 208 value = calculator(**suppress_pd(pars)) 209 # Now run the timing loop 210 toc = tic() 211 for _ in range(max(Nevals, 1)): # force at least one eval 212 value = calculator(**pars) 213 average_time = toc()*1000./Nevals 214 return value, average_time 215 216 217 def make_data(qmax, is2D, Nq=128, resolution=0.0, accuracy='Low', view='log'): 218 if is2D: 219 data = empty_data2D(np.linspace(-qmax, qmax, Nq), resolution=resolution) 220 data.accuracy = accuracy 228 def make_data(opts): 229 qmax, nq, res = opts['qmax'], opts['nq'], opts['res'] 230 if opts['is2d']: 231 data = empty_data2D(np.linspace(-qmax, qmax, nq), resolution=res) 232 data.accuracy = opts['accuracy'] 221 233 set_beam_stop(data, 0.004) 222 234 index = ~data.mask 223 235 else: 224 if view== 'log':236 if opts['view'] == 'log': 225 237 qmax = math.log10(qmax) 226 q = np.logspace(qmax-3, qmax, Nq)238 q = np.logspace(qmax-3, qmax, nq) 227 239 else: 228 q = np.linspace(0.001*qmax, qmax, Nq)229 data = empty_data1D(q, resolution=res olution)240 q = np.linspace(0.001*qmax, qmax, nq) 241 data = empty_data1D(q, resolution=res) 230 242 index = slice(None, None) 231 243 return data, index 232 244 233 def compare(name, pars, Ncomp, Nbase, opts, set_pars): 234 model_definition = core.load_model_definition(name) 235 236 view = ('linear' if '-linear' in opts 237 else 'log' if '-log' in opts 238 else 'q4' if '-q4' in opts 239 else 'log') 240 241 opt_values = dict(split 242 for s in opts for split in ((s.split('='),)) 243 if len(split) == 2) 244 # Sort out data 245 qmax = (10.0 if '-exq' in opts 246 else 1.0 if '-highq' in opts 247 else 0.2 if '-midq' in opts 248 else 0.05) 249 Nq = int(opt_values.get('-Nq', '128')) 250 res = float(opt_values.get('-res', '0')) 251 accuracy = opt_values.get('-accuracy', 'Low') 252 is2D = "-2d" in opts 253 data, index = make_data(qmax, is2D, Nq, res, accuracy, view=view) 254 255 256 # modelling accuracy is determined by dtype and cutoff 257 dtype = ('longdouble' if '-quad' in opts 258 else 'double' if '-double' in opts 259 else 'half' if '-half' in opts 260 else 'single') 261 cutoff = float(opt_values.get('-cutoff','1e-5')) 262 fast = "-fast" in opts and dtype is 'single' 263 264 # randomize parameters 265 #pars.update(set_pars) # set value before random to control range 266 if '-random' in opts or '-random' in opt_values: 267 seed = int(opt_values['-random']) if '-random' in opt_values else None 268 pars, seed = randomize_model(pars, seed=seed) 269 print("Randomize using -random=%i"%seed) 270 pars.update(set_pars) # set value after random to control value 271 constrain_pars(model_definition, pars) 272 constrain_new_to_old(model_definition, pars) 273 274 # parameter selection 275 if '-mono' in opts: 276 pars = suppress_pd(pars) 277 if '-pars' in opts: 278 print("pars "+str(parlist(pars))) 245 def make_engine(model_definition, data, dtype, cutoff): 246 if dtype == 'sasview': 247 return eval_sasview(model_definition, data) 248 elif dtype.endswith('!'): 249 return eval_ctypes(model_definition, data, dtype=dtype[:-1], 250 cutoff=cutoff) 251 else: 252 return eval_opencl(model_definition, data, dtype=dtype, 253 cutoff=cutoff) 254 255 def compare(opts): 256 Nbase, Ncomp = opts['N1'], opts['N2'] 257 pars = opts['pars'] 258 data = opts['data'] 279 259 280 260 # Base calculation 281 if 0: 282 from sasmodels.models import sphere as target 283 base_name = target.name 284 base, base_time = eval_ctypes(target, pars, data, 285 dtype='longdouble', cutoff=0., Nevals=Ncomp) 286 elif Nbase > 0 and "-ctypes" in opts and "-sasview" in opts: 261 if Nbase > 0: 262 base = opts['engines'][0] 287 263 try: 288 base, base_time = eval_sasview(model_definition, pars, data, Ncomp) 289 base_name = "sasview" 290 #print("base/sasview", (base-pars['background'])/(comp-pars['background'])) 291 print("sasview t=%.1f ms, intensity=%.0f"%(base_time, sum(base))) 292 #print("sasview",comp) 264 base_value, base_time = time_calculation(base, pars, Nbase) 265 print("%s t=%.1f ms, intensity=%.0f"%(base.engine, base_time, sum(base_value))) 293 266 except ImportError: 294 267 traceback.print_exc() 295 268 Nbase = 0 296 elif Nbase > 0:297 base, base_time = eval_opencl(model_definition, pars, data,298 dtype=dtype, cutoff=cutoff, Nevals=Nbase, fast=fast)299 base_name = "ocl"300 print("opencl t=%.1f ms, intensity=%.0f"%(base_time, sum(base)))301 #print("base " + base)302 #print(max(base), min(base))303 269 304 270 # Comparison calculation 305 if Ncomp > 0 and "-ctypes" in opts: 306 comp, comp_time = eval_ctypes(model_definition, pars, data, 307 dtype=dtype, cutoff=cutoff, Nevals=Ncomp) 308 comp_name = "ctypes" 309 print("ctypes t=%.1f ms, intensity=%.0f"%(comp_time, sum(comp))) 310 elif Ncomp > 0: 271 if Ncomp > 0: 272 comp = opts['engines'][1] 311 273 try: 312 comp, comp_time = eval_sasview(model_definition, pars, data, Ncomp) 313 comp_name = "sasview" 314 #print("base/sasview", (base-pars['background'])/(comp-pars['background'])) 315 print("sasview t=%.1f ms, intensity=%.0f"%(comp_time, sum(comp))) 316 #print("sasview",comp) 274 comp_value, comp_time = time_calculation(comp, pars, Ncomp) 275 print("%s t=%.1f ms, intensity=%.0f"%(comp.engine, comp_time, sum(comp_value))) 317 276 except ImportError: 318 277 traceback.print_exc() … … 322 281 if Nbase > 0 and Ncomp > 0: 323 282 #print("speedup %.2g"%(comp_time/base_time)) 324 #print("max |base/comp|", max(abs(base/comp)), "%.15g"%max(abs(base)), "%.15g"%max(abs(comp))) 325 #comp *= max(base/comp) 326 resid = (base - comp) 327 relerr = resid/comp 328 #bad = (relerr>1e-4) 329 #print(relerr[bad],comp[bad],base[bad],data.qx_data[bad],data.qy_data[bad]) 330 _print_stats("|%s-%s|"%(base_name,comp_name)+(" "*(3+len(comp_name))), resid) 331 _print_stats("|(%s-%s)/%s|"%(base_name,comp_name,comp_name), relerr) 283 #print("max |base/comp|", max(abs(base_value/comp_value)), "%.15g"%max(abs(base_value)), "%.15g"%max(abs(comp_value))) 284 #comp *= max(base_value/comp_value) 285 resid = (base_value - comp_value) 286 relerr = resid/comp_value 287 _print_stats("|%s - %s|"%(base.engine,comp.engine)+(" "*(3+len(comp.engine))), resid) 288 _print_stats("|(%s - %s) / %s|"%(base.engine,comp.engine,comp.engine), relerr) 332 289 333 290 # Plot if requested 334 if '-noplot' in opts: return 291 if not opts['plot'] and not opts['explore']: return 292 view = opts['view'] 335 293 import matplotlib.pyplot as plt 294 if Nbase > 0: 295 if Ncomp > 0: plt.subplot(131) 296 plot_theory(data, base_value, view=view, plot_data=False) 297 plt.title("%s t=%.1f ms"%(base.engine, base_time)) 298 #cbar_title = "log I" 336 299 if Ncomp > 0: 337 if Nbase > 0: plt.subplot(131) 338 plot_theory(data, comp, view=view, plot_data=False) 339 plt.title("%s t=%.1f ms"%(comp_name,comp_time)) 340 #cbar_title = "log I" 341 if Nbase > 0: 342 if Ncomp > 0: plt.subplot(132) 343 plot_theory(data, base, view=view, plot_data=False) 344 plt.title("%s t=%.1f ms"%(base_name,base_time)) 300 if Nbase > 0: plt.subplot(132) 301 plot_theory(data, comp_value, view=view, plot_data=False) 302 plt.title("%s t=%.1f ms"%(comp.engine,comp_time)) 345 303 #cbar_title = "log I" 346 304 if Ncomp > 0 and Nbase > 0: … … 363 321 v[v==0] = 0.5*np.min(np.abs(v[v!=0])) 364 322 plt.hist(np.log10(np.abs(v)), normed=1, bins=50); 365 plt.xlabel('log10(err), err = | F(q) single - F(q) double| / | F(q) double |');323 plt.xlabel('log10(err), err = |(%s - %s) / %s|'%(base.engine, comp.engine, comp.engine)); 366 324 plt.ylabel('P(err)') 367 plt.title('Comparison of single and double precision models for %s'%name) 368 369 plt.show() 325 plt.title('Distribution of relative error between calculation engines') 326 327 if not opts['explore']: 328 plt.show() 370 329 371 330 def _print_stats(label, err): … … 387 346 # 388 347 USAGE=""" 389 usage: compare.py model [Nopencl] [Nsasview][options...] [key=val]348 usage: compare.py model N1 N2 [options...] [key=val] 390 349 391 350 Compare the speed and value for a model between the SasView original and the 392 OpenCLrewrite.351 sasmodels rewrite. 393 352 394 353 model is the name of the model to compare (see below). 395 N opencl is the number of times to run the OpenCL model (default=5)396 N sasview is the number of times to run the Sasview model (default=1)354 N1 is the number of times to run sasmodels (default=1). 355 N2 is the number times to run sasview (default=1). 397 356 398 357 Options (* for default): 399 358 400 359 -plot*/-noplot plots or suppress the plot of the model 401 -half/-single*/-double/-quad/-fast sets the calculation precision402 360 -lowq*/-midq/-highq/-exq use q values up to 0.05, 0.2, 1.0, 10.0 403 - Nq=128 sets the number of Q points in the data set361 -nq=128 sets the number of Q points in the data set 404 362 -1d*/-2d computes 1d or 2d data 405 363 -preset*/-random[=seed] preset or random parameters 406 364 -mono/-poly* force monodisperse/polydisperse 407 -ctypes/-sasview* selects gpu:cpu, gpu:sasview, or sasview:cpu if both408 365 -cutoff=1e-5* cutoff value for including a point in polydispersity 409 366 -pars/-nopars* prints the parameter set or not 410 367 -abs/-rel* plot relative or absolute error 411 -linear/-log /-q4 intensity scaling368 -linear/-log*/-q4 intensity scaling 412 369 -hist/-nohist* plot histogram of relative error 413 370 -res=0 sets the resolution width dQ/Q if calculating with resolution 414 371 -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh 415 416 Key=value pairs allow you to set specific values to any of the model 417 parameters. 372 -edit starts the parameter explorer 373 374 Any two calculation engines can be selected for comparison: 375 376 -single/-double/-half/-fast sets an OpenCL calculation engine 377 -single!/-double!/-quad! sets an OpenMP calculation engine 378 -sasview sets the sasview calculation engine 379 380 The default is -single -sasview 381 382 Key=value pairs allow you to set specific values for the model parameters. 418 383 419 384 Available models: … … 423 388 NAME_OPTIONS = set([ 424 389 'plot', 'noplot', 425 'half', 'single', 'double', 'quad', 'fast', 390 'half', 'fast', 'single', 'double', 391 'single!', 'double!', 'quad!', 'sasview', 426 392 'lowq', 'midq', 'highq', 'exq', 427 393 '2d', '1d', 428 394 'preset', 'random', 429 395 'poly', 'mono', 430 'sasview', 'ctypes',431 396 'nopars', 'pars', 432 397 'rel', 'abs', 433 398 'linear', 'log', 'q4', 434 399 'hist', 'nohist', 400 'edit', 435 401 ]) 436 402 VALUE_OPTIONS = [ 437 403 # Note: random is both a name option and a value option 438 'cutoff', 'random', ' Nq', 'res', 'accuracy',404 'cutoff', 'random', 'nq', 'res', 'accuracy', 439 405 ] 440 406 … … 453 419 def get_demo_pars(model_definition): 454 420 info = generate.make_info(model_definition) 421 # Get the default values for the parameters 455 422 pars = dict((p[0],p[2]) for p in info['parameters']) 423 424 # Fill in default values for the polydispersity parameters 425 for p in info['parameters']: 426 if p[4] in ('volume', 'orientation'): 427 pars[p[0]+'_pd'] = 0.0 428 pars[p[0]+'_pd_n'] = 0 429 pars[p[0]+'_pd_nsigma'] = 3.0 430 pars[p[0]+'_pd_type'] = "gaussian" 431 432 # Plug in values given in demo 456 433 pars.update(info['demo']) 457 434 return pars 458 435 459 def main():460 opts = [arg for arg in sys.argv[1:] if arg.startswith('-')]461 popts = [arg for arg in sys.argv[1:] if not arg.startswith('-') and '=' in arg]436 def parse_opts(): 437 flags = [arg for arg in sys.argv[1:] if arg.startswith('-')] 438 values = [arg for arg in sys.argv[1:] if not arg.startswith('-') and '=' in arg] 462 439 args = [arg for arg in sys.argv[1:] if not arg.startswith('-') and '=' not in arg] 463 440 models = "\n ".join("%-15s"%v for v in MODELS) … … 472 449 print("expected parameters: model Nopencl Nsasview") 473 450 474 invalid = [o[1:] for o in opts451 invalid = [o[1:] for o in flags 475 452 if o[1:] not in NAME_OPTIONS 476 453 and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)] … … 479 456 sys.exit(1) 480 457 458 459 # Interpret the flags 460 opts = { 461 'plot' : True, 462 'view' : 'log', 463 'is2d' : False, 464 'qmax' : 0.05, 465 'nq' : 128, 466 'res' : 0.0, 467 'accuracy' : 'Low', 468 'cutoff' : 1e-5, 469 'seed' : -1, # default to preset 470 'mono' : False, 471 'show_pars' : False, 472 'show_hist' : False, 473 'rel_err' : True, 474 'explore' : False, 475 } 476 engines = [] 477 for arg in flags: 478 if arg == '-noplot': opts['plot'] = False 479 elif arg == '-plot': opts['plot'] = True 480 elif arg == '-linear': opts['view'] = 'linear' 481 elif arg == '-log': opts['view'] = 'log' 482 elif arg == '-q4': opts['view'] = 'q4' 483 elif arg == '-1d': opts['is2d'] = False 484 elif arg == '-2d': opts['is2d'] = True 485 elif arg == '-exq': opts['qmax'] = 10.0 486 elif arg == '-highq': opts['qmax'] = 1.0 487 elif arg == '-midq': opts['qmax'] = 0.2 488 elif arg == '-loq': opts['qmax'] = 0.05 489 elif arg.startswith('-nq='): opts['nq'] = int(arg[4:]) 490 elif arg.startswith('-res='): opts['res'] = float(arg[5:]) 491 elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:] 492 elif arg.startswith('-cutoff='): opts['cutoff'] = float(arg[8:]) 493 elif arg.startswith('-random='): opts['seed'] = int(arg[8:]) 494 elif arg == '-random': opts['seed'] = np.random.randint(1e6) 495 elif arg == '-preset': opts['seed'] = -1 496 elif arg == '-mono': opts['mono'] = True 497 elif arg == '-poly': opts['mono'] = False 498 elif arg == '-pars': opts['show_pars'] = True 499 elif arg == '-nopars': opts['show_pars'] = False 500 elif arg == '-hist': opts['show_hist'] = True 501 elif arg == '-nohist': opts['show_hist'] = False 502 elif arg == '-rel': opts['rel_err'] = True 503 elif arg == '-abs': opts['rel_err'] = False 504 elif arg == '-half': engines.append(arg[1:]) 505 elif arg == '-fast': engines.append(arg[1:]) 506 elif arg == '-single': engines.append(arg[1:]) 507 elif arg == '-double': engines.append(arg[1:]) 508 elif arg == '-single!': engines.append(arg[1:]) 509 elif arg == '-double!': engines.append(arg[1:]) 510 elif arg == '-quad!': engines.append(arg[1:]) 511 elif arg == '-sasview': engines.append(arg[1:]) 512 elif arg == '-edit': opts['explore'] = True 513 514 if len(engines) == 0: 515 engines.extend(['single','sasview']) 516 elif len(engines) == 1: 517 if engines[0][0] != 'sasview': 518 engines.append('sasview') 519 else: 520 engines.append('single') 521 elif len(engines) > 2: 522 del engines[2:] 523 524 name = args[0] 525 model_definition = core.load_model_definition(name) 526 527 N1 = int(args[1]) if len(args) > 1 else 1 528 N2 = int(args[2]) if len(args) > 2 else 1 529 481 530 # Get demo parameters from model definition, or use default parameters 482 531 # if model does not define demo parameters 483 name = args[0]484 model_definition = core.load_model_definition(name)485 532 pars = get_demo_pars(model_definition) 486 533 487 Ncomp = int(args[1]) if len(args) > 1 else 5488 Nbase = int(args[2]) if len(args) > 2 else 1489 490 # Fill in default polydispersity parameters491 pds = set(p.split('_pd')[0] for p in pars if p.endswith('_pd'))492 for p in pds:493 if p+"_pd_nsigma" not in pars: pars[p+"_pd_nsigma"] = 3494 if p+"_pd_type" not in pars: pars[p+"_pd_type"] = "gaussian"495 496 534 # Fill in parameters given on the command line 497 set_pars = {}498 for arg in popts:535 presets = {} 536 for arg in values: 499 537 k,v = arg.split('=',1) 500 538 if k not in pars: 501 # extract base name without distribution539 # extract base name without polydispersity info 502 540 s = set(p.split('_pd')[0] for p in pars) 503 541 print("%r invalid; parameters are: %s"%(k,", ".join(sorted(s)))) 504 542 sys.exit(1) 505 set_pars[k] = float(v) if not v.endswith('type') else v 506 507 compare(name, pars, Ncomp, Nbase, opts, set_pars) 543 presets[k] = float(v) if not k.endswith('type') else v 544 545 # randomize parameters 546 #pars.update(set_pars) # set value before random to control range 547 if opts['seed'] > -1: 548 pars = randomize_pars(pars, seed=opts['seed']) 549 print("Randomize using -random=%i"%opts['seed']) 550 pars.update(presets) # set value after random to control value 551 constrain_pars(model_definition, pars) 552 constrain_new_to_old(model_definition, pars) 553 if opts['mono']: 554 pars = suppress_pd(pars) 555 if opts['show_pars']: 556 print("pars " + str(parlist(pars))) 557 558 # Create the computational engines 559 data, _index = make_data(opts) 560 if N1: 561 base = make_engine(model_definition, data, engines[0], opts['cutoff']) 562 else: 563 base = None 564 if N2: 565 comp = make_engine(model_definition, data, engines[1], opts['cutoff']) 566 else: 567 comp = None 568 569 # Remember it all 570 opts.update({ 571 'name' : name, 572 'def' : model_definition, 573 'N1' : N1, 574 'N2' : N2, 575 'presets' : presets, 576 'pars' : pars, 577 'data' : data, 578 'engines' : [base, comp], 579 }) 580 581 return opts 582 583 def main(): 584 opts = parse_opts() 585 if opts['explore']: 586 explore(opts) 587 else: 588 compare(opts) 589 590 def explore(opts): 591 import wx 592 from bumps.names import FitProblem 593 from bumps.gui.app_frame import AppFrame 594 595 problem = FitProblem(Explore(opts)) 596 isMac = "cocoa" in wx.version() 597 app = wx.App() 598 frame = AppFrame(parent=None, title="explore") 599 if not isMac: frame.Show() 600 frame.panel.set_model(model=problem) 601 frame.panel.Layout() 602 frame.panel.aui.Split(0, wx.TOP) 603 if isMac: frame.Show() 604 app.MainLoop() 605 606 class Explore(object): 607 """ 608 Return a bumps wrapper for a SAS model comparison. 609 """ 610 def __init__(self, opts): 611 from bumps.cli import config_matplotlib 612 import bumps_model 613 config_matplotlib() 614 self.opts = opts 615 info = generate.make_info(opts['def']) 616 pars, pd_types = bumps_model.create_parameters(info, **opts['pars']) 617 if not opts['is2d']: 618 active = [base + ext 619 for base in info['partype']['pd-1d'] 620 for ext in ['','_pd','_pd_n','_pd_nsigma']] 621 active.extend(info['partype']['fixed-1d']) 622 for k in active: 623 v = pars[k] 624 v.range(*parameter_range(k, v.value)) 625 else: 626 for k, v in self.pars.items(): 627 v.range(*parameter_range(k, v.value)) 628 629 self.pars = pars 630 self.pd_types = pd_types 631 632 def numpoints(self): 633 """ 634 Return the number of points 635 """ 636 return len(self.pars) + 1 # so dof is 1 637 638 def parameters(self): 639 """ 640 Return a dictionary of parameters 641 """ 642 return self.pars 643 644 def nllf(self): 645 return 0. # No nllf 646 647 def plot(self, view='log'): 648 """ 649 Plot the data and residuals. 650 """ 651 pars = dict((k, v.value) for k,v in self.pars.items()) 652 pars.update(self.pd_types) 653 self.opts['pars'] = pars 654 compare(self.opts) 655 508 656 509 657 if __name__ == "__main__": -
sasmodels/compare_many.py
rf4f3919 rec7e360 7 7 from . import core 8 8 from .kernelcl import environment 9 from .compare import (MODELS, randomize_ model, suppress_pd, eval_sasview,9 from .compare import (MODELS, randomize_pars, suppress_pd, eval_sasview, 10 10 eval_opencl, eval_ctypes, make_data, get_demo_pars, 11 columnize, constrain_pars, constrain_new_to_old) 11 columnize, constrain_pars, constrain_new_to_old, 12 make_engine) 12 13 13 14 def calc_stats(target, value, index): … … 34 35 print(','.join('"%s"'%c for c in columns)) 35 36 37 PRECISION = { 38 'fast': 1e-3, 39 'half': 1e-3, 40 'single': 5e-5, 41 'double': 5e-14, 42 'single!': 5e-5, 43 'double!': 5e-14, 44 'quad!': 5e-18, 45 'sasview': 5e-14, 46 } 36 47 def compare_instance(name, data, index, N=1, mono=True, cutoff=1e-5, 37 precision='double'):48 base='sasview', comp='double'): 38 49 model_definition = core.load_model_definition(name) 39 50 pars = get_demo_pars(model_definition) … … 47 58 # to allow them to update values in the current scope since nonlocal 48 59 # declarations are not available in python 2.7. 49 def try_model(fn, *args, **kw):60 def try_model(fn, pars): 50 61 try: 51 result , _ = fn(model_definition, pars_i, data, *args, **kw)62 result = fn(**pars) 52 63 except KeyboardInterrupt: 53 64 raise … … 60 71 result = np.NaN*data.x 61 72 return result 62 def check_model( label, target, value, acceptable):63 stats = calc_stats(target, value, index)64 co lumns.extend(stats)65 labels.append('GPU single')73 def check_model(pars): 74 base_value = try_model(calc_base, pars) 75 comp_value = try_model(calc_comp, pars) 76 stats = calc_stats(base_value, comp_value, index) 66 77 max_diff[0] = max(max_diff[0], stats[0]) 67 good[0] = good[0] and (stats[0] < acceptable) 78 good[0] = good[0] and (stats[0] < expected) 79 return list(stats) 80 81 82 calc_base = make_engine(model_definition, data, base, cutoff) 83 calc_comp = make_engine(model_definition, data, comp, cutoff) 84 expected = max(PRECISION[base], PRECISION[comp]) 68 85 69 86 num_good = 0 … … 72 89 for k in range(N): 73 90 print("%s %d"%(name, k)) 74 pars_i, seed = randomize_model(pars) 91 seed = np.random.randint(1e6) 92 pars_i = randomize_pars(pars, seed) 75 93 constrain_pars(model_definition, pars_i) 76 94 constrain_new_to_old(model_definition, pars_i) … … 79 97 80 98 good = [True] 81 labels = [] 82 columns = [] 83 target = try_model(eval_sasview) 84 #target = try_model(eval_ctypes, dtype='double', cutoff=0.) 85 #target = try_model(eval_ctypes, dtype='longdouble', cutoff=0.) 86 if precision == 'single': 87 value = try_model(eval_opencl, dtype='single', cutoff=cutoff) 88 check_model('GPU single', target, value, 5e-5) 89 single_value = value # remember for single/double comparison 90 elif precision == 'double': 91 if environment().has_type('double'): 92 label = 'GPU double' 93 value = try_model(eval_opencl, dtype='double', cutoff=cutoff) 94 else: 95 label = 'CPU double' 96 value = try_model(eval_ctypes, dtype='double', cutoff=cutoff) 97 check_model(label, target, value, 5e-14) 98 double_value = value # remember for single/double comparison 99 elif precision == 'quad': 100 value = try_model(eval_opencl, dtype='longdouble', cutoff=cutoff) 101 check_model('CPU quad', target, value, 5e-14) 102 if 0: 103 check_model('single/double', double_value, single_value, 5e-5) 104 99 columns = check_model(pars_i) 105 100 columns += [v for _,v in sorted(pars_i.items())] 106 101 if first: 102 labels = [" vs. ".join((calc_base.engine, calc_comp.engine))] 107 103 print_column_headers(pars_i, labels) 108 104 first = False … … 110 106 num_good += 1 111 107 else: 112 print(("%d,"%seed)+','.join("% g"%v for v in columns))108 print(("%d,"%seed)+','.join("%s"%v for v in columns)) 113 109 print('"good","%d/%d","max diff",%g'%(num_good, N, max_diff[0])) 114 110 … … 144 140 is set in compare.py defaults for each model. 145 141 146 PRECISION is the floating point precision to use for comparisons. 142 PRECISION is the floating point precision to use for comparisons. If two 143 precisions are given, then compare one to the other, ignoring sasview. 147 144 148 145 Available models: … … 151 148 152 149 def main(): 153 if len(sys.argv) != 6:150 if len(sys.argv) not in (6,7): 154 151 print_help() 155 152 sys.exit(1) … … 167 164 mono = sys.argv[4] == 'mono' 168 165 cutoff = float(sys.argv[4]) if not mono else 0 169 precision = sys.argv[5] 166 base = sys.argv[5] 167 comp = sys.argv[6] if len(sys.argv) > 6 else "sasview" 170 168 except: 171 169 traceback.print_exc() … … 173 171 sys.exit(1) 174 172 175 data, index = make_data(qmax=1.0, is2D=is2D, Nq=Nq) 173 data, index = make_data({'qmax':1.0, 'is2d':is2D, 'nq':Nq, 'res':0., 174 'accuracy': 'Low', 'view':'log'}) 176 175 model_list = [model] if model != "all" else MODELS 177 176 for model in model_list: 178 177 compare_instance(model, data, index, N=count, mono=mono, 179 cutoff=cutoff, precision=precision)178 cutoff=cutoff, base=base, comp=comp) 180 179 181 180 if __name__ == "__main__":
Note: See TracChangeset
for help on using the changeset viewer.