Simple adaptive integration

explore/precision.py

@@ -803,6 +803,21 @@ def np_2J1x_x(x):
     ocl_function=make_ocl("return sas_2J1x_x(q);", "sas_2J1x_x", ["lib/polevl.c", "lib/sas_J1.c"]),
 )
 
+add_function(
+    name="sin_arccos",
+    mp_function=lambda x: (mp.sin(mp.acos(x))),
+    np_function=lambda x: (np.sin(np.arccos(x))),
+    ocl_function=make_ocl("return sin(acos(q));", "sas_sin_arccos"),
+    limits=(0, 1),
+)
+add_function(
+    name="sin_from_cos",
+    mp_function=lambda x: (mp.sqrt(1 - x*x)),
+    np_function=lambda x: (np.sqrt(1-x*x)),
+    ocl_function=make_ocl("return sqrt(1.-q*q);", "sas_sin_from_cos"),
+    limits=(0, 1),
+)
+
 ALL_FUNCTIONS = set(FUNCTIONS.keys())
 ALL_FUNCTIONS.discard("loggamma")  # use cephes-based gammaln instead
 ALL_FUNCTIONS.discard("3j1/x:taylor")

sasmodels/compare.py

@@ -75,7 +75,7 @@ def __call__(self, **par: float) -> np.ndarray: ...
 
     === model parameters ===
     -preset*/-random[=seed] preset or random parameters
-    -sets=n generates n random datasets with the seed given by -random=seed
+    -sets=n generates n random datasets using -seed=seed
     -pars/-nopars* prints the parameter set or not
     -sphere[=150] set up spherical integration over theta/phi using n points
     -mono*/-poly suppress or allow polydispersity on generated parameters
@@ -1049,7 +1049,7 @@ def plot_models(opts, result, limits=None, setnum=0):
     '2d', '1d', 'sesans',
 
     # Parameter set
-    'preset', 'random', 'random=', 'sets=',
+    'preset', 'random', 'random=', 'sets=', 'seed=',
     'nopars', 'pars',
     'sphere', 'sphere=', # integrate over a sphere in 2d with n points
     'poly', 'mono',
@@ -1253,6 +1253,7 @@ def parse_opts(argv):
         elif arg.startswith('-res='):      opts['res'] = arg[5:]
         elif arg.startswith('-noise='):    opts['noise'] = float(arg[7:])
         elif arg.startswith('-sets='):     opts['sets'] = int(arg[6:])
+        elif arg.startswith('-seed='):     opts['seed'] = int(arg[6:])
         elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:]
         elif arg.startswith('-cutoff='):   opts['cutoff'] = arg[8:]
         elif arg.startswith('-title='):    opts['title'] = arg[7:]

sasmodels/direct_model.py

@@ -453,14 +453,18 @@ def test_reparameterize():
     except Exception:
         pass
 
-def _direct_calculate(model, data, pars):
+def _direct_calculate(model, data, pars, ngauss=0):
     from .core import build_model, load_model_info
+    from .generate import set_integration_size
+
     model_info = load_model_info(model)
+    if ngauss != 0:
+        set_integration_size(model_info, ngauss)
     kernel = build_model(model_info)
     calculator = DirectModel(data, kernel)
     return calculator(**pars)
 
-def Iq(model, q, dq=None, ql=None, qw=None, **pars):
+def Iq(model, q, dq=None, ql=None, qw=None, ngauss=0, **pars):
     """
     Compute I(q) for *model*. Resolution is *dq* for pinhole or *ql* and *qw*
     for slit geometry. Use 0 or None for infinite slits.
@@ -498,16 +502,16 @@ def broadcast(v):
             else np.full(len(q), v) if np.isscalar(v)
             else _as_numpy(v))
     data.dxl, data.dxw = broadcast(ql), broadcast(qw)
-    return _direct_calculate(model, data, pars)
+    return _direct_calculate(model, data, pars, ngauss=ngauss)
 
-def Iqxy(model, qx, qy, dqx=None, dqy=None, **pars):
+def Iqxy(model, qx, qy, dqx=None, dqy=None, ngauss=0, **pars):
     """
     Compute I(qx, qy) for *model*. Resolution is *dqx* and *dqy*.
     See :func:`Iq` for details on model and parameters.
     """
     from .data import Data2D
     data = Data2D(x=qx, y=qy, dx=dqx, dy=dqy)
-    return _direct_calculate(model, data, pars)
+    return _direct_calculate(model, data, pars, ngauss=ngauss)
 
 def Gxi(model, xi, **pars):
     """
@@ -528,6 +532,8 @@ def main():
     if len(sys.argv) < 3:
         print("usage: python -m sasmodels.direct_model modelname (q|qx,qy) par=val ...")
         sys.exit(1)
+
+    ngauss = 0
     model = sys.argv[1]
     call = sys.argv[2].upper()
     pars = dict((k, (float(v) if not k.endswith("_pd_type") else v))
@@ -542,13 +548,13 @@ def main():
         dq = dqw = dql = None
         #dq = [q*0.05] # 5% pinhole resolution
         dqw, dql = [q*0.05], [1.0] # 5% horizontal slit resolution
-        print(Iq(model, [q], dq=dq, qw=dqw, ql=dql, **pars)[0])
+        print(Iq(model, [q], dq=dq, qw=dqw, ql=dql, ngauss=ngauss, **pars)[0])
         #print(Gxi(model, [q], **pars)[0])
     elif len(values) == 2:
         qx, qy = values
         dq = None
         #dq = [0.005] # 5% pinhole resolution at q = 0.1
-        print(Iqxy(model, [qx], [qy], dqx=dq, dqy=dq, **pars)[0])
+        print(Iqxy(model, [qx], [qy], dqx=dq, dqy=dq, ngauss=ngauss, **pars)[0])
     else:
         print("use q or qx,qy")
         sys.exit(1)

sasmodels/generate.py

@@ -291,13 +291,27 @@ def set_integration_size(info, n):
     Note: this really ought to be a method in modelinfo, but that leads to
     import loops.
     """
-    if info.source and any(lib.startswith('lib/gauss') for lib in info.source):
-        from .gengauss import gengauss
-        path = joinpath(MODEL_PATH, "lib", "gauss%d.c"%n)
-        if not exists(path):
-            gengauss(n, path)
-        info.source = ["lib/gauss%d.c"%n if lib.startswith('lib/gauss')
-                       else lib for lib in info.source]
+    from .gengauss import gengauss
+
+    if not info.source:
+        return
+
+    # Generate the integration points
+    path = joinpath(MODEL_PATH, "lib", f"gauss{n}.c")
+    if not exists(path):
+        # print(f"building Gaussian integration points of size {n} in {str(path)}")
+        gengauss(n, path)
+
+    # Replace adaptive.c or lib/gauss<n>.c
+    try:
+        index = info.source.index("lib/adaptive.c")
+        info.source[index:index+1] = [f"lib/gauss{n}.c", "lib/nonadaptive.c"]
+    except ValueError:
+        for index in range(len(info.source)-1, -1, -1):
+            if info.source[index].startswith("lib/gauss"):
+                info.source[index] = f"lib/gauss{n}.c"
+                break
+    # print("info.source is now", info.source)
 
 def format_units(units):
     # type: (str) -> str

sasmodels/gengauss.py

@@ -21,19 +21,18 @@ def gengauss(n, path):
         array_size = n
 
     with open(path, "w") as fid:
-        fid.write("""\
-// Generated by sasmodels.gengauss.gengauss(%d)
+        fid.write(f"""\
+// Generated by sasmodels.gengauss.gengauss({n})
 
 #ifdef GAUSS_N
 # undef GAUSS_N
 # undef GAUSS_Z
 # undef GAUSS_W
 #endif
-#define GAUSS_N %d
-#define GAUSS_Z Gauss%dZ
-#define GAUSS_W Gauss%dWt
-
-"""%(n, n, n, n))
+#define GAUSS_N {n}
+#define GAUSS_Z Gauss{n}Z
+#define GAUSS_W Gauss{n}Wt
+""")
 
         if array_size != n:
             fid.write("// Note: using array size %d so that it is a multiple of 4\n\n"%array_size)

sasmodels/model_test.py

@@ -579,6 +579,31 @@ def _build_test(test):
     for test in tests:
         yield _build_test(test)
 
+def _generate_target_values(modelname, ngauss=0):
+    from .generate import set_integration_size
+
+    model_info = load_model_info(modelname)
+    if ngauss != 0:
+        set_integration_size(model_info, ngauss)
+    model = build_model(model_info, platform="dll", dtype="d")
+
+    for pars, q, Iq in model_info.tests:
+        qin = q
+        if isinstance(Iq, float):
+            q, Iq = [q], [Iq]
+        if isinstance(q[0], tuple):
+            qx, qy = zip(*q)
+            q_vectors = [np.array(qx), np.array(qy)]
+        else:
+            q_vectors = [np.array(q)]
+        kernel = model.make_kernel(q_vectors)
+        target = np.array(Iq)
+        actual = call_kernel(kernel, pars)
+        if True or (actual != target).any():
+            print("Test:", modelname, pars)
+            print("  q = ", qin)
+            print(f"  current => [{', '.join(f'{v:.15g}' for v in target)}]")
+            print(f"  ngauss={ngauss} => [{', '.join(f'{v:.15g}' for v in actual)}]")
 
 def main():
     # type: () -> int
@@ -601,6 +626,11 @@ def main():
                         help="Engines on which to run the test.  "
                         "Valid values are opencl, cuda, dll, and all. "
                         "Defaults to all if no value is given")
+    parser.add_argument("-t", "--targets", action="store_true",
+                        help="Generate target values for test.")
+    parser.add_argument("--ngauss", type=int, default=10000,
+                        help="Number of gauss points to use in integration for "
+                        "target values. Warning: this is very slow the first time.")
     parser.add_argument("models", nargs="*",
                         help="The names of the models to be tested.  "
                         "If the first model is 'all', then all but the listed "
@@ -630,9 +660,13 @@ def main():
         print("unknown engine " + opts.engine)
         return 1
 
-    runner = TestRunner(verbosity=opts.verbose, **test_args)
-    result = runner.run(make_suite(loaders, opts.models))
-    return 1 if result.failures or result.errors else 0
+    if opts.targets:
+        for model in opts.models:
+            _generate_target_values(model, ngauss=opts.ngauss)
+    else:
+        runner = TestRunner(verbosity=opts.verbose, **test_args)
+        result = runner.run(make_suite(loaders, opts.models))
+        return 1 if result.failures or result.errors else 0
 
 
 if __name__ == "__main__":

sasmodels/models/barbell.c

@@ -19,13 +19,18 @@ _bell_kernel(double qab, double qc, double h, double radius_bell,
     const double m = radius_bell*qc; // cos argument slope
     const double b = (half_length+h)*qc; // cos argument intercept
     const double qab_r = radius_bell*qab; // Q*R*sin(theta)
+
+    const double qr_max = fmax(qab_r, m+b);
+    constant double *w, *z;
+    int n = gauss_weights(qr_max, &w, &z);
+
     double total = 0.0;
-    for (int i = 0; i < GAUSS_N; i++){
-        const double t = GAUSS_Z[i]*zm + zb;
+    for (int i = 0; i < n; i++){
+        const double t = z[i]*zm + zb;
         const double radical = 1.0 - t*t;
         const double bj = sas_2J1x_x(qab_r*sqrt(radical));
         const double Fq = cos(m*t + b) * radical * bj;
-        total += GAUSS_W[i] * Fq;
+        total += w[i] * Fq;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
     const double integral = total*zm;
@@ -110,21 +115,30 @@ Fq(double q,double *F1, double *F2, double sld, double solvent_sld,
     const double h = sqrt(square(radius_bell) - square(radius));
     const double half_length = 0.5*length;
 
+    // The term h comes from solving the right triangle with diagonal
+    // equal to the bell radius and horizontal equal to the bar radius.
+    // The result is the height of the equator above the end of the rod.
+    // To get the total length of bar+bell use bar length + 2*(bell radius + h).
+    // We want the radius, so divide that by two.
+    const double qr_max = q*(half_length + radius_bell + h);
+    constant double *w, *z;
+    int n = gauss_weights(qr_max, &w, &z);
+
     // translate a point in [-1,1] to a point in [0, pi/2]
     const double zm = M_PI_4;
     const double zb = M_PI_4;
     double total_F1 = 0.0;
     double total_F2 = 0.0;
-    for (int i = 0; i < GAUSS_N; i++){
-        const double theta = GAUSS_Z[i]*zm + zb;
+    for (int i = 0; i < n; i++){
+        const double theta = z[i]*zm + zb;
         double sin_theta, cos_theta; // slots to hold sincos function output
         SINCOS(theta, sin_theta, cos_theta);
         const double qab = q*sin_theta;
         const double qc = q*cos_theta;
         const double Aq = _fq(qab, qc, h, radius_bell, radius, half_length);
         // scale by sin_theta for spherical coord integration
-        total_F1 += GAUSS_W[i] * Aq * sin_theta;
-        total_F2 += GAUSS_W[i] * Aq * Aq * sin_theta;
+        total_F1 += w[i] * Aq * sin_theta;
+        total_F2 += w[i] * Aq * Aq * sin_theta;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
     const double form_avg = total_F1 * zm;

sasmodels/models/barbell.py

@@ -117,7 +117,7 @@
     ]
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "barbell.c"]
+source = ["lib/polevl.c", "lib/sas_J1.c", "lib/adaptive.c", "barbell.c"]
 valid = "radius_bell >= radius"
 have_Fq = True
 radius_effective_modes = [

sasmodels/models/capped_cylinder.c

@@ -7,7 +7,7 @@
 //   length is the cylinder length, or the separation between the lens halves
 //   theta is the angle of the cylinder wrt q.
 static double
-_cap_kernel(double qab, double qc, double h, double radius_cap,
+_cap_kernel(double qab, double qc, double h, double radius_cap, double radius,
     double half_length)
 {
     // translate a point in [-1,1] to a point in [lower,upper]
@@ -26,13 +26,25 @@ _cap_kernel(double qab, double qc, double h, double radius_cap,
     const double m = radius_cap*qc; // cos argument slope
     const double b = (half_length+h)*qc; // cos argument intercept
     const double qab_r = radius_cap*qab; // Q*R*sin(theta)
+
+    // m+b = qc*(half_length + radius_cap + h). With h in [-radius_cap, 0] depending
+    // on cylinder radius, that means m+b is in qc*[length/2, length_2 + radius_cap].
+    // The qab_r term will be very large for mostly flat caps. Since the bj term will
+    // oscillate at this frequency, it seems like we should increase the number of
+    // gauss points to accomodate. However, if we use the radius of the cylinder
+    // we seem to get good results, so use that to set the number of integration points.
+    //const double qr_max = fmax(qab_r, m+b);
+    const double qr_max = fmax(qab*radius, m+b);
+    constant double *w, *z;
+    int n = gauss_weights(qr_max, &w, &z);
+
     double total = 0.0;
-    for (int i=0; i<GAUSS_N; i++) {
-        const double t = GAUSS_Z[i]*zm + zb;
+    for (int i=0; i<n; i++) {
+        const double t = z[i]*zm + zb;
         const double radical = 1.0 - t*t;
         const double bj = sas_2J1x_x(qab_r*sqrt(radical));
         const double Fq = cos(m*t + b) * radical * bj;
-        total += GAUSS_W[i] * Fq;
+        total += w[i] * Fq;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
     const double integral = total*zm;
@@ -43,7 +55,7 @@ _cap_kernel(double qab, double qc, double h, double radius_cap,
 static double
 _fq(double qab, double qc, double h, double radius_cap, double radius, double half_length)
 {
-    const double cap_Fq = _cap_kernel(qab, qc, h, radius_cap, half_length);
+    const double cap_Fq = _cap_kernel(qab, qc, h, radius_cap, radius, half_length);
     const double bj = sas_2J1x_x(radius*qab);
     const double si = sas_sinx_x(half_length*qc);
     const double cyl_Fq = 2.0*M_PI*radius*radius*half_length*bj*si;
@@ -114,21 +126,33 @@ Fq(double q,double *F1, double *F2, double sld, double solvent_sld,
     const double h = -sqrt(square(radius_cap) - square(radius));
     const double half_length = 0.5*length;
 
+    // The term h comes from solving the right triangle with diagonal
+    // equal to the cap radius and horizontal equal to the bar radius.
+    // The result is the (negative) height of the equator above the end of the rod.
+    // To get the total length of bar+cap use bar length + 2*(cap radius + h).
+    // We want the radius, so divide that by two.
+    // For a lentil with length=0, the radius will be the dominant term, hence
+    // the fmax in the calculation below. This isn't needed for the barbell shape
+    // since the bell length is always greater than the bar radius.
+    const double qr_max = q*fmax(half_length + radius_cap + h, radius);
+    constant double *w, *z;
+    int n = gauss_weights(qr_max, &w, &z);
+
     // translate a point in [-1,1] to a point in [0, pi/2]
     const double zm = M_PI_4;
     const double zb = M_PI_4;
     double total_F1 = 0.0;
     double total_F2 = 0.0;
-    for (int i=0; i<GAUSS_N ;i++) {
-        const double theta = GAUSS_Z[i]*zm + zb;
+    for (int i=0; i<n ;i++) {
+        const double theta = z[i]*zm + zb;
         double sin_theta, cos_theta; // slots to hold sincos function output
         SINCOS(theta, sin_theta, cos_theta);
         const double qab = q*sin_theta;
         const double qc = q*cos_theta;
         const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length);
         // scale by sin_theta for spherical coord integration
-        total_F1 += GAUSS_W[i] * Aq * sin_theta;
-        total_F2 += GAUSS_W[i] * Aq * Aq * sin_theta;
+        total_F1 += w[i] * Aq * sin_theta;
+        total_F2 += w[i] * Aq * Aq * sin_theta;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
     const double form_avg = total_F1 * zm;

sasmodels/models/capped_cylinder.py

@@ -137,7 +137,7 @@
              ]
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "capped_cylinder.c"]
+source = ["lib/polevl.c", "lib/sas_J1.c", "lib/adaptive.c", "capped_cylinder.c"]
 valid = "radius_cap >= radius"
 have_Fq = True
 radius_effective_modes = [