Changes in explore/precision.py [fba9ca0:ee60aa7] in sasmodels

File:

• explore/precision.py (modified) (13 diffs)

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"
 
@@ -105 +102 @@
         linear = not xrange.startswith("log")
         if xrange == "zoom":
-            start, stop, steps = 1000, 1010, 2000
+            lin_min, lin_max, lin_steps = 1000, 1010, 2000
         elif xrange == "neg":
-            start, stop, steps = -100.1, 100.1, 2000
+            lin_min, lin_max, lin_steps = -100.1, 100.1, 2000
         elif xrange == "linear":
-            start, stop, steps = 1, 1000, 2000
-            start, stop, steps = 0.001, 2, 2000
+            lin_min, lin_max, lin_steps = 1, 1000, 2000
+            lin_min, lin_max, lin_steps = 0.001, 2, 2000
         elif xrange == "log":
-            start, stop, steps = -3, 5, 400
+            log_min, log_max, log_steps = -3, 5, 400
         elif xrange == "logq":
-            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])
-
+            log_min, log_max, log_steps = -4, 1, 400
         else:
             raise ValueError("unknown range "+xrange)
@@ -131 +121 @@
             # value to x in the given precision.
             if linear:
-                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')
+                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')
                 qr = [mp.mpf(float(v)) for v in qrf]
-                #qr = mp.linspace(start, stop, steps)
+                #qr = mp.linspace(lin_min, lin_max, lin_steps)
             else:
-                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')
+                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')
                 qr = [mp.mpf(float(v)) for v in qrf]
-                #qr = [10**v for v in mp.linspace(start, stop, steps)]
+                #qr = [10**v for v in mp.linspace(log_min, log_max, log_steps)]
 
         target = self.call_mpmath(qr, bits=500)
@@ -186 +176 @@
     """
     if diff == "relative":
-        err = np.array([(abs((t-a)/t) if t != 0 else a) for t, a in zip(target, actual)], 'd')
+        err = np.array([abs((t-a)/t) for t, a in zip(target, actual)], 'd')
         #err = np.clip(err, 0, 1)
         pylab.loglog(x, err, '-', label=label)
@@ -207 +197 @@
     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 ================
@@ -326 +299 @@
 )
 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(
     name="erf(x)",
     mp_function=mp.erf,
@@ -403 +357 @@
 )
 add_function(
-    name="debye",
+    name="gauss_coil",
     mp_function=lambda x: 2*(mp.exp(-x**2) + x**2 - 1)/x**4,
     np_function=lambda x: 2*(np.expm1(-x**2) + x**2)/x**4,
     ocl_function=make_ocl("""
     const double qsq = q*q;
-    if (qsq < 1.0) { // Pade approximation
+    // For double: use O(5) Pade with 0.5 cutoff (10 mad + 1 divide)
+    // For single: use O(7) Taylor with 0.8 cutoff (7 mad)
+    if (qsq < 0.0) {
         const double x = qsq;
         if (0) { // 0.36 single
@@ -418 +374 @@
             const double B1=3./8., B2=3./56., B3=1./336.;
             return (((A3*x + A2)*x + A1)*x + 1.)/(((B3*x + B2)*x + B1)*x + 1.);
-        } else if (1) { // 1.0 for single, 0.25 for double
+        } else if (0) { // 1.0 for single, 0.25 for double
             // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}]
             const double A1=1./15., A2=1./60, A3=0., A4=1./75600.;
@@ -431 +387 @@
                   /(((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.);
         }
-    } else if (qsq < 1.) { // Taylor series; 0.9 for single, 0.25 for double
+    } else if (qsq < 0.8) {
         const double x = qsq;
         const double C0 = +1.;
@@ -509 +465 @@
 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
             };
@@ -521 +477 @@
 """
 add_function(
-    name="loggamma(x)",
+    name="log gamma(x)",
     mp_function=mp.loggamma,
     np_function=scipy.special.gammaln,
@@ -645 +601 @@
 
 ALL_FUNCTIONS = set(FUNCTIONS.keys())
-ALL_FUNCTIONS.discard("loggamma")  # use cephes-based gammaln instead
+ALL_FUNCTIONS.discard("loggamma")  # OCL version not ready yet
 ALL_FUNCTIONS.discard("3j1/x:taylor")
 ALL_FUNCTIONS.discard("3j1/x:trig")
@@ -661 +617 @@
     -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]
-        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:
+      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:
     """+names)
     sys.exit(1)