Simple adaptive integration

[pkienzle]

Alternative to the #608 using a simple heuristic based on qr.

Implements adaptive integration for all shapes except superball. The paracrystal models (bcc, fcc, sc) need a different approach.

Accuracy is usually comparable to a 10000 point gaussian integration for every qr. The target is 0.1% difference, though it isn't always achieved. For example:

$ python -m sasmodels.compare background=0 core_shell_cylinder -ngauss=0,10000 -engine=single,single! -random=83174 -nq=1000 -pars -neval=10
Randomize using -random=83174
scale: 0.291548
background: 0
sld_core: 8.92736
sld_shell: 11.2834
sld_solvent: 10.7719
radius: 291.224
thickness: 33528.1
length: 1.53983
GPU[32] t=12.22 ms, intensity=36122148864
DLL[32] t=293.91 ms, intensity=36122148864
|GPU[32]-DLL[32]|            max:3.648e+03  median:2.930e-03  98%:6.400e+02  rms:2.012e+02  zero-offset:+3.928e+01
|(GPU[32]-DLL[32])/DLL[32]|  max:3.443e-02  median:2.903e-06  98%:2.554e-03  rms:1.932e-03  zero-offset:+2.716e-04

Because we include a 20 point gaussian integration scheme, speed is frequently faster than the fixed 76 point gaussian integration in master, at least for small shapes. For large shapes it can be several times slower than the fixed scheme, though the increase in accuracy easily justifies the cost.

Shapes with nested integrals (e.g., triaxial ellipsoid) can be very slow. For example:

$ python -m sasmodels.compare background=0 triaxial_ellipsoid -random=709580 -nq=100 -pars
Randomize using -random=709580
scale: 0.0515113
background: 0
sld: 4.27038
sld_solvent: 8.05591
radius_equat_minor: 30.0115
radius_equat_major: 28083
radius_polar: 42754.1
GPU[32] t=11001.53 ms, intensity=69025144

Because the cost for a 10000 point gaussian with nested integration is so high these models have only be checked for accuracy at a few Q points.

Refs #248

[pkienzle]

Example of bad triaxial ellipsoid (20% error):

$ python -m sasmodels.compare background=0 triaxial_ellipsoid -ngauss=0,10000 -engine=single,single! -nq=30 -random=716856 -pars
Randomize using -random=716856
scale: 0.00343363
background: 0
sld: 11.2141
sld_solvent: 10.9297
radius_equat_minor: 41.5349
radius_equat_major: 9142.92
radius_polar: 74.0436
GPU[32] t=58.31 ms, intensity=33
DLL[32] t=12646.33 ms, intensity=33
|GPU[32]-DLL[32]|            max:1.941e-03  median:1.907e-06  98%:1.849e-03  rms:6.002e-04  zero-offset:+2.745e-04
|(GPU[32]-DLL[32])/DLL[32]|  max:1.884e-01  median:1.179e-06  98%:1.810e-01  rms:5.891e-02  zero-offset:+2.319e-02

The fixed 76 point integration scheme works better for this example (0.3% error).

Maybe it is worth exploring Lebedev and other surface quadrature schemes for these nested integrals. It is messy, though, because not all of them are of the form ∫∫ F(q) sin(θ) dφ dθ.

[butlerpd]

This was briefly discussed at today's fortnightly call and tagged as of interest to the upcoming camp. Question is whether it provides a minimal change to provide a reasonable speedup. It is noted that this PR not only adds the new adaptive integation it changes all the model files that currently use the GaussXX methods with this one. Probably would have been cleaner as two separate PRs?

Also at issue is what to do with the integration speedup already proposed a few years earlier and sitting in #608
Michael Wagner agreed to look at this.

NOTE: there are conflicts that will need to be resolved before this can be merged

[pkienzle]

This works well for rotationally symmetric shapes that only use 1D integrals.

Performance is unsatisfactory on shapes such as triaxial ellipsoid that need 2D integrals.

I could revert changes for those models until we've had a chance to explore other schemes such as Lebedev or Fibonacci.