Changes in / [31ae428:a85a569] in sasmodels


Ignore:
Location:
sasmodels
Files:
80 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/compare.py

    r630156b rd9ec8f9  
    5656 
    5757USAGE = """ 
    58 usage: compare.py model N1 N2 [options...] [key=val] 
    59  
    60 Compare the speed and value for a model between the SasView original and the 
    61 sasmodels rewrite. 
     58usage: sascomp model [options...] [key=val] 
     59 
     60Generate and compare SAS models.  If a single model is specified it shows 
     61a plot of that model.  Different models can be compared, or the same model 
     62with different parameters.  The same model with the same parameters can 
     63be compared with different calculation engines to see the effects of precision 
     64on the resultant values. 
    6265 
    6366model or model1,model2 are the names of the models to compare (see below). 
    64 N1 is the number of times to run sasmodels (default=1). 
    65 N2 is the number times to run sasview (default=1). 
    6667 
    6768Options (* for default): 
    6869 
    69     -plot*/-noplot plots or suppress the plot of the model 
     70    === data generation === 
     71    -data="path" uses q, dq from the data file 
     72    -noise=0 sets the measurement error dI/I 
     73    -res=0 sets the resolution width dQ/Q if calculating with resolution 
    7074    -lowq*/-midq/-highq/-exq use q values up to 0.05, 0.2, 1.0, 10.0 
    7175    -nq=128 sets the number of Q points in the data set 
     76    -1d*/-2d computes 1d or 2d data 
    7277    -zero indicates that q=0 should be included 
    73     -1d*/-2d computes 1d or 2d data 
     78 
     79    === model parameters === 
    7480    -preset*/-random[=seed] preset or random parameters 
     81    -sets=n generates n random datasets with the seed given by -random=seed 
     82    -pars/-nopars* prints the parameter set or not 
     83    -default/-demo* use demo vs default parameters 
     84 
     85    === calculation options === 
    7586    -mono*/-poly force monodisperse or allow polydisperse demo parameters 
     87    -cutoff=1e-5* cutoff value for including a point in polydispersity 
    7688    -magnetic/-nonmagnetic* suppress magnetism 
    77     -cutoff=1e-5* cutoff value for including a point in polydispersity 
    78     -pars/-nopars* prints the parameter set or not 
    79     -abs/-rel* plot relative or absolute error 
    80     -linear/-log*/-q4 intensity scaling 
    81     -hist/-nohist* plot histogram of relative error 
    82     -res=0 sets the resolution width dQ/Q if calculating with resolution 
    8389    -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh 
    84     -edit starts the parameter explorer 
    85     -default/-demo* use demo vs default parameters 
    86     -help/-html shows the model docs instead of running the model 
    87     -title="note" adds note to the plot title, after the model name 
    88     -data="path" uses q, dq from the data file 
    89  
    90 Any two calculation engines can be selected for comparison: 
    91  
     90 
     91    === precision options === 
     92    -calc=default uses the default calcution precision 
    9293    -single/-double/-half/-fast sets an OpenCL calculation engine 
    9394    -single!/-double!/-quad! sets an OpenMP calculation engine 
    9495    -sasview sets the sasview calculation engine 
    9596 
    96 The default is -single -double.  Note that the interpretation of quad 
    97 precision depends on architecture, and may vary from 64-bit to 128-bit, 
    98 with 80-bit floats being common (1e-19 precision). 
     97    === plotting === 
     98    -plot*/-noplot plots or suppress the plot of the model 
     99    -linear/-log*/-q4 intensity scaling on plots 
     100    -hist/-nohist* plot histogram of relative error 
     101    -abs/-rel* plot relative or absolute error 
     102    -title="note" adds note to the plot title, after the model name 
     103 
     104    === output options === 
     105    -edit starts the parameter explorer 
     106    -help/-html shows the model docs instead of running the model 
     107 
     108The interpretation of quad precision depends on architecture, and may 
     109vary from 64-bit to 128-bit, with 80-bit floats being common (1e-19 precision). 
     110On unix and mac you may need single quotes around the DLL computation 
     111engines, such as -calc='single!,double!' since !, is treated as a history 
     112expansion request in the shell. 
    99113 
    100114Key=value pairs allow you to set specific values for the model parameters. 
    101 Key=value1,value2 to compare different values of the same parameter. 
    102 value can be an expression including other parameters 
     115Key=value1,value2 to compare different values of the same parameter. The 
     116value can be an expression including other parameters. 
     117 
     118Items later on the command line override those that appear earlier. 
     119 
     120Examples: 
     121 
     122    # compare single and double precision calculation for a barbell 
     123    sascomp barbell -calc=single,double 
     124 
     125    # generate 10 random lorentz models, with seed=27 
     126    sascomp lorentz -sets=10 -seed=27 
     127 
     128    # compare ellipsoid with R = R_polar = R_equatorial to sphere of radius R 
     129    sascomp sphere,ellipsoid radius_polar=radius radius_equatorial=radius 
     130 
     131    # model timing test requires multiple evals to perform the estimate 
     132    sascomp pringle -calc=single,double -timing=100,100 -noplot 
    103133""" 
    104134 
    … …  
    111141------------------- 
    112142 
    113 """ 
    114            + USAGE) 
     143""" + USAGE) 
    115144 
    116145kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True 
    … …  
    264293 
    265294 
    266 def _randomize_one(model_info, p, v): 
     295def _randomize_one(model_info, name, value): 
    267296    # type: (ModelInfo, str, float) -> float 
    268297    # type: (ModelInfo, str, str) -> str 
    … …  
    270299    Randomize a single parameter. 
    271300    """ 
    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  
     301    # Set the amount of polydispersity/angular dispersion, but by default pd_n 
     302    # is zero so there is no polydispersity.  This allows us to turn on/off 
     303    # pd by setting pd_n, and still have randomly generated values 
     304    if name.endswith('_pd'): 
     305        par = model_info.parameters[name[:-3]] 
     306        if par.type == 'orientation': 
     307            # Let oriention variation peak around 13 degrees; 95% < 42 degrees 
     308            return 180*np.random.beta(2.5, 20) 
     309        else: 
     310            # Let polydispersity peak around 15%; 95% < 0.4; max=100% 
     311            return np.random.beta(1.5, 7) 
     312 
     313    # pd is selected globally rather than per parameter, so set to 0 for no pd 
     314    # In particular, when multiple pd dimensions, want to decrease the number 
     315    # of points per dimension for faster computation 
     316    if name.endswith('_pd_n'): 
     317        return 0 
     318 
     319    # Don't mess with distribution type for now 
     320    if name.endswith('_pd_type'): 
     321        return 'gaussian' 
     322 
     323    # type-dependent value of number of sigmas; for gaussian use 3. 
     324    if name.endswith('_pd_nsigma'): 
     325        return 3. 
     326 
     327    # background in the range [0.01, 1] 
     328    if name == 'background': 
     329        return 10**np.random.uniform(-2, 0) 
     330 
     331    # scale defaults to 0.1% to 30% volume fraction 
     332    if name == 'scale': 
     333        return 10**np.random.uniform(-3, -0.5) 
     334 
     335    # If it is a list of choices, pick one at random with equal probability 
     336    # In practice, the model specific random generator will override. 
     337    par = model_info.parameters[name] 
    282338    if len(par.limits) > 2:  # choice list 
    283339        return np.random.randint(len(par.limits)) 
    284340 
    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  
     341    # If it is a fixed range, pick from it with equal probability. 
     342    # For logarithmic ranges, the model will have to override. 
     343    if np.isfinite(par.limits).all(): 
     344        return np.random.uniform(*par.limits) 
     345 
     346    # If the paramter is marked as an sld use the range of neutron slds 
     347    # TODO: ought to randomly contrast match a pair of SLDs 
     348    if par.type == 'sld': 
     349        return np.random.uniform(-0.5, 12) 
     350 
     351    # Limit magnetic SLDs to a smaller range, from zero to iron=5/A^2 
     352    if par.name.startswith('M0:'): 
     353        return np.random.uniform(0, 5) 
     354 
     355    # Guess at the random length/radius/thickness.  In practice, all models 
     356    # are going to set their own reasonable ranges. 
     357    if par.type == 'volume': 
     358        if ('length' in par.name or 
     359                'radius' in par.name or 
     360                'thick' in par.name): 
     361            return 10**np.random.uniform(2, 4) 
     362 
     363    # In the absence of any other info, select a value in [0, 2v], or 
     364    # [-2|v|, 2|v|] if v is negative, or [0, 1] if v is zero.  Mostly the 
     365    # model random parameter generators will override this default. 
     366    low, high = parameter_range(par.name, value) 
     367    limits = (max(par.limits[0], low), min(par.limits[1], high)) 
    290368    return np.random.uniform(*limits) 
    291369 
    292  
    293 def randomize_pars(model_info, pars, seed=None): 
    294     # type: (ModelInfo, ParameterSet, int) -> ParameterSet 
     370def _random_pd(model_info, pars): 
     371    pd = [p for p in model_info.parameters.kernel_parameters if p.polydisperse] 
     372    pd_volume = [] 
     373    pd_oriented = [] 
     374    for p in pd: 
     375        if p.type == 'orientation': 
     376            pd_oriented.append(p.name) 
     377        elif p.length_control is not None: 
     378            n = int(pars.get(p.length_control, 1) + 0.5) 
     379            pd_volume.extend(p.name+str(k+1) for k in range(n)) 
     380        elif p.length > 1: 
     381            pd_volume.extend(p.name+str(k+1) for k in range(p.length)) 
     382        else: 
     383            pd_volume.append(p.name) 
     384    u = np.random.rand() 
     385    n = len(pd_volume) 
     386    if u < 0.01 or n < 1: 
     387        pass  # 1% chance of no polydispersity 
     388    elif u < 0.86 or n < 2: 
     389        pars[np.random.choice(pd_volume)+"_pd_n"] = 35 
     390    elif u < 0.99 or n < 3: 
     391        choices = np.random.choice(len(pd_volume), size=2) 
     392        pars[pd_volume[choices[0]]+"_pd_n"] = 25 
     393        pars[pd_volume[choices[1]]+"_pd_n"] = 10 
     394    else: 
     395        choices = np.random.choice(len(pd_volume), size=3) 
     396        pars[pd_volume[choices[0]]+"_pd_n"] = 25 
     397        pars[pd_volume[choices[1]]+"_pd_n"] = 10 
     398        pars[pd_volume[choices[2]]+"_pd_n"] = 5 
     399    if pd_oriented: 
     400        pars['theta_pd_n'] = 20 
     401        if np.random.rand() < 0.1: 
     402            pars['phi_pd_n'] = 5 
     403        if np.random.rand() < 0.1: 
     404            if any(p.name == 'psi' for p in model_info.parameters.kernel_parameters): 
     405                #print("generating psi_pd_n") 
     406                pars['psi_pd_n'] = 5 
     407 
     408    ## Show selected polydispersity 
     409    #for name, value in pars.items(): 
     410    #    if name.endswith('_pd_n') and value > 0: 
     411    #        print(name, value, pars.get(name[:-5], 0), pars.get(name[:-2], 0)) 
     412 
     413 
     414def randomize_pars(model_info, pars): 
     415    # type: (ModelInfo, ParameterSet) -> ParameterSet 
    295416    """ 
    296417    Generate random values for all of the parameters. 
    … …  
    301422    :func:`constrain_pars` needs to be called afterward.. 
    302423    """ 
    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())) 
     424    # Note: the sort guarantees order of calls to random number generator 
     425    random_pars = dict((p, _randomize_one(model_info, p, v)) 
     426                       for p, v in sorted(pars.items())) 
     427    if model_info.random is not None: 
     428        random_pars.update(model_info.random()) 
     429    _random_pd(model_info, random_pars) 
    307430    return random_pars 
     431 
    308432 
    309433def constrain_pars(model_info, pars): 
    … …  
    318442    Warning: this updates the *pars* dictionary in place. 
    319443    """ 
     444    # TODO: move the model specific code to the individual models 
    320445    name = model_info.id 
    321446    # if it is a product model, then just look at the form factor since 
    … …  
    338463    elif name == 'guinier': 
    339464        # Limit guinier to an Rg such that Iq > 1e-30 (single precision cutoff) 
     465        # I(q) = A e^-(Rg^2 q^2/3) > e^-(30 ln 10) 
     466        # => ln A - (Rg^2 q^2/3) > -30 ln 10 
     467        # => Rg^2 q^2/3 < 30 ln 10 + ln A 
     468        # => Rg < sqrt(90 ln 10 + 3 ln A)/q 
    340469        #q_max = 0.2  # mid q maximum 
    341470        q_max = 1.0  # high q maximum 
    … …  
    367496    parameters = model_info.parameters 
    368497    magnetic = False 
     498    magnetic_pars = [] 
    369499    for p in parameters.user_parameters(pars, is2d): 
    370500        if any(p.id.startswith(x) for x in ('M0:', 'mtheta:', 'mphi:')): 
    371501            continue 
    372         if p.id.startswith('up:') and not magnetic: 
     502        if p.id.startswith('up:'): 
     503            magnetic_pars.append("%s=%s"%(p.id, pars.get(p.id, p.default))) 
    373504            continue 
    374505        fields = dict( 
    … …  
    385516        lines.append(_format_par(p.name, **fields)) 
    386517        magnetic = magnetic or fields['M0'] != 0. 
     518    if magnetic and magnetic_pars: 
     519        lines.append(" ".join(magnetic_pars)) 
    387520    return "\n".join(lines) 
    388521 
    … …  
    402535    return line 
    403536 
    404 def suppress_pd(pars): 
     537def suppress_pd(pars, suppress=True): 
    405538    # type: (ParameterSet) -> ParameterSet 
    406539    """ 
    407     Suppress theta_pd for now until the normalization is resolved. 
    408  
    409     May also suppress complete polydispersity of the model to test 
    410     models more quickly. 
     540    If suppress is True complete eliminate polydispersity of the model to test 
     541    models more quickly.  If suppress is False, make sure at least one 
     542    parameter is polydisperse, setting the first polydispersity parameter to 
     543    15% if no polydispersity is given (with no explicit demo parameters given 
     544    in the model, there will be no default polydispersity). 
    411545    """ 
    412546    pars = pars.copy() 
    413     for p in pars: 
    414         if p.endswith("_pd_n"): pars[p] = 0 
     547    #print("pars=", pars) 
     548    if suppress: 
     549        for p in pars: 
     550            if p.endswith("_pd_n"): 
     551                pars[p] = 0 
     552    else: 
     553        any_pd = False 
     554        first_pd = None 
     555        for p in pars: 
     556            if p.endswith("_pd_n"): 
     557                pd = pars.get(p[:-2], 0.) 
     558                any_pd |= (pars[p] != 0 and pd != 0.) 
     559                if first_pd is None: 
     560                    first_pd = p 
     561        if not any_pd and first_pd is not None: 
     562            if pars[first_pd] == 0: 
     563                pars[first_pd] = 35 
     564            if first_pd[:-2] not in pars or pars[first_pd[:-2]] == 0: 
     565                pars[first_pd[:-2]] = 0.15 
    415566    return pars 
    416567 
    417 def suppress_magnetism(pars): 
     568def suppress_magnetism(pars, suppress=True): 
    418569    # type: (ParameterSet) -> ParameterSet 
    419570    """ 
    420     Suppress theta_pd for now until the normalization is resolved. 
    421  
    422     May also suppress complete polydispersity of the model to test 
    423     models more quickly. 
     571    If suppress is True complete eliminate magnetism of the model to test 
     572    models more quickly.  If suppress is False, make sure at least one sld 
     573    parameter is magnetic, setting the first parameter to have a strong 
     574    magnetic sld (8/A^2) at 60 degrees (with no explicit demo parameters given 
     575    in the model, there will be no default magnetism). 
    424576    """ 
    425577    pars = pars.copy() 
    426     for p in pars: 
    427         if p.startswith("M0:"): pars[p] = 0 
     578    if suppress: 
     579        for p in pars: 
     580            if p.startswith("M0:"): 
     581                pars[p] = 0 
     582    else: 
     583        any_mag = False 
     584        first_mag = None 
     585        for p in pars: 
     586            if p.startswith("M0:"): 
     587                any_mag |= (pars[p] != 0) 
     588                if first_mag is None: 
     589                    first_mag = p 
     590        if not any_mag and first_mag is not None: 
     591            pars[first_mag] = 8. 
    428592    return pars 
    429593 
    … …  
    561725    model = core.build_model(model_info, dtype=dtype, platform="ocl") 
    562726    calculator = DirectModel(data, model, cutoff=cutoff) 
    563     calculator.engine = "OCL%s"%DTYPE_MAP[dtype] 
     727    calculator.engine = "OCL%s"%DTYPE_MAP[str(model.dtype)] 
    564728    return calculator 
    565729 
    … …  
    571735    model = core.build_model(model_info, dtype=dtype, platform="dll") 
    572736    calculator = DirectModel(data, model, cutoff=cutoff) 
    573     calculator.engine = "OMP%s"%DTYPE_MAP[dtype] 
     737    calculator.engine = "OMP%s"%DTYPE_MAP[str(model.dtype)] 
    574738    return calculator 
    575739 
    … …  
    632796    if dtype == 'sasview': 
    633797        return eval_sasview(model_info, data) 
    634     elif dtype.endswith('!'): 
     798    elif dtype is None or not dtype.endswith('!'): 
     799        return eval_opencl(model_info, data, dtype=dtype, cutoff=cutoff) 
     800    else: 
    635801        return eval_ctypes(model_info, data, dtype=dtype[:-1], cutoff=cutoff) 
    636     else: 
    637         return eval_opencl(model_info, data, dtype=dtype, cutoff=cutoff) 
    638802 
    639803def _show_invalid(data, theory): 
    … …  
    662826    parameters. 
    663827    """ 
    664     result = run_models(opts, verbose=True) 
    665     if opts['plot']:  # Note: never called from explore 
    666         plot_models(opts, result, limits=limits) 
     828    limits = np.Inf, -np.Inf 
     829    for k in range(opts['sets']): 
     830        opts['pars'] = parse_pars(opts) 
     831        if opts['pars'] is None: 
     832            return 
     833        result = run_models(opts, verbose=True) 
     834        if opts['plot']: 
     835            limits = plot_models(opts, result, limits=limits, setnum=k) 
     836    if opts['plot']: 
     837        import matplotlib.pyplot as plt 
     838        plt.show() 
    667839 
    668840def run_models(opts, verbose=False): 
    669841    # type: (Dict[str, Any]) -> Dict[str, Any] 
    670842 
    671     n_base, n_comp = opts['count'] 
    672     pars, pars2 = opts['pars'] 
     843    base, comp = opts['engines'] 
     844    base_n, comp_n = opts['count'] 
     845    base_pars, comp_pars = opts['pars'] 
    673846    data = opts['data'] 
    674847 
    675     # silence the linter 
    676     base = opts['engines'][0] if n_base else None 
    677     comp = opts['engines'][1] if n_comp else None 
     848    comparison = comp is not None 
    678849 
    679850    base_time = comp_time = None 
    … …  
    681852 
    682853    # Base calculation 
    683     if n_base > 0: 
     854    try: 
     855        base_raw, base_time = time_calculation(base, base_pars, base_n) 
     856        base_value = np.ma.masked_invalid(base_raw) 
     857        if verbose: 
     858            print("%s t=%.2f ms, intensity=%.0f" 
     859                  % (base.engine, base_time, base_value.sum())) 
     860        _show_invalid(data, base_value) 
     861    except ImportError: 
     862        traceback.print_exc() 
     863 
     864    # Comparison calculation 
     865    if comparison: 
    684866        try: 
    685             base_raw, base_time = time_calculation(base, pars, n_base) 
    686             base_value = np.ma.masked_invalid(base_raw) 
    687             if verbose: 
    688                 print("%s t=%.2f ms, intensity=%.0f" 
    689                       % (base.engine, base_time, base_value.sum())) 
    690             _show_invalid(data, base_value) 
    691         except ImportError: 
    692             traceback.print_exc() 
    693             n_base = 0 
    694  
    695     # Comparison calculation 
    696     if n_comp > 0: 
    697         try: 
    698             comp_raw, comp_time = time_calculation(comp, pars2, n_comp) 
     867            comp_raw, comp_time = time_calculation(comp, comp_pars, comp_n) 
    699868            comp_value = np.ma.masked_invalid(comp_raw) 
    700869            if verbose: 
    … …  
    704873        except ImportError: 
    705874            traceback.print_exc() 
    706             n_comp = 0 
    707875 
    708876    # Compare, but only if computing both forms 
    709     if n_base > 0 and n_comp > 0: 
     877    if comparison: 
    710878        resid = (base_value - comp_value) 
    711879        relerr = resid/np.where(comp_value != 0., abs(comp_value), 1.0) 
    … …  
    743911 
    744912 
    745 def plot_models(opts, result, limits=None): 
     913def plot_models(opts, result, limits=(np.Inf, -np.Inf), setnum=0): 
    746914    # type: (Dict[str, Any], Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float] 
    747     base_value, comp_value= result['base_value'], result['comp_value'] 
     915    base_value, comp_value = result['base_value'], result['comp_value'] 
    748916    base_time, comp_time = result['base_time'], result['comp_time'] 
    749917    resid, relerr = result['resid'], result['relerr'] 
    750918 
    751919    have_base, have_comp = (base_value is not None), (comp_value is not None) 
    752     base = opts['engines'][0] if have_base else None 
    753     comp = opts['engines'][1] if have_comp else None 
     920    base, comp = opts['engines'] 
    754921    data = opts['data'] 
    755922    use_data = (opts['datafile'] is not None) and (have_base ^ have_comp) 
    … …  
    758925    view = opts['view'] 
    759926    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()) 
     927    vmin, vmax = limits 
     928    if have_base: 
     929        vmin = min(vmin, base_value.min()) 
     930        vmax = max(vmax, base_value.max()) 
     931    if have_comp: 
     932        vmin = min(vmin, comp_value.min()) 
     933        vmax = max(vmax, comp_value.max()) 
     934    limits = vmin, vmax 
     935 
     936    if have_base: 
    765937        if have_comp: 
    766             vmin = min(vmin, comp_value.min()) 
    767             vmax = max(vmax, comp_value.max()) 
    768         limits = vmin, vmax 
    769  
    770     if have_base: 
    771         if have_comp: plt.subplot(131) 
     938            plt.subplot(131) 
    772939        plot_theory(data, base_value, view=view, use_data=use_data, limits=limits) 
    773940        plt.title("%s t=%.2f ms"%(base.engine, base_time)) 
    774941        #cbar_title = "log I" 
    775942    if have_comp: 
    776         if have_base: plt.subplot(132) 
     943        if have_base: 
     944            plt.subplot(132) 
    777945        if not opts['is2d'] and have_base: 
    778946            plot_theory(data, base_value, view=view, use_data=use_data, limits=limits) 
    … …  
    789957            sorted = np.sort(err.flatten()) 
    790958            cutoff = sorted[int(sorted.size*0.95)] 
    791             err[err>cutoff] = cutoff 
     959            err[err > cutoff] = cutoff 
    792960        #err,errstr = base/comp,"ratio" 
    793961        plot_theory(data, None, resid=err, view=errview, use_data=use_data) 
    … …  
    813981        plt.title('Distribution of relative error between calculation engines') 
    814982 
    815     if not opts['explore']: 
    816         plt.show() 
    817  
    818983    return limits 
    819  
    820  
    821984 
    822985 
    823986# =========================================================================== 
    824987# 
    825 NAME_OPTIONS = set([ 
     988 
     989# Set of command line options. 
     990# Normal options such as -plot/-noplot are specified as 'name'. 
     991# For options such as -nq=500 which require a value use 'name='. 
     992# 
     993OPTIONS = [ 
     994    # Plotting 
    826995    'plot', 'noplot', 
    827     'half', 'fast', 'single', 'double', 
    828     'single!', 'double!', 'quad!', 'sasview', 
    829     'lowq', 'midq', 'highq', 'exq', 'zero', 
     996    'linear', 'log', 'q4', 
     997    'rel', 'abs', 
     998    'hist', 'nohist', 
     999    'title=', 
     1000 
     1001    # Data generation 
     1002    'data=', 'noise=', 'res=', 
     1003    'nq=', 'lowq', 'midq', 'highq', 'exq', 'zero', 
    8301004    '2d', '1d', 
    831     'preset', 'random', 
    832     'poly', 'mono', 
     1005 
     1006    # Parameter set 
     1007    'preset', 'random', 'random=', 'sets=', 
     1008    'demo', 'default',  # TODO: remove demo/default 
     1009    'nopars', 'pars', 
     1010 
     1011    # Calculation options 
     1012    'poly', 'mono', 'cutoff=', 
    8331013    'magnetic', 'nonmagnetic', 
    834     'nopars', 'pars', 
    835     'rel', 'abs', 
    836     'linear', 'log', 'q4', 
    837     'hist', 'nohist', 
    838     'edit', 'html', 'help', 
    839     'demo', 'default', 
    840     ]) 
    841 VALUE_OPTIONS = [ 
    842     # Note: random is both a name option and a value option 
    843     'cutoff', 'random', 'nq', 'res', 'accuracy', 'title', 'data', 
     1014    'accuracy=', 
     1015 
     1016    # Precision options 
     1017    'calc=', 
     1018    'half', 'fast', 'single', 'double', 'single!', 'double!', 'quad!', 
     1019    'sasview',  # TODO: remove sasview 3.x support 
     1020    'timing=', 
     1021 
     1022    # Output options 
     1023    'help', 'html', 'edit', 
    8441024    ] 
     1025 
     1026NAME_OPTIONS = set(k for k in OPTIONS if not k.endswith('=')) 
     1027VALUE_OPTIONS = [k[:-1] for k in OPTIONS if k.endswith('=')] 
     1028 
    8451029 
    8461030def columnize(items, indent="", width=79): 
    … …  
    9021086# not sure about brackets [] {} 
    9031087# maybe one of the following @ ? ^ ! , 
    904 MODEL_SPLIT = ',' 
     1088PAR_SPLIT = ',' 
    9051089def parse_opts(argv): 
    9061090    # type: (List[str]) -> Dict[str, Any] 
    … …  
    9141098              if not arg.startswith('-') and '=' in arg] 
    9151099    positional_args = [arg for arg in argv 
    916             if not arg.startswith('-') and '=' not in arg] 
     1100                       if not arg.startswith('-') and '=' not in arg] 
    9171101    models = "\n    ".join("%-15s"%v for v in MODELS) 
    9181102    if len(positional_args) == 0: 
    … …  
    9211105        print(columnize(MODELS, indent="  ")) 
    9221106        return None 
    923     if len(positional_args) > 3: 
    924         print("expected parameters: model N1 N2") 
    9251107 
    9261108    invalid = [o[1:] for o in flags 
    … …  
    9311113        return None 
    9321114 
    933     name = positional_args[0] 
    934     n1 = int(positional_args[1]) if len(positional_args) > 1 else 1 
    935     n2 = int(positional_args[2]) if len(positional_args) > 2 else 1 
     1115    name = positional_args[-1] 
    9361116 
    9371117    # pylint: disable=bad-whitespace 
    … …  
    9441124        'nq'        : 128, 
    9451125        'res'       : 0.0, 
     1126        'noise'     : 0.0, 
    9461127        'accuracy'  : 'Low', 
    947         'cutoff'    : 0.0, 
     1128        'cutoff'    : '0.0', 
    9481129        'seed'      : -1,  # default to preset 
    9491130        'mono'      : True, 
    … …  
    9591140        'title'     : None, 
    9601141        'datafile'  : None, 
     1142        'sets'      : 0, 
     1143        'engine'    : 'default', 
     1144        'evals'     : '1', 
    9611145    } 
    962     engines = [] 
    9631146    for arg in flags: 
    9641147        if arg == '-noplot':    opts['plot'] = False 
    … …  
    9761159        elif arg.startswith('-nq='):       opts['nq'] = int(arg[4:]) 
    9771160        elif arg.startswith('-res='):      opts['res'] = float(arg[5:]) 
     1161        elif arg.startswith('-noise='):    opts['noise'] = float(arg[7:]) 
     1162        elif arg.startswith('-sets='):     opts['sets'] = int(arg[6:]) 
    9781163        elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:] 
    979         elif arg.startswith('-cutoff='):   opts['cutoff'] = float(arg[8:]) 
     1164        elif arg.startswith('-cutoff='):   opts['cutoff'] = arg[8:] 
    9801165        elif arg.startswith('-random='):   opts['seed'] = int(arg[8:]) 
    9811166        elif arg.startswith('-title='):    opts['title'] = arg[7:] 
    9821167        elif arg.startswith('-data='):     opts['datafile'] = arg[6:] 
     1168        elif arg.startswith('-calc='):     opts['engine'] = arg[6:] 
     1169        elif arg.startswith('-neval='):    opts['evals'] = arg[7:] 
    9831170        elif arg == '-random':  opts['seed'] = np.random.randint(1000000) 
    9841171        elif arg == '-preset':  opts['seed'] = -1 
    … …  
    9931180        elif arg == '-rel':     opts['rel_err'] = True 
    9941181        elif arg == '-abs':     opts['rel_err'] = False 
    995         elif arg == '-half':    engines.append(arg[1:]) 
    996         elif arg == '-fast':    engines.append(arg[1:]) 
    997         elif arg == '-single':  engines.append(arg[1:]) 
    998         elif arg == '-double':  engines.append(arg[1:]) 
    999         elif arg == '-single!': engines.append(arg[1:]) 
    1000         elif arg == '-double!': engines.append(arg[1:]) 
    1001         elif arg == '-quad!':   engines.append(arg[1:]) 
    1002         elif arg == '-sasview': engines.append(arg[1:]) 
     1182        elif arg == '-half':    opts['engine'] = 'half' 
     1183        elif arg == '-fast':    opts['engine'] = 'fast' 
     1184        elif arg == '-single':  opts['engine'] = 'single' 
     1185        elif arg == '-double':  opts['engine'] = 'double' 
     1186        elif arg == '-single!': opts['engine'] = 'single!' 
     1187        elif arg == '-double!': opts['engine'] = 'double!' 
     1188        elif arg == '-quad!':   opts['engine'] = 'quad!' 
     1189        elif arg == '-sasview': opts['engine'] = 'sasview' 
    10031190        elif arg == '-edit':    opts['explore'] = True 
    10041191        elif arg == '-demo':    opts['use_demo'] = True 
    1005         elif arg == '-default':    opts['use_demo'] = False 
     1192        elif arg == '-default': opts['use_demo'] = False 
    10061193        elif arg == '-html':    opts['html'] = True 
    10071194        elif arg == '-help':    opts['html'] = True 
    10081195    # pylint: enable=bad-whitespace 
    10091196 
    1010     if MODEL_SPLIT in name: 
    1011         name, name2 = name.split(MODEL_SPLIT, 2) 
     1197    # Magnetism forces 2D for now 
     1198    if opts['magnetic']: 
     1199        opts['is2d'] = True 
     1200 
     1201    # Force random if sets is used 
     1202    if opts['sets'] >= 1 and opts['seed'] < 0: 
     1203        opts['seed'] = np.random.randint(1000000) 
     1204    if opts['sets'] == 0: 
     1205        opts['sets'] = 1 
     1206 
     1207    # Create the computational engines 
     1208    if opts['datafile'] is not None: 
     1209        data = load_data(os.path.expanduser(opts['datafile'])) 
    10121210    else: 
    1013         name2 = name 
     1211        data, _ = make_data(opts) 
     1212 
     1213    comparison = any(PAR_SPLIT in v for v in values) 
     1214    if PAR_SPLIT in name: 
     1215        names = name.split(PAR_SPLIT, 2) 
     1216        comparison = True 
     1217    else: 
     1218        names = [name]*2 
    10141219    try: 
    1015         model_info = core.load_model_info(name) 
    1016         model_info2 = core.load_model_info(name2) if name2 != name else model_info 
     1220        model_info = [core.load_model_info(k) for k in names] 
    10171221    except ImportError as exc: 
    10181222        print(str(exc)) 
    10191223        print("Could not find model; use one of:\n    " + models) 
    10201224        return None 
     1225 
     1226    if PAR_SPLIT in opts['engine']: 
     1227        engine_types = opts['engine'].split(PAR_SPLIT, 2) 
     1228        comparison = True 
     1229    else: 
     1230        engine_types = [opts['engine']]*2 
     1231 
     1232    if PAR_SPLIT in opts['evals']: 
     1233        evals = [int(k) for k in opts['evals'].split(PAR_SPLIT, 2)] 
     1234        comparison = True 
     1235    else: 
     1236        evals = [int(opts['evals'])]*2 
     1237 
     1238    if PAR_SPLIT in opts['cutoff']: 
     1239        cutoff = [float(k) for k in opts['cutoff'].split(PAR_SPLIT, 2)] 
     1240        comparison = True 
     1241    else: 
     1242        cutoff = [float(opts['cutoff'])]*2 
     1243 
     1244    base = make_engine(model_info[0], data, engine_types[0], cutoff[0]) 
     1245    if comparison: 
     1246        comp = make_engine(model_info[1], data, engine_types[1], cutoff[1]) 
     1247    else: 
     1248        comp = None 
     1249 
     1250    # pylint: disable=bad-whitespace 
     1251    # Remember it all 
     1252    opts.update({ 
     1253        'data'      : data, 
     1254        'name'      : names, 
     1255        'def'       : model_info, 
     1256        'count'     : evals, 
     1257        'engines'   : [base, comp], 
     1258        'values'    : values, 
     1259    }) 
     1260    # pylint: enable=bad-whitespace 
     1261 
     1262    return opts 
     1263 
     1264def parse_pars(opts): 
     1265    model_info, model_info2 = opts['def'] 
    10211266 
    10221267    # Get demo parameters from model definition, or use default parameters 
    … …  
    10281273    #pars.update(set_pars)  # set value before random to control range 
    10291274    if opts['seed'] > -1: 
    1030         pars = randomize_pars(model_info, pars, seed=opts['seed']) 
     1275        pars = randomize_pars(model_info, pars) 
    10311276        if model_info != model_info2: 
    1032             pars2 = randomize_pars(model_info2, pars2, seed=opts['seed']) 
     1277            pars2 = randomize_pars(model_info2, pars2) 
    10331278            # Share values for parameters with the same name 
    10341279            for k, v in pars.items(): 
    … …  
    10391284        constrain_pars(model_info, pars) 
    10401285        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) 
     1286    pars = suppress_pd(pars, opts['mono']) 
     1287    pars2 = suppress_pd(pars2, opts['mono']) 
     1288    pars = suppress_magnetism(pars, not opts['magnetic']) 
     1289    pars2 = suppress_magnetism(pars2, not opts['magnetic']) 
    10481290 
    10491291    # Fill in parameters given on the command line 
    10501292    presets = {} 
    10511293    presets2 = {} 
    1052     for arg in values: 
     1294    for arg in opts['values']: 
    10531295        k, v = arg.split('=', 1) 
    10541296        if k not in pars and k not in pars2: 
    … …  
    10571299            print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s)))) 
    10581300            return None 
    1059         v1, v2 = v.split(MODEL_SPLIT, 2) if MODEL_SPLIT in v else (v,v) 
     1301        v1, v2 = v.split(PAR_SPLIT, 2) if PAR_SPLIT in v else (v,v) 
    10601302        if v1 and k in pars: 
    10611303            presets[k] = float(v1) if isnumber(v1) else v1 
    … …  
    10751317    context['np'] = np 
    10761318    context.update(pars) 
    1077     context.update((k,v) for k,v in presets.items() if isinstance(v, float)) 
     1319    context.update((k, v) for k, v in presets.items() if isinstance(v, float)) 
    10781320    for k, v in presets.items(): 
    10791321        if not isinstance(v, float) and not k.endswith('_type'): 
    10801322            presets[k] = eval(v, context) 
    10811323    context.update(presets) 
    1082     context.update((k,v) for k,v in presets2.items() if isinstance(v, float)) 
     1324    context.update((k, v) for k, v in presets2.items() if isinstance(v, float)) 
    10831325    for k, v in presets2.items(): 
    10841326        if not isinstance(v, float) and not k.endswith('_type'): 
    … …  
    10901332    #import pprint; pprint.pprint(model_info) 
    10911333 
    1092     same_model = name == name2 and pars == pars 
    1093     if len(engines) == 0: 
    1094         if same_model: 
    1095             engines.extend(['single', 'double']) 
    1096         else: 
    1097             engines.extend(['single', 'single']) 
    1098     elif len(engines) == 1: 
    1099         if not same_model: 
    1100             engines.append(engines[0]) 
    1101         elif engines[0] == 'double': 
    1102             engines.append('single') 
    1103         else: 
    1104             engines.append('double') 
    1105     elif len(engines) > 2: 
    1106         del engines[2:] 
    1107  
    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  
    11141334    if opts['show_pars']: 
    1115         if not same_model: 
     1335        if model_info.name != model_info2.name or pars != pars2: 
    11161336            print("==== %s ====="%model_info.name) 
    11171337            print(str(parlist(model_info, pars, opts['is2d']))) 
    … …  
    11211341            print(str(parlist(model_info, pars, opts['is2d']))) 
    11221342 
    1123     # Create the computational engines 
    1124     if opts['datafile'] is not None: 
    1125         data = load_data(os.path.expanduser(opts['datafile'])) 
    1126     else: 
    1127         data, _ = make_data(opts) 
    1128     if n1: 
    1129         base = make_engine(model_info, data, engines[0], opts['cutoff']) 
    1130     else: 
    1131         base = None 
    1132     if n2: 
    1133         comp = make_engine(model_info2, data, engines[1], opts['cutoff']) 
    1134     else: 
    1135         comp = None 
    1136  
    1137     # pylint: disable=bad-whitespace 
    1138     # Remember it all 
    1139     opts.update({ 
    1140         'data'      : data, 
    1141         'name'      : [name, name2], 
    1142         'def'       : [model_info, model_info2], 
    1143         'count'     : [n1, n2], 
    1144         'presets'   : [presets, presets2], 
    1145         'pars'      : [pars, pars2], 
    1146         'engines'   : [base, comp], 
    1147     }) 
    1148     # pylint: enable=bad-whitespace 
    1149  
    1150     return opts 
     1343    return pars, pars2 
    11511344 
    11521345def show_docs(opts): 
    … …  
    11801373    model = Explore(opts) 
    11811374    problem = FitProblem(model) 
    1182     frame = AppFrame(parent=None, title="explore", size=(1000,700)) 
    1183     if not is_mac: frame.Show() 
     1375    frame = AppFrame(parent=None, title="explore", size=(1000, 700)) 
     1376    if not is_mac: 
     1377        frame.Show() 
    11841378    frame.panel.set_model(model=problem) 
    11851379    frame.panel.Layout() 
    … …  
    11911385    if is_mac: frame.Show() 
    11921386    # If running withing an app, start the main loop 
    1193     if app: app.MainLoop() 
     1387    if app: 
     1388        app.MainLoop() 
    11941389 
    11951390class Explore(object): 
    … …  
    12061401        config_matplotlib() 
    12071402        self.opts = opts 
     1403        opts['pars'] = list(opts['pars']) 
    12081404        p1, p2 = opts['pars'] 
    12091405        m1, m2 = opts['def'] 
    … …  
    12261422        self.starting_values = dict((k, v.value) for k, v in pars.items()) 
    12271423        self.pd_types = pd_types 
    1228         self.limits = None 
     1424        self.limits = np.Inf, -np.Inf 
    12291425 
    12301426    def revert_values(self): 
    … …  
    12831479    opts = parse_opts(argv) 
    12841480    if opts is not None: 
     1481        if opts['seed'] > -1: 
     1482            print("Randomize using -random=%i"%opts['seed']) 
     1483            np.random.seed(opts['seed']) 
    12851484        if opts['html']: 
    12861485            show_docs(opts) 
    12871486        elif opts['explore']: 
     1487            opts['pars'] = parse_pars(opts) 
     1488            if opts['pars'] is None: 
     1489                return 
    12881490            explore(opts) 
    12891491        else: 

  • sasmodels/core.py

    r60335cc r60335cc  
    308308 
    309309    # Convert dtype string to numpy dtype. 
    310     if dtype is None: 
     310    if dtype is None or dtype == "default": 
    311311        numpy_dtype = (generate.F32 if platform == "ocl" and model_info.single 
    312312                       else generate.F64) 

  • sasmodels/direct_model.py

    ra769b54 rd1ff3a5  
    293293        self.Iq = y 
    294294        if self.data_type in ('Iq', 'Iq-oriented'): 
     295            if self._data.y is None: 
     296                self._data.y = np.empty(len(self._data.x), 'd') 
     297            if self._data.dy is None: 
     298                self._data.dy = np.empty(len(self._data.x), 'd') 
    295299            self._data.dy[self.index] = dy 
    296300            self._data.y[self.index] = y 
    297301        elif self.data_type == 'Iqxy': 
     302            if self._data.data is None: 
     303                self._data.data = np.empty_like(self._data.qx_data, 'd') 
     304            if self._data.err_data is None: 
     305                self._data.err_data = np.empty_like(self._data.qx_data, 'd') 
    298306            self._data.data[self.index] = y 
     307            self._data.err_data[self.index] = dy 
    299308        elif self.data_type == 'sesans': 
     309            if self._data.y is None: 
     310                self._data.y = np.empty(len(self._data.x), 'd') 
    300311            self._data.y[self.index] = y 
    301312        else: 
    … …  
    315326        # TODO: extend plotting of calculate Iq to other measurement types 
    316327        # TODO: refactor so we don't store the result in the model 
    317         self.Iq_calc = None 
     328        self.Iq_calc = Iq_calc 
    318329        if self.data_type == 'sesans': 
    319330            Iq_mono = (call_kernel(self._kernel_mono, pars, mono=True) 

  • sasmodels/kernelpy.py

    rdbfd471 r883ecf4  
    99from __future__ import division, print_function 
    1010 
     11import logging 
     12 
    1113import numpy as np  # type: ignore 
    1214from numpy import pi, sin, cos  #type: ignore 
    … …  
    1820try: 
    1921    from typing import Union, Callable 
    20 except: 
     22except ImportError: 
    2123    pass 
    2224else: 
    … …  
    3133        _create_default_functions(model_info) 
    3234        self.info = model_info 
     35        self.dtype = np.dtype('d') 
    3336 
    3437    def make_kernel(self, q_vectors): 
     38        logging.info("creating python kernel " + self.info.name) 
    3539        q_input = PyInput(q_vectors, dtype=F64) 
    3640        kernel = self.info.Iqxy if q_input.is_2d else self.info.Iq 
    … …  
    236240            # INVALID expression like the C models, but that is too expensive. 
    237241            Iq = np.asarray(form(), 'd') 
    238             if np.isnan(Iq).any(): continue 
     242            if np.isnan(Iq).any(): 
     243                continue 
    239244 
    240245            # update value and norm 
    … …  
    300305        default_Iqxy.vectorized = True 
    301306        model_info.Iqxy = default_Iqxy 
    302  

  • sasmodels/list_pars.py

    r3330bb4 r72be531  
    1212 
    1313import sys 
     14import argparse 
    1415 
    1516from .core import load_model_info, list_models 
    1617from .compare import columnize 
    1718 
    18 def find_pars(): 
     19def find_pars(type=None): 
    1920    """ 
    2021    Find all parameters in all models. 
    … …  
    2627        model_info = load_model_info(name) 
    2728        for p in model_info.parameters.kernel_parameters: 
    28             partable.setdefault(p.name, []) 
    29             partable[p.name].append(name) 
     29            if type is None or p.type == type: 
     30                partable.setdefault(p.name, []) 
     31                partable[p.name].append(name) 
    3032    return partable 
    3133 
    32 def list_pars(names_only=True): 
     34def list_pars(names_only=True, type=None): 
    3335    """ 
    3436    Print all parameters in all models. 
    … …  
    3739    occurs in. 
    3840    """ 
    39     partable = find_pars() 
     41    partable = find_pars(type) 
    4042    if names_only: 
    4143        print(columnize(list(sorted(partable.keys())))) 
    … …  
    4850    Program to list the parameters used across all models. 
    4951    """ 
    50     if len(sys.argv) == 2 and sys.argv[1] == '-v': 
    51         verbose = True 
    52     elif len(sys.argv) == 1: 
    53         verbose = False 
    54     else: 
    55         print(__doc__) 
    56         sys.exit(1) 
    57     list_pars(names_only=not verbose) 
     52    parser = argparse.ArgumentParser( 
     53        description="Find all parameters in all models", 
     54        ) 
     55    parser.add_argument( 
     56        '-v', '--verbose', 
     57        action='store_true', 
     58        help="list models which use this argument") 
     59    parser.add_argument( 
     60        'type', default="any", nargs='?', 
     61        metavar="volume|orientation|sld|none|any", 
     62        choices=['volume', 'orientation', 'sld', None, 'any'], 
     63        type=lambda v: None if v == 'any' else '' if v == 'none' else v, 
     64        help="only list arguments of the given type") 
     65    args = parser.parse_args() 
     66 
     67    list_pars(names_only=not args.verbose, type=args.type) 
    5868 
    5969if __name__ == "__main__": 

  • sasmodels/modelinfo.py

    r65314f7 r65314f7  
    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): 
    … …  
    757766    info.opencl = getattr(kernel_module, 'opencl', not callable(info.Iq)) 
    758767    info.single = getattr(kernel_module, 'single', not callable(info.Iq)) 
     768    info.random = getattr(kernel_module, 'random', None) 
    759769 
    760770    # multiplicity info 

  • sasmodels/models/adsorbed_layer.py

    rb0c4271 r8f04da4  
    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 as done by default in compare.py 
     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/barbell.py

    r9802ab3 r31df0c9  
    115115source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "barbell.c"] 
    116116 
     117def random(): 
     118    import numpy as np 
     119    # TODO: increase volume range once problem with bell radius is fixed 
     120    # The issue is that bell radii of more than about 200 fail at high q 
     121    V = 10**np.random.uniform(7, 9) 
     122    bar_volume = 10**np.random.uniform(-4, -1)*V 
     123    bell_volume = V - bar_volume 
     124    bell_radius = (bell_volume/6)**0.3333  # approximate 
     125    min_bar = bar_volume/np.pi/bell_radius**2 
     126    bar_length = 10**np.random.uniform(0, 3)*min_bar 
     127    bar_radius = np.sqrt(bar_volume/bar_length/np.pi) 
     128    if bar_radius > bell_radius: 
     129        bell_radius, bar_radius = bar_radius, bell_radius 
     130    pars = dict( 
     131        #background=0, 
     132        radius_bell=bell_radius, 
     133        radius=bar_radius, 
     134        length=bar_length, 
     135    ) 
     136    return pars 
     137 
    117138# parameters for demo 
    118139demo = dict(scale=1, background=0, 

  • sasmodels/models/bcc_paracrystal.py

    r69e1afc r8f04da4  
    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 
    133 # parameters for demo 
    134 demo = dict( 
    135     scale=1, background=0, 
    136     dnn=220, d_factor=0.06, sld=4, sld_solvent=1, 
    137     radius=40, 
    138     theta=60, phi=60, psi=60, 
    139     radius_pd=.2, radius_pd_n=2, 
    140     theta_pd=15, theta_pd_n=0, 
    141     phi_pd=15, phi_pd_n=0, 
    142     psi_pd=15, psi_pd_n=0, 
     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    radius = 10**np.random.uniform(1.3, 4) 
     143    d_factor = 10**np.random.uniform(-2, -0.7)  # sigma_d in 0.01-0.7 
     144    dnn_fraction = np.random.beta(a=10, b=1) 
     145    dnn = radius*4/np.sqrt(3)/dnn_fraction 
     146    pars = dict( 
     147        #sld=1, sld_solvent=0, scale=1, background=1e-32, 
     148        dnn=dnn, 
     149        d_factor=d_factor, 
     150        radius=radius, 
    143151    ) 
     152    return pars 
     153 
    144154# april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
    145155# add 2d test later 
    146 q =4.*pi/220. 
     156q = 4.*pi/220. 
    147157tests = [ 
    148     [{ }, 
    149      [0.001, q, 0.215268], [1.46601394721, 2.85851284174, 0.00866710287078]], 
    150     [{'theta':20.0,'phi':30,'psi':40.0},(-0.017,0.035),2082.20264399 ], 
    151     [{'theta':20.0,'phi':30,'psi':40.0},(-0.081,0.011),0.436323144781 ] 
     158    [{}, [0.001, q, 0.215268], [1.46601394721, 2.85851284174, 0.00866710287078]], 
     159    [{'theta': 20.0, 'phi': 30, 'psi': 40.0}, (-0.017, 0.035), 2082.20264399], 
     160    [{'theta': 20.0, 'phi': 30, 'psi': 40.0}, (-0.081, 0.011), 0.436323144781], 
    152161    ] 

  • sasmodels/models/be_polyelectrolyte.py

    r3330bb4 r8f04da4  
    1717    r_{0}^2 = \frac{1}{\alpha \sqrt{C_a} \left( b/\sqrt{48\pi L_b}\right)} 
    1818 
    19 where  
     19where 
    2020 
    2121$K$ is the contrast factor for the polymer which is defined differently than in 
    … …  
    2828 
    2929    a = b_p - (v_p/v_s) b_s 
    30      
     30 
    3131where $b_p$ and $b_s$ are sum of the scattering lengths of the atoms 
    3232constituting the monomer of the polymer and the sum of the scattering lengths 
    … …  
    3535 
    3636$L_b$ is the Bjerrum length(|Ang|) - **Note:** This parameter needs to be 
    37 kept constant for a given solvent and temperature!  
     37kept constant for a given solvent and temperature! 
    3838 
    3939$h$ is the virial parameter (|Ang^3|/mol) - **Note:** See [#Borue]_ for the 
    … …  
    6767* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    6868* **Last Modified by:** Paul Kienzle **Date:** July 24, 2016 
    69 * **Last Reviewed by:** Paul Butler and Richard Heenan **Date:**  
     69* **Last Reviewed by:** Paul Butler and Richard Heenan **Date:** 
    7070  October 07, 2016 
    7171""" 
    … …  
    139139Iq.vectorized = True  # Iq accepts an array of q values 
    140140 
     141def random(): 
     142    import numpy as np 
     143    # TODO: review random be_polyelectrolyte model generation 
     144    pars = dict( 
     145        scale=10000, #background=0, 
     146        #polymer_concentration=0.7, 
     147        polymer_concentration=np.random.beta(5, 3), # around 70% 
     148        #salt_concentration=0.0, 
     149        # keep salt concentration extremely low 
     150        # and use explicit molar to match polymer concentration 
     151        salt_concentration=np.random.beta(1, 100)*6.022136e-4, 
     152        #contrast_factor=10.0, 
     153        contrast_fact=np.random.uniform(1, 100), 
     154        #bjerrum_length=7.1, 
     155        bjerrum_length=np.random.uniform(1, 10), 
     156        #virial_param=12.0, 
     157        virial_param=np.random.uniform(-1000, 30), 
     158        #monomer_length=10.0, 
     159        monomer_length=10.0**(4*np.random.beta(1.5, 3)), 
     160        #ionization_degree=0.05, 
     161        ionization_degree=np.random.beta(1.5, 4), 
     162    ) 
     163    return pars 
    141164 
    142165demo = dict(scale=1, background=0.1, 

  • sasmodels/models/binary_hard_sphere.py

    r925ad6e r8f04da4  
    110110source = ["lib/sas_3j1x_x.c", "binary_hard_sphere.c"] 
    111111 
     112def random(): 
     113    import numpy as np 
     114    # TODO: binary_hard_sphere fails at low qr 
     115    radius_lg = 10**np.random.uniform(2, 4.7) 
     116    radius_sm = 10**np.random.uniform(2, 4.7) 
     117    volfraction_lg = 10**np.random.uniform(-3, -0.3) 
     118    volfraction_sm = 10**np.random.uniform(-3, -0.3) 
     119    # TODO: Get slightly different results if large and small are swapped 
     120    # modify the model so it doesn't care which is which 
     121    if radius_lg < radius_sm: 
     122        radius_lg, radius_sm = radius_sm, radius_lg 
     123        volfraction_lg, volfraction_sm = volfraction_sm, volfraction_lg 
     124    pars = dict( 
     125        radius_lg=radius_lg, 
     126        radius_sm=radius_sm, 
     127        volfraction_lg=volfraction_lg, 
     128        volfraction_sm=volfraction_sm, 
     129    ) 
     130    return pars 
     131 
    112132# parameters for demo and documentation 
    113133demo = dict(sld_lg=3.5, sld_sm=0.5, sld_solvent=6.36, 

  • 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, 

  • sasmodels/models/capped_cylinder.py

    r9802ab3 r31df0c9  
    8080 
    8181.. [#] H Kaya, *J. Appl. Cryst.*, 37 (2004) 223-230 
    82 .. [#] H Kaya and N-R deSouza, *J. Appl. Cryst.*, 37 (2004) 508-509 (addenda  
     82.. [#] H Kaya and N-R deSouza, *J. Appl. Cryst.*, 37 (2004) 508-509 (addenda 
    8383   and errata) 
    8484 
    … …  
    136136source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "capped_cylinder.c"] 
    137137 
     138def random(): 
     139    import numpy as np 
     140    # TODO: increase volume range once problem with bell radius is fixed 
     141    # The issue is that bell radii of more than about 200 fail at high q 
     142    V = 10**np.random.uniform(7, 9) 
     143    bar_volume = 10**np.random.uniform(-4, -1)*V 
     144    bell_volume = V - bar_volume 
     145    bell_radius = (bell_volume/6)**0.3333  # approximate 
     146    min_bar = bar_volume/np.pi/bell_radius**2 
     147    bar_length = 10**np.random.uniform(0, 3)*min_bar 
     148    bar_radius = np.sqrt(bar_volume/bar_length/np.pi) 
     149    if bar_radius > bell_radius: 
     150        bell_radius, bar_radius = bar_radius, bell_radius 
     151    pars = dict( 
     152        #background=0, 
     153        radius_cap=bell_radius, 
     154        radius=bar_radius, 
     155        length=bar_length, 
     156    ) 
     157    return pars 
     158 
     159 
    138160demo = dict(scale=1, background=0, 
    139161            sld=6, sld_solvent=1, 

  • sasmodels/models/core_multi_shell.py

    r5a0b3d7 r1511c37c  
    7070        sld_shell: the SLD of the shell# 
    7171        A_shell#: the coefficient in the exponential function 
    72          
    73          
     72 
     73 
    7474    scale: 1.0 if data is on absolute scale 
    7575    volfraction: volume fraction of spheres 
    … …  
    102102 
    103103source = ["lib/sas_3j1x_x.c", "core_multi_shell.c"] 
     104 
     105def random(): 
     106    import numpy as np 
     107    num_shells = np.minimum(np.random.poisson(3)+1, 10) 
     108    total_radius = 10**np.random.uniform(1.7, 4) 
     109    thickness = np.random.exponential(size=num_shells+1) 
     110    thickness *= total_radius/np.sum(thickness) 
     111    pars = dict( 
     112        #background=0, 
     113        n=num_shells, 
     114        radius=thickness[0], 
     115    ) 
     116    for k, v in enumerate(thickness[1:]): 
     117        pars['thickness%d'%(k+1)] = v 
     118    return pars 
    104119 
    105120def profile(sld_core, radius, sld_solvent, n, sld, thickness): 

  • sasmodels/models/core_shell_bicelle.py

    r9802ab3 ra151caa  
    1717.. figure:: img/core_shell_bicelle_parameters.png 
    1818 
    19    Cross section of cylindrical symmetry model used here. Users will have  
    20    to decide how to distribute "heads" and "tails" between the rim, face  
     19   Cross section of cylindrical symmetry model used here. Users will have 
     20   to decide how to distribute "heads" and "tails" between the rim, face 
    2121   and core regions in order to estimate appropriate starting parameters. 
    2222 
    … …  
    2727.. math:: 
    2828 
    29     \rho(r) =  
    30       \begin{cases}  
     29    \rho(r) = 
     30      \begin{cases} 
    3131      &\rho_c \text{ for } 0 \lt r \lt R; -L \lt z\lt L \\[1.5ex] 
    3232      &\rho_f \text{ for } 0 \lt r \lt R; -(L+2t) \lt z\lt -L; 
    … …  
    4747.. math:: 
    4848 
    49     \begin{align}     
    50     F(Q,\alpha) = &\bigg[  
     49    \begin{align} 
     50    F(Q,\alpha) = &\bigg[ 
    5151    (\rho_c - \rho_f) V_c \frac{2J_1(QRsin \alpha)}{QRsin\alpha}\frac{sin(QLcos\alpha/2)}{Q(L/2)cos\alpha} \\ 
    5252    &+(\rho_f - \rho_r) V_{c+f} \frac{2J_1(QRsin\alpha)}{QRsin\alpha}\frac{sin(Q(L/2+t_f)cos\alpha)}{Q(L/2+t_f)cos\alpha} \\ 
    5353    &+(\rho_r - \rho_s) V_t \frac{2J_1(Q(R+t_r)sin\alpha)}{Q(R+t_r)sin\alpha}\frac{sin(Q(L/2+t_f)cos\alpha)}{Q(L/2+t_f)cos\alpha} 
    5454    \bigg] 
    55     \end{align}  
     55    \end{align} 
    5656 
    5757where $V_t$ is the total volume of the bicelle, $V_c$ the volume of the core, 
    … …  
    6363cylinders is then given by integrating over all possible $\theta$ and $\phi$. 
    6464 
    65 For oriented bicelles the *theta*, and *phi* orientation parameters will appear when fitting 2D data,  
     65For oriented bicelles the *theta*, and *phi* orientation parameters will appear when fitting 2D data, 
    6666see the :ref:`cylinder` model for further information. 
    6767Our implementation of the scattering kernel and the 1D scattering intensity 
    … …  
    9696title = "Circular cylinder with a core-shell scattering length density profile.." 
    9797description = """ 
    98     P(q,alpha)= (scale/Vs)*f(q)^(2) + bkg,  where:  
     98    P(q,alpha)= (scale/Vs)*f(q)^(2) + bkg,  where: 
    9999    f(q)= Vt(sld_rim - sld_solvent)* sin[qLt.cos(alpha)/2] 
    100100    /[qLt.cos(alpha)/2]*J1(qRout.sin(alpha)) 
    … …  
    147147          "core_shell_bicelle.c"] 
    148148 
     149def random(): 
     150    import numpy as np 
     151    pars = dict( 
     152        radius=10**np.random.uniform(1.3, 3), 
     153        length=10**np.random.uniform(1.3, 4), 
     154        thick_rim=10**np.random.uniform(0, 1.7), 
     155        thick_face=10**np.random.uniform(0, 1.7), 
     156    ) 
     157    return pars 
     158 
    149159demo = dict(scale=1, background=0, 
    150160            radius=20.0, 

  • sasmodels/models/core_shell_bicelle_elliptical.py

    r9802ab3 r8f04da4  
    66core-shell scattering length density profile. Thus this is a variation 
    77of the core-shell bicelle model, but with an elliptical cylinder for the core. 
    8 Outer shells on the rims and flat ends may be of different thicknesses and  
     8Outer shells on the rims and flat ends may be of different thicknesses and 
    99scattering length densities. The form factor is normalized by the total particle volume. 
    1010 
    … …  
    1717.. figure:: img/core_shell_bicelle_parameters.png 
    1818 
    19    Cross section of model used here. Users will have  
    20    to decide how to distribute "heads" and "tails" between the rim, face  
     19   Cross section of model used here. Users will have 
     20   to decide how to distribute "heads" and "tails" between the rim, face 
    2121   and core regions in order to estimate appropriate starting parameters. 
    2222 
    … …  
    2727.. math:: 
    2828 
    29     \rho(r) =  
    30       \begin{cases}  
     29    \rho(r) = 
     30      \begin{cases} 
    3131      &\rho_c \text{ for } 0 \lt r \lt R; -L \lt z\lt L \\[1.5ex] 
    3232      &\rho_f \text{ for } 0 \lt r \lt R; -(L+2t) \lt z\lt -L; 
    … …  
    4848.. math:: 
    4949 
    50         \begin{align}     
    51     F(Q,\alpha,\psi) = &\bigg[  
     50        \begin{align} 
     51    F(Q,\alpha,\psi) = &\bigg[ 
    5252    (\rho_c - \rho_f) V_c \frac{2J_1(QR'sin \alpha)}{QR'sin\alpha}\frac{sin(QLcos\alpha/2)}{Q(L/2)cos\alpha} \\ 
    5353    &+(\rho_f - \rho_r) V_{c+f} \frac{2J_1(QR'sin\alpha)}{QR'sin\alpha}\frac{sin(Q(L/2+t_f)cos\alpha)}{Q(L/2+t_f)cos\alpha} \\ 
    5454    &+(\rho_r - \rho_s) V_t \frac{2J_1(Q(R'+t_r)sin\alpha)}{Q(R'+t_r)sin\alpha}\frac{sin(Q(L/2+t_f)cos\alpha)}{Q(L/2+t_f)cos\alpha} 
    5555    \bigg] 
    56     \end{align}  
     56    \end{align} 
    5757 
    5858where 
    … …  
    6161 
    6262    R'=\frac{R}{\sqrt{2}}\sqrt{(1+X_{core}^{2}) + (1-X_{core}^{2})cos(\psi)} 
    63      
    64      
    65 and $V_t = \pi.(R+t_r)(Xcore.R+t_r)^2.(L+2.t_f)$ is the total volume of the bicelle,  
    66 $V_c = \pi.Xcore.R^2.L$ the volume of the core, $V_{c+f} = \pi.Xcore.R^2.(L+2.t_f)$  
     63 
     64 
     65and $V_t = \pi.(R+t_r)(Xcore.R+t_r)^2.(L+2.t_f)$ is the total volume of the bicelle, 
     66$V_c = \pi.Xcore.R^2.L$ the volume of the core, $V_{c+f} = \pi.Xcore.R^2.(L+2.t_f)$ 
    6767the volume of the core plus the volume of the faces, $R$ is the radius 
    68 of the core, $Xcore$ is the axial ratio of the core, $L$ the length of the core,  
    69 $t_f$ the thickness of the face, $t_r$ the thickness of the rim and $J_1$ the usual  
    70 first order bessel function. The core has radii $R$ and $Xcore.R$ so is circular,  
    71 as for the core_shell_bicelle model, for $Xcore$ =1. Note that you may need to  
    72 limit the range of $Xcore$, especially if using the Monte-Carlo algorithm, as  
     68of the core, $Xcore$ is the axial ratio of the core, $L$ the length of the core, 
     69$t_f$ the thickness of the face, $t_r$ the thickness of the rim and $J_1$ the usual 
     70first order bessel function. The core has radii $R$ and $Xcore.R$ so is circular, 
     71as for the core_shell_bicelle model, for $Xcore$ =1. Note that you may need to 
     72limit the range of $Xcore$, especially if using the Monte-Carlo algorithm, as 
    7373setting radius to $R/Xcore$ and axial ratio to $1/Xcore$ gives an equivalent solution! 
    7474 
    … …  
    7676bicelles is then given by integrating over all possible $\alpha$ and $\psi$. 
    7777 
    78 For oriented bicelles the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,  
     78For oriented bicelles the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data, 
    7979see the :ref:`elliptical-cylinder` model for further information. 
    8080 
    … …  
    8282.. figure:: img/elliptical_cylinder_angle_definition.png 
    8383 
    84     Definition of the angles for the oriented core_shell_bicelle_elliptical particles.    
     84    Definition of the angles for the oriented core_shell_bicelle_elliptical particles. 
    8585 
    8686 
    … …  
    105105description = """ 
    106106    core_shell_bicelle_elliptical 
    107     Elliptical cylinder core, optional shell on the two flat faces, and shell of  
    108     uniform thickness on its rim (extending around the end faces).     
     107    Elliptical cylinder core, optional shell on the two flat faces, and shell of 
     108    uniform thickness on its rim (extending around the end faces). 
    109109    Please see full documentation for equations and further details. 
    110110    Involves a double numerical integral around the ellipsoid diameter 
    … …  
    136136          "core_shell_bicelle_elliptical.c"] 
    137137 
    138 demo = dict(scale=1, background=0, 
    139             radius=30.0, 
    140             x_core=3.0, 
    141             thick_rim=8.0, 
    142             thick_face=14.0, 
    143             length=50.0, 
    144             sld_core=4.0, 
    145             sld_face=7.0, 
    146             sld_rim=1.0, 
    147             sld_solvent=6.0, 
    148             theta=90, 
    149             phi=0, 
    150             psi=0) 
     138def random(): 
     139    import numpy as np 
     140    outer_major = 10**np.random.uniform(1, 4.7) 
     141    outer_minor = 10**np.random.uniform(1, 4.7) 
     142    # Use a distribution with a preference for thin shell or thin core, 
     143    # limited by the minimum radius. Avoid core,shell radii < 1 
     144    min_radius = min(outer_major, outer_minor) 
     145    thick_rim = np.random.beta(0.5, 0.5)*(min_radius-2) + 1 
     146    radius_major = outer_major - thick_rim 
     147    radius_minor = outer_minor - thick_rim 
     148    radius = radius_major 
     149    x_core = radius_minor/radius_major 
     150    outer_length = 10**np.random.uniform(1, 4.7) 
     151    # Caps should be a small percentage of the total length, but at least one 
     152    # angstrom long.  Since outer length >= 10, the following will suffice 
     153    thick_face = 10**np.random.uniform(-np.log10(outer_length), -1)*outer_length 
     154    length = outer_length - thick_face 
     155    pars = dict( 
     156        radius=radius, 
     157        x_core=x_core, 
     158        thick_rim=thick_rim, 
     159        thick_face=thick_face, 
     160        length=length 
     161    ) 
     162    return pars 
     163 
    151164 
    152165q = 0.1 

  • sasmodels/models/core_shell_cylinder.py

    r9b79f29 r9f6823b  
    142142    return whole, whole - core 
    143143 
     144def random(): 
     145    import numpy as np 
     146    outer_radius = 10**np.random.uniform(1, 4.7) 
     147    # Use a distribution with a preference for thin shell or thin core 
     148    # Avoid core,shell radii < 1 
     149    radius = np.random.beta(0.5, 0.5)*(outer_radius-2) + 1 
     150    thickness = outer_radius - radius 
     151    length = np.random.uniform(1, 4.7) 
     152    pars = dict( 
     153        radius=radius, 
     154        thickness=thickness, 
     155        length=length, 
     156    ) 
     157    return pars 
     158 
    144159demo = dict(scale=1, background=0, 
    145160            sld_core=6, sld_shell=8, sld_solvent=1, 

  • sasmodels/models/core_shell_ellipsoid.py

    r9802ab3 r9f6823b  
    2525ellipsoid. This may have some undesirable effects if the aspect ratio of the 
    2626ellipsoid is large (ie, if $X << 1$ or $X >> 1$ ), when the $S(q)$ 
    27 - which assumes spheres - will not in any case be valid.  Generating a  
    28 custom product model will enable separate effective volume fraction and effective  
     27- which assumes spheres - will not in any case be valid.  Generating a 
     28custom product model will enable separate effective volume fraction and effective 
    2929radius in the $S(q)$. 
    3030 
    … …  
    4444 
    4545.. math:: 
    46     \begin{align}     
     46    \begin{align} 
    4747    F(q,\alpha) = &f(q,radius\_equat\_core,radius\_equat\_core.x\_core,\alpha) \\ 
    4848    &+ f(q,radius\_equat\_core + thick\_shell,radius\_equat\_core.x\_core + thick\_shell.x\_polar\_shell,\alpha) 
    49     \end{align}  
     49    \end{align} 
    5050 
    5151where 
    52   
     52 
    5353.. math:: 
    5454 
    … …  
    7777   F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha} 
    7878 
    79 For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,  
     79For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data, 
    8080see the :ref:`elliptical-cylinder` model for further information. 
    8181 
    … …  
    151151    return ellipsoid_ER(polar_outer, equat_outer) 
    152152 
    153  
    154 demo = dict(scale=0.05, background=0.001, 
    155             radius_equat_core=20.0, 
    156             x_core=3.0, 
    157             thick_shell=30.0, 
    158             x_polar_shell=1.0, 
    159             sld_core=2.0, 
    160             sld_shell=1.0, 
    161             sld_solvent=6.3, 
    162             theta=0, 
    163             phi=0) 
     153def random(): 
     154    import numpy as np 
     155    V = 10**np.random.uniform(5, 12) 
     156    outer_polar = 10**np.random.uniform(1.3, 4) 
     157    outer_equatorial = np.sqrt(V/outer_polar) # ignore 4/3 pi 
     158    # Use a distribution with a preference for thin shell or thin core 
     159    # Avoid core,shell radii < 1 
     160    thickness_polar = np.random.beta(0.5, 0.5)*(outer_polar-2) + 1 
     161    thickness_equatorial = np.random.beta(0.5, 0.5)*(outer_equatorial-2) + 1 
     162    radius_polar = outer_polar - thickness_polar 
     163    radius_equatorial = outer_equatorial - thickness_equatorial 
     164    x_core = radius_polar/radius_equatorial 
     165    x_polar_shell = thickness_polar/thickness_equatorial 
     166    pars = dict( 
     167        #background=0, sld=0, sld_solvent=1, 
     168        radius_equat_core=radius_equatorial, 
     169        x_core=x_core, 
     170        thick_shell=thickness_equatorial, 
     171        x_polar_shell=x_polar_shell, 
     172    ) 
     173    return pars 
    164174 
    165175q = 0.1 

  • sasmodels/models/core_shell_parallelepiped.py

    r9b79f29 r8f04da4  
    44 
    55Calculates the form factor for a rectangular solid with a core-shell structure. 
    6 The thickness and the scattering length density of the shell or  
     6The thickness and the scattering length density of the shell or 
    77"rim" can be different on each (pair) of faces. However at this time 
    88the 1D calculation does **NOT** actually calculate a c face rim despite the presence of 
    … …  
    4242 
    4343**meaning that there are "gaps" at the corners of the solid.**  Again note that 
    44 $t_C = 0$ currently.  
     44$t_C = 0$ currently. 
    4545 
    4646The intensity calculated follows the :ref:`parallelepiped` model, with the 
    … …  
    8282based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
    8383and length $(C+2t_C)$ values, after appropriately 
    84 sorting the three dimensions to give an oblate or prolate particle, to give an  
     84sorting the three dimensions to give an oblate or prolate particle, to give an 
    8585effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
    8686 
    … …  
    126126title = "Rectangular solid with a core-shell structure." 
    127127description = """ 
    128      P(q)=  
     128     P(q)= 
    129129""" 
    130130category = "shape:parallelepiped" 
    … …  
    179179# VR defaults to 1.0 
    180180 
     181def random(): 
     182    import numpy as np 
     183    outer = 10**np.random.uniform(1, 4.7, size=3) 
     184    thick = np.random.beta(0.5, 0.5, size=3)*(outer-2) + 1 
     185    length = outer - thick 
     186    pars = dict( 
     187        length_a=length[0], 
     188        length_b=length[1], 
     189        length_c=length[2], 
     190        thick_rim_a=thick[0], 
     191        thick_rim_b=thick[1], 
     192        thick_rim_c=thick[2], 
     193    ) 
     194    return pars 
     195 
    181196# parameters for demo 
    182197demo = dict(scale=1, background=0.0, 

  • sasmodels/models/core_shell_sphere.py

    r925ad6e r9f6823b  
    5757title = "Form factor for a monodisperse spherical particle with particle with a core-shell structure." 
    5858description = """ 
    59     F^2(q) = 3/V_s [V_c (sld_core-sld_shell) (sin(q*radius)-q*radius*cos(q*radius))/(q*radius)^3  
     59    F^2(q) = 3/V_s [V_c (sld_core-sld_shell) (sin(q*radius)-q*radius*cos(q*radius))/(q*radius)^3 
    6060                   + V_s (sld_shell-sld_solvent) (sin(q*r_s)-q*r_s*cos(q*r_s))/(q*r_s)^3] 
    6161 
    … …  
    9999    return whole, whole - core 
    100100 
     101def random(): 
     102    import numpy as np 
     103    outer_radius = 10**np.random.uniform(1.3, 4.3) 
     104    # Use a distribution with a preference for thin shell or thin core 
     105    # Avoid core,shell radii < 1 
     106    radius = np.random.beta(0.5, 0.5)*(outer_radius-2) + 1 
     107    thickness = outer_radius - radius 
     108    pars = dict( 
     109        radius=radius, 
     110        thickness=thickness, 
     111    ) 
     112    return pars 
     113 
    101114tests = [ 
    102115    [{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0], 
    103      # TODO: VR test suppressed until we sort out new product model 
    104      # and determine what to do with volume ratio. 
    105      #[{'radius': 20.0, 'thickness': 10.0}, 'VR', 0.703703704], 
     116    # TODO: VR test suppressed until we sort out new product model 
     117    # and determine what to do with volume ratio. 
     118    #[{'radius': 20.0, 'thickness': 10.0}, 'VR', 0.703703704], 
    106119 
    107      # The SasView test result was 0.00169, with a background of 0.001 
    108      [{'radius': 60.0, 'thickness': 10.0, 'sld_core': 1.0, 'sld_shell':2.0, 
    109        'sld_solvent':3.0, 'background':0.0}, 
    110       0.4, 0.000698838], 
     120    # The SasView test result was 0.00169, with a background of 0.001 
     121    [{'radius': 60.0, 'thickness': 10.0, 'sld_core': 1.0, 'sld_shell': 2.0, 
     122      'sld_solvent': 3.0, 'background': 0.0}, 0.4, 0.000698838], 
    111123] 

  • sasmodels/models/cylinder.py

    r9802ab3 r31df0c9  
    6363.. figure:: img/cylinder_angle_definition.png 
    6464 
    65     Definition of the $\theta$ and $\phi$ orientation angles for a cylinder relative  
    66     to the beam line coordinates, plus an indication of their orientation distributions  
    67     which are described as rotations about each of the perpendicular axes $\delta_1$ and $\delta_2$  
     65    Definition of the $\theta$ and $\phi$ orientation angles for a cylinder relative 
     66    to the beam line coordinates, plus an indication of their orientation distributions 
     67    which are described as rotations about each of the perpendicular axes $\delta_1$ and $\delta_2$ 
    6868    in the frame of the cylinder itself, which when $\theta = \phi = 0$ are parallel to the $Y$ and $X$ axes. 
    6969 
    … …  
    7272    Examples for oriented cylinders. 
    7373 
    74 The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data.  
     74The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data. 
    7575On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will 
    76 appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, which when $\theta = \phi = 0$ are parallel  
     76appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, which when $\theta = \phi = 0$ are parallel 
    7777to the $Y$ and $X$ axes of the instrument respectively. Some experimentation may be required to understand the 2d patterns fully. 
    78 (Earlier implementations had numerical integration issues in some circumstances when orientation distributions passed through 90 degrees, such  
    79 situations, with very broad distributions, should still be approached with care.)  
     78(Earlier implementations had numerical integration issues in some circumstances when orientation distributions passed through 90 degrees, such 
     79situations, with very broad distributions, should still be approached with care.) 
    8080 
    8181Validation 
    … …  
    150150    return 0.5 * (ddd) ** (1. / 3.) 
    151151 
     152def random(): 
     153    import numpy as np 
     154    V = 10**np.random.uniform(5, 12) 
     155    length = 10**np.random.uniform(-2, 2)*V**0.333 
     156    radius = np.sqrt(V/length/np.pi) 
     157    pars = dict( 
     158        #scale=1, 
     159        #background=0, 
     160        length=length, 
     161        radius=radius, 
     162    ) 
     163    return pars 
     164 
     165 
    152166# parameters for demo 
    153167demo = dict(scale=1, background=0, 

  • sasmodels/models/dab.py

    r4962519 r404ebbd  
    6161    double numerator   = cube(cor_length); 
    6262    double denominator = square(1 + square(q*cor_length)); 
    63      
     63 
    6464    return numerator / denominator ; 
    6565    """ 
     66 
     67def random(): 
     68    import numpy as np 
     69    pars = dict( 
     70        scale=10**np.random.uniform(1, 4), 
     71        cor_length=10**np.random.uniform(0.3, 3), 
     72        #background = 0, 
     73    ) 
     74    pars['scale'] /= pars['cor_length']**3 
     75    return pars 
    6676 
    6777# ER defaults to 1.0 

  • sasmodels/models/ellipsoid.py

    r9b79f29 r92708d8  
    22# Note: model title and parameter table are inserted automatically 
    33r""" 
    4 The form factor is normalized by the particle volume  
     4The form factor is normalized by the particle volume 
    55 
    66Definition 
    … …  
    112112Plenum Press, New York, 1987. 
    113113 
     114A. Isihara. J. Chem. Phys. 18(1950) 1446-1449 
     115 
    114116Authorship and Verification 
    115117---------------------------- 
    … …  
    119121* **Last Modified by:** Paul Kienzle **Date:** March 22, 2017 
    120122""" 
     123from __future__ import division 
    121124 
    122125from numpy import inf, sin, cos, pi 
    … …  
    161164def ER(radius_polar, radius_equatorial): 
    162165    import numpy as np 
    163 # see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
     166    # see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
    164167    ee = np.empty_like(radius_polar) 
    165168    idx = radius_polar > radius_equatorial 
    … …  
    181184    return 0.5 * ddd ** (1.0 / 3.0) 
    182185 
     186def random(): 
     187    import numpy as np 
     188    V = 10**np.random.uniform(5, 12) 
     189    radius_polar = 10**np.random.uniform(1.3, 4) 
     190    radius_equatorial = np.sqrt(V/radius_polar) # ignore 4/3 pi 
     191    pars = dict( 
     192        #background=0, sld=0, sld_solvent=1, 
     193        radius_polar=radius_polar, 
     194        radius_equatorial=radius_equatorial, 
     195    ) 
     196    return pars 
    183197 
    184198demo = dict(scale=1, background=0, 

  • sasmodels/models/elliptical_cylinder.py

    r9802ab3 rd9ec8f9  
    3333 
    3434    a = qr'\sin(\alpha) 
    35      
     35 
    3636    b = q\frac{L}{2}\cos(\alpha) 
    37      
     37 
    3838    r'=\frac{r_{minor}}{\sqrt{2}}\sqrt{(1+\nu^{2}) + (1-\nu^{2})cos(\psi)} 
    3939 
    … …  
    5757define the axis of the cylinder using two angles $\theta$, $\phi$ and $\Psi$ 
    5858(see :ref:`cylinder orientation <cylinder-angle-definition>`). The angle 
    59 $\Psi$ is the rotational angle around its own long_c axis.  
     59$\Psi$ is the rotational angle around its own long_c axis. 
    6060 
    6161All angle parameters are valid and given only for 2D calculation; ie, an 
    … …  
    7272    detector plane, with $\Psi$ = 0. 
    7373 
    74 The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data.  
     74The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data. 
    7575On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will 
    76 appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, the $b$ and $a$ axes of the  
    77 cylinder cross section. (When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.)  
    78 The third orientation distribution, in $\psi$, is about the $c$ axis of the particle. Some experimentation may be required to  
    79 understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation  
    80 distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.)  
     76appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, the $b$ and $a$ axes of the 
     77cylinder cross section. (When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.) 
     78The third orientation distribution, in $\psi$, is about the $c$ axis of the particle. Some experimentation may be required to 
     79understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation 
     80distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.) 
    8181 
    8282NB: The 2nd virial coefficient of the cylinder is calculated based on the 
    … …  
    110110 
    111111* **Author:** 
    112 * **Last Modified by:**  
     112* **Last Modified by:** 
    113113* **Last Reviewed by:**  Richard Heenan - corrected equation in docs **Date:** December 21, 2016 
    114114 
    … …  
    156156                           + (length + radius) * (length + pi * radius)) 
    157157    return 0.5 * (ddd) ** (1. / 3.) 
     158 
     159def random(): 
     160    import numpy as np 
     161    # V = pi * radius_major * radius_minor * length; 
     162    V = 10**np.random.uniform(3, 9) 
     163    length = 10**np.random.uniform(1, 3) 
     164    axis_ratio = 10**np.random.uniform(0, 2) 
     165    radius_minor = np.sqrt(V/length/axis_ratio) 
     166    Vf = 10**np.random.uniform(-4, -2) 
     167    pars = dict( 
     168        #background=0, sld=0, sld_solvent=1, 
     169        scale=1e9*Vf/V, 
     170        length=length, 
     171        radius_minor=radius_minor, 
     172        axis_ratio=axis_ratio, 
     173    ) 
     174    return pars 
     175 
    158176q = 0.1 
    159177# april 6 2017, rkh added a 2d unit test, NOT READY YET pull #890 branch assume correct! 

  • sasmodels/models/fcc_paracrystal.py

    r69e1afc r8f04da4  
    7878.. figure:: img/parallelepiped_angle_definition.png 
    7979 
    80     Orientation of the crystal with respect to the scattering plane, when  
     80    Orientation of the crystal with respect to the scattering plane, when 
    8181    $\theta = \phi = 0$ the $c$ axis is along the beam direction (the $z$ axis). 
    8282 
    … …  
    119119source = ["lib/sas_3j1x_x.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc_paracrystal.c"] 
    120120 
    121 # parameters for demo 
    122 demo = dict(scale=1, background=0, 
    123             dnn=220, d_factor=0.06, sld=4, sld_solvent=1, 
    124             radius=40, 
    125             theta=60, phi=60, psi=60, 
    126             radius_pd=.2, radius_pd_n=0.2, 
    127             theta_pd=15, theta_pd_n=0, 
    128             phi_pd=15, phi_pd_n=0, 
    129             psi_pd=15, psi_pd_n=0, 
    130            ) 
     121def random(): 
     122    import numpy as np 
     123    # copied from bcc_paracrystal 
     124    radius = 10**np.random.uniform(1.3, 4) 
     125    d_factor = 10**np.random.uniform(-2, -0.7)  # sigma_d in 0.01-0.7 
     126    dnn_fraction = np.random.beta(a=10, b=1) 
     127    dnn = radius*4/np.sqrt(2)/dnn_fraction 
     128    pars = dict( 
     129        #sld=1, sld_solvent=0, scale=1, background=1e-32, 
     130        dnn=dnn, 
     131        d_factor=d_factor, 
     132        radius=radius, 
     133    ) 
     134    return pars 
     135 
    131136# april 10 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 
    132 q =4.*pi/220. 
     137q = 4.*pi/220. 
    133138tests = [ 
    134     [{ }, 
    135      [0.001, q, 0.215268], [0.275164706668, 5.7776842567, 0.00958167119232]], 
    136      [{}, (-0.047,-0.007), 238.103096286], 
    137      [{}, (0.053,0.063), 0.863609587796 ], 
     139    [{}, [0.001, q, 0.215268], [0.275164706668, 5.7776842567, 0.00958167119232]], 
     140    [{}, (-0.047, -0.007), 238.103096286], 
     141    [{}, (0.053, 0.063), 0.863609587796], 
    138142] 

  • sasmodels/models/flexible_cylinder.py

    r42356c8 r31df0c9  
    8686source = ["lib/polevl.c", "lib/sas_J1.c", "lib/wrc_cyl.c", "flexible_cylinder.c"] 
    8787 
    88 demo = dict(scale=1.0, background=0.0001, 
    89             length=1000.0, 
    90             kuhn_length=100.0, 
    91             radius=20.0, 
    92             sld=1.0, 
    93             sld_solvent=6.3) 
     88def random(): 
     89    import numpy as np 
     90    length = 10**np.random.uniform(2, 6) 
     91    radius = 10**np.random.uniform(1, 3) 
     92    kuhn_length = 10**np.random.uniform(-2, -0.7)*length  # at least 10 segments 
     93    pars = dict( 
     94        length=length, 
     95        radius=radius, 
     96        kuhn_length=kuhn_length, 
     97    ) 
     98    return pars 
    9499 
    95100tests = [ 

  • sasmodels/models/flexible_cylinder_elliptical.py

    r40a87fa r31df0c9  
    112112          "flexible_cylinder_elliptical.c"] 
    113113 
    114 demo = dict(scale=1.0, background=0.0001, 
    115             length=1000.0, 
    116             kuhn_length=100.0, 
    117             radius=20.0, 
    118             axis_ratio=1.5, 
    119             sld=1.0, 
    120             sld_solvent=6.3) 
     114def random(): 
     115    import numpy as np 
     116    length = 10**np.random.uniform(2, 6) 
     117    radius = 10**np.random.uniform(1, 3) 
     118    axis_ratio = 10**np.random.uniform(-1, 1) 
     119    kuhn_length = 10**np.random.uniform(-2, -0.7)*length  # at least 10 segments 
     120    pars = dict( 
     121        length=length, 
     122        radius=radius, 
     123        axis_ratio=axis_ratio, 
     124        kuhn_length=kuhn_length, 
     125    ) 
     126    return pars 
    121127 
    122128tests = [ 

  • sasmodels/models/fractal.py

    rdf89d77 r1511c37c  
    2727 
    2828where $\xi$ is the correlation length representing the cluster size and $D_f$ 
    29 is the fractal dimension, representing the self similarity of the structure.  
    30 Note that S(q) here goes negative if $D_f$ is too large, and the Gamma function  
     29is the fractal dimension, representing the self similarity of the structure. 
     30Note that S(q) here goes negative if $D_f$ is too large, and the Gamma function 
    3131diverges at $D_f=0$ and $D_f=1$. 
    3232 
    … …  
    5555 
    5656""" 
     57from __future__ import division 
    5758 
    5859from numpy import inf 
    … …  
    9899source = ["lib/sas_3j1x_x.c", "lib/sas_gamma.c", "lib/fractal_sq.c", "fractal.c"] 
    99100 
     101def random(): 
     102    import numpy as np 
     103    radius = 10**np.random.uniform(0.7, 4) 
     104    #radius = 5 
     105    cor_length = 10**np.random.uniform(0.7, 2)*radius 
     106    #cor_length = 20*radius 
     107    volfraction = 10**np.random.uniform(-3, -1) 
     108    #volfraction = 0.05 
     109    fractal_dim = 2*np.random.beta(3, 4) + 1 
     110    #fractal_dim = 2 
     111    pars = dict( 
     112        #background=0, sld_block=1, sld_solvent=0, 
     113        volfraction=volfraction, 
     114        radius=radius, 
     115        cor_length=cor_length, 
     116        fractal_dim=fractal_dim, 
     117    ) 
     118    return pars 
     119 
    100120demo = dict(volfraction=0.05, 
    101121            radius=5.0, 

  • sasmodels/models/fractal_core_shell.py

    r925ad6e r8f04da4  
    2424    \frac{\sin(qr_s)-qr_s\cos(qr_s)}{(qr_s)^3}\right]^2 
    2525 
    26     S(q) &= 1 + \frac{D_f\ \Gamma\!(D_f-1)}{[1+1/(q\xi)^2]^{(D_f-1)/2}}  
     26    S(q) &= 1 + \frac{D_f\ \Gamma\!(D_f-1)}{[1+1/(q\xi)^2]^{(D_f-1)/2}} 
    2727    \frac{\sin[(D_f-1)\tan^{-1}(q\xi)]}{(qr_s)^{D_f}} 
    2828 
    … …  
    3333of the whole particle respectively, $D_f$ is the fractal dimension, and |xi| the 
    3434correlation length. 
    35   
     35 
    3636Polydispersity of radius and thickness are also provided for. 
    3737 
    … …  
    9898          "lib/fractal_sq.c", "fractal_core_shell.c"] 
    9999 
     100def random(): 
     101    import numpy as np 
     102    outer_radius = 10**np.random.uniform(0.7, 4) 
     103    # Use a distribution with a preference for thin shell or thin core 
     104    # Avoid core,shell radii < 1 
     105    thickness = np.random.beta(0.5, 0.5)*(outer_radius-2) + 1 
     106    radius = outer_radius - thickness 
     107    cor_length = 10**np.random.uniform(0.7, 2)*outer_radius 
     108    volfraction = 10**np.random.uniform(-3, -1) 
     109    fractal_dim = 2*np.random.beta(3, 4) + 1 
     110    pars = dict( 
     111        #background=0, sld_block=1, sld_solvent=0, 
     112        volfraction=volfraction, 
     113        radius=radius, 
     114        cor_length=cor_length, 
     115        fractal_dim=fractal_dim, 
     116    ) 
     117    return pars 
     118 
    100119demo = dict(scale=0.05, 
    101120            background=0, 

  • sasmodels/models/fuzzy_sphere.py

    r925ad6e r31df0c9  
    105105# VR defaults to 1.0 
    106106 
     107def random(): 
     108    import numpy as np 
     109    radius = 10**np.random.uniform(1, 4.7) 
     110    fuzziness = 10**np.random.uniform(-2, -0.5)*radius  # 1% to 31% fuzziness 
     111    pars = dict( 
     112        radius=radius, 
     113        fuzziness=fuzziness, 
     114    ) 
     115    return pars 
     116 
    107117demo = dict(scale=1, background=0.001, 
    108118            sld=1, sld_solvent=3, 

  • sasmodels/models/gauss_lorentz_gel.py

    ra807206 r48462b0  
    33but typically a physical rather than chemical network. 
    44It is modeled as a sum of a low-q exponential decay (which happens to 
    5 give a functional form similar to Guinier scattering, so interpret with  
     5give a functional form similar to Guinier scattering, so interpret with 
    66care) plus a Lorentzian at higher-q values. See also the gel_fit model. 
    77 
    … …  
    8888 
    8989 
     90def random(): 
     91    import numpy as np 
     92    gauss_scale = 10**np.random.uniform(1, 3) 
     93    lorentz_scale = 10**np.random.uniform(1, 3) 
     94    cor_length_static = 10**np.random.uniform(0, 3) 
     95    cor_length_dynamic = 10**np.random.uniform(0, 3) 
     96    pars = dict( 
     97        #background=0, 
     98        scale=1, 
     99        gauss_scale=gauss_scale, 
     100        lorentz_scale=lorentz_scale, 
     101        cor_length_static=cor_length_static, 
     102        cor_length_dynamic=cor_length_dynamic, 
     103    ) 
     104    return pars 
     105 
     106 
    90107demo = dict(scale=1, background=0.1, 
    91108            gauss_scale=100.0, 

  • sasmodels/models/gaussian_peak.py

    ra807206 r48462b0  
    5151# VR defaults to 1.0 
    5252 
    53 demo = dict(scale=1, background=0, peak_pos=0.05, sigma=0.005) 
     53def random(): 
     54    import numpy as np 
     55    peak_pos = 10**np.random.uniform(-3, -1) 
     56    sigma = 10**np.random.uniform(-1.3, -0.3)*peak_pos 
     57    scale = 10**np.random.uniform(0, 4) 
     58    pars = dict( 
     59        #background=1e-8, 
     60        scale=scale, 
     61        peak_pos=peak_pos, 
     62        sigam=sigma, 
     63    ) 
     64    return pars 

  • sasmodels/models/gel_fit.c

    ra807206 r48462b0  
    1010    double lorentzian_term = square(q*cor_length); 
    1111    lorentzian_term = 1.0 + ((fractal_dim + 1.0)/3.0)*lorentzian_term; 
    12     lorentzian_term = pow(lorentzian_term, (fractal_dim/2.0) ); 
     12    lorentzian_term = pow(lorentzian_term, fractal_dim/2.0 ); 
    1313 
    1414    // Exponential Term 
    1515    ////////////////////////double d(x[i]*x[i]*rg*rg); 
    1616    double exp_term = square(q*rg); 
    17     exp_term = exp(-1.0*(exp_term/3.0)); 
     17    exp_term = exp(-exp_term/3.0); 
    1818 
    1919    // Scattering Law 

  • sasmodels/models/gel_fit.py

    ra807206 r48462b0  
    7171source = ["gel_fit.c"] 
    7272 
     73def random(): 
     74    import numpy as np 
     75    guinier_scale = 10**np.random.uniform(1, 3) 
     76    lorentz_scale = 10**np.random.uniform(1, 3) 
     77    rg = 10**np.random.uniform(1, 5) 
     78    fractal_dim = np.random.uniform(0, 6) 
     79    cor_length = 10**np.random.uniform(0, 3) 
     80    pars = dict( 
     81        #background=0, 
     82        scale=1, 
     83        guinier_scale=guinier_scale, 
     84        lorentz_scale=lorentz_scale, 
     85        rg=rg, 
     86        fractal_dim=fractal_dim, 
     87        cor_length=cor_length 
     88    ) 
     89    return pars 
     90 
    7391demo = dict(background=0.01, 
    7492            guinier_scale=1.7, 

  • sasmodels/models/guinier.py

    r3330bb4 r48462b0  
    3333description = """ 
    3434 I(q) = scale.exp ( - rg^2 q^2 / 3.0 ) 
    35   
     35 
    3636    List of default parameters: 
    3737    scale = scale 
    … …  
    4949""" 
    5050 
     51def random(): 
     52    import numpy as np 
     53    scale = 10**np.random.uniform(1, 4) 
     54    # Note: compare.py has Rg cutoff of 1e-30 at q=1 for guinier, so use that 
     55    # log_10 Ae^(-(q Rg)^2/3) = log_10(A) - (q Rg)^2/ (3 ln 10) > -30 
     56    #   => log_10(A) > Rg^2/(3 ln 10) - 30 
     57    q_max = 1.0 
     58    rg_max = np.sqrt(90*np.log(10) + 3*np.log(scale))/q_max 
     59    rg = 10**np.random.uniform(0, np.log10(rg_max)) 
     60    pars = dict( 
     61        #background=0, 
     62        scale=scale, 
     63        rg=rg, 
     64    ) 
     65    return pars 
     66 
    5167# parameters for demo 
    5268demo = dict(scale=1.0, rg=60.0) 

  • sasmodels/models/guinier_porod.py

    rcdcebf1 r48462b0  
    44and dimensionality of scattering objects, including asymmetric objects 
    55such as rods or platelets, and shapes intermediate between spheres 
    6 and rods or between rods and platelets, and overcomes some of the  
     6and rods or between rods and platelets, and overcomes some of the 
    77deficiencies of the (Beaucage) Unified_Power_Rg model (see Hammouda, 2010). 
    88 
    … …  
    7777                        scale = Guinier Scale 
    7878                        s = Dimension Variable 
    79                         Rg = Radius of Gyration [A]  
     79                        Rg = Radius of Gyration [A] 
    8080                        porod_exp = Porod Exponent 
    8181                        background  = Background [1/cm]""" 
    … …  
    114114Iq.vectorized = True # Iq accepts an array of q values 
    115115 
     116def random(): 
     117    import numpy as np 
     118    rg = 10**np.random.uniform(1, 5) 
     119    s = np.random.uniform(0, 3) 
     120    porod_exp = s + np.random.uniform(0, 3) 
     121    pars = dict( 
     122        #scale=1, background=0, 
     123        rg=rg, 
     124        s=s, 
     125        porod_exp=porod_exp, 
     126    ) 
     127    return pars 
     128 
    116129demo = dict(scale=1.5, background=0.5, rg=60, s=1.0, porod_exp=3.0) 
    117130 

  • sasmodels/models/hardsphere.py

    r40a87fa r8f04da4  
    7575Iq = r""" 
    7676      double D,A,B,G,X,X2,X4,S,C,FF,HARDSPH; 
    77       // these are c compiler instructions, can also put normal code inside the "if else" structure  
     77      // these are c compiler instructions, can also put normal code inside the "if else" structure 
    7878      #if FLOAT_SIZE > 4 
    7979      // double precision    orig had 0.2, don't call the variable cutoff as PAK already has one called that! Must use UPPERCASE name please. 
    … …  
    8282      #else 
    8383      // 0.1 bad, 0.2 OK, 0.3 good, 0.4 better, 0.8 no good 
    84       #define CUTOFFHS 0.4   
     84      #define CUTOFFHS 0.4 
    8585      #endif 
    8686 
    … …  
    110110      if(X < CUTOFFHS) { 
    111111      // RKH Feb 2016, use Taylor series expansion for small X 
    112       // else no obvious way to rearrange the equations to avoid needing a very high number of significant figures.  
    113       // Series expansion found using Mathematica software. Numerical test in .xls showed terms to X^2 are sufficient  
    114       // for 5 or 6 significant figures, but I put the X^4 one in anyway  
     112      // else no obvious way to rearrange the equations to avoid needing a very high number of significant figures. 
     113      // Series expansion found using Mathematica software. Numerical test in .xls showed terms to X^2 are sufficient 
     114      // for 5 or 6 significant figures, but I put the X^4 one in anyway 
    115115            //FF = 8*A +6*B + 4*G - (0.8*A +2.0*B/3.0 +0.5*G)*X2 +(A/35. +B/40. +G/50.)*X4; 
    116116            // refactoring the polynomial makes it very slightly faster (0.5 not 0.6 msec) 
    … …  
    153153   """ 
    154154 
     155def random(): 
     156    import numpy as np 
     157    pars = dict( 
     158        scale=1, background=0, 
     159        radius_effective=10**np.random.uniform(1, 4), 
     160        volfraction=10**np.random.uniform(-2, 0),  # high volume fraction 
     161    ) 
     162    return pars 
     163 
    155164# ER defaults to 0.0 
    156165# VR defaults to 1.0 

  • sasmodels/models/hayter_msa.py

    r5702f52 r8f04da4  
    6464        Routine takes absolute value of charge, use HardSphere if charge 
    6565        goes to zero. 
    66         In sasview the effective radius and volume fraction may be calculated  
     66        In sasview the effective radius and volume fraction may be calculated 
    6767        from the parameters used in P(Q). 
    6868""" 
    … …  
    8989# ER defaults to 0.0 
    9090# VR defaults to 1.0 
     91 
     92def random(): 
     93    import numpy as np 
     94    # TODO: too many failures for random hayter_msa parameters 
     95    pars = dict( 
     96        scale=1, background=0, 
     97        radius_effective=10**np.random.uniform(1, 4.7), 
     98        volfraction=10**np.random.uniform(-2, 0),  # high volume fraction 
     99        charge=min(int(10**np.random.uniform(0, 1.3)+0.5), 200), 
     100        temperature=10**np.random.uniform(0, np.log10(450)), # max T = 450 
     101        #concentration_salt=10**np.random.uniform(-3, 1), 
     102        dialectconst=10**np.random.uniform(0, 6), 
     103        #charge=10, 
     104        #temperature=318.16, 
     105        concentration_salt=0.0, 
     106        #dielectconst=71.08, 
     107    ) 
     108    return pars 
    91109 
    92110# default parameter set,  use  compare.sh -midQ -linear 

  • sasmodels/models/hollow_cylinder.py

    rf102a96 r8f04da4  
    5454 
    5555* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    56 * **Last Modified by:** Richard Heenan **Date:** October 06, 2016  
     56* **Last Modified by:** Richard Heenan **Date:** October 06, 2016 
    5757   (reparametrised to use thickness, not outer radius) 
    5858* **Last Reviewed by:** Richard Heenan **Date:** October 06, 2016 
    … …  
    121121    return vol_shell, vol_total 
    122122 
     123def random(): 
     124    import numpy as np 
     125    length = 10**np.random.uniform(1, 4.7) 
     126    outer_radius = 10**np.random.uniform(1, 4.7) 
     127    # Use a distribution with a preference for thin shell or thin core 
     128    # Avoid core,shell radii < 1 
     129    thickness = np.random.beta(0.5, 0.5)*(outer_radius-2) + 1 
     130    radius = outer_radius - thickness 
     131    pars = dict( 
     132        length=length, 
     133        radius=radius, 
     134        thickness=thickness, 
     135    ) 
     136    return pars 
     137 
    123138# parameters for demo 
    124139demo = dict(scale=1.0, background=0.0, length=400.0, radius=20.0, 

  • sasmodels/models/hollow_rectangular_prism.py

    rab2aea8 r8f04da4  
    146146 
    147147 
     148def random(): 
     149    import numpy as np 
     150    a, b, c = 10**np.random.uniform(1, 4.7, size=3) 
     151    # Thickness is limited to 1/2 the smallest dimension 
     152    # Use a distribution with a preference for thin shell or thin core 
     153    # Avoid core,shell radii < 1 
     154    min_dim = 0.5*min(a, b, c) 
     155    thickness = np.random.beta(0.5, 0.5)*(min_dim-2) + 1 
     156    #print(a, b, c, thickness, thickness/min_dim) 
     157    pars = dict( 
     158        length_a=a, 
     159        b2a_ratio=b/a, 
     160        c2a_ratio=c/a, 
     161        thickness=thickness, 
     162    ) 
     163    return pars 
     164 
     165 
    148166# parameters for demo 
    149167demo = dict(scale=1, background=0, 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.py

    rab2aea8 r31df0c9  
    127127 
    128128 
     129def random(): 
     130    import numpy as np 
     131    a, b, c = 10**np.random.uniform(1, 4.7, size=3) 
     132    pars = dict( 
     133        length_a=a, 
     134        b2a_ratio=b/a, 
     135        c2a_ratio=c/a, 
     136    ) 
     137    return pars 
     138 
     139 
    129140# parameters for demo 
    130141demo = dict(scale=1, background=0, 

  • sasmodels/models/lamellar.py

    r40a87fa r1511c37c  
    8989    """ 
    9090 
     91def random(): 
     92    import numpy as np 
     93    thickness = 10**np.random.uniform(1, 4) 
     94    pars = dict( 
     95        thickness=thickness, 
     96    ) 
     97    return pars 
     98 
    9199# ER defaults to 0.0 
    92100# VR defaults to 1.0 

  • sasmodels/models/lamellar_hg.py

    ra807206 r1511c37c  
    9898    """ 
    9999 
     100def random(): 
     101    import numpy as np 
     102    thickness = 10**np.random.uniform(1, 4) 
     103    length_head = thickness * np.random.uniform(0, 1) 
     104    length_tail = thickness - length_head 
     105    pars = dict( 
     106        length_head=length_head, 
     107        length_tail=length_tail, 
     108    ) 
     109    return pars 
     110 
    100111# ER defaults to 0.0 
    101112# VR defaults to 1.0 

  • sasmodels/models/lamellar_hg_stack_caille.py

    ra57b31d r1511c37c  
    123123# VR defaults to 1.0 
    124124 
     125def random(): 
     126    import numpy as np 
     127    total_thickness = 10**np.random.uniform(2, 4.7) 
     128    Nlayers = np.random.randint(2, 200) 
     129    d_spacing = total_thickness / Nlayers 
     130    thickness = d_spacing * np.random.uniform(0, 1) 
     131    length_head = thickness * np.random.uniform(0, 1) 
     132    length_tail = thickness - length_head 
     133    Caille_parameter = np.random.uniform(0, 0.8) 
     134    pars = dict( 
     135        length_head=length_head, 
     136        length_tail=length_tail, 
     137        Nlayers=Nlayers, 
     138        d_spacing=d_spacing, 
     139        Caille_parameter=Caille_parameter, 
     140    ) 
     141    return pars 
     142 
    125143demo = dict( 
    126144    scale=1, background=0, 

  • sasmodels/models/lamellar_stack_caille.py

    ra57b31d r1511c37c  
    9898source = ["lamellar_stack_caille.c"] 
    9999 
     100def random(): 
     101    import numpy as np 
     102    total_thickness = 10**np.random.uniform(2, 4.7) 
     103    Nlayers = np.random.randint(2, 200) 
     104    d_spacing = total_thickness / Nlayers 
     105    thickness = d_spacing * np.random.uniform(0, 1) 
     106    Caille_parameter = np.random.uniform(0, 0.8) 
     107    pars = dict( 
     108        thickness=thickness, 
     109        Nlayers=Nlayers, 
     110        d_spacing=d_spacing, 
     111        Caille_parameter=Caille_parameter, 
     112    ) 
     113    return pars 
     114 
    100115# No volume normalization despite having a volume parameter 
    101116# This should perhaps be volume normalized? 

  • sasmodels/models/lamellar_stack_paracrystal.py

    ra0168e8 r8f04da4  
    130130form_volume = """ 
    131131    return 1.0; 
    132     """ 
     132""" 
     133 
     134def random(): 
     135    import numpy as np 
     136    total_thickness = 10**np.random.uniform(2, 4.7) 
     137    Nlayers = np.random.randint(2, 200) 
     138    d_spacing = total_thickness / Nlayers 
     139    thickness = d_spacing * np.random.uniform(0, 1) 
     140    # Let polydispersity peak around 15%; 95% < 0.4; max=100% 
     141    sigma_d = np.random.beta(1.5, 7) 
     142    pars = dict( 
     143        thickness=thickness, 
     144        Nlayers=Nlayers, 
     145        d_spacing=d_spacing, 
     146        sigma_d=sigma_d, 
     147    ) 
     148    return pars 
    133149 
    134150# ER defaults to 0.0 

  • sasmodels/models/line.py

    r3330bb4 r48462b0  
    6868Iqxy.vectorized = True  # Iqxy accepts an array of qx qy values 
    6969 
     70def random(): 
     71    import numpy as np 
     72    scale = 10**np.random.uniform(0, 3) 
     73    slope = np.random.uniform(-1, 1)*1e2 
     74    offset = 1e-5 + (0 if slope > 0 else -slope) 
     75    intercept = 10**np.random.uniform(0, 1) + offset 
     76    pars = dict( 
     77        #background=0, 
     78        scale=scale, 
     79        slope=slope, 
     80        intercept=intercept, 
     81    ) 
     82    return pars 
    7083 
    7184tests = [ 

  • sasmodels/models/linear_pearls.py

    r925ad6e r8f04da4  
    6565source = ["lib/sas_3j1x_x.c", "linear_pearls.c"] 
    6666 
    67 demo = dict(scale=1.0, background=0.0, 
    68             radius=80.0, 
    69             edge_sep=350.0, 
    70             num_pearls=3, 
    71             sld=1.0, 
    72             sld_solvent=6.3) 
     67def random(): 
     68    import numpy as np 
     69    radius = 10**np.random.uniform(1, 3) # 1 - 1000 
     70    edge_sep = 10**np.random.uniform(0, 3)  # 1 - 1000 
     71    num_pearls = np.round(10**np.random.uniform(0.3, 3)) # 2 - 1000 
     72    pars = dict( 
     73        radius=radius, 
     74        edge_sep=edge_sep, 
     75        num_pearls=num_pearls, 
     76    ) 
     77    return pars 
    7378 
    7479""" 

  • sasmodels/models/lorentz.py

    r2c74c11 r404ebbd  
    3030description = """ 
    3131Model that evaluates a Lorentz (Ornstein-Zernicke) model. 
    32          
     32 
    3333I(q) = scale/( 1 + (q*L)^2 ) + bkd 
    34          
    35 The model has three parameters:  
    36     length     = screening Length 
    37     scale  = scale factor 
    38     background    = incoherent background 
     34 
     35The model has three parameters: 
     36    length = screening Length 
     37    scale = scale factor 
     38    background = incoherent background 
    3939""" 
    4040category = "shape-independent" 
    … …  
    4848""" 
    4949 
     50def random(): 
     51    import numpy as np 
     52    pars = dict( 
     53        #background=0, 
     54        scale=10**np.random.uniform(1, 4), 
     55        cor_length=10**np.random.uniform(0, 3), 
     56    ) 
     57    return pars 
     58 
    5059# parameters for demo 
    5160demo = dict(scale=1.0, background=0.0, cor_length=50.0) 

  • sasmodels/models/mass_fractal.py

    r925ad6e r4553dae  
    2828.. math:: 
    2929 
    30     scale = scale\_factor \times NV^2(\rho_{particle} - \rho_{solvent})^2 
     30    scale = scale\_factor \times NV^2(\rho_\text{particle} - \rho_\text{solvent})^2 
    3131 
    3232.. math:: 
    … …  
    3535 
    3636where $R$ is the radius of the building block, $D_m$ is the **mass** fractal 
    37 dimension, | \zeta\|  is the cut-off length, $\rho_{solvent}$ is the scattering 
    38 length density of the solvent, 
    39 and $\rho_{particle}$ is the scattering length density of particles. 
     37dimension, $\zeta$  is the cut-off length, $\rho_\text{solvent}$ is the scattering 
     38length density of the solvent, and $\rho_\text{particle}$ is the scattering 
     39length density of particles. 
    4040 
    4141.. note:: 
    4242 
    4343    The mass fractal dimension ( $D_m$ ) is only 
    44     valid if $0 < mass\_dim < 6$. It is also only valid over a limited 
     44    valid if $1 < mass\_dim < 6$. It is also only valid over a limited 
    4545    $q$ range (see the reference for details). 
    4646 
    … …  
    8888source = ["lib/sas_3j1x_x.c", "lib/sas_gamma.c", "mass_fractal.c"] 
    8989 
     90def random(): 
     91    import numpy as np 
     92    radius = 10**np.random.uniform(0.7, 4) 
     93    cutoff_length = 10**np.random.uniform(0.7, 2)*radius 
     94    # TODO: fractal dimension should range from 1 to 6 
     95    fractal_dim_mass = 2*np.random.beta(3, 4) + 1 
     96    Vf = 10**np.random.uniform(-4, -1) 
     97    pars = dict( 
     98        #background=0, 
     99        scale=1, #1e5*Vf/radius**(fractal_dim_mass), 
     100        radius=radius, 
     101        cutoff_length=cutoff_length, 
     102        fractal_dim_mass=fractal_dim_mass, 
     103    ) 
     104    return pars 
     105 
    90106demo = dict(scale=1, background=0, 
    91107            radius=10.0, 

  • sasmodels/models/mass_surface_fractal.py

    r6d96b66 r232bb12  
    8989source = ["mass_surface_fractal.c"] 
    9090 
     91def random(): 
     92    import numpy as np 
     93    fractal_dim = np.random.uniform(0, 6) 
     94    surface_portion = np.random.uniform(0, 1) 
     95    fractal_dim_surf = fractal_dim*surface_portion 
     96    fractal_dim_mass = fractal_dim - fractal_dim_surf 
     97    rg_cluster = 10**np.random.uniform(1, 5) 
     98    rg_primary = rg_cluster*10**np.random.uniform(-4, -1) 
     99    scale = 10**np.random.uniform(2, 5) 
     100    pars = dict( 
     101        #background=0, 
     102        scale=scale, 
     103        fractal_dim_mass=fractal_dim_mass, 
     104        fractal_dim_surf=fractal_dim_surf, 
     105        rg_cluster=rg_cluster, 
     106        rg_primary=rg_primary, 
     107    ) 
     108    return pars 
     109 
     110 
    91111demo = dict(scale=1, background=0, 
    92112            fractal_dim_mass=1.8, 

  • sasmodels/models/mono_gauss_coil.py

    r3330bb4 r404ebbd  
    6262 
    6363description = """ 
    64     Evaluates the scattering from  
     64    Evaluates the scattering from 
    6565    monodisperse polymer chains. 
    6666    """ 
    … …  
    8686Iq.vectorized = True # Iq accepts an array of q values 
    8787 
     88def random(): 
     89    import numpy as np 
     90    rg = 10**np.random.uniform(0, 4) 
     91    #rg = 1e3 
     92    pars = dict( 
     93        #scale=1, background=0, 
     94        i_zero=1e7, # i_zero is a simple scale 
     95        rg=rg, 
     96    ) 
     97    return pars 
     98 
    8899demo = dict(scale=1.0, i_zero=70.0, rg=75.0, background=0.0) 
    89100 

  • sasmodels/models/multilayer_vesicle.py

    r870a2f4 r870a2f4  
    150150    return radius + n_shells * (thick_shell + thick_solvent) - thick_solvent 
    151151 
    152 demo = dict(scale=1, background=0, 
    153             volfraction=0.05, 
    154             radius=60.0, 
    155             thick_shell=10.0, 
    156             thick_solvent=10.0, 
    157             sld_solvent=6.4, 
    158             sld=0.4, 
    159             n_shells=2.0) 
     152def random(): 
     153    import numpy as np 
     154    volfraction = 10**np.random.uniform(-3, -0.5)  # scale from 0.1% to 30% 
     155    radius = 10**np.random.uniform(0, 2.5) # core less than 300 A 
     156    total_thick = 10**np.random.uniform(2, 4) # up to 10000 A of shells 
     157    # random number of shells, with shell+solvent thickness > 10 A 
     158    n_shells = int(10**np.random.uniform(0, np.log10(total_thick)-1)+0.5) 
     159    # split total shell thickness with preference for shell over solvent; 
     160    # make sure that shell thickness is at least 1 A 
     161    one_thick = total_thick/n_shells 
     162    thick_solvent = 10**np.random.uniform(-2, 0)*(one_thick - 1) 
     163    thick_shell = one_thick - thick_solvent 
     164    pars = dict( 
     165        scale=1, 
     166        volfraction=volfraction, 
     167        radius=radius, 
     168        thick_shell=thick_shell, 
     169        thick_solvent=thick_solvent, 
     170        n_shells=n_shells, 
     171    ) 
     172    return pars 
    160173 
    161174tests = [ 

  • sasmodels/models/parallelepiped.py

    r34a9e4e r8f04da4  
    2222.. note:: 
    2323 
    24 The three dimensions of the parallelepiped (strictly here a cuboid) may be given in  
     24The three dimensions of the parallelepiped (strictly here a cuboid) may be given in 
    2525$any$ size order. To avoid multiple fit solutions, especially 
    26 with Monte-Carlo fit methods, it may be advisable to restrict their ranges. There may  
    27 be a number of closely similar "best fits", so some trial and error, or fixing of some  
     26with Monte-Carlo fit methods, it may be advisable to restrict their ranges. There may 
     27be a number of closely similar "best fits", so some trial and error, or fixing of some 
    2828dimensions at expected values, may help. 
    2929 
    … …  
    115115 
    116116On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will 
    117 appear. These are actually rotations about axes $\delta_1$ and $\delta_2$ of the parallelepiped, perpendicular to the $a$ x $c$ and $b$ x $c$ faces.  
    118 (When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.) The third orientation distribution, in $\psi$, is  
    119 about the $c$ axis of the particle, perpendicular to the $a$ x $b$ face. Some experimentation may be required to  
    120 understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation  
    121 distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.)  
    122  
    123      
     117appear. These are actually rotations about axes $\delta_1$ and $\delta_2$ of the parallelepiped, perpendicular to the $a$ x $c$ and $b$ x $c$ faces. 
     118(When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.) The third orientation distribution, in $\psi$, is 
     119about the $c$ axis of the particle, perpendicular to the $a$ x $b$ face. Some experimentation may be required to 
     120understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation 
     121distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.) 
     122 
     123 
    124124For a given orientation of the parallelepiped, the 2D form factor is 
    125125calculated as 
    … …  
    241241 
    242242# VR defaults to 1.0 
     243 
     244 
     245def random(): 
     246    import numpy as np 
     247    length = 10**np.random.uniform(1, 4.7, size=3) 
     248    pars = dict( 
     249        length_a=length[0], 
     250        length_b=length[1], 
     251        length_c=length[2], 
     252    ) 
     253    return pars 
     254 
    243255 
    244256# parameters for demo 

  • sasmodels/models/peak_lorentz.py

    r2c74c11 r404ebbd  
    5959Iq.vectorized = True  # Iq accepts an array of q values 
    6060 
     61def random(): 
     62    import numpy as np 
     63    peak_pos = 10**np.random.uniform(-3, -1) 
     64    peak_hwhm = peak_pos * 10**np.random.uniform(-3, 0) 
     65    pars = dict( 
     66        #background=0, 
     67        scale=10**np.random.uniform(2, 6), 
     68        peak_pos=peak_pos, 
     69        peak_hwhm=peak_hwhm, 
     70    ) 
     71    return pars 
     72 
    6173demo = dict(scale=100, background=1.0, 
    6274            peak_pos=0.05, peak_hwhm=0.005) 

  • sasmodels/models/pearl_necklace.py

    r1bd1ea2 r8f04da4  
    117117    return rad_out 
    118118 
     119def random(): 
     120    import numpy as np 
     121    radius = 10**np.random.uniform(1, 3) # 1 - 1000 
     122    thick_string = 10**np.random.uniform(0, np.log10(radius)-1) # 1 - radius/10 
     123    edge_sep = 10**np.random.uniform(0, 3)  # 1 - 1000 
     124    num_pearls = np.round(10**np.random.uniform(0.3, 3)) # 2 - 1000 
     125    pars = dict( 
     126        radius=radius, 
     127        edge_sep=edge_sep, 
     128        thick_string=thick_string, 
     129        num_pearls=num_pearls, 
     130    ) 
     131    return pars 
     132 
    119133# parameters for demo 
    120134demo = dict(scale=1, background=0, radius=80.0, edge_sep=350.0, 

  • sasmodels/models/poly_gauss_coil.py

    r3330bb4 r404ebbd  
    6666 
    6767description = """ 
    68     Evaluates the scattering from  
     68    Evaluates the scattering from 
    6969    polydisperse polymer chains. 
    7070    """ 
    … …  
    9696        p = [ 
    9797            #(-1 - 20*u - 155*u**2 - 580*u**3 - 1044*u**4 - 720*u**5) / 2520., 
    98             #( 1 + 14*u + 71*u**2 + 154*u**3 + 120*u**4) / 360., 
     98            #(+1 + 14*u + 71*u**2 + 154*u**3 + 120*u**4) / 360., 
    9999            #(-1 - 9*u - 26*u**2 - 24*u**3) / 60., 
    100             ( 1 + 5*u + 6*u**2) / 12., 
     100            (+1 + 5*u + 6*u**2) / 12., 
    101101            (-1 - 2*u) / 3., 
    102             ( 1 ), 
     102            (+1), 
    103103            ] 
    104104        result = 2.0 * (power(1.0 + u*z, -1.0/u) + z - 1.0) / (1.0 + u) 
    … …  
    108108    return i_zero * result 
    109109Iq.vectorized = True  # Iq accepts an array of q values 
     110 
     111def random(): 
     112    import numpy as np 
     113    rg = 10**np.random.uniform(0, 4) 
     114    #rg = 1e3 
     115    polydispersity = 10**np.random.uniform(0, 3) 
     116    pars = dict( 
     117        #scale=1, background=0, 
     118        i_zero=1e7, # i_zero is a simple scale 
     119        rg=rg, 
     120        polydispersity=polydispersity, 
     121    ) 
     122    return pars 
    110123 
    111124demo = dict(scale=1.0, 

  • sasmodels/models/polymer_excl_volume.py

    r40a87fa r404ebbd  
    108108#   ["name", "units", default, [lower, upper], "type", "description"], 
    109109parameters = [ 
    110     ["rg",        "Ang", 60.0, [0, inf],    "", "Radius of Gyration"], 
    111     ["porod_exp", "",     3.0, [-inf, inf], "", "Porod exponent"], 
     110    ["rg",        "Ang", 60.0, [0, inf], "", "Radius of Gyration"], 
     111    ["porod_exp", "",     3.0, [0, inf], "", "Porod exponent"], 
    112112] 
    113113# pylint: enable=bad-whitespace, line-too-long 
    … …  
    133133Iq.vectorized = True  # Iq accepts an array of q values 
    134134 
     135def random(): 
     136    import numpy as np 
     137    rg = 10**np.random.uniform(0, 4) 
     138    porod_exp = np.random.uniform(1e-3, 6) 
     139    scale = 10**np.random.uniform(1, 5) 
     140    pars = dict( 
     141        #background=0, 
     142        scale=scale, 
     143        rg=rg, 
     144        porod_exp=porod_exp, 
     145    ) 
     146    return pars 
     147 
    135148tests = [ 
    136149    # Accuracy tests based on content in test/polyexclvol_default_igor.txt 

  • sasmodels/models/polymer_micelle.py

    r925ad6e r404ebbd  
    1616which then defines a micelle core of $radius\_core$, which is a separate parameter 
    1717even though it could be directly determined. 
    18 The Gaussian random coil tails, of gyration radius $rg$, are imagined uniformly  
     18The Gaussian random coil tails, of gyration radius $rg$, are imagined uniformly 
    1919distributed around the spherical core, centred at a distance $radius\_core + d\_penetration.rg$ 
    2020from the micelle centre, where $d\_penetration$ is of order unity. 
    … …  
    2727.. math:: 
    2828    P(q) = N^2\beta^2_s\Phi(qR)^2+N\beta^2_cP_c(q)+2N^2\beta_s\beta_cS_{sc}s_c(q)+N(N-1)\beta_c^2S_{cc}(q) 
    29      
     29 
    3030    \beta_s = v\_core(sld\_core - sld\_solvent) 
    31      
     31 
    3232    \beta_c = v\_corona(sld\_corona - sld\_solvent) 
    3333 
    34 where $N = n\_aggreg$, and for the spherical core of radius $R$  
     34where $N = n\_aggreg$, and for the spherical core of radius $R$ 
    3535 
    36 .. math::    
     36.. math:: 
    3737   \Phi(qR)= \frac{\sin(qr) - qr\cos(qr)}{(qr)^3} 
    3838 
    … …  
    4949 
    5050.. math:: 
    51     
     51 
    5252   S_{sc}(q)=\Phi(qR)\psi(Z)\frac{sin(q(R+d.R_g))}{q(R+d.R_g)} 
    53     
     53 
    5454   S_{cc}(q)=\psi(Z)^2\left[\frac{sin(q(R+d.R_g))}{q(R+d.R_g)} \right ]^2 
    55     
     55 
    5656   \psi(Z)=\frac{[1-exp^{-Z}]}{Z} 
    5757 
    … …  
    6060 
    6161$P(q)$ above is multiplied by $ndensity$, and a units conversion of 10^{-13}, so $scale$ 
    62 is likely 1.0 if the scattering data is in absolute units. This model has not yet been  
     62is likely 1.0 if the scattering data is in absolute units. This model has not yet been 
    6363independently validated. 
    6464 
    … …  
    7171""" 
    7272 
    73 from numpy import inf 
     73from numpy import inf, pi 
    7474 
    7575name = "polymer_micelle" 
    … …  
    8080to block copolymer micelles. To work well the Gaussian chains must be much 
    8181smaller than the core, which is often not the case.  Please study the 
    82 reference to Pedersen and full documentation carefully.  
     82reference to Pedersen and full documentation carefully. 
    8383    """ 
    8484 
    … …  
    106106source = ["lib/sas_3j1x_x.c", "polymer_micelle.c"] 
    107107 
    108 demo = dict(scale=1, background=0, 
    109             ndensity=8.94, 
    110             v_core=62624.0, 
    111             v_corona=61940.0, 
    112             sld_solvent=6.4, 
    113             sld_core=0.34, 
    114             sld_corona=0.8, 
    115             radius_core=45.0, 
    116             rg=20.0, 
    117             d_penetration=1.0, 
    118             n_aggreg=6.0) 
    119  
     108def random(): 
     109    import numpy as np 
     110    radius_core = 10**np.random.uniform(1, 3) 
     111    rg = radius_core * 10**np.random.uniform(-2, -0.3) 
     112    d_penetration = np.random.randn()*0.05 + 1 
     113    n_aggreg = np.random.randint(3, 30) 
     114    # volume of head groups is the core volume over the number of groups, 
     115    # with a correction for packing fraction of the head groups. 
     116    v_core = 4*pi/3*radius_core**3/n_aggreg * 0.68 
     117    # Rg^2 for gaussian coil is a^2n/6 => a^2 = 6 Rg^2/n 
     118    # a=2r => r = Rg sqrt(3/2n) 
     119    # v = 4/3 pi r^3 n => v = 4/3 pi Rg^3 (3/2n)^(3/2) n = pi Rg^3 sqrt(6/n) 
     120    tail_segments = np.random.randint(6, 30) 
     121    v_corona = pi * rg**3 * np.sqrt(6/tail_segments) 
     122    V = 4*pi/3*(radius_core + rg)**3 
     123    pars = dict( 
     124        background=0, 
     125        scale=1e7/V, 
     126        ndensity=8.94, 
     127        v_core=v_core, 
     128        v_corona=v_corona, 
     129        radius_core=radius_core, 
     130        rg=rg, 
     131        d_penetration=d_penetration, 
     132        n_aggreg=n_aggreg, 
     133    ) 
     134    return pars 
    120135 
    121136tests = [ 

  • sasmodels/models/porod.py

    r4962519 r232bb12  
    4545Iq.vectorized = True  # Iq accepts an array of q values 
    4646 
     47def random(): 
     48    import numpy as np 
     49    sld, solvent = np.random.uniform(-0.5, 12, size=2) 
     50    radius = 10**np.random.uniform(1, 4.7) 
     51    Vf = 10**np.random.uniform(-3, -1) 
     52    scale = 1e-4 * Vf * 2*np.pi*(sld-solvent)**2/(3*radius) 
     53    pars = dict( 
     54        scale=scale, 
     55    ) 
     56    return pars 
     57 
    4758demo = dict(scale=1.5, background=0.5) 
    4859 

  • sasmodels/models/power_law.py

    r40a87fa r404ebbd  
    5050Iq.vectorized = True  # Iq accepts an array of q values 
    5151 
     52def random(): 
     53    import numpy as np 
     54    power = np.random.uniform(1, 6) 
     55    pars = dict( 
     56        scale=0.1**power*10**np.random.uniform(-4, 2), 
     57        power=power, 
     58    ) 
     59    return pars 
     60 
    5261demo = dict(scale=1.0, power=4.0, background=0.0) 
    5362 

  • sasmodels/models/pringle.py

    r30fbe2e r8f04da4  
    8585    return 0.5 * (ddd) ** (1. / 3.) 
    8686 
    87 demo = dict(background=0.0, 
    88             scale=1.0, 
    89             radius=60.0, 
    90             thickness=10.0, 
    91             alpha=0.001, 
    92             beta=0.02, 
    93             sld=1.0, 
    94             sld_solvent=6.35) 
     87def random(): 
     88    import numpy as np 
     89    alpha, beta = 10**np.random.uniform(-1, 1, size=2) 
     90    radius = 10**np.random.uniform(1, 3) 
     91    thickness = 10**np.random.uniform(0.7, 2) 
     92    pars = dict( 
     93        radius=radius, 
     94        thickness=thickness, 
     95        alpha=alpha, 
     96        beta=beta, 
     97    ) 
     98    return pars 
    9599 
    96100tests = [ 

  • sasmodels/models/raspberry.py

    r8e68ea0 r8f04da4  
    154154source = ["lib/sas_3j1x_x.c", "raspberry.c"] 
    155155 
     156def random(): 
     157    import numpy as np 
     158    # Limit volume fraction to 20% each 
     159    volfraction_lg = 10**np.random.uniform(-3, -0.3) 
     160    volfraction_sm = 10**np.random.uniform(-3, -0.3) 
     161    # Prefer most particles attached (peak near 60%), but not all or none 
     162    surface_fraction = np.random.beta(6, 4) 
     163    radius_lg = 10**np.random.uniform(1.7, 4.7)  # 500 - 50000 A 
     164    radius_sm = 10**np.random.uniform(-3, -0.3)*radius_lg  # 0.1% - 20% 
     165    penetration = np.random.beta(1, 10) # up to 20% pen. for 90% of examples 
     166    pars = dict( 
     167        volfraction_lg=volfraction_lg, 
     168        volfraction_sm=volfraction_sm, 
     169        surface_fraction=surface_fraction, 
     170        radius_lg=radius_lg, 
     171        radius_sm=radius_sm, 
     172        penetration=penetration, 
     173    ) 
     174    return pars 
     175 
    156176# parameters for demo 
    157177demo = dict(scale=1, background=0.001, 

  • sasmodels/models/rectangular_prism.py

    rab2aea8 r31df0c9  
    104104              ["length_a", "Ang", 35, [0, inf], "volume", 
    105105               "Shorter side of the parallelepiped"], 
    106               ["b2a_ratio", "Ang", 1, [0, inf], "volume", 
     106              ["b2a_ratio", "", 1, [0, inf], "volume", 
    107107               "Ratio sides b/a"], 
    108               ["c2a_ratio", "Ang", 1, [0, inf], "volume", 
     108              ["c2a_ratio", "", 1, [0, inf], "volume", 
    109109               "Ratio sides c/a"], 
    110110             ] 
    … …  
    125125    return 0.5 * (ddd) ** (1. / 3.) 
    126126 
     127def random(): 
     128    import numpy as np 
     129    a, b, c = 10**np.random.uniform(1, 4.7, size=3) 
     130    pars = dict( 
     131        length_a=a, 
     132        b2a_ratio=b/a, 
     133        c2a_ratio=c/a, 
     134    ) 
     135    return pars 
    127136 
    128137# parameters for demo 

  • sasmodels/models/sc_paracrystal.py

    r69e1afc r8f04da4  
    138138source = ["lib/sas_3j1x_x.c", "lib/sphere_form.c", "lib/gauss150.c", "sc_paracrystal.c"] 
    139139 
    140 demo = dict(scale=1, background=0, 
    141             dnn=220.0, 
    142             d_factor=0.06, 
    143             radius=40.0, 
    144             sld=3.0, 
    145             sld_solvent=6.3, 
    146             theta=0.0, 
    147             phi=0.0, 
    148             psi=0.0) 
     140def random(): 
     141    import numpy as np 
     142    # copied from bcc_paracrystal 
     143    radius = 10**np.random.uniform(1.3, 4) 
     144    d_factor = 10**np.random.uniform(-2, -0.7)  # sigma_d in 0.01-0.7 
     145    dnn_fraction = np.random.beta(a=10, b=1) 
     146    dnn = radius*4/np.sqrt(4)/dnn_fraction 
     147    pars = dict( 
     148        #sld=1, sld_solvent=0, scale=1, background=1e-32, 
     149        dnn=dnn, 
     150        d_factor=d_factor, 
     151        radius=radius, 
     152    ) 
     153    return pars 
    149154 
    150155tests = [ 
    … …  
    152157    [{}, 0.001, 10.3048], 
    153158    [{}, 0.215268, 0.00814889], 
    154     [{}, (0.414467), 0.001313289], 
    155     [{'theta':10.0,'phi':20,'psi':30.0},(0.045,-0.035),18.0397138402 ], 
    156     [{'theta':10.0,'phi':20,'psi':30.0},(0.023,0.045),0.0177333171285 ] 
     159    [{}, 0.414467, 0.001313289], 
     160    [{'theta': 10.0, 'phi': 20, 'psi': 30.0}, (0.045, -0.035), 18.0397138402], 
     161    [{'theta': 10.0, 'phi': 20, 'psi': 30.0}, (0.023, 0.045), 0.0177333171285], 
    157162    ] 
    158163 

  • sasmodels/models/sphere.py

    r925ad6e r97d89af  
    8686# VR defaults to 1.0 
    8787 
    88 demo = dict(scale=1, background=0, 
    89             sld=6, sld_solvent=1, 
    90             radius=120, 
    91             radius_pd=.2, radius_pd_n=45) 
     88def random(): 
     89    import numpy as np 
     90    radius = 10**np.random.uniform(1.3, 4) 
     91    pars = dict( 
     92        radius=radius, 
     93    ) 
     94    return pars 
    9295 
    9396tests = [ 

  • sasmodels/models/spinodal.py

    rbba9361 r48462b0  
    33---------- 
    44 
    5 This model calculates the SAS signal of a phase separating solution under spinodal decomposition.  
     5This model calculates the SAS signal of a phase separating solution under spinodal decomposition. 
    66The scattering intensity $I(q)$ is calculated as 
    77 
    8 .. math::  
     8.. math:: 
    99    I(q) = I_{max}\frac{(1+\gamma/2)x^2}{\gamma/2+x^{2+\gamma}}+B 
    1010 
    11 where $x=q/q_0$ and $B$ is a flat background. The characteristic structure length  
    12 scales with the correlation peak at $q_0$. The exponent $\gamma$ is equal to  
     11where $x=q/q_0$ and $B$ is a flat background. The characteristic structure length 
     12scales with the correlation peak at $q_0$. The exponent $\gamma$ is equal to 
    1313$d+1$ with d the dimensionality of the off-critical concentration mixtures. 
    14 A transition to $\gamma=2d$ is seen near the percolation threshold into the  
     14A transition to $\gamma=2d$ is seen near the percolation threshold into the 
    1515critical concentration regime. 
    1616 
    … …  
    1818---------- 
    1919 
    20 H. Furukawa. Dynamics-scaling theory for phase-separating unmixing mixtures:  
     20H. Furukawa. Dynamics-scaling theory for phase-separating unmixing mixtures: 
    2121Growth rates of droplets and scaling properties of autocorrelation functions. Physica A 123,497 (1984). 
    2222 
    … …  
    2525 
    2626* **Author:** Dirk Honecker **Date:** Oct 7, 2016 
    27 * **Last Modified by:**  
    28 * **Last Reviewed by:**  
     27* **Last Modified by:** 
     28* **Last Reviewed by:** 
    2929""" 
    3030 
    … …  
    4646# pylint: disable=bad-whitespace, line-too-long 
    4747#             ["name", "units", default, [lower, upper], "type", "description"], 
    48 parameters = [["scale",    "",      1.0, [-inf, inf], "", "Scale factor"], 
    49               ["gamma",      "",    3.0, [-inf, inf], "", "Exponent"], 
     48parameters = [["gamma",      "",    3.0, [-inf, inf], "", "Exponent"], 
    5049              ["q_0",  "1/Ang",     0.1, [-inf, inf], "", "Correlation peak position"] 
    5150             ] 
    … …  
    5352 
    5453def Iq(q, 
    55        scale=1.0, 
    5654       gamma=3.0, 
    5755       q_0=0.1): 
    5856    """ 
    5957    :param q:              Input q-value 
    60     :param scale:          Scale factor 
    6158    :param gamma:          Exponent 
    6259    :param q_0:            Correlation peak position 
    6360    :return:               Calculated intensity 
    6461    """ 
    65      
     62 
    6663    with errstate(divide='ignore'): 
    6764        x = q/q_0 
    68         inten = scale * ((1 + gamma / 2) * x ** 2) / (gamma / 2 + x ** (2 + gamma))  
     65        inten = ((1 + gamma / 2) * x ** 2) / (gamma / 2 + x ** (2 + gamma)) 
    6966    return inten 
    7067Iq.vectorized = True  # Iq accepts an array of q values 
    7168 
     69def random(): 
     70    import numpy as np 
     71    pars = dict( 
     72        scale=10**np.random.uniform(1, 3), 
     73        gamma=np.random.uniform(0, 6), 
     74        q_0=10**np.random.uniform(-3, -1), 
     75    ) 
     76    return pars 
     77 
    7278demo = dict(scale=1, background=0, 
    7379            gamma=1, q_0=0.1) 

  • sasmodels/models/squarewell.py

    r3330bb4 r8f04da4  
    9191        double S,C,SL,CL; 
    9292        x= q; 
    93          
     93 
    9494        req = radius_effective; 
    9595        phis = volfraction; 
    9696        edibkb = welldepth; 
    9797        lambda = wellwidth; 
    98          
     98 
    9999        sigma = req*2.; 
    100100        eta = phis; 
    … …  
    110110        beta = -(eta/3.0) * ( 18. + 20.*eta - 12.*eta2 + eta4 )/etam14; 
    111111        gamm = 0.5*eta*( sk + eta3*(eta-4.) )/etam14; 
    112          
     112 
    113113        //  calculate the structure factor 
    114          
     114 
    115115        sk = x*sigma; 
    116116        sk2 = sk*sk; 
    … …  
    125125        ck =  -24.0*eta*( t1 + t2 + t3 + t4 )/sk3/sk3; 
    126126        struc  = 1./(1.-ck); 
    127          
     127 
    128128        return(struc); 
    129129""" 
    … …  
    132132# VR defaults to 1.0 
    133133 
     134def random(): 
     135    import numpy as np 
     136    pars = dict( 
     137        scale=1, background=0, 
     138        radius_effective=10**np.random.uniform(1, 4.7), 
     139        volfraction=np.random.uniform(0.00001, 0.08), 
     140        welldepth=np.random.uniform(0, 1.5), 
     141        wellwidth=np.random.uniform(1, 1.2), 
     142    ) 
     143    return pars 
     144 
    134145demo = dict(radius_effective=50, volfraction=0.04, welldepth=1.5, 
    135146            wellwidth=1.2, radius_effective_pd=0, radius_effective_pd_n=0) 
    136147# 
    137148tests = [ 
    138     [{'scale': 1.0, 'background' : 0.0, 'radius_effective' : 50.0, 
    139       'volfraction' : 0.04,'welldepth' : 1.5, 'wellwidth' : 1.2, 
    140       'radius_effective_pd' : 0}, 
    141      [0.001], [0.97665742]], 
     149    [{'scale': 1.0, 'background': 0.0, 'radius_effective': 50.0, 
     150      'volfraction': 0.04,'welldepth': 1.5, 'wellwidth': 1.2, 
     151      'radius_effective_pd': 0}, [0.001], [0.97665742]], 
    142152    ] 
    143153# ADDED by: converting from sasview RKH  ON: 16Mar2016 

  • sasmodels/models/stacked_disks.py

    r9d50a1e r8f04da4  
    137137 
    138138source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "stacked_disks.c"] 
     139 
     140def random(): 
     141    import numpy as np 
     142    radius = 10**np.random.uniform(1, 4.7) 
     143    total_stack = 10**np.random.uniform(1, 4.7) 
     144    n_stacking = int(10**np.random.uniform(0, np.log10(total_stack)-1) + 0.5) 
     145    d = total_stack/n_stacking 
     146    thick_core = np.random.uniform(0, d-2)  # at least 1 A for each layer 
     147    thick_layer = (d - thick_core)/2 
     148    # Let polydispersity peak around 15%; 95% < 0.4; max=100% 
     149    sigma_d = d * np.random.beta(1.5, 7) 
     150    pars = dict( 
     151        thick_core=thick_core, 
     152        thick_layer=thick_layer, 
     153        radius=radius, 
     154        n_stacking=n_stacking, 
     155        sigma_d=sigma_d, 
     156    ) 
     157    return pars 
    139158 
    140159demo = dict(background=0.001, 

  • sasmodels/models/star_polymer.py

    rd439007 rd439007  
    8383source = ["star_polymer.c"] 
    8484 
    85 demo = dict(scale=1, background=0, 
    86             rg_squared=100.0, 
    87             arms=3.0) 
     85def random(): 
     86    import numpy as np 
     87    pars = dict( 
     88        #background=0, 
     89        scale=10**np.random.uniform(1, 4), 
     90        rg_squared=10**np.random.uniform(1, 8), 
     91        arms=np.random.uniform(1, 6), 
     92    ) 
     93    return pars 
    8894 
    8995tests = [[{'rg_squared': 2.0, 

  • sasmodels/models/stickyhardsphere.py

    r40a87fa r8f04da4  
    9898    ] 
    9999 
     100def random(): 
     101    import numpy as np 
     102    pars = dict( 
     103        scale=1, background=0, 
     104        radius_effective=10**np.random.uniform(1, 4.7), 
     105        volfraction=np.random.uniform(0.00001, 0.74), 
     106        perturb=10**np.random.uniform(-2, -1), 
     107        stickiness=np.random.uniform(0, 1), 
     108    ) 
     109    return pars 
     110 
    100111# No volume normalization despite having a volume parameter 
    101112# This should perhaps be volume normalized? 

  • sasmodels/models/surface_fractal.py

    r925ad6e r48462b0  
    7474source = ["lib/sas_3j1x_x.c", "lib/sas_gamma.c", "surface_fractal.c"] 
    7575 
    76 demo = dict(scale=1, background=1e-5, 
    77             radius=10, fractal_dim_surf=2.0, cutoff_length=500) 
     76def random(): 
     77    import numpy as np 
     78    radius = 10**np.random.uniform(1, 4) 
     79    fractal_dim_surf = np.random.uniform(1, 3-1e-6) 
     80    cutoff_length = 1e6  # Sets the low q limit; keep it big for sim 
     81    pars = dict( 
     82        #background=0, 
     83        scale=1, 
     84        radius=radius, 
     85        fractal_dim_surf=fractal_dim_surf, 
     86        cutoff_length=cutoff_length, 
     87    ) 
     88    return pars 
    7889 
    7990tests = [ 

  • sasmodels/models/teubner_strey.py

    r8393c74 r48462b0  
    102102Iq.vectorized = True  # Iq accepts an array of q values 
    103103 
     104def random(): 
     105    import numpy as np 
     106    d = 10**np.random.uniform(1, 4) 
     107    xi = 10**np.random.uniform(-0.3, 2)*d 
     108    pars = dict( 
     109        #background=0, 
     110        scale=100, 
     111        volfraction_a=10**np.random.uniform(-3, 0), 
     112        sld_a=np.random.uniform(-0.5, 12), 
     113        sld_b=np.random.uniform(-0.5, 12), 
     114        d=d, 
     115        xi=xi, 
     116    ) 
     117    return pars 
     118 
    104119demo = dict(scale=1, background=0, volfraction_a=0.5, 
    105                      sld_a=0.3, sld_b=6.3, 
    106                      d=100.0, xi=30.0) 
     120            sld_a=0.3, sld_b=6.3, 
     121            d=100.0, xi=30.0) 
    107122tests = [[{}, 0.06, 41.5918888453]] 

  • sasmodels/models/triaxial_ellipsoid.py

    r34a9e4e r31df0c9  
    7171small angle diffraction situations there may be a number of closely similar "best fits", 
    7272so some trial and error, or fixing of some radii at expected values, may help. 
    73      
     73 
    7474To provide easy access to the orientation of the triaxial ellipsoid, 
    7575we define the axis of the cylinder using the angles $\theta$, $\phi$ 
    … …  
    7979.. figure:: img/elliptical_cylinder_angle_definition.png 
    8080 
    81     Definition of angles for oriented triaxial ellipsoid, where radii are for illustration here  
     81    Definition of angles for oriented triaxial ellipsoid, where radii are for illustration here 
    8282    $a < b << c$ and angle $\Psi$ is a rotation around the axis of the particle. 
    8383 
    84 For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,  
     84For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data, 
    8585see the :ref:`elliptical-cylinder` model for further information. 
    8686 
    … …  
    173173    return ellipsoid_ER(polar, equatorial) 
    174174 
     175def random(): 
     176    import numpy as np 
     177    a, b, c = 10**np.random.uniform(1, 4.7, size=3) 
     178    pars = dict( 
     179        radius_equat_minor=a, 
     180        radius_equat_major=b, 
     181        radius_polar=c, 
     182    ) 
     183    return pars 
     184 
     185 
    175186demo = dict(scale=1, background=0, 
    176187            sld=6, sld_solvent=1, 

  • sasmodels/models/two_lorentzian.py

    r2c74c11 r48462b0  
    9393Iq.vectorized = True  # Iq accepts an array of q values 
    9494 
     95def random(): 
     96    import numpy as np 
     97    scale = 10**np.random.uniform(0, 4, 2) 
     98    length = 10**np.random.uniform(1, 4, 2) 
     99    expon = np.random.uniform(1, 6, 2) 
     100 
     101    pars = dict( 
     102        #background=0, 
     103        scale=1, # scale provided in model 
     104        lorentz_scale_1=scale[0], 
     105        lorentz_length_1=length[0], 
     106        lorentz_exp_1=expon[0], 
     107        lorentz_scale_2=scale[1], 
     108        lorentz_length_2=length[1], 
     109        lorentz_exp_2=expon[1], 
     110    ) 
     111    return pars 
     112 
    95113 
    96114demo = dict(scale=1, background=0.1, 

  • sasmodels/models/two_power_law.py

    rbb4b509 r48462b0  
    9898Iq.vectorized = True  # Iq accepts an array of q values 
    9999 
     100def random(): 
     101    import numpy as np 
     102    coefficient_1 = 1 
     103    crossover = 10**np.random.uniform(-3, -1) 
     104    power_1 = np.random.uniform(1, 6) 
     105    power_2 = np.random.uniform(1, 6) 
     106    pars = dict( 
     107        scale=1, #background=0, 
     108        coefficient_1=coefficient_1, 
     109        crossover=crossover, 
     110        power_1=power_1, 
     111        power_2=power_2, 
     112    ) 
     113    return pars 
     114 
    100115demo = dict(scale=1, background=0.0, 
    101116            coefficent_1=1.0, 

  • sasmodels/models/unified_power_Rg.py

    r66ca2a6 r48462b0  
    33---------- 
    44 
    5 This model employs the empirical multiple level unified Exponential/Power-law  
    6 fit method developed by Beaucage. Four functions are included so that 1, 2, 3,  
    7 or 4 levels can be used. In addition a 0 level has been added which simply  
     5This model employs the empirical multiple level unified Exponential/Power-law 
     6fit method developed by Beaucage. Four functions are included so that 1, 2, 3, 
     7or 4 levels can be used. In addition a 0 level has been added which simply 
    88calculates 
    99 
    … …  
    1616(Debye equation), ellipsoidal particles, etc. 
    1717 
    18 The model works best for mass fractal systems characterized by Porod exponents  
    19 between 5/3 and 3. It should not be used for surface fractal systems. Hammouda  
    20 (2010) has pointed out a deficiency in the way this model handles the  
    21 transitioning between the Guinier and Porod regimes and which can create  
     18The model works best for mass fractal systems characterized by Porod exponents 
     19between 5/3 and 3. It should not be used for surface fractal systems. Hammouda 
     20(2010) has pointed out a deficiency in the way this model handles the 
     21transitioning between the Guinier and Porod regimes and which can create 
    2222artefacts that appear as kinks in the fitted model function. 
    2323 
    … …  
    8888# pylint: disable=bad-whitespace, line-too-long 
    8989parameters = [ 
    90     ["level",     "",     1,      [0, 6], "", "Level number"], 
     90    ["level",     "",     1,      [1, 6], "", "Level number"], 
    9191    ["rg[level]", "Ang",  15.8,   [0, inf], "", "Radius of gyration"], 
    9292    ["power[level]", "",  4,      [-inf, inf], "", "Power"], 
    … …  
    118118Iq.vectorized = True 
    119119 
     120def random(): 
     121    import numpy as np 
     122    from scipy.special import gamma 
     123    level = np.minimum(np.random.poisson(0.5) + 1, 6) 
     124    n = level 
     125    power = np.random.uniform(1.6, 3, n) 
     126    rg = 10**np.random.uniform(1, 5, n) 
     127    G = np.random.uniform(0.1, 10, n)**2 * 10**np.random.uniform(0.3, 3, n) 
     128    B = G * power / rg**power * gamma(power/2) 
     129    scale = 10**np.random.uniform(1, 4) 
     130    pars = dict( 
     131        #background=0, 
     132        scale=scale, 
     133        level=level, 
     134    ) 
     135    pars.update(("power%d"%(k+1), v) for k, v in enumerate(power)) 
     136    pars.update(("rg%d"%(k+1), v) for k, v in enumerate(rg)) 
     137    pars.update(("B%d"%(k+1), v) for k, v in enumerate(B)) 
     138    pars.update(("G%d"%(k+1), v) for k, v in enumerate(G)) 
     139    return pars 
     140 
    120141demo = dict( 
    121142    level=2, 

  • sasmodels/models/vesicle.py

    r925ad6e r48462b0  
    120120    return whole, whole - core 
    121121 
     122def random(): 
     123    import numpy as np 
     124    total_radius = 10**np.random.uniform(1.3, 5) 
     125    radius = total_radius * np.random.uniform(0, 1) 
     126    thickness = total_radius - radius 
     127    volfraction = 10**np.random.uniform(-3, -1) 
     128    pars = dict( 
     129        #background=0, 
     130        scale=1,  # volfraction is part of the model, so scale=1 
     131        radius=radius, 
     132        thickness=thickness, 
     133        volfraction=volfraction, 
     134    ) 
     135    return pars 
    122136 
    123137# parameters for demo 

  • sasmodels/product.py

    r6a5ccfb r6a5ccfb  
    9393    # Remember the component info blocks so we can build the model 
    9494    model_info.composition = ('product', [p_info, s_info]) 
     95    # TODO: delegate random to p_info, s_info 
     96    #model_info.random = lambda: {} 
    9597    model_info.demo = demo 
    9698 

  • sasmodels/sesans.py

    r94d13f1 r9f91afe  
    4141    _H0 = None # type: np.ndarray 
    4242 
    43     def __init__(self, z, SElength, zaccept, Rmax): 
     43    def __init__(self, z, SElength, lam, zaccept, Rmax): 
    4444        # type: (np.ndarray, float, float) -> None 
    4545        #import logging; logging.info("creating SESANS transform") 
    4646        self.q = z 
    47         self._set_hankel(SElength, zaccept, Rmax) 
     47        self._set_hankel(SElength, lam, zaccept, Rmax) 
    4848 
    4949    def apply(self, Iq): 
    … …  
    5454        return P 
    5555 
    56     def _set_hankel(self, SElength, zaccept, Rmax): 
     56    def _set_hankel(self, SElength, lam, zaccept, Rmax): 
    5757        # type: (np.ndarray, float, float) -> None 
    5858        # Force float32 arrays, otherwise run into memory problems on some machines 
    … …  
    7171        H = np.float32(dq/(2*pi)) * j0(repSE*repq) * repq 
    7272 
     73        replam = np.tile(lam, (q.size, 1)) 
     74        reptheta = np.arcsin(repq*replam/2*np.pi) 
     75        mask = reptheta > zaccept 
     76        H[mask] = 0 
     77 
    7378        self.q_calc = q 
    7479        self._H, self._H0 = H, H0 
Note: See TracChangeset for help on using the changeset viewer.