Changeset 0a9fcab in sasmodels


Ignore:
Timestamp:
Dec 4, 2017 8:13:55 AM (6 years ago)
Author:
Paul Kienzle <pkienzle@…>
Children:
2db9fe4
Parents:
ef6a512
Message:

simplify handling of gauss.z in py2c

Files:
3 edited

Legend:

Unmodified
Added
Removed

  • doc/guide/plugin.rst

    r3048ec6 r0a9fcab  
    543543    M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E: 
    544544        $\pi$, $\pi/2$, $\pi/4$, $1/\sqrt{2}$ and Euler's constant $e$ 
    545     exp, log, pow(x,y), expm1, sqrt: 
    546         Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\sqrt{x}$. 
    547         The function expm1(x) is accurate across all $x$, including $x$ 
    548         very close to zero. 
     545    exp, log, pow(x,y), expm1, log1p, sqrt, cbrt: 
     546        Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\ln 1 + x$, 
     547        $\sqrt{x}$, $\sqrt[3]{x}$. The functions expm1(x) and log1p(x) 
     548        are accurate across all $x$, including $x$ very close to zero. 
    549549    sin, cos, tan, asin, acos, atan: 
    550550        Trigonometry functions and inverses, operating on radians. 
    … …  
    557557        atan(y/x) would return a value in quadrant I. Similarly for 
    558558        quadrants II and IV when $x$ and $y$ have opposite sign. 
    559     fmin(x,y), fmax(x,y), trunc, rint: 
     559    fabs(x), fmin(x,y), fmax(x,y), trunc, rint: 
    560560        Floating point functions.  rint(x) returns the nearest integer. 
    561561    NAN: 
    562562        NaN, Not a Number, $0/0$.  Use isnan(x) to test for NaN.  Note that 
    563563        you cannot use :code:`x == NAN` to test for NaN values since that 
    564         will always return false.  NAN does not equal NAN! 
     564        will always return false.  NAN does not equal NAN!  The alternative, 
     565        :code:`x != x` may fail if the compiler optimizes the test away. 
    565566    INFINITY: 
    566567        $\infty, 1/0$.  Use isinf(x) to test for infinity, or isfinite(x) 
    … …  
    734735        Similar arrays are available in :code:`gauss20.c` for 20-point 
    735736        quadrature and in :code:`gauss150.c` for 150-point quadrature. 
     737        The macros :code:`GAUSS_N`, :code:`GAUSS_Z` and :code:`GAUSS_W` are 
     738        defined so that you can change the order of the integration by 
     739        selecting an different source without touching the C code. 
    736740 
    737741        :code:`source = ["lib/gauss76.c", ...]` 

  • sasmodels/autoc.py

    ref6a512 r0a9fcab  
    6767    while translate: 
    6868        function_name, function = translate.pop(0) 
    69         source = inspect.getsource(function) 
    70         tree = ast.parse(source) 
    7169        filename = function.__code__.co_filename 
    7270        offset = function.__code__.co_firstlineno 
    7371        refs = function.__code__.co_names 
    74         snippet = codegen.to_source(tree) #, filename, offset) 
    75         code[function_name] = snippet 
    7672        depends[function_name] = set(refs) 
     73        source = inspect.getsource(function) 
    7774        for name in refs: 
    7875            if name in tagged or name in DEFINES: 
    … …  
    9996            elif isinstance(obj, special.Gauss): 
    10097                libs.append('lib/gauss%d.c'%obj.n) 
     98                source = (source.replace(name+'.n', 'GAUSS_N') 
     99                          .replace(name+'.z', 'GAUSS_Z') 
     100                          .replace(name+'.w', 'GAUSS_W')) 
    101101            else: 
    102102                raise TypeError("Could not convert global %s of type %s" 
    103103                                % (name, str(type(obj)))) 
    104104 
     105        tree = ast.parse(source) 
     106        snippet = codegen.to_source(tree) #, filename, offset) 
     107        code[function_name] = snippet 
     108 
    105109    # remove duplicates from the dependecy list 
    106     unique_libs= [] 
     110    unique_libs = [] 
    107111    for filename in libs: 
    108112        if filename not in unique_libs: 

  • sasmodels/special.py

    ref6a512 r0a9fcab  
    33................. 
    44 
    5 The C code follows the C99 standard, with the usual math functions, 
    6 as defined in 
    7 `OpenCL <https://www.khronos.org/registry/cl/sdk/1.1/docs/man/xhtml/mathFunctions.html>`_. 
    8 This includes the following: 
     5This following standard C99 math functions are available: 
    96 
    107    M_PI, M_PI_2, M_PI_4, M_SQRT1_2, M_E: 
    118        $\pi$, $\pi/2$, $\pi/4$, $1/\sqrt{2}$ and Euler's constant $e$ 
    129 
    13     exp, log, pow(x,y), expm1, sqrt: 
    14         Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\sqrt{x}$. 
    15         The function expm1(x) is accurate across all $x$, including $x$ 
    16         very close to zero. 
     10    exp, log, pow(x,y), expm1, log1p, sqrt, cbrt: 
     11        Power functions $e^x$, $\ln x$, $x^y$, $e^x - 1$, $\ln 1 + x$, 
     12        $\sqrt{x}$, $\sqrt[3]{x}$. The functions expm1(x) and log1p(x) 
     13        are accurate across all $x$, including $x$ very close to zero. 
    1714 
    1815    sin, cos, tan, asin, acos, atan: 
    … …  
    2926        quadrants II and IV when $x$ and $y$ have opposite sign. 
    3027 
    31     fmin(x,y), fmax(x,y), trunc, rint: 
     28    fabs(x), fmin(x,y), fmax(x,y), trunc, rint: 
    3229        Floating point functions.  rint(x) returns the nearest integer. 
    3330 
    … …  
    3532        NaN, Not a Number, $0/0$.  Use isnan(x) to test for NaN.  Note that 
    3633        you cannot use :code:`x == NAN` to test for NaN values since that 
    37         will always return false.  NAN does not equal NAN! 
     34        will always return false.  NAN does not equal NAN!  The alternative, 
     35        :code:`x != x` may fail if the compiler optimizes the test away. 
    3836 
    3937    INFINITY: 
    … …  
    8987        for forcing a constant to stay double precision. 
    9088 
    91 The following special functions and scattering calculations are defined in 
    92 `sasmodels/models/lib <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib>`_. 
     89The following special functions and scattering calculations are defined. 
    9390These functions have been tuned to be fast and numerically stable down 
    9491to $q=0$ even in single precision.  In some cases they work around bugs 
    … …  
    184181 
    185182 
    186     Gauss76Z[i], Gauss76Wt[i]: 
     183    gauss76.n, gauss76.z[i], gauss76.w[i]: 
    187184        Points $z_i$ and weights $w_i$ for 76-point Gaussian quadrature, respectively, 
    188185        computing $\int_{-1}^1 f(z)\,dz \approx \sum_{i=1}^{76} w_i\,f(z_i)$. 
    189  
    190         Similar arrays are available in :code:`gauss20.c` for 20-point 
    191         quadrature and in :code:`gauss150.c` for 150-point quadrature. 
    192  
     186        When translating the model to C, include 'lib/gauss76.c' in the source 
     187        and use :code:`GAUSS_N`, :code:`GAUSS_Z`, and :code:`GAUSS_W`. 
     188 
     189        Similar arrays are available in :code:`gauss20` for 20-point quadrature 
     190        and :code:`gauss150.c` for 150-point quadrature. By using 
     191        :code:`import gauss76 as gauss` it is easy to change the number of 
     192        points in the integration. 
    193193""" 
    194194# pylint: disable=unused-import 
    … …  
    200200 
    201201# C99 standard math library functions 
    202 from numpy import exp, log, power as pow, expm1, sqrt 
     202from numpy import exp, log, power as pow, expm1, logp1, sqrt, cbrt 
    203203from numpy import sin, cos, tan, arcsin as asin, arccos as acos, arctan as atan 
    204204from numpy import sinh, cosh, tanh, arcsinh as asinh, arccosh as acosh, arctanh as atanh 
    205205from numpy import arctan2 as atan2 
    206 from numpy import fmin, fmax, trunc, rint 
     206from numpy import fabs, fmin, fmax, trunc, rint 
    207207from numpy import pi, nan, inf 
    208208from scipy.special import gamma as sas_gamma 
Note: See TracChangeset for help on using the changeset viewer.