Changes in / [aa8c6e0:153f8f6] in sasmodels


Ignore:
Files:
1 deleted
6 edited

Legend:

Unmodified
Added
Removed

  • doc/guide/plugin.rst

    r81751c2 r81751c2  
    744744    erf, erfc, tgamma, lgamma:  **do not use** 
    745745        Special functions that should be part of the standard, but are missing 
    746         or inaccurate on some platforms. Use sas_erf, sas_erfc, sas_gamma 
    747         and sas_lgamma instead (see below). 
     746        or inaccurate on some platforms. Use sas_erf, sas_erfc and sas_gamma 
     747        instead (see below). Note: lgamma(x) has not yet been tested. 
    748748 
    749749Some non-standard constants and functions are also provided: 
    … …  
    812812        Gamma function sas_gamma\ $(x) = \Gamma(x)$. 
    813813 
    814         The standard math function, tgamma(x), is unstable for $x < 1$ 
     814        The standard math function, tgamma(x) is unstable for $x < 1$ 
    815815        on some platforms. 
    816816 
    817817        :code:`source = ["lib/sas_gamma.c", ...]` 
    818818        (`sas_gamma.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gamma.c>`_) 
    819  
    820     sas_gammaln(x): 
    821         log gamma function sas_gammaln\ $(x) = \log \Gamma(|x|)$. 
    822  
    823         The standard math function, lgamma(x), is incorrect for single 
    824         precision on some platforms. 
    825  
    826         :code:`source = ["lib/sas_gammainc.c", ...]` 
    827         (`sas_gammainc.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gammainc.c>`_) 
    828  
    829     sas_gammainc(a, x), sas_gammaincc(a, x): 
    830         Incomplete gamma function 
    831         sas_gammainc\ $(a, x) = \int_0^x t^{a-1}e^{-t}\,dt / \Gamma(a)$ 
    832         and complementary incomplete gamma function 
    833         sas_gammaincc\ $(a, x) = \int_x^\infty t^{a-1}e^{-t}\,dt / \Gamma(a)$ 
    834  
    835         :code:`source = ["lib/sas_gammainc.c", ...]` 
    836         (`sas_gammainc.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gammainc.c>`_) 
    837819 
    838820    sas_erf(x), sas_erfc(x): 
    … …  
    872854        If $n$ = 0 or 1, it uses sas_J0($x$) or sas_J1($x$), respectively. 
    873855 
    874         Warning: JN(n,x) can be very inaccurate (0.1%) for x not in [0.1, 100]. 
    875  
    876856        The standard math function jn(n, x) is not available on all platforms. 
    877857 
    … …  
    882862        Sine integral Si\ $(x) = \int_0^x \tfrac{\sin t}{t}\,dt$. 
    883863 
    884         Warning: Si(x) can be very inaccurate (0.1%) for x in [0.1, 100]. 
    885  
    886864        This function uses Taylor series for small and large arguments: 
    887865 
    888         For large arguments use the following Taylor series, 
     866        For large arguments, 
    889867 
    890868        .. math:: 
    … …  
    894872             - \frac{\sin(x)}{x}\left(\frac{1}{x} - \frac{3!}{x^3} + \frac{5!}{x^5} - \frac{7!}{x^7}\right) 
    895873 
    896         For small arguments , 
     874        For small arguments, 
    897875 
    898876        .. math:: 

  • explore/precision.py

    rfba9ca0 ree60aa7  
    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  
    10097        *diff* is "relative", "absolute" or "none" 
    10198 
    … …  
    105102        linear = not xrange.startswith("log") 
    106103        if xrange == "zoom": 
    107             start, stop, steps = 1000, 1010, 2000 
     104            lin_min, lin_max, lin_steps = 1000, 1010, 2000 
    108105        elif xrange == "neg": 
    109             start, stop, steps = -100.1, 100.1, 2000 
     106            lin_min, lin_max, lin_steps = -100.1, 100.1, 2000 
    110107        elif xrange == "linear": 
    111             start, stop, steps = 1, 1000, 2000 
    112             start, stop, steps = 0.001, 2, 2000 
     108            lin_min, lin_max, lin_steps = 1, 1000, 2000 
     109            lin_min, lin_max, lin_steps = 0.001, 2, 2000 
    113110        elif xrange == "log": 
    114             start, stop, steps = -3, 5, 400 
     111            log_min, log_max, log_steps = -3, 5, 400 
    115112        elif xrange == "logq": 
    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  
     113            log_min, log_max, log_steps = -4, 1, 400 
    124114        else: 
    125115            raise ValueError("unknown range "+xrange) 
    … …  
    131121            # value to x in the given precision. 
    132122            if linear: 
    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') 
     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') 
    137127                qr = [mp.mpf(float(v)) for v in qrf] 
    138                 #qr = mp.linspace(start, stop, steps) 
     128                #qr = mp.linspace(lin_min, lin_max, lin_steps) 
    139129            else: 
    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') 
     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') 
    144134                qr = [mp.mpf(float(v)) for v in qrf] 
    145                 #qr = [10**v for v in mp.linspace(start, stop, steps)] 
     135                #qr = [10**v for v in mp.linspace(log_min, log_max, log_steps)] 
    146136 
    147137        target = self.call_mpmath(qr, bits=500) 
    … …  
    186176    """ 
    187177    if diff == "relative": 
    188         err = np.array([(abs((t-a)/t) if t != 0 else a) for t, a in zip(target, actual)], 'd') 
     178        err = np.array([abs((t-a)/t) for t, a in zip(target, actual)], 'd') 
    189179        #err = np.clip(err, 0, 1) 
    190180        pylab.loglog(x, err, '-', label=label) 
    … …  
    207197    return model_info 
    208198 
    209 # Hack to allow second parameter A in two parameter functions 
    210 A = 1 
    211 def 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  
    224 parse_extra_pars() 
    225  
    226199 
    227200# =============== FUNCTION DEFINITIONS ================ 
    … …  
    324297    ocl_function=make_ocl("return sas_gamma(q);", "sas_gamma", ["lib/sas_gamma.c"]), 
    325298    limits=(-3.1, 10), 
    326 ) 
    327 add_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 ) 
    334 add_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 ) 
    340 add_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"]), 
    345299) 
    346300add_function( 
    … …  
    511465lanczos_gamma = """\ 
    512466    const double coeff[] = { 
    513             76.18009172947146, -86.50532032941677, 
    514             24.01409824083091, -1.231739572450155, 
     467            76.18009172947146,     -86.50532032941677, 
     468            24.01409824083091,     -1.231739572450155, 
    515469            0.1208650973866179e-2,-0.5395239384953e-5 
    516470            }; 
    … …  
    523477""" 
    524478add_function( 
    525     name="loggamma(x)", 
     479    name="log gamma(x)", 
    526480    mp_function=mp.loggamma, 
    527481    np_function=scipy.special.gammaln, 
    … …  
    647601 
    648602ALL_FUNCTIONS = set(FUNCTIONS.keys()) 
    649 ALL_FUNCTIONS.discard("loggamma")  # use cephes-based gammaln instead 
     603ALL_FUNCTIONS.discard("loggamma")  # OCL version not ready yet 
    650604ALL_FUNCTIONS.discard("3j1/x:taylor") 
    651605ALL_FUNCTIONS.discard("3j1/x:trig") 
    … …  
    663617    -r indicates that the relative error should be plotted (default), 
    664618    -x<range> indicates the steps in x, where <range> is one of the following 
    665         log indicates log stepping in [10^-3, 10^5] (default) 
    666         logq indicates log stepping in [10^-4, 10^1] 
    667         linear indicates linear stepping in [1, 1000] 
    668         zoom indicates linear stepping in [1000, 1010] 
    669         neg indicates linear stepping in [-100.1, 100.1] 
    670         start:stop:n[:stepping] indicates an n-step plot in [start, stop] 
    671             or [10^start, 10^stop] if stepping is "log" (default n=400) 
    672 Some functions (notably gammainc/gammaincc) have an additional parameter A 
    673 which can be set from the command line as A=value.  Default is A=1. 
    674  
    675 Name is one of: 
     619      log indicates log stepping in [10^-3, 10^5] (default) 
     620      logq indicates log stepping in [10^-4, 10^1] 
     621      linear indicates linear stepping in [1, 1000] 
     622      zoom indicates linear stepping in [1000, 1010] 
     623      neg indicates linear stepping in [-100.1, 100.1] 
     624and name is "all" or one of: 
    676625    """+names) 
    677626    sys.exit(1) 

  • sasmodels/kernelpy.py

    re44432d re44432d  
    4242        self.info = model_info 
    4343        self.dtype = np.dtype('d') 
    44         logger.info("make python model " + self.info.name) 
    4544 
    4645    def make_kernel(self, q_vectors): 

  • sasmodels/model_test.py

    r81751c2 r81751c2  
    4747import sys 
    4848import unittest 
    49 import traceback 
    5049 
    5150try: 
    … …  
    7776# pylint: enable=unused-import 
    7877 
     78 
    7979def make_suite(loaders, models): 
    8080    # type: (List[str], List[str]) -> unittest.TestSuite 
    … …  
    8787    *models* is the list of models to test, or *["all"]* to test all models. 
    8888    """ 
     89    ModelTestCase = _hide_model_case_from_nose() 
    8990    suite = unittest.TestSuite() 
    9091 
    … …  
    9596        skip = [] 
    9697    for model_name in models: 
    97         if model_name not in skip: 
    98             model_info = load_model_info(model_name) 
    99             _add_model_to_suite(loaders, suite, model_info) 
     98        if model_name in skip: 
     99            continue 
     100        model_info = load_model_info(model_name) 
     101 
     102        #print('------') 
     103        #print('found tests in', model_name) 
     104        #print('------') 
     105 
     106        # if ispy then use the dll loader to call pykernel 
     107        # don't try to call cl kernel since it will not be 
     108        # available in some environmentes. 
     109        is_py = callable(model_info.Iq) 
     110 
     111        # Some OpenCL drivers seem to be flaky, and are not producing the 
     112        # expected result.  Since we don't have known test values yet for 
     113        # all of our models, we are instead going to compare the results 
     114        # for the 'smoke test' (that is, evaluation at q=0.1 for the default 
     115        # parameters just to see that the model runs to completion) between 
     116        # the OpenCL and the DLL.  To do this, we define a 'stash' which is 
     117        # shared between OpenCL and DLL tests.  This is just a list.  If the 
     118        # list is empty (which it will be when DLL runs, if the DLL runs 
     119        # first), then the results are appended to the list.  If the list 
     120        # is not empty (which it will be when OpenCL runs second), the results 
     121        # are compared to the results stored in the first element of the list. 
     122        # This is a horrible stateful hack which only makes sense because the 
     123        # test suite is thrown away after being run once. 
     124        stash = [] 
     125 
     126        if is_py:  # kernel implemented in python 
     127            test_name = "%s-python"%model_name 
     128            test_method_name = "test_%s_python" % model_info.id 
     129            test = ModelTestCase(test_name, model_info, 
     130                                 test_method_name, 
     131                                 platform="dll",  # so that 
     132                                 dtype="double", 
     133                                 stash=stash) 
     134            suite.addTest(test) 
     135        else:   # kernel implemented in C 
     136 
     137            # test using dll if desired 
     138            if 'dll' in loaders: 
     139                test_name = "%s-dll"%model_name 
     140                test_method_name = "test_%s_dll" % model_info.id 
     141                test = ModelTestCase(test_name, model_info, 
     142                                     test_method_name, 
     143                                     platform="dll", 
     144                                     dtype="double", 
     145                                     stash=stash) 
     146                suite.addTest(test) 
     147 
     148            # test using opencl if desired and available 
     149            if 'opencl' in loaders and use_opencl(): 
     150                test_name = "%s-opencl"%model_name 
     151                test_method_name = "test_%s_opencl" % model_info.id 
     152                # Using dtype=None so that the models that are only 
     153                # correct for double precision are not tested using 
     154                # single precision.  The choice is determined by the 
     155                # presence of *single=False* in the model file. 
     156                test = ModelTestCase(test_name, model_info, 
     157                                     test_method_name, 
     158                                     platform="ocl", dtype=None, 
     159                                     stash=stash) 
     160                #print("defining", test_name) 
     161                suite.addTest(test) 
     162 
     163            # test using cuda if desired and available 
     164            if 'cuda' in loaders and use_cuda(): 
     165                test_name = "%s-cuda"%model_name 
     166                test_method_name = "test_%s_cuda" % model_info.id 
     167                # Using dtype=None so that the models that are only 
     168                # correct for double precision are not tested using 
     169                # single precision.  The choice is determined by the 
     170                # presence of *single=False* in the model file. 
     171                test = ModelTestCase(test_name, model_info, 
     172                                     test_method_name, 
     173                                     platform="cuda", dtype=None, 
     174                                     stash=stash) 
     175                #print("defining", test_name) 
     176                suite.addTest(test) 
    100177 
    101178    return suite 
    102  
    103 def _add_model_to_suite(loaders, suite, model_info): 
    104     ModelTestCase = _hide_model_case_from_nose() 
    105  
    106     #print('------') 
    107     #print('found tests in', model_name) 
    108     #print('------') 
    109  
    110     # if ispy then use the dll loader to call pykernel 
    111     # don't try to call cl kernel since it will not be 
    112     # available in some environmentes. 
    113     is_py = callable(model_info.Iq) 
    114  
    115     # Some OpenCL drivers seem to be flaky, and are not producing the 
    116     # expected result.  Since we don't have known test values yet for 
    117     # all of our models, we are instead going to compare the results 
    118     # for the 'smoke test' (that is, evaluation at q=0.1 for the default 
    119     # parameters just to see that the model runs to completion) between 
    120     # the OpenCL and the DLL.  To do this, we define a 'stash' which is 
    121     # shared between OpenCL and DLL tests.  This is just a list.  If the 
    122     # list is empty (which it will be when DLL runs, if the DLL runs 
    123     # first), then the results are appended to the list.  If the list 
    124     # is not empty (which it will be when OpenCL runs second), the results 
    125     # are compared to the results stored in the first element of the list. 
    126     # This is a horrible stateful hack which only makes sense because the 
    127     # test suite is thrown away after being run once. 
    128     stash = [] 
    129  
    130     if is_py:  # kernel implemented in python 
    131         test_name = "%s-python"%model_info.name 
    132         test_method_name = "test_%s_python" % model_info.id 
    133         test = ModelTestCase(test_name, model_info, 
    134                                 test_method_name, 
    135                                 platform="dll",  # so that 
    136                                 dtype="double", 
    137                                 stash=stash) 
    138         suite.addTest(test) 
    139     else:   # kernel implemented in C 
    140  
    141         # test using dll if desired 
    142         if 'dll' in loaders or not use_opencl(): 
    143             test_name = "%s-dll"%model_info.name 
    144             test_method_name = "test_%s_dll" % model_info.id 
    145             test = ModelTestCase(test_name, model_info, 
    146                                     test_method_name, 
    147                                     platform="dll", 
    148                                     dtype="double", 
    149                                     stash=stash) 
    150             suite.addTest(test) 
    151  
    152         # test using opencl if desired and available 
    153         if 'opencl' in loaders and use_opencl(): 
    154             test_name = "%s-opencl"%model_info.name 
    155             test_method_name = "test_%s_opencl" % model_info.id 
    156             # Using dtype=None so that the models that are only 
    157             # correct for double precision are not tested using 
    158             # single precision.  The choice is determined by the 
    159             # presence of *single=False* in the model file. 
    160             test = ModelTestCase(test_name, model_info, 
    161                                     test_method_name, 
    162                                     platform="ocl", dtype=None, 
    163                                     stash=stash) 
    164             #print("defining", test_name) 
    165             suite.addTest(test) 
    166  
    167         # test using cuda if desired and available 
    168         if 'cuda' in loaders and use_cuda(): 
    169             test_name = "%s-cuda"%model_name 
    170             test_method_name = "test_%s_cuda" % model_info.id 
    171             # Using dtype=None so that the models that are only 
    172             # correct for double precision are not tested using 
    173             # single precision.  The choice is determined by the 
    174             # presence of *single=False* in the model file. 
    175             test = ModelTestCase(test_name, model_info, 
    176                                     test_method_name, 
    177                                     platform="cuda", dtype=None, 
    178                                     stash=stash) 
    179             #print("defining", test_name) 
    180             suite.addTest(test) 
    181  
    182179 
    183180def _hide_model_case_from_nose(): 
    … …  
    406403    return abs(target-actual)/shift < 1.5*10**-digits 
    407404 
    408 # CRUFT: old interface; should be deprecated and removed 
    409 def run_one(model_name): 
    410     # msg = "use check_model(model_info) rather than run_one(model_name)" 
    411     # warnings.warn(msg, category=DeprecationWarning, stacklevel=2) 
    412     try: 
    413         model_info = load_model_info(model_name) 
    414     except Exception: 
    415         output = traceback.format_exc() 
    416         return output 
    417  
    418     success, output = check_model(model_info) 
    419     return output 
    420  
    421 def check_model(model_info): 
    422     # type: (ModelInfo) -> str 
    423     """ 
    424     Run the tests for a single model, capturing the output. 
    425  
    426     Returns success status and the output string. 
     405def run_one(model): 
     406    # type: (str) -> str 
     407    """ 
     408    Run the tests for a single model, printing the results to stdout. 
     409 
     410    *model* can by a python file, which is handy for checking user defined 
     411    plugin models. 
    427412    """ 
    428413    # Note that running main() directly did not work from within the 
    … …  
    438423 
    439424    # Build a test suite containing just the model 
    440     loaders = ['opencl' if use_opencl() else 'cuda' if use_cuda() else 'dll'] 
    441     suite = unittest.TestSuite() 
    442     _add_model_to_suite(loaders, suite, model_info) 
     425    loader = 'opencl' if use_opencl() else 'cuda' if use_cuda() else 'dll' 
     426    models = [model] 
     427    try: 
     428        suite = make_suite([loader], models) 
     429    except Exception: 
     430        import traceback 
     431        stream.writeln(traceback.format_exc()) 
     432        return 
    443433 
    444434    # Warn if there are no user defined tests. 
    … …  
    455445    for test in suite: 
    456446        if not test.info.tests: 
    457             stream.writeln("Note: %s has no user defined tests."%model_info.name) 
     447            stream.writeln("Note: %s has no user defined tests."%model) 
    458448        break 
    459449    else: 
    … …  
    471461    output = stream.getvalue() 
    472462    stream.close() 
    473     return result.wasSuccessful(), output 
     463    return output 
    474464 
    475465 

  • sasmodels/sasview_model.py

    r39a06c9 r39a06c9  
    819819            return value, [value], [1.0] 
    820820 
    821     @classmethod 
    822     def runTests(cls): 
    823         """ 
    824         Run any tests built into the model and captures the test output. 
    825  
    826         Returns success flag and output 
    827         """ 
    828         from .model_test import check_model 
    829         return check_model(cls._model_info) 
    830  
    831821def test_cylinder(): 
    832822    # type: () -> float 

  • sasmodels/special.py

    rfba9ca0 rdf69efa  
    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)$ 
    127115 
    128116    sas_erf(x), sas_erfc(x): 
    … …  
    219207from numpy import pi, nan, inf 
    220208from scipy.special import gamma as sas_gamma 
    221 from scipy.special import gammaln as sas_gammaln 
    222 from scipy.special import gammainc as sas_gammainc 
    223 from scipy.special import gammaincc as sas_gammaincc 
    224209from scipy.special import erf as sas_erf 
    225210from scipy.special import erfc as sas_erfc 
Note: See TracChangeset for help on using the changeset viewer.