Changes in / [df87acf:275b07dc] in sasmodels


Ignore:
Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed

  • doc/guide/plugin.rst

    r2015f02 r57c609b  
    428428        def random(): 
    429429        ... 
    430          
    431 This function provides a model-specific random parameter set which shows model  
    432 features in the USANS to SANS range.  For example, core-shell sphere sets the  
    433 outer radius of the sphere logarithmically in `[20, 20,000]`, which sets the Q  
    434 value for the transition from flat to falling.  It then uses a beta distribution  
    435 to set the percentage of the shape which is shell, giving a preference for very  
    436 thin or very thick shells (but never 0% or 100%).  Using `-sets=10` in sascomp  
    437 should show a reasonable variety of curves over the default sascomp q range.   
    438 The parameter set is returned as a dictionary of `{parameter: value, ...}`.   
    439 Any model parameters not included in the dictionary will default according to  
     430 
     431This function provides a model-specific random parameter set which shows model 
     432features in the USANS to SANS range.  For example, core-shell sphere sets the 
     433outer radius of the sphere logarithmically in `[20, 20,000]`, which sets the Q 
     434value for the transition from flat to falling.  It then uses a beta distribution 
     435to set the percentage of the shape which is shell, giving a preference for very 
     436thin or very thick shells (but never 0% or 100%).  Using `-sets=10` in sascomp 
     437should show a reasonable variety of curves over the default sascomp q range. 
     438The parameter set is returned as a dictionary of `{parameter: value, ...}`. 
     439Any model parameters not included in the dictionary will default according to 
    440440the code in the `_randomize_one()` function from sasmodels/compare.py. 
    441441 
    … …  
    701701    erf, erfc, tgamma, lgamma:  **do not use** 
    702702        Special functions that should be part of the standard, but are missing 
    703         or inaccurate on some platforms. Use sas_erf, sas_erfc and sas_gamma 
    704         instead (see below). Note: lgamma(x) has not yet been tested. 
     703        or inaccurate on some platforms. Use sas_erf, sas_erfc, sas_gamma 
     704        and sas_lgamma instead (see below). 
    705705 
    706706Some non-standard constants and functions are also provided: 
    … …  
    769769        Gamma function sas_gamma\ $(x) = \Gamma(x)$. 
    770770 
    771         The standard math function, tgamma(x) is unstable for $x < 1$ 
     771        The standard math function, tgamma(x), is unstable for $x < 1$ 
    772772        on some platforms. 
    773773 
    774774        :code:`source = ["lib/sas_gamma.c", ...]` 
    775775        (`sas_gamma.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gamma.c>`_) 
     776 
     777    sas_gammaln(x): 
     778        log gamma function sas_gammaln\ $(x) = \log \Gamma(|x|)$. 
     779 
     780        The standard math function, lgamma(x), is incorrect for single 
     781        precision on some platforms. 
     782 
     783        :code:`source = ["lib/sas_gammainc.c", ...]` 
     784        (`sas_gammainc.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gammainc.c>`_) 
     785 
     786    sas_gammainc(a, x), sas_gammaincc(a, x): 
     787        Incomplete gamma function 
     788        sas_gammainc\ $(a, x) = \int_0^x t^{a-1}e^{-t}\,dt / \Gamma(a)$ 
     789        and complementary incomplete gamma function 
     790        sas_gammaincc\ $(a, x) = \int_x^\infty t^{a-1}e^{-t}\,dt / \Gamma(a)$ 
     791 
     792        :code:`source = ["lib/sas_gammainc.c", ...]` 
     793        (`sas_gammainc.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gammainc.c>`_) 
    776794 
    777795    sas_erf(x), sas_erfc(x): 
    … …  
    811829        If $n$ = 0 or 1, it uses sas_J0($x$) or sas_J1($x$), respectively. 
    812830 
     831        Warning: JN(n,x) can be very inaccurate (0.1%) for x not in [0.1, 100]. 
     832 
    813833        The standard math function jn(n, x) is not available on all platforms. 
    814834 
    … …  
    819839        Sine integral Si\ $(x) = \int_0^x \tfrac{\sin t}{t}\,dt$. 
    820840 
     841        Warning: Si(x) can be very inaccurate (0.1%) for x in [0.1, 100]. 
     842 
    821843        This function uses Taylor series for small and large arguments: 
    822844 
    823         For large arguments, 
     845        For large arguments use the following Taylor series, 
    824846 
    825847        .. math:: 
    … …  
    829851             - \frac{\sin(x)}{x}\left(\frac{1}{x} - \frac{3!}{x^3} + \frac{5!}{x^5} - \frac{7!}{x^7}\right) 
    830852 
    831         For small arguments, 
     853        For small arguments , 
    832854 
    833855        .. math:: 

  • explore/precision.py

    r2a7e20e rfba9ca0  
    9595            neg:    [-100,100] 
    9696 
     97        For arbitrary range use "start:stop:steps:scale" where scale is 
     98        one of log, lin, or linear. 
     99 
    97100        *diff* is "relative", "absolute" or "none" 
    98101 
    … …  
    102105        linear = not xrange.startswith("log") 
    103106        if xrange == "zoom": 
    104             lin_min, lin_max, lin_steps = 1000, 1010, 2000 
     107            start, stop, steps = 1000, 1010, 2000 
    105108        elif xrange == "neg": 
    106             lin_min, lin_max, lin_steps = -100.1, 100.1, 2000 
     109            start, stop, steps = -100.1, 100.1, 2000 
    107110        elif xrange == "linear": 
    108             lin_min, lin_max, lin_steps = 1, 1000, 2000 
    109             lin_min, lin_max, lin_steps = 0.001, 2, 2000 
     111            start, stop, steps = 1, 1000, 2000 
     112            start, stop, steps = 0.001, 2, 2000 
    110113        elif xrange == "log": 
    111             log_min, log_max, log_steps = -3, 5, 400 
     114            start, stop, steps = -3, 5, 400 
    112115        elif xrange == "logq": 
    113             log_min, log_max, log_steps = -4, 1, 400 
     116            start, stop, steps = -4, 1, 400 
     117        elif ':' in xrange: 
     118            parts = xrange.split(':') 
     119            linear = parts[3] != "log" if len(parts) == 4 else True 
     120            steps = int(parts[2]) if len(parts) > 2 else 400 
     121            start = float(parts[0]) 
     122            stop = float(parts[1]) 
     123 
    114124        else: 
    115125            raise ValueError("unknown range "+xrange) 
    … …  
    121131            # value to x in the given precision. 
    122132            if linear: 
    123                 lin_min = max(lin_min, self.limits[0]) 
    124                 lin_max = min(lin_max, self.limits[1]) 
    125                 qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='single') 
    126                 #qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='double') 
     133                start = max(start, self.limits[0]) 
     134                stop = min(stop, self.limits[1]) 
     135                qrf = np.linspace(start, stop, steps, dtype='single') 
     136                #qrf = np.linspace(start, stop, steps, dtype='double') 
    127137                qr = [mp.mpf(float(v)) for v in qrf] 
    128                 #qr = mp.linspace(lin_min, lin_max, lin_steps) 
     138                #qr = mp.linspace(start, stop, steps) 
    129139            else: 
    130                 log_min = np.log10(max(10**log_min, self.limits[0])) 
    131                 log_max = np.log10(min(10**log_max, self.limits[1])) 
    132                 qrf = np.logspace(log_min, log_max, log_steps, dtype='single') 
    133                 #qrf = np.logspace(log_min, log_max, log_steps, dtype='double') 
     140                start = np.log10(max(10**start, self.limits[0])) 
     141                stop = np.log10(min(10**stop, self.limits[1])) 
     142                qrf = np.logspace(start, stop, steps, dtype='single') 
     143                #qrf = np.logspace(start, stop, steps, dtype='double') 
    134144                qr = [mp.mpf(float(v)) for v in qrf] 
    135                 #qr = [10**v for v in mp.linspace(log_min, log_max, log_steps)] 
     145                #qr = [10**v for v in mp.linspace(start, stop, steps)] 
    136146 
    137147        target = self.call_mpmath(qr, bits=500) 
    … …  
    176186    """ 
    177187    if diff == "relative": 
    178         err = np.array([abs((t-a)/t) for t, a in zip(target, actual)], 'd') 
     188        err = np.array([(abs((t-a)/t) if t != 0 else a) for t, a in zip(target, actual)], 'd') 
    179189        #err = np.clip(err, 0, 1) 
    180190        pylab.loglog(x, err, '-', label=label) 
    … …  
    197207    return model_info 
    198208 
     209# Hack to allow second parameter A in two parameter functions 
     210A = 1 
     211def parse_extra_pars(): 
     212    global A 
     213 
     214    A_str = str(A) 
     215    pop = [] 
     216    for k, v in enumerate(sys.argv[1:]): 
     217        if v.startswith("A="): 
     218            A_str = v[2:] 
     219            pop.append(k+1) 
     220    if pop: 
     221        sys.argv = [v for k, v in enumerate(sys.argv) if k not in pop] 
     222        A = float(A_str) 
     223 
     224parse_extra_pars() 
     225 
    199226 
    200227# =============== FUNCTION DEFINITIONS ================ 
    … …  
    297324    ocl_function=make_ocl("return sas_gamma(q);", "sas_gamma", ["lib/sas_gamma.c"]), 
    298325    limits=(-3.1, 10), 
     326) 
     327add_function( 
     328    name="gammaln(x)", 
     329    mp_function=mp.loggamma, 
     330    np_function=scipy.special.gammaln, 
     331    ocl_function=make_ocl("return sas_gammaln(q);", "sas_gammaln", ["lib/sas_gammainc.c"]), 
     332    #ocl_function=make_ocl("return lgamma(q);", "sas_gammaln"), 
     333) 
     334add_function( 
     335    name="gammainc(x)", 
     336    mp_function=lambda x, a=A: mp.gammainc(a, a=0, b=x)/mp.gamma(a), 
     337    np_function=lambda x, a=A: scipy.special.gammainc(a, x), 
     338    ocl_function=make_ocl("return sas_gammainc(%.15g,q);"%A, "sas_gammainc", ["lib/sas_gammainc.c"]), 
     339) 
     340add_function( 
     341    name="gammaincc(x)", 
     342    mp_function=lambda x, a=A: mp.gammainc(a, a=x, b=mp.inf)/mp.gamma(a), 
     343    np_function=lambda x, a=A: scipy.special.gammaincc(a, x), 
     344    ocl_function=make_ocl("return sas_gammaincc(%.15g,q);"%A, "sas_gammaincc", ["lib/sas_gammainc.c"]), 
    299345) 
    300346add_function( 
    … …  
    463509lanczos_gamma = """\ 
    464510    const double coeff[] = { 
    465             76.18009172947146,     -86.50532032941677, 
    466             24.01409824083091,     -1.231739572450155, 
     511            76.18009172947146, -86.50532032941677, 
     512            24.01409824083091, -1.231739572450155, 
    467513            0.1208650973866179e-2,-0.5395239384953e-5 
    468514            }; 
    … …  
    475521""" 
    476522add_function( 
    477     name="log gamma(x)", 
     523    name="loggamma(x)", 
    478524    mp_function=mp.loggamma, 
    479525    np_function=scipy.special.gammaln, 
    … …  
    599645 
    600646ALL_FUNCTIONS = set(FUNCTIONS.keys()) 
    601 ALL_FUNCTIONS.discard("loggamma")  # OCL version not ready yet 
     647ALL_FUNCTIONS.discard("loggamma")  # use cephes-based gammaln instead 
    602648ALL_FUNCTIONS.discard("3j1/x:taylor") 
    603649ALL_FUNCTIONS.discard("3j1/x:trig") 
    … …  
    615661    -r indicates that the relative error should be plotted (default), 
    616662    -x<range> indicates the steps in x, where <range> is one of the following 
    617       log indicates log stepping in [10^-3, 10^5] (default) 
    618       logq indicates log stepping in [10^-4, 10^1] 
    619       linear indicates linear stepping in [1, 1000] 
    620       zoom indicates linear stepping in [1000, 1010] 
    621       neg indicates linear stepping in [-100.1, 100.1] 
    622 and name is "all" or one of: 
     663        log indicates log stepping in [10^-3, 10^5] (default) 
     664        logq indicates log stepping in [10^-4, 10^1] 
     665        linear indicates linear stepping in [1, 1000] 
     666        zoom indicates linear stepping in [1000, 1010] 
     667        neg indicates linear stepping in [-100.1, 100.1] 
     668        start:stop:n[:stepping] indicates an n-step plot in [start, stop] 
     669            or [10^start, 10^stop] if stepping is "log" (default n=400) 
     670Some functions (notably gammainc/gammaincc) have an additional parameter A 
     671which can be set from the command line as A=value.  Default is A=1. 
     672 
     673Name is one of: 
    623674    """+names) 
    624675    sys.exit(1) 

  • sasmodels/special.py

    rdf69efa rfba9ca0  
    113113        The standard math function, tgamma(x) is unstable for $x < 1$ 
    114114        on some platforms. 
     115 
     116    sas_gammaln(x): 
     117        log gamma function sas_gammaln\ $(x) = \log \Gamma(|x|)$. 
     118 
     119        The standard math function, lgamma(x), is incorrect for single 
     120        precision on some platforms. 
     121 
     122    sas_gammainc(a, x), sas_gammaincc(a, x): 
     123        Incomplete gamma function 
     124        sas_gammainc\ $(a, x) = \int_0^x t^{a-1}e^{-t}\,dt / \Gamma(a)$ 
     125        and complementary incomplete gamma function 
     126        sas_gammaincc\ $(a, x) = \int_x^\infty t^{a-1}e^{-t}\,dt / \Gamma(a)$ 
    115127 
    116128    sas_erf(x), sas_erfc(x): 
    … …  
    207219from numpy import pi, nan, inf 
    208220from scipy.special import gamma as sas_gamma 
     221from scipy.special import gammaln as sas_gammaln 
     222from scipy.special import gammainc as sas_gammainc 
     223from scipy.special import gammaincc as sas_gammaincc 
    209224from scipy.special import erf as sas_erf 
    210225from scipy.special import erfc as sas_erfc 
Note: See TracChangeset for help on using the changeset viewer.