use adaptive integration in sasmodels

[pkienzle]

Adaptive integration can get much more accurate results for the same number of function evaluations. This will require a rewrite of the polydispersity loops, the resolution calculation, and the spherical integration used to evaluate 1D models from oriented shapes. Our simple use of opencl which evaluates the form factor over a fixed set of q will become much more complicated.

See ticket #312 for slit_test.py, an implementation of slit smearing with adaptive and non-adaptive integration using scipy.

[pkienzle]

At least we should move to using trapezoid rule rather than mid-point rule to make better use of our calculated q points.

For both Q_parallel and Q_perpendicular resolution, the trapezoid rule implementation in resolution.romberg_slit_1d (yes, the function is poorly named…) give a better result than the current implementation in resolution.slit_resolution. We should be able to incorporate the interpolation into the weight matrix calculation without too much trouble, and give a better result for no additional cost.

[pkienzle]

Comments pulled in from #312 so we can close that ticket:

I compared both Igor and sasview to a direct implementation of the slit smearing integral using adaptive integration. Although neither is correct, the Igor implementation is much closer.

Looking at the code in sasview, I do not see where the residual term is added. This is the power law approximation to the tail of the distribution from q_max to q_max + Delta q_vertical. Maybe that is enough to explain the difference.

I implemented replacement resolution functions in sasmodels, though I do not yet use them in the fit. The basic idea is to continue to use a weight matrix, but allow it to be non-square, with calculated q values extrapolated to the full range of the resolution function. At least for the test problem (spheres of 6 micro radius), rectangular integration gets with 2.5% relative error for about the same number of evaluations. Trapezoid and Simpson’s rule don’t do any better. [Note: still hard to believe that trapz and simpson don't out-perform rectangular; may want to try across a broader range of shapes and sizes.]

The attached file (slit_test.py) will let you explore the various implementations. It contains a python implementation of the sphere model so that it depends only on scipy, numpy and matplotlib.

Made slit_test.py into a GIST so that we can edit it.

​https://gist.github.com/pkienzle/2faa18591b7d82d4d62d

Things to try:

use power law interpolation between points rather than linear interpolation.

use sampling around each point to average over the oscillations

use logarithmic rather than linear bin edges

These can be tried independently or together.

[pkienzle]

Note that the resolution matrix is fixed throughout the fit; it would be easy enough to through this matrix onto the GPU and apply the resolution before returning the q values.

[pkienzle]

In lieu of adaptive integration, we could define an oversampling factor based on the maximum length scale of the form, so that we have at least N points in each 2 pi/d interval.