Changeset 7609046 in sasmodels for sasmodels/compare.py


Ignore:
Timestamp:
Jul 19, 2018 12:47:44 PM (11 months ago)
Author:
GitHub <noreply@…>
Branches:
master
Children:
f41027b
Parents:
9ce5bcb (diff), 1d9998c (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
git-author:
Paul Kienzle <pkienzle@…> (07/19/18 12:47:44)
git-committer:
GitHub <noreply@…> (07/19/18 12:47:44)
Message:

Merge branch 'master' into ticket-608-user-defined-weights

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/compare.py

    raa25fc7 r7609046  
    107107    -title="note" adds note to the plot title, after the model name 
    108108    -weights shows weights plots for the polydisperse parameters 
     109    -profile shows the sld profile if the model has a plottable sld profile 
    109110 
    110111    === output options === 
    111112    -edit starts the parameter explorer 
    112113    -help/-html shows the model docs instead of running the model 
     114 
     115    === environment variables === 
     116    -DSAS_MODELPATH=path sets directory containing custom models 
     117    -DSAS_OPENCL=vendor:device|none sets the target OpenCL device 
     118    -DXDG_CACHE_HOME=~/.cache sets the pyopencl cache root (linux only) 
     119    -DSAS_COMPILER=tinycc|msvc|mingw|unix sets the DLL compiler 
     120    -DSAS_OPENMP=1 turns on OpenMP for the DLLs 
     121    -DSAS_DLL_PATH=path sets the path to the compiled modules 
    113122 
    114123The interpretation of quad precision depends on architecture, and may 
     
    645654 
    646655def make_data(opts): 
    647     # type: (Dict[str, Any]) -> Tuple[Data, np.ndarray] 
     656    # type: (Dict[str, Any], float) -> Tuple[Data, np.ndarray] 
    648657    """ 
    649658    Generate an empty dataset, used with the model to set Q points 
     
    667676        if opts['zero']: 
    668677            q = np.hstack((0, q)) 
    669         data = empty_data1D(q, resolution=res) 
     678        # TODO: provide command line control of lambda and Delta lambda/lambda 
     679        #L, dLoL = 5, 0.14/np.sqrt(6)  # wavelength and 14% triangular FWHM 
     680        L, dLoL = 0, 0 
     681        data = empty_data1D(q, resolution=res, L=L, dL=L*dLoL) 
    670682        index = slice(None, None) 
    671683    return data, index 
     
    772784            dim = base._kernel.dim 
    773785            weights.plot_weights(model_info, get_mesh(model_info, base_pars, dim=dim)) 
     786        if opts['show_profile']: 
     787            import pylab 
     788            base, comp = opts['engines'] 
     789            base_pars, comp_pars = opts['pars'] 
     790            have_base = base._kernel.info.profile is not None 
     791            have_comp = ( 
     792                comp is not None 
     793                and comp._kernel.info.profile is not None 
     794                and base_pars != comp_pars 
     795            ) 
     796            if have_base or have_comp: 
     797                pylab.figure() 
     798            if have_base: 
     799                plot_profile(base._kernel.info, **base_pars) 
     800            if have_comp: 
     801                plot_profile(comp._kernel.info, label='comp', **comp_pars) 
     802                pylab.legend() 
    774803    if opts['plot']: 
    775804        import matplotlib.pyplot as plt 
     
    777806    return limits 
    778807 
     808def plot_profile(model_info, label='base', **args): 
     809    # type: (ModelInfo, List[Tuple[float, np.ndarray, np.ndarray]]) -> None 
     810    """ 
     811    Plot the profile returned by the model profile method. 
     812 
     813    *model_info* defines model parameters, etc. 
     814 
     815    *mesh* is a list of tuples containing (*value*, *dispersity*, *weights*) 
     816    for each parameter, where (*dispersity*, *weights*) pairs are the 
     817    distributions to be plotted. 
     818    """ 
     819    import pylab 
     820 
     821    args = dict((k, v) for k, v in args.items() 
     822                if "_pd" not in k 
     823                   and ":" not in k 
     824                   and k not in ("background", "scale", "theta", "phi", "psi")) 
     825    args = args.copy() 
     826 
     827    args.pop('scale', 1.) 
     828    args.pop('background', 0.) 
     829    z, rho = model_info.profile(**args) 
     830    #pylab.interactive(True) 
     831    pylab.plot(z, rho, '-', label=label) 
     832    pylab.grid(True) 
     833    #pylab.show() 
     834 
     835 
     836 
    779837def run_models(opts, verbose=False): 
    780838    # type: (Dict[str, Any]) -> Dict[str, Any] 
     
    786844    base_n, comp_n = opts['count'] 
    787845    base_pars, comp_pars = opts['pars'] 
    788     data = opts['data'] 
     846    base_data, comp_data = opts['data'] 
    789847 
    790848    comparison = comp is not None 
     
    800858            print("%s t=%.2f ms, intensity=%.0f" 
    801859                  % (base.engine, base_time, base_value.sum())) 
    802         _show_invalid(data, base_value) 
     860        _show_invalid(base_data, base_value) 
    803861    except ImportError: 
    804862        traceback.print_exc() 
     
    812870                print("%s t=%.2f ms, intensity=%.0f" 
    813871                      % (comp.engine, comp_time, comp_value.sum())) 
    814             _show_invalid(data, comp_value) 
     872            _show_invalid(base_data, comp_value) 
    815873        except ImportError: 
    816874            traceback.print_exc() 
     
    866924    have_base, have_comp = (base_value is not None), (comp_value is not None) 
    867925    base, comp = opts['engines'] 
    868     data = opts['data'] 
     926    base_data, comp_data = opts['data'] 
    869927    use_data = (opts['datafile'] is not None) and (have_base ^ have_comp) 
    870928 
    871929    # Plot if requested 
    872930    view = opts['view'] 
     931    #view = 'log' 
    873932    if limits is None: 
    874933        vmin, vmax = np.inf, -np.inf 
     
    884943        if have_comp: 
    885944            plt.subplot(131) 
    886         plot_theory(data, base_value, view=view, use_data=use_data, limits=limits) 
     945        plot_theory(base_data, base_value, view=view, use_data=use_data, limits=limits) 
    887946        plt.title("%s t=%.2f ms"%(base.engine, base_time)) 
    888947        #cbar_title = "log I" 
     
    891950            plt.subplot(132) 
    892951        if not opts['is2d'] and have_base: 
    893             plot_theory(data, base_value, view=view, use_data=use_data, limits=limits) 
    894         plot_theory(data, comp_value, view=view, use_data=use_data, limits=limits) 
     952            plot_theory(comp_data, base_value, view=view, use_data=use_data, limits=limits) 
     953        plot_theory(comp_data, comp_value, view=view, use_data=use_data, limits=limits) 
    895954        plt.title("%s t=%.2f ms"%(comp.engine, comp_time)) 
    896955        #cbar_title = "log I" 
     
    908967            err[err > cutoff] = cutoff 
    909968        #err,errstr = base/comp,"ratio" 
    910         plot_theory(data, None, resid=err, view=errview, use_data=use_data) 
     969        # Note: base_data only since base and comp have same q values (though 
     970        # perhaps different resolution), and we are plotting the difference 
     971        # at each q 
     972        plot_theory(base_data, None, resid=err, view=errview, use_data=use_data) 
    911973        plt.xscale('log' if view == 'log' and not opts['is2d'] else 'linear') 
    912974        plt.legend(['P%d'%(k+1) for k in range(setnum+1)], loc='best') 
     
    9421004OPTIONS = [ 
    9431005    # Plotting 
    944     'plot', 'noplot', 'weights', 
     1006    'plot', 'noplot', 
     1007    'weights', 'profile', 
    9451008    'linear', 'log', 'q4', 
    9461009    'rel', 'abs', 
     
    10581121 
    10591122    invalid = [o[1:] for o in flags 
    1060                if o[1:] not in NAME_OPTIONS 
    1061                and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)] 
     1123               if not (o[1:] in NAME_OPTIONS 
     1124                       or any(o.startswith('-%s='%t) for t in VALUE_OPTIONS) 
     1125                       or o.startswith('-D'))] 
    10621126    if invalid: 
    10631127        print("Invalid options: %s"%(", ".join(invalid))) 
     
    10751139        'qmax'      : 0.05, 
    10761140        'nq'        : 128, 
    1077         'res'       : 0.0, 
     1141        'res'       : '0.0', 
    10781142        'noise'     : 0.0, 
    10791143        'accuracy'  : 'Low', 
     
    10961160        'count'     : '1', 
    10971161        'show_weights' : False, 
     1162        'show_profile' : False, 
    10981163        'sphere'    : 0, 
    10991164        'ngauss'    : '0', 
     
    11151180        elif arg.startswith('-q='): 
    11161181            opts['qmin'], opts['qmax'] = [float(v) for v in arg[3:].split(':')] 
    1117         elif arg.startswith('-res='):      opts['res'] = float(arg[5:]) 
     1182        elif arg.startswith('-res='):      opts['res'] = arg[5:] 
    11181183        elif arg.startswith('-noise='):    opts['noise'] = float(arg[7:]) 
    11191184        elif arg.startswith('-sets='):     opts['sets'] = int(arg[6:]) 
     
    11561221        elif arg == '-default': opts['use_demo'] = False 
    11571222        elif arg == '-weights': opts['show_weights'] = True 
     1223        elif arg == '-profile': opts['show_profile'] = True 
    11581224        elif arg == '-html':    opts['html'] = True 
    11591225        elif arg == '-help':    opts['html'] = True 
     1226        elif arg.startswith('-D'): 
     1227            var, val = arg[2:].split('=') 
     1228            os.environ[var] = val 
    11601229    # pylint: enable=bad-whitespace,C0321 
    11611230 
     
    11731242    if opts['qmin'] is None: 
    11741243        opts['qmin'] = 0.001*opts['qmax'] 
    1175     if opts['datafile'] is not None: 
    1176         data = load_data(os.path.expanduser(opts['datafile'])) 
    1177     else: 
    1178         data, _ = make_data(opts) 
    11791244 
    11801245    comparison = any(PAR_SPLIT in v for v in values) 
     
    12161281        opts['cutoff'] = [float(opts['cutoff'])]*2 
    12171282 
    1218     base = make_engine(model_info[0], data, opts['engine'][0], 
     1283    if PAR_SPLIT in opts['res']: 
     1284        opts['res'] = [float(k) for k in opts['res'].split(PAR_SPLIT, 2)] 
     1285        comparison = True 
     1286    else: 
     1287        opts['res'] = [float(opts['res'])]*2 
     1288 
     1289    if opts['datafile'] is not None: 
     1290        data = load_data(os.path.expanduser(opts['datafile'])) 
     1291    else: 
     1292        # Hack around the fact that make_data doesn't take a pair of resolutions 
     1293        res = opts['res'] 
     1294        opts['res'] = res[0] 
     1295        data0, _ = make_data(opts) 
     1296        if res[0] != res[1]: 
     1297            opts['res'] = res[1] 
     1298            data1, _ = make_data(opts) 
     1299        else: 
     1300            data1 = data0 
     1301        opts['res'] = res 
     1302        data = data0, data1 
     1303 
     1304    base = make_engine(model_info[0], data[0], opts['engine'][0], 
    12191305                       opts['cutoff'][0], opts['ngauss'][0]) 
    12201306    if comparison: 
    1221         comp = make_engine(model_info[1], data, opts['engine'][1], 
     1307        comp = make_engine(model_info[1], data[1], opts['engine'][1], 
    12221308                           opts['cutoff'][1], opts['ngauss'][1]) 
    12231309    else: 
Note: See TracChangeset for help on using the changeset viewer.