Changes in / [153f8f6:aa8c6e0] in sasmodels

Files:

• doc/guide/plugin.rst (modified) (5 diffs)
• explore/precision.py (modified) (10 diffs)
• sasmodels/kernelpy.py (modified) (1 diff)
• sasmodels/model_test.py (modified) (8 diffs)
• sasmodels/models/lib/sas_gammainc.c (added)
• sasmodels/sasview_model.py (modified) (1 diff)
• sasmodels/special.py (modified) (2 diffs)

doc/guide/plugin.rst

@@ -744 +744 @@
     erf, erfc, tgamma, lgamma:  **do not use**
         Special functions that should be part of the standard, but are missing
-        or inaccurate on some platforms. Use sas_erf, sas_erfc and sas_gamma
-        instead (see below). Note: lgamma(x) has not yet been tested.
+        or inaccurate on some platforms. Use sas_erf, sas_erfc, sas_gamma
+        and sas_lgamma instead (see below).
 
 Some non-standard constants and functions are also provided:
@@ -812 +812 @@
         Gamma function sas_gamma\ $(x) = \Gamma(x)$.
 
-        The standard math function, tgamma(x) is unstable for $x < 1$
+        The standard math function, tgamma(x), is unstable for $x < 1$
         on some platforms.
 
         :code:`source = ["lib/sas_gamma.c", ...]`
         (`sas_gamma.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gamma.c>`_)
+
+    sas_gammaln(x):
+        log gamma function sas_gammaln\ $(x) = \log \Gamma(|x|)$.
+
+        The standard math function, lgamma(x), is incorrect for single
+        precision on some platforms.
+
+        :code:`source = ["lib/sas_gammainc.c", ...]`
+        (`sas_gammainc.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gammainc.c>`_)
+
+    sas_gammainc(a, x), sas_gammaincc(a, x):
+        Incomplete gamma function
+        sas_gammainc\ $(a, x) = \int_0^x t^{a-1}e^{-t}\,dt / \Gamma(a)$
+        and complementary incomplete gamma function
+        sas_gammaincc\ $(a, x) = \int_x^\infty t^{a-1}e^{-t}\,dt / \Gamma(a)$
+
+        :code:`source = ["lib/sas_gammainc.c", ...]`
+        (`sas_gammainc.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gammainc.c>`_)
 
     sas_erf(x), sas_erfc(x):
@@ -854 +872 @@
         If $n$ = 0 or 1, it uses sas_J0($x$) or sas_J1($x$), respectively.
 
+        Warning: JN(n,x) can be very inaccurate (0.1%) for x not in [0.1, 100].
+
         The standard math function jn(n, x) is not available on all platforms.
 
@@ -862 +882 @@
         Sine integral Si\ $(x) = \int_0^x \tfrac{\sin t}{t}\,dt$.
 
+        Warning: Si(x) can be very inaccurate (0.1%) for x in [0.1, 100].
+
         This function uses Taylor series for small and large arguments:
 
-        For large arguments,
+        For large arguments use the following Taylor series,
 
         .. math::
@@ -872 +894 @@
              - \frac{\sin(x)}{x}\left(\frac{1}{x} - \frac{3!}{x^3} + \frac{5!}{x^5} - \frac{7!}{x^7}\right)
 
-        For small arguments,
+        For small arguments ,
 
         .. math::

explore/precision.py

@@ -95 +95 @@
             neg:    [-100,100]
 
+        For arbitrary range use "start:stop:steps:scale" where scale is
+        one of log, lin, or linear.
+
         *diff* is "relative", "absolute" or "none"
 
@@ -102 +105 @@
         linear = not xrange.startswith("log")
         if xrange == "zoom":
-            lin_min, lin_max, lin_steps = 1000, 1010, 2000
+            start, stop, steps = 1000, 1010, 2000
         elif xrange == "neg":
-            lin_min, lin_max, lin_steps = -100.1, 100.1, 2000
+            start, stop, steps = -100.1, 100.1, 2000
         elif xrange == "linear":
-            lin_min, lin_max, lin_steps = 1, 1000, 2000
-            lin_min, lin_max, lin_steps = 0.001, 2, 2000
+            start, stop, steps = 1, 1000, 2000
+            start, stop, steps = 0.001, 2, 2000
         elif xrange == "log":
-            log_min, log_max, log_steps = -3, 5, 400
+            start, stop, steps = -3, 5, 400
         elif xrange == "logq":
-            log_min, log_max, log_steps = -4, 1, 400
+            start, stop, steps = -4, 1, 400
+        elif ':' in xrange:
+            parts = xrange.split(':')
+            linear = parts[3] != "log" if len(parts) == 4 else True
+            steps = int(parts[2]) if len(parts) > 2 else 400
+            start = float(parts[0])
+            stop = float(parts[1])
+
         else:
             raise ValueError("unknown range "+xrange)
@@ -121 +131 @@
             # value to x in the given precision.
             if linear:
-                lin_min = max(lin_min, self.limits[0])
-                lin_max = min(lin_max, self.limits[1])
-                qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='single')
-                #qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='double')
+                start = max(start, self.limits[0])
+                stop = min(stop, self.limits[1])
+                qrf = np.linspace(start, stop, steps, dtype='single')
+                #qrf = np.linspace(start, stop, steps, dtype='double')
                 qr = [mp.mpf(float(v)) for v in qrf]
-                #qr = mp.linspace(lin_min, lin_max, lin_steps)
+                #qr = mp.linspace(start, stop, steps)
             else:
-                log_min = np.log10(max(10**log_min, self.limits[0]))
-                log_max = np.log10(min(10**log_max, self.limits[1]))
-                qrf = np.logspace(log_min, log_max, log_steps, dtype='single')
-                #qrf = np.logspace(log_min, log_max, log_steps, dtype='double')
+                start = np.log10(max(10**start, self.limits[0]))
+                stop = np.log10(min(10**stop, self.limits[1]))
+                qrf = np.logspace(start, stop, steps, dtype='single')
+                #qrf = np.logspace(start, stop, steps, dtype='double')
                 qr = [mp.mpf(float(v)) for v in qrf]
-                #qr = [10**v for v in mp.linspace(log_min, log_max, log_steps)]
+                #qr = [10**v for v in mp.linspace(start, stop, steps)]
 
         target = self.call_mpmath(qr, bits=500)
@@ -176 +186 @@
     """
     if diff == "relative":
-        err = np.array([abs((t-a)/t) for t, a in zip(target, actual)], 'd')
+        err = np.array([(abs((t-a)/t) if t != 0 else a) for t, a in zip(target, actual)], 'd')
         #err = np.clip(err, 0, 1)
         pylab.loglog(x, err, '-', label=label)
@@ -197 +207 @@
     return model_info
 
+# Hack to allow second parameter A in two parameter functions
+A = 1
+def parse_extra_pars():
+    global A
+
+    A_str = str(A)
+    pop = []
+    for k, v in enumerate(sys.argv[1:]):
+        if v.startswith("A="):
+            A_str = v[2:]
+            pop.append(k+1)
+    if pop:
+        sys.argv = [v for k, v in enumerate(sys.argv) if k not in pop]
+        A = float(A_str)
+
+parse_extra_pars()
+
 
 # =============== FUNCTION DEFINITIONS ================
@@ -297 +324 @@
     ocl_function=make_ocl("return sas_gamma(q);", "sas_gamma", ["lib/sas_gamma.c"]),
     limits=(-3.1, 10),
+)
+add_function(
+    name="gammaln(x)",
+    mp_function=mp.loggamma,
+    np_function=scipy.special.gammaln,
+    ocl_function=make_ocl("return sas_gammaln(q);", "sas_gammaln", ["lib/sas_gammainc.c"]),
+    #ocl_function=make_ocl("return lgamma(q);", "sas_gammaln"),
+)
+add_function(
+    name="gammainc(x)",
+    mp_function=lambda x, a=A: mp.gammainc(a, a=0, b=x)/mp.gamma(a),
+    np_function=lambda x, a=A: scipy.special.gammainc(a, x),
+    ocl_function=make_ocl("return sas_gammainc(%.15g,q);"%A, "sas_gammainc", ["lib/sas_gammainc.c"]),
+)
+add_function(
+    name="gammaincc(x)",
+    mp_function=lambda x, a=A: mp.gammainc(a, a=x, b=mp.inf)/mp.gamma(a),
+    np_function=lambda x, a=A: scipy.special.gammaincc(a, x),
+    ocl_function=make_ocl("return sas_gammaincc(%.15g,q);"%A, "sas_gammaincc", ["lib/sas_gammainc.c"]),
 )
 add_function(
@@ -465 +511 @@
 lanczos_gamma = """\
     const double coeff[] = {
-            76.18009172947146,     -86.50532032941677,
-            24.01409824083091,     -1.231739572450155,
+            76.18009172947146, -86.50532032941677,
+            24.01409824083091, -1.231739572450155,
             0.1208650973866179e-2,-0.5395239384953e-5
             };
@@ -477 +523 @@
 """
 add_function(
-    name="log gamma(x)",
+    name="loggamma(x)",
     mp_function=mp.loggamma,
     np_function=scipy.special.gammaln,
@@ -601 +647 @@
 
 ALL_FUNCTIONS = set(FUNCTIONS.keys())
-ALL_FUNCTIONS.discard("loggamma")  # OCL version not ready yet
+ALL_FUNCTIONS.discard("loggamma")  # use cephes-based gammaln instead
 ALL_FUNCTIONS.discard("3j1/x:taylor")
 ALL_FUNCTIONS.discard("3j1/x:trig")
@@ -617 +663 @@
     -r indicates that the relative error should be plotted (default),
     -x<range> indicates the steps in x, where <range> is one of the following
-      log indicates log stepping in [10^-3, 10^5] (default)
-      logq indicates log stepping in [10^-4, 10^1]
-      linear indicates linear stepping in [1, 1000]
-      zoom indicates linear stepping in [1000, 1010]
-      neg indicates linear stepping in [-100.1, 100.1]
-and name is "all" or one of:
+        log indicates log stepping in [10^-3, 10^5] (default)
+        logq indicates log stepping in [10^-4, 10^1]
+        linear indicates linear stepping in [1, 1000]
+        zoom indicates linear stepping in [1000, 1010]
+        neg indicates linear stepping in [-100.1, 100.1]
+        start:stop:n[:stepping] indicates an n-step plot in [start, stop]
+            or [10^start, 10^stop] if stepping is "log" (default n=400)
+Some functions (notably gammainc/gammaincc) have an additional parameter A
+which can be set from the command line as A=value.  Default is A=1.
+
+Name is one of:
     """+names)
     sys.exit(1)

sasmodels/kernelpy.py

@@ -42 +42 @@
         self.info = model_info
         self.dtype = np.dtype('d')
+        logger.info("make python model " + self.info.name)
 
     def make_kernel(self, q_vectors):

sasmodels/model_test.py

@@ -47 +47 @@
 import sys
 import unittest
+import traceback
 
 try:
@@ -76 +77 @@
 # pylint: enable=unused-import
 
-
 def make_suite(loaders, models):
     # type: (List[str], List[str]) -> unittest.TestSuite
@@ -87 +87 @@
     *models* is the list of models to test, or *["all"]* to test all models.
     """
-    ModelTestCase = _hide_model_case_from_nose()
     suite = unittest.TestSuite()
 
@@ -96 +95 @@
         skip = []
     for model_name in models:
-        if model_name in skip:
-            continue
-        model_info = load_model_info(model_name)
-
-        #print('------')
-        #print('found tests in', model_name)
-        #print('------')
-
-        # if ispy then use the dll loader to call pykernel
-        # don't try to call cl kernel since it will not be
-        # available in some environmentes.
-        is_py = callable(model_info.Iq)
-
-        # Some OpenCL drivers seem to be flaky, and are not producing the
-        # expected result.  Since we don't have known test values yet for
-        # all of our models, we are instead going to compare the results
-        # for the 'smoke test' (that is, evaluation at q=0.1 for the default
-        # parameters just to see that the model runs to completion) between
-        # the OpenCL and the DLL.  To do this, we define a 'stash' which is
-        # shared between OpenCL and DLL tests.  This is just a list.  If the
-        # list is empty (which it will be when DLL runs, if the DLL runs
-        # first), then the results are appended to the list.  If the list
-        # is not empty (which it will be when OpenCL runs second), the results
-        # are compared to the results stored in the first element of the list.
-        # This is a horrible stateful hack which only makes sense because the
-        # test suite is thrown away after being run once.
-        stash = []
-
-        if is_py:  # kernel implemented in python
-            test_name = "%s-python"%model_name
-            test_method_name = "test_%s_python" % model_info.id
+        if model_name not in skip:
+            model_info = load_model_info(model_name)
+            _add_model_to_suite(loaders, suite, model_info)
+
+    return suite
+
+def _add_model_to_suite(loaders, suite, model_info):
+    ModelTestCase = _hide_model_case_from_nose()
+
+    #print('------')
+    #print('found tests in', model_name)
+    #print('------')
+
+    # if ispy then use the dll loader to call pykernel
+    # don't try to call cl kernel since it will not be
+    # available in some environmentes.
+    is_py = callable(model_info.Iq)
+
+    # Some OpenCL drivers seem to be flaky, and are not producing the
+    # expected result.  Since we don't have known test values yet for
+    # all of our models, we are instead going to compare the results
+    # for the 'smoke test' (that is, evaluation at q=0.1 for the default
+    # parameters just to see that the model runs to completion) between
+    # the OpenCL and the DLL.  To do this, we define a 'stash' which is
+    # shared between OpenCL and DLL tests.  This is just a list.  If the
+    # list is empty (which it will be when DLL runs, if the DLL runs
+    # first), then the results are appended to the list.  If the list
+    # is not empty (which it will be when OpenCL runs second), the results
+    # are compared to the results stored in the first element of the list.
+    # This is a horrible stateful hack which only makes sense because the
+    # test suite is thrown away after being run once.
+    stash = []
+
+    if is_py:  # kernel implemented in python
+        test_name = "%s-python"%model_info.name
+        test_method_name = "test_%s_python" % model_info.id
+        test = ModelTestCase(test_name, model_info,
+                                test_method_name,
+                                platform="dll",  # so that
+                                dtype="double",
+                                stash=stash)
+        suite.addTest(test)
+    else:   # kernel implemented in C
+
+        # test using dll if desired
+        if 'dll' in loaders or not use_opencl():
+            test_name = "%s-dll"%model_info.name
+            test_method_name = "test_%s_dll" % model_info.id
             test = ModelTestCase(test_name, model_info,
-                                 test_method_name,
-                                 platform="dll",  # so that
-                                 dtype="double",
-                                 stash=stash)
+                                    test_method_name,
+                                    platform="dll",
+                                    dtype="double",
+                                    stash=stash)
             suite.addTest(test)
-        else:   # kernel implemented in C
-
-            # test using dll if desired
-            if 'dll' in loaders:
-                test_name = "%s-dll"%model_name
-                test_method_name = "test_%s_dll" % model_info.id
-                test = ModelTestCase(test_name, model_info,
-                                     test_method_name,
-                                     platform="dll",
-                                     dtype="double",
-                                     stash=stash)
-                suite.addTest(test)
-
-            # test using opencl if desired and available
-            if 'opencl' in loaders and use_opencl():
-                test_name = "%s-opencl"%model_name
-                test_method_name = "test_%s_opencl" % model_info.id
-                # Using dtype=None so that the models that are only
-                # correct for double precision are not tested using
-                # single precision.  The choice is determined by the
-                # presence of *single=False* in the model file.
-                test = ModelTestCase(test_name, model_info,
-                                     test_method_name,
-                                     platform="ocl", dtype=None,
-                                     stash=stash)
-                #print("defining", test_name)
-                suite.addTest(test)
-
-            # test using cuda if desired and available
-            if 'cuda' in loaders and use_cuda():
-                test_name = "%s-cuda"%model_name
-                test_method_name = "test_%s_cuda" % model_info.id
-                # Using dtype=None so that the models that are only
-                # correct for double precision are not tested using
-                # single precision.  The choice is determined by the
-                # presence of *single=False* in the model file.
-                test = ModelTestCase(test_name, model_info,
-                                     test_method_name,
-                                     platform="cuda", dtype=None,
-                                     stash=stash)
-                #print("defining", test_name)
-                suite.addTest(test)
-
-    return suite
+
+        # test using opencl if desired and available
+        if 'opencl' in loaders and use_opencl():
+            test_name = "%s-opencl"%model_info.name
+            test_method_name = "test_%s_opencl" % model_info.id
+            # Using dtype=None so that the models that are only
+            # correct for double precision are not tested using
+            # single precision.  The choice is determined by the
+            # presence of *single=False* in the model file.
+            test = ModelTestCase(test_name, model_info,
+                                    test_method_name,
+                                    platform="ocl", dtype=None,
+                                    stash=stash)
+            #print("defining", test_name)
+            suite.addTest(test)
+
+        # test using cuda if desired and available
+        if 'cuda' in loaders and use_cuda():
+            test_name = "%s-cuda"%model_name
+            test_method_name = "test_%s_cuda" % model_info.id
+            # Using dtype=None so that the models that are only
+            # correct for double precision are not tested using
+            # single precision.  The choice is determined by the
+            # presence of *single=False* in the model file.
+            test = ModelTestCase(test_name, model_info,
+                                    test_method_name,
+                                    platform="cuda", dtype=None,
+                                    stash=stash)
+            #print("defining", test_name)
+            suite.addTest(test)
+
 
 def _hide_model_case_from_nose():
@@ -403 +406 @@
     return abs(target-actual)/shift < 1.5*10**-digits
 
-def run_one(model):
-    # type: (str) -> str
-    """
-    Run the tests for a single model, printing the results to stdout.
-
-    *model* can by a python file, which is handy for checking user defined
-    plugin models.
+# CRUFT: old interface; should be deprecated and removed
+def run_one(model_name):
+    # msg = "use check_model(model_info) rather than run_one(model_name)"
+    # warnings.warn(msg, category=DeprecationWarning, stacklevel=2)
+    try:
+        model_info = load_model_info(model_name)
+    except Exception:
+        output = traceback.format_exc()
+        return output
+
+    success, output = check_model(model_info)
+    return output
+
+def check_model(model_info):
+    # type: (ModelInfo) -> str
+    """
+    Run the tests for a single model, capturing the output.
+
+    Returns success status and the output string.
     """
     # Note that running main() directly did not work from within the
@@ -423 +438 @@
 
     # Build a test suite containing just the model
-    loader = 'opencl' if use_opencl() else 'cuda' if use_cuda() else 'dll'
-    models = [model]
-    try:
-        suite = make_suite([loader], models)
-    except Exception:
-        import traceback
-        stream.writeln(traceback.format_exc())
-        return
+    loaders = ['opencl' if use_opencl() else 'cuda' if use_cuda() else 'dll']
+    suite = unittest.TestSuite()
+    _add_model_to_suite(loaders, suite, model_info)
 
     # Warn if there are no user defined tests.
@@ -445 +455 @@
     for test in suite:
         if not test.info.tests:
-            stream.writeln("Note: %s has no user defined tests."%model)
+            stream.writeln("Note: %s has no user defined tests."%model_info.name)
         break
     else:
@@ -461 +471 @@
     output = stream.getvalue()
     stream.close()
-    return output
+    return result.wasSuccessful(), output
 
 

sasmodels/sasview_model.py

@@ -819 +819 @@
             return value, [value], [1.0]
 
+    @classmethod
+    def runTests(cls):
+        """
+        Run any tests built into the model and captures the test output.
+
+        Returns success flag and output
+        """
+        from .model_test import check_model
+        return check_model(cls._model_info)
+
 def test_cylinder():
     # type: () -> float

sasmodels/special.py

@@ -113 +113 @@
         The standard math function, tgamma(x) is unstable for $x < 1$
         on some platforms.
+
+    sas_gammaln(x):
+        log gamma function sas_gammaln\ $(x) = \log \Gamma(|x|)$.
+
+        The standard math function, lgamma(x), is incorrect for single
+        precision on some platforms.
+
+    sas_gammainc(a, x), sas_gammaincc(a, x):
+        Incomplete gamma function
+        sas_gammainc\ $(a, x) = \int_0^x t^{a-1}e^{-t}\,dt / \Gamma(a)$
+        and complementary incomplete gamma function
+        sas_gammaincc\ $(a, x) = \int_x^\infty t^{a-1}e^{-t}\,dt / \Gamma(a)$
 
     sas_erf(x), sas_erfc(x):
@@ -207 +219 @@
 from numpy import pi, nan, inf
 from scipy.special import gamma as sas_gamma
+from scipy.special import gammaln as sas_gammaln
+from scipy.special import gammainc as sas_gammainc
+from scipy.special import gammaincc as sas_gammaincc
 from scipy.special import erf as sas_erf
 from scipy.special import erfc as sas_erfc