Changeset 0bdddc2 in sasmodels


Ignore:
Timestamp:
Jul 28, 2017 8:59:19 PM (7 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
a151caa
Parents:
72be531
Message:

extend sascomp so it can display sets of random models; extend model def to allow random parameter generation

Location:
sasmodels
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/compare.py

    r630156b r0bdddc2  
    264264 
    265265 
    266 def _randomize_one(model_info, p, v): 
     266def _randomize_one(model_info, name, value): 
    267267    # type: (ModelInfo, str, float) -> float 
    268268    # type: (ModelInfo, str, str) -> str 
    … …  
    270270    Randomize a single parameter. 
    271271    """ 
    272     if any(p.endswith(s) for s in ('_pd', '_pd_n', '_pd_nsigma', '_pd_type')): 
    273         return v 
    274  
    275     # Find the parameter definition 
    276     for par in model_info.parameters.call_parameters: 
    277         if par.name == p: 
    278             break 
    279     else: 
    280         raise ValueError("unknown parameter %r"%p) 
    281  
     272    if name.endswith('_pd'): 
     273        par = model_info.parameters[name[:-3]] 
     274        if par.type == 'orientation': 
     275            # Let oriention variation peak around 13 degrees; 95% < 42 degrees 
     276            return 180*np.random.beta(2.5, 20) 
     277        else: 
     278            # Let polydispersity peak around 15%; 95% < 0.4; max=100% 
     279            return np.random.beta(1.5, 7) 
     280 
     281    if name.endswith('_pd_n'): 
     282        # let pd be selected globally rather than per parameter 
     283        return 0 
     284 
     285    if name.endswith('_pd_type'): 
     286        # Don't mess with distribution type for now 
     287        return 'gaussian' 
     288 
     289    if name.endswith('_pd_nsigma'): 
     290        # type-dependent value; for gaussian use 3. 
     291        return 3. 
     292 
     293    if name == 'background': 
     294        return np.random.uniform(0, 1) 
     295 
     296    if name == 'scale': 
     297        return 10**np.random.uniform(-5,0) 
     298 
     299    par = model_info.parameters[name] 
    282300    if len(par.limits) > 2:  # choice list 
    283301        return np.random.randint(len(par.limits)) 
    284302 
    285     limits = par.limits 
    286     if not np.isfinite(limits).all(): 
    287         low, high = parameter_range(p, v) 
    288         limits = (max(limits[0], low), min(limits[1], high)) 
    289  
     303    if np.isfinite(par.limits).all(): 
     304        return np.random.uniform(*par.limits) 
     305 
     306    if par.type == 'sld': 
     307        # Range of neutron SLDs 
     308        return np.random.uniform(-0.5, 12) 
     309 
     310    if par.type == 'volume': 
     311        if ('length' in par.name or 
     312                'radius' in par.name or 
     313                'thick' in par.name): 
     314            return 10**np.random.uniform(2,4) 
     315 
     316    low, high = parameter_range(par.name, value) 
     317    limits = (max(par.limits[0], low), min(par.limits[1], high)) 
    290318    return np.random.uniform(*limits) 
    291319 
    … …  
    301329    :func:`constrain_pars` needs to be called afterward.. 
    302330    """ 
    303     with push_seed(seed): 
    304         # Note: the sort guarantees order `of calls to random number generator 
    305         random_pars = dict((p, _randomize_one(model_info, p, v)) 
    306                            for p, v in sorted(pars.items())) 
     331    # Note: the sort guarantees order of calls to random number generator 
     332    random_pars = dict((p, _randomize_one(model_info, p, v)) 
     333                       for p, v in sorted(pars.items())) 
     334    if model_info.random is not None: 
     335        random_pars.update(model_info.random()) 
     336 
    307337    return random_pars 
    308338 
    … …  
    662692    parameters. 
    663693    """ 
    664     result = run_models(opts, verbose=True) 
    665     if opts['plot']:  # Note: never called from explore 
    666         plot_models(opts, result, limits=limits) 
     694    limits = np.Inf, -np.Inf 
     695    for k in range(opts['sets']): 
     696        opts['pars'] = parse_pars(opts) 
     697        result = run_models(opts, verbose=True) 
     698        if opts['plot']: 
     699            limits = plot_models(opts, result, limits=limits, setnum=k) 
     700    if opts['plot']: 
     701        import matplotlib.pyplot as plt 
     702        plt.show() 
    667703 
    668704def run_models(opts, verbose=False): 
    … …  
    743779 
    744780 
    745 def plot_models(opts, result, limits=None): 
     781def plot_models(opts, result, limits=(np.Inf, -np.Inf), setnum=0): 
    746782    # type: (Dict[str, Any], Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float] 
    747783    base_value, comp_value= result['base_value'], result['comp_value'] 
    … …  
    758794    view = opts['view'] 
    759795    import matplotlib.pyplot as plt 
    760     if limits is None and not use_data: 
    761         vmin, vmax = np.Inf, -np.Inf 
    762         if have_base: 
    763             vmin = min(vmin, base_value.min()) 
    764             vmax = max(vmax, base_value.max()) 
    765         if have_comp: 
    766             vmin = min(vmin, comp_value.min()) 
    767             vmax = max(vmax, comp_value.max()) 
    768         limits = vmin, vmax 
     796    vmin, vmax = limits 
     797    if have_base: 
     798        vmin = min(vmin, base_value.min()) 
     799        vmax = max(vmax, base_value.max()) 
     800    if have_comp: 
     801        vmin = min(vmin, comp_value.min()) 
     802        vmax = max(vmax, comp_value.max()) 
     803    limits = vmin, vmax 
    769804 
    770805    if have_base: 
    … …  
    813848        plt.title('Distribution of relative error between calculation engines') 
    814849 
    815     if not opts['explore']: 
    816         plt.show() 
    817  
    818850    return limits 
    819  
    820  
    821851 
    822852 
    … …  
    841871VALUE_OPTIONS = [ 
    842872    # Note: random is both a name option and a value option 
    843     'cutoff', 'random', 'nq', 'res', 'accuracy', 'title', 'data', 
     873    'cutoff', 'random', 'nq', 'res', 'accuracy', 'title', 'data', 'sets' 
    844874    ] 
    845875 
    … …  
    914944              if not arg.startswith('-') and '=' in arg] 
    915945    positional_args = [arg for arg in argv 
    916             if not arg.startswith('-') and '=' not in arg] 
     946                       if not arg.startswith('-') and '=' not in arg] 
    917947    models = "\n    ".join("%-15s"%v for v in MODELS) 
    918948    if len(positional_args) == 0: 
    … …  
    959989        'title'     : None, 
    960990        'datafile'  : None, 
     991        'sets'      : 1, 
    961992    } 
    962993    engines = [] 
    … …  
    9761007        elif arg.startswith('-nq='):       opts['nq'] = int(arg[4:]) 
    9771008        elif arg.startswith('-res='):      opts['res'] = float(arg[5:]) 
     1009        elif arg.startswith('-sets='):     opts['sets'] = int(arg[6:]) 
    9781010        elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:] 
    9791011        elif arg.startswith('-cutoff='):   opts['cutoff'] = float(arg[8:]) 
    … …  
    10081040    # pylint: enable=bad-whitespace 
    10091041 
     1042    # Force random if more than one set 
     1043    if opts['sets'] > 1 and opts['seed'] < 0: 
     1044        opts['seed'] = np.random.randint(1000000) 
     1045 
    10101046    if MODEL_SPLIT in name: 
    10111047        name, name2 = name.split(MODEL_SPLIT, 2) 
    … …  
    10201056        return None 
    10211057 
    1022     # Get demo parameters from model definition, or use default parameters 
    1023     # if model does not define demo parameters 
    1024     pars = get_pars(model_info, opts['use_demo']) 
    1025     pars2 = get_pars(model_info2, opts['use_demo']) 
    1026     pars2.update((k, v) for k, v in pars.items() if k in pars2) 
    1027     # randomize parameters 
    1028     #pars.update(set_pars)  # set value before random to control range 
    1029     if opts['seed'] > -1: 
    1030         pars = randomize_pars(model_info, pars, seed=opts['seed']) 
    1031         if model_info != model_info2: 
    1032             pars2 = randomize_pars(model_info2, pars2, seed=opts['seed']) 
    1033             # Share values for parameters with the same name 
    1034             for k, v in pars.items(): 
    1035                 if k in pars2: 
    1036                     pars2[k] = v 
    1037         else: 
    1038             pars2 = pars.copy() 
    1039         constrain_pars(model_info, pars) 
    1040         constrain_pars(model_info2, pars2) 
    1041         print("Randomize using -random=%i"%opts['seed']) 
    1042     if opts['mono']: 
    1043         pars = suppress_pd(pars) 
    1044         pars2 = suppress_pd(pars2) 
    1045     if not opts['magnetic']: 
    1046         pars = suppress_magnetism(pars) 
    1047         pars2 = suppress_magnetism(pars2) 
    1048  
    1049     # Fill in parameters given on the command line 
    1050     presets = {} 
    1051     presets2 = {} 
    1052     for arg in values: 
    1053         k, v = arg.split('=', 1) 
    1054         if k not in pars and k not in pars2: 
    1055             # extract base name without polydispersity info 
    1056             s = set(p.split('_pd')[0] for p in pars) 
    1057             print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s)))) 
    1058             return None 
    1059         v1, v2 = v.split(MODEL_SPLIT, 2) if MODEL_SPLIT in v else (v,v) 
    1060         if v1 and k in pars: 
    1061             presets[k] = float(v1) if isnumber(v1) else v1 
    1062         if v2 and k in pars2: 
    1063             presets2[k] = float(v2) if isnumber(v2) else v2 
    1064  
    1065     # If pd given on the command line, default pd_n to 35 
    1066     for k, v in list(presets.items()): 
    1067         if k.endswith('_pd'): 
    1068             presets.setdefault(k+'_n', 35.) 
    1069     for k, v in list(presets2.items()): 
    1070         if k.endswith('_pd'): 
    1071             presets2.setdefault(k+'_n', 35.) 
    1072  
    1073     # Evaluate preset parameter expressions 
    1074     context = MATH.copy() 
    1075     context['np'] = np 
    1076     context.update(pars) 
    1077     context.update((k,v) for k,v in presets.items() if isinstance(v, float)) 
    1078     for k, v in presets.items(): 
    1079         if not isinstance(v, float) and not k.endswith('_type'): 
    1080             presets[k] = eval(v, context) 
    1081     context.update(presets) 
    1082     context.update((k,v) for k,v in presets2.items() if isinstance(v, float)) 
    1083     for k, v in presets2.items(): 
    1084         if not isinstance(v, float) and not k.endswith('_type'): 
    1085             presets2[k] = eval(v, context) 
    1086  
    1087     # update parameters with presets 
    1088     pars.update(presets)  # set value after random to control value 
    1089     pars2.update(presets2)  # set value after random to control value 
    1090     #import pprint; pprint.pprint(model_info) 
    1091  
    1092     same_model = name == name2 and pars == pars 
     1058    # TODO: check if presets are different when deciding if models are same 
     1059    same_model = name == name2 
    10931060    if len(engines) == 0: 
    10941061        if same_model: 
    … …  
    11061073        del engines[2:] 
    11071074 
    1108     use_sasview = any(engine == 'sasview' and count > 0 
    1109                       for engine, count in zip(engines, [n1, n2])) 
    1110     if use_sasview: 
    1111         constrain_new_to_old(model_info, pars) 
    1112         constrain_new_to_old(model_info2, pars2) 
    1113  
    1114     if opts['show_pars']: 
    1115         if not same_model: 
    1116             print("==== %s ====="%model_info.name) 
    1117             print(str(parlist(model_info, pars, opts['is2d']))) 
    1118             print("==== %s ====="%model_info2.name) 
    1119             print(str(parlist(model_info2, pars2, opts['is2d']))) 
    1120         else: 
    1121             print(str(parlist(model_info, pars, opts['is2d']))) 
    1122  
    11231075    # Create the computational engines 
    11241076    if opts['datafile'] is not None: 
    … …  
    11421094        'def'       : [model_info, model_info2], 
    11431095        'count'     : [n1, n2], 
    1144         'presets'   : [presets, presets2], 
    1145         'pars'      : [pars, pars2], 
    11461096        'engines'   : [base, comp], 
     1097        'values'    : values, 
    11471098    }) 
    11481099    # pylint: enable=bad-whitespace 
    11491100 
    11501101    return opts 
     1102 
     1103def parse_pars(opts): 
     1104    model_info, model_info2 = opts['def'] 
     1105 
     1106    # Get demo parameters from model definition, or use default parameters 
     1107    # if model does not define demo parameters 
     1108    pars = get_pars(model_info, opts['use_demo']) 
     1109    pars2 = get_pars(model_info2, opts['use_demo']) 
     1110    pars2.update((k, v) for k, v in pars.items() if k in pars2) 
     1111    # randomize parameters 
     1112    #pars.update(set_pars)  # set value before random to control range 
     1113    if opts['seed'] > -1: 
     1114        pars = randomize_pars(model_info, pars) 
     1115        if model_info != model_info2: 
     1116            pars2 = randomize_pars(model_info2, pars2) 
     1117            # Share values for parameters with the same name 
     1118            for k, v in pars.items(): 
     1119                if k in pars2: 
     1120                    pars2[k] = v 
     1121        else: 
     1122            pars2 = pars.copy() 
     1123        constrain_pars(model_info, pars) 
     1124        constrain_pars(model_info2, pars2) 
     1125    if opts['mono']: 
     1126        pars = suppress_pd(pars) 
     1127        pars2 = suppress_pd(pars2) 
     1128    if not opts['magnetic']: 
     1129        pars = suppress_magnetism(pars) 
     1130        pars2 = suppress_magnetism(pars2) 
     1131 
     1132    # Fill in parameters given on the command line 
     1133    presets = {} 
     1134    presets2 = {} 
     1135    for arg in opts['values']: 
     1136        k, v = arg.split('=', 1) 
     1137        if k not in pars and k not in pars2: 
     1138            # extract base name without polydispersity info 
     1139            s = set(p.split('_pd')[0] for p in pars) 
     1140            print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s)))) 
     1141            return None 
     1142        v1, v2 = v.split(MODEL_SPLIT, 2) if MODEL_SPLIT in v else (v,v) 
     1143        if v1 and k in pars: 
     1144            presets[k] = float(v1) if isnumber(v1) else v1 
     1145        if v2 and k in pars2: 
     1146            presets2[k] = float(v2) if isnumber(v2) else v2 
     1147 
     1148    # If pd given on the command line, default pd_n to 35 
     1149    for k, v in list(presets.items()): 
     1150        if k.endswith('_pd'): 
     1151            presets.setdefault(k+'_n', 35.) 
     1152    for k, v in list(presets2.items()): 
     1153        if k.endswith('_pd'): 
     1154            presets2.setdefault(k+'_n', 35.) 
     1155 
     1156    # Evaluate preset parameter expressions 
     1157    context = MATH.copy() 
     1158    context['np'] = np 
     1159    context.update(pars) 
     1160    context.update((k, v) for k, v in presets.items() if isinstance(v, float)) 
     1161    for k, v in presets.items(): 
     1162        if not isinstance(v, float) and not k.endswith('_type'): 
     1163            presets[k] = eval(v, context) 
     1164    context.update(presets) 
     1165    context.update((k, v) for k, v in presets2.items() if isinstance(v, float)) 
     1166    for k, v in presets2.items(): 
     1167        if not isinstance(v, float) and not k.endswith('_type'): 
     1168            presets2[k] = eval(v, context) 
     1169 
     1170    # update parameters with presets 
     1171    pars.update(presets)  # set value after random to control value 
     1172    pars2.update(presets2)  # set value after random to control value 
     1173    #import pprint; pprint.pprint(model_info) 
     1174 
     1175    if opts['show_pars']: 
     1176        if model_info.name != model_info2.name or pars != pars2: 
     1177            print("==== %s ====="%model_info.name) 
     1178            print(str(parlist(model_info, pars, opts['is2d']))) 
     1179            print("==== %s ====="%model_info2.name) 
     1180            print(str(parlist(model_info2, pars2, opts['is2d']))) 
     1181        else: 
     1182            print(str(parlist(model_info, pars, opts['is2d']))) 
     1183 
     1184    return pars, pars2 
    11511185 
    11521186def show_docs(opts): 
    … …  
    11801214    model = Explore(opts) 
    11811215    problem = FitProblem(model) 
    1182     frame = AppFrame(parent=None, title="explore", size=(1000,700)) 
    1183     if not is_mac: frame.Show() 
     1216    frame = AppFrame(parent=None, title="explore", size=(1000, 700)) 
     1217    if not is_mac: 
     1218        frame.Show() 
    11841219    frame.panel.set_model(model=problem) 
    11851220    frame.panel.Layout() 
    … …  
    11911226    if is_mac: frame.Show() 
    11921227    # If running withing an app, start the main loop 
    1193     if app: app.MainLoop() 
     1228    if app: 
     1229        app.MainLoop() 
    11941230 
    11951231class Explore(object): 
    … …  
    12061242        config_matplotlib() 
    12071243        self.opts = opts 
     1244        opts['pars'] = list(opts['pars']) 
    12081245        p1, p2 = opts['pars'] 
    12091246        m1, m2 = opts['def'] 
    … …  
    12261263        self.starting_values = dict((k, v.value) for k, v in pars.items()) 
    12271264        self.pd_types = pd_types 
    1228         self.limits = None 
     1265        self.limits = np.Inf, -np.Inf 
    12291266 
    12301267    def revert_values(self): 
    … …  
    12821319    """ 
    12831320    opts = parse_opts(argv) 
     1321    if opts['seed'] > -1: 
     1322        print("Randomize using -random=%i"%opts['seed']) 
     1323        np.random.seed(opts['seed']) 
    12841324    if opts is not None: 
    12851325        if opts['html']: 
    12861326            show_docs(opts) 
    12871327        elif opts['explore']: 
     1328            opts['pars'] = parse_pars(opts) 
    12881329            explore(opts) 
    12891330        else: 

  • sasmodels/modelinfo.py

    r724257c r0bdddc2  
    467467                         if p.polydisperse and p.type not in ('orientation', 'magnetic')) 
    468468        self.pd_2d = set(p.name for p in self.call_parameters if p.polydisperse) 
     469 
     470    def __getitem__(self, key): 
     471        # Find the parameter definition 
     472        for par in self.call_parameters: 
     473            if par.name == key: 
     474                break 
     475        else: 
     476            raise KeyError("unknown parameter %r"%key) 
     477        return par 
    469478 
    470479    def _set_vector_lengths(self): 
    … …  
    754763    info.opencl = getattr(kernel_module, 'opencl', not callable(info.Iq)) 
    755764    info.single = getattr(kernel_module, 'single', not callable(info.Iq)) 
     765    info.random = getattr(kernel_module, 'random', None) 
    756766 
    757767    # multiplicity info 

  • sasmodels/models/adsorbed_layer.py

    rb0c4271 r0bdddc2  
    9494Iq.vectorized = True  # Iq accepts an array of q values 
    9595 
     96def random(): 
     97    # only care about the value of second_moment: 
     98    #    curve = scale * e**(-second_moment^2 q^2)/q^2 
     99    #    scale = 6 pi/100 (contrast/density*absorbed_amount)^2 * Vf/radius 
     100    # the remaining parameters can be randomly generated from zero to 
     101    # twice the default value. 
     102    import numpy as np 
     103    pars = dict( 
     104        scale=1, 
     105        second_moment=10**np.random.uniform(1, 3), 
     106    ) 
     107    return pars 
     108 
    96109# unit test values taken from SasView 3.1.2 
    97110tests = [ 

  • sasmodels/models/bcc_paracrystal.py

    r69e1afc r0bdddc2  
    8181.. figure:: img/parallelepiped_angle_definition.png 
    8282 
    83     Orientation of the crystal with respect to the scattering plane, when  
     83    Orientation of the crystal with respect to the scattering plane, when 
    8484    $\theta = \phi = 0$ the $c$ axis is along the beam direction (the $z$ axis). 
    8585 
    … …  
    131131source = ["lib/sas_3j1x_x.c", "lib/gauss150.c", "lib/sphere_form.c", "bcc_paracrystal.c"] 
    132132 
     133def random(): 
     134    import numpy as np 
     135    # Define lattice spacing as a multiple of the particle radius 
     136    # using the formulat a = 4 r/sqrt(3).  Systems which are ordered 
     137    # are probably mostly filled, so use a distribution which goes from 
     138    # zero to one, but leaving 90% of them within 80% of the 
     139    # maximum bcc packing.  Lattice distortion values are empirically 
     140    # useful between 0.01 and 0.7.  Use an exponential distribution 
     141    # in this range 'cuz its easy. 
     142    dnn_fraction = np.random.beta(a=10, b=1) 
     143    pars = dict( 
     144        #sld=1, sld_solvent=0, scale=1, background=1e-32, 
     145        radius=10**np.random.uniform(1.3, 4), 
     146        d_factor=10**np.random.uniform(-2, -0.7),  # sigma_d in 0.01-0.7 
     147    ) 
     148    pars['dnn'] = pars['radius']*4/np.sqrt(3)/dnn_fraction 
     149    #pars['scale'] = 1/(0.68*dnn_fraction**3)  # bcc packing fraction is 0.68 
     150    pars['scale'] = 1 
     151    return pars 
     152 
    133153# parameters for demo 
    134154demo = dict( 

  • sasmodels/models/broad_peak.py

    r43fe34b r0bdddc2  
    9494Iq.vectorized = True  # Iq accepts an array of q values 
    9595 
     96def random(): 
     97    import numpy as np 
     98    pars = dict( 
     99        scale=1, 
     100        porod_scale=10**np.random.uniform(-8, -5), 
     101        porod_exp=np.random.uniform(1, 6), 
     102        lorentz_scale=10**np.random.uniform(0.3, 6), 
     103        lorentz_length=10**np.random.uniform(0, 2), 
     104        peak_pos=10**np.random.uniform(-3, -1), 
     105        lorentz_exp=np.random.uniform(1, 4), 
     106    ) 
     107    pars['lorentz_length'] /= pars['peak_pos'] 
     108    pars['lorentz_scale'] *= pars['porod_scale'] / pars['peak_pos']**pars['porod_exp'] 
     109    #pars['porod_scale'] = 0. 
     110    return pars 
     111 
    96112demo = dict(scale=1, background=0, 
    97113            porod_scale=1.0e-05, porod_exp=3, 
Note: See TracChangeset for help on using the changeset viewer.