Changeset d15a908 in sasmodels
- Timestamp:
- Jan 27, 2016 5:35:15 PM (9 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:
- fdc538a
- Parents:
- 37a7252
- Location:
- sasmodels
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/compare.py
r608e31e rd15a908 72 72 # doc string so that we can display it at run time if there is an error. 73 73 # lin 74 __doc__ = __doc__ + """ 74 __doc__ = (__doc__ # pylint: disable=redefined-builtin 75 + """ 75 76 Program description 76 77 ------------------- 77 78 78 """ + USAGE 79 """ 80 + USAGE) 79 81 80 82 … … 281 283 theory = lambda: smearer.get_value() 282 284 else: 283 theory = lambda: model.evalDistribution([data.qx_data[index], data.qy_data[index]]) 285 theory = lambda: model.evalDistribution([data.qx_data[index], 286 data.qy_data[index]]) 284 287 elif smearer is not None: 285 288 theory = lambda: smearer(model.evalDistribution(data.x)) … … 416 419 try: 417 420 base_value, base_time = time_calculation(base, pars, Nbase) 418 print("%s t=%.1f ms, intensity=%.0f"%(base.engine, base_time, sum(base_value))) 421 print("%s t=%.1f ms, intensity=%.0f" 422 % (base.engine, base_time, sum(base_value))) 419 423 except ImportError: 420 424 traceback.print_exc() … … 426 430 try: 427 431 comp_value, comp_time = time_calculation(comp, pars, Ncomp) 428 print("%s t=%.1f ms, intensity=%.0f"%(comp.engine, comp_time, sum(comp_value))) 432 print("%s t=%.1f ms, intensity=%.0f" 433 % (comp.engine, comp_time, sum(comp_value))) 429 434 except ImportError: 430 435 traceback.print_exc() … … 433 438 # Compare, but only if computing both forms 434 439 if Nbase > 0 and Ncomp > 0: 435 #print("speedup %.2g"%(comp_time/base_time))436 #print("max |base/comp|", max(abs(base_value/comp_value)), "%.15g"%max(abs(base_value)), "%.15g"%max(abs(comp_value)))437 #comp *= max(base_value/comp_value)438 440 resid = (base_value - comp_value) 439 441 relerr = resid/comp_value 440 _print_stats("|%s-%s|"%(base.engine, comp.engine) + (" "*(3+len(comp.engine))), 442 _print_stats("|%s-%s|" 443 % (base.engine, comp.engine) + (" "*(3+len(comp.engine))), 441 444 resid) 442 _print_stats("|(%s-%s)/%s|"%(base.engine, comp.engine, comp.engine), 445 _print_stats("|(%s-%s)/%s|" 446 % (base.engine, comp.engine, comp.engine), 443 447 relerr) 444 448 … … 591 595 invalid = [o[1:] for o in flags 592 596 if o[1:] not in NAME_OPTIONS 593 597 and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)] 594 598 if invalid: 595 599 print("Invalid options: %s"%(", ".join(invalid))) … … 597 601 598 602 603 # pylint: disable=bad-whitespace 599 604 # Interpret the flags 600 605 opts = { … … 651 656 elif arg == '-sasview': engines.append(arg[1:]) 652 657 elif arg == '-edit': opts['explore'] = True 658 # pylint: enable=bad-whitespace 653 659 654 660 if len(engines) == 0: … … 675 681 presets = {} 676 682 for arg in values: 677 k, v = arg.split('=',1)683 k, v = arg.split('=', 1) 678 684 if k not in pars: 679 685 # extract base name without polydispersity info 680 686 s = set(p.split('_pd')[0] for p in pars) 681 print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s))))687 print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s)))) 682 688 sys.exit(1) 683 689 presets[k] = float(v) if not k.endswith('type') else v … … 697 703 698 704 # Create the computational engines 699 data, _ index= make_data(opts)705 data, _ = make_data(opts) 700 706 if n1: 701 707 base = make_engine(model_definition, data, engines[0], opts['cutoff']) … … 707 713 comp = None 708 714 715 # pylint: disable=bad-whitespace 709 716 # Remember it all 710 717 opts.update({ … … 718 725 'engines' : [base, comp], 719 726 }) 727 # pylint: enable=bad-whitespace 720 728 721 729 return opts 722 730 723 def main():724 opts = parse_opts()725 if opts['explore']:726 explore(opts)727 else:728 compare(opts)729 730 731 def explore(opts): 732 """ 733 Explore the model using the Bumps GUI. 734 """ 731 735 import wx 732 736 from bumps.names import FitProblem … … 734 738 735 739 problem = FitProblem(Explore(opts)) 736 is Mac = "cocoa" in wx.version()740 is_mac = "cocoa" in wx.version() 737 741 app = wx.App() 738 742 frame = AppFrame(parent=None, title="explore") 739 if not is Mac: frame.Show()743 if not is_mac: frame.Show() 740 744 frame.panel.set_model(model=problem) 741 745 frame.panel.Layout() 742 746 frame.panel.aui.Split(0, wx.TOP) 743 if is Mac: frame.Show()747 if is_mac: frame.Show() 744 748 app.MainLoop() 745 749 746 750 class Explore(object): 747 751 """ 748 Return a bumps wrapper for a SAS model comparison. 752 Bumps wrapper for a SAS model comparison. 753 754 The resulting object can be used as a Bumps fit problem so that 755 parameters can be adjusted in the GUI, with plots updated on the fly. 749 756 """ 750 757 def __init__(self, opts): … … 787 794 Return cost. 788 795 """ 796 # pylint: disable=no-self-use 789 797 return 0. # No nllf 790 798 … … 804 812 805 813 814 def main(): 815 """ 816 Main program. 817 """ 818 opts = parse_opts() 819 if opts['explore']: 820 explore(opts) 821 else: 822 compare(opts) 823 806 824 if __name__ == "__main__": 807 825 main() -
sasmodels/compare_many.py
r6458608 rd15a908 1 1 #!/usr/bin/env python 2 """ 3 Program to compare results from many random parameter sets for a given model. 4 5 The result is a comma separated value (CSV) table that can be redirected 6 from standard output into a file and loaded into a spreadsheet. 7 8 The models are compared for each parameter set and if the difference is 9 greater than expected for that precision, the parameter set is labeled 10 as bad and written to the output, along with the random seed used to 11 generate that parameter value. This seed can be used with :mod:`compare` 12 to reload and display the details of the model. 13 """ 2 14 from __future__ import print_function 3 15 … … 14 26 15 27 def calc_stats(target, value, index): 28 """ 29 Calculate statistics between the target value and the computed value. 30 31 *target* and *value* are the vectors being compared, with the 32 difference normalized by *target* to get relative error. Only 33 the elements listed in *index* are used, though index may be 34 and empty slice defined by *slice(None, None)*. 35 36 Returns: 37 38 *maxrel* the maximum relative difference 39 40 *rel95* the relative difference with the 5% biggest differences ignored 41 42 *maxabs* the maximum absolute difference for the 5% biggest differences 43 44 *maxval* the maximum value for the 5% biggest differences 45 """ 16 46 resid = abs(value-target)[index] 17 47 relerr = resid/target[index] 18 s rel= np.argsort(relerr)48 sorted_rel_index = np.argsort(relerr) 19 49 #p90 = int(len(relerr)*0.90) 20 50 p95 = int(len(relerr)*0.95) 21 51 maxrel = np.max(relerr) 22 rel95 = relerr[s rel[p95]]23 maxabs = np.max(resid[s rel[p95:]])24 maxval = np.max(value[s rel[p95:]])25 return maxrel, rel95,maxabs,maxval52 rel95 = relerr[sorted_rel_index[p95]] 53 maxabs = np.max(resid[sorted_rel_index[p95:]]) 54 maxval = np.max(value[sorted_rel_index[p95:]]) 55 return maxrel, rel95, maxabs, maxval 26 56 27 57 def print_column_headers(pars, parts): 58 """ 59 Generate column headers for the differences and for the parameters, 60 and print them to standard output. 61 """ 28 62 stats = list('Max rel err|95% rel err|Max abs err above 90% rel|Max value above 90% rel'.split('|')) 29 63 groups = [''] … … 36 70 print(','.join('"%s"'%c for c in columns)) 37 71 72 # Target 'good' value for various precision levels. 38 73 PRECISION = { 39 74 'fast': 1e-3, … … 48 83 def compare_instance(name, data, index, N=1, mono=True, cutoff=1e-5, 49 84 base='sasview', comp='double'): 50 is2D = hasattr(data, 'qx_data') 85 r""" 86 Compare the model under different calculation engines. 87 88 *name* is the name of the model. 89 90 *data* is the data object giving $q, \Delta q$ calculation points. 91 92 *index* is the active set of points. 93 94 *N* is the number of comparisons to make. 95 96 *cutoff* is the polydispersity weight cutoff to make the calculation 97 a little bit faster. 98 99 *base* and *comp* are the names of the calculation engines to compare. 100 """ 101 102 is_2d = hasattr(data, 'qx_data') 51 103 model_definition = core.load_model_definition(name) 52 104 pars = get_demo_pars(model_definition) 53 105 header = ('\n"Model","%s","Count","%d","Dimension","%s"' 54 % (name, N, "2D" if is 2Delse "1D"))106 % (name, N, "2D" if is_2d else "1D")) 55 107 if not mono: header += ',"Cutoff",%g'%(cutoff,) 56 108 print(header) 57 109 58 if is 2D:110 if is_2d: 59 111 info = generate.make_info(model_definition) 60 112 partype = info['partype'] … … 69 121 # declarations are not available in python 2.7. 70 122 def try_model(fn, pars): 123 """ 124 Return the model evaluated at *pars*. If there is an exception, 125 print it and return NaN of the right shape. 126 """ 71 127 try: 72 128 result = fn(**pars) … … 82 138 return result 83 139 def check_model(pars): 140 """ 141 Run the two calculators against *pars*, returning statistics 142 on the differences. See :func:`calc_stats` for the list of stats. 143 """ 84 144 base_value = try_model(calc_base, pars) 85 145 comp_value = try_model(calc_comp, pars) … … 108 168 good = [True] 109 169 columns = check_model(pars_i) 110 columns += [v for _, v in sorted(pars_i.items())]170 columns += [v for _, v in sorted(pars_i.items())] 111 171 if first: 112 172 labels = [" vs. ".join((calc_base.engine, calc_comp.engine))] … … 121 181 122 182 def print_usage(): 183 """ 184 Print the command usage string. 185 """ 123 186 print("usage: compare_many.py MODEL COUNT (1dNQ|2dNQ) (CUTOFF|mono) (single|double|quad)") 124 187 125 188 126 189 def print_models(): 190 """ 191 Print the list of available models in columns. 192 """ 127 193 print(columnize(MODELS, indent=" ")) 128 194 129 195 130 196 def print_help(): 197 """ 198 Print usage string, the option description and the list of available models. 199 """ 131 200 print_usage() 132 201 print("""\ … … 158 227 159 228 def main(): 160 if len(sys.argv) not in (6,7): 229 """ 230 Main program. 231 """ 232 if len(sys.argv) not in (6, 7): 161 233 print_help() 162 234 sys.exit(1) … … 182 254 183 255 data, index = make_data({'qmax':1.0, 'is2d':is2D, 'nq':Nq, 'res':0., 184 256 'accuracy': 'Low', 'view':'log'}) 185 257 model_list = [model] if model != "all" else MODELS 186 258 for model in model_list: -
sasmodels/convert.py
r29da213 rd15a908 49 49 """ 50 50 return dict((p, (v*1e6 if p.endswith('sld') else v)) 51 for p, v in pars.items())51 for p, v in pars.items()) 52 52 53 53 def convert_model(name, pars): … … 55 55 Convert model from old style parameter names to new style. 56 56 """ 57 _, _ = name,pars # lint57 _, _ = name, pars # lint 58 58 raise NotImplementedError 59 59 # need to load all new models in order to determine old=>new … … 67 67 """ 68 68 return dict((p, (v*1e-6 if p.endswith('sld') else v)) 69 for p, v in pars.items())69 for p, v in pars.items()) 70 70 71 71 def _remove_pd(pars, key, name): -
sasmodels/core.py
rcde11f0f rd15a908 2 2 Core model handling routines. 3 3 """ 4 __all__ = ["list_models", "load_model_definition", "precompile_dll", 5 "load_model", "make_kernel", "call_kernel", "call_ER", "call_VR" ] 4 __all__ = [ 5 "list_models", "load_model_definition", "precompile_dll", 6 "load_model", "make_kernel", "call_kernel", "call_ER", "call_VR", 7 ] 6 8 7 9 from os.path import basename, dirname, join as joinpath … … 61 63 62 64 def isstr(s): 65 """ 66 Return True if *s* is a string-like object. 67 """ 63 68 try: s + '' 64 69 except: return False … … 99 104 # source = open(info['name']+'.cl','r').read() 100 105 101 if (platform =="dll"106 if (platform == "dll" 102 107 or not HAVE_OPENCL 103 108 or not kernelcl.environment().has_type(dtype)): … … 128 133 width = pars.get(name+'_pd', 0.0) 129 134 nsigma = pars.get(name+'_pd_nsigma', 3.0) 130 value, weight = weights.get_weights(135 value, weight = weights.get_weights( 131 136 disperser, npts, width, nsigma, value, limits, relative) 132 137 return value, weight / np.sum(weight) … … 195 200 for name in info['partype']['volume']] 196 201 value, weight = dispersion_mesh(vol_pars) 197 whole, part = VR(*value)202 whole, part = VR(*value) 198 203 return np.sum(weight*part)/np.sum(weight*whole) 199 204 -
sasmodels/data.py
r013adb7 rd15a908 284 284 mdata[mdata <= 0] = masked 285 285 plt.errorbar(data.x/10, scale*mdata, yerr=data.dy, fmt='.') 286 all_positive = all_positive and (mdata >0).all()286 all_positive = all_positive and (mdata > 0).all() 287 287 some_present = some_present or (mdata.count() > 0) 288 288 … … 292 292 mtheory[~np.isfinite(mtheory)] = masked 293 293 if view is 'log': 294 mtheory[mtheory <=0] = masked294 mtheory[mtheory <= 0] = masked 295 295 plt.plot(data.x/10, scale*mtheory, '-', hold=True) 296 all_positive = all_positive and (mtheory >0).all()296 all_positive = all_positive and (mtheory > 0).all() 297 297 some_present = some_present or (mtheory.count() > 0) 298 298 … … 396 396 plt.title('theory') 397 397 h = plt.colorbar() 398 h.set_label(r'$\log_{10}I(q)$' if view =='log'398 h.set_label(r'$\log_{10}I(q)$' if view == 'log' 399 399 else r'$q^4 I(q)$' if view == 'q4' 400 400 else '$I(q)$') … … 411 411 plt.title('residuals') 412 412 h = plt.colorbar() 413 h.set_label( '$\Delta I(q)$')413 h.set_label(r'$\Delta I(q)$') 414 414 415 415
Note: See TracChangeset
for help on using the changeset viewer.