Changeset bc873053 in sasview

Timestamp:

Nov 10, 2015 3:39:44 PM (9 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload

Children:

80ba1a2

Parents:

40c69b3

Message:

Fix slit height resolution problem. Still issues with slit width. Refs #472

Files:

• .gitignore (modified) (1 diff)
• src/sas/models/resolution.py (modified) (20 diffs)
• test/sasfit/test/utest_fit_line.py (modified) (2 diffs)
• test/sasfit/test/utest_fit_smeared.py (modified) (2 diffs)

.gitignore

@@ -23 +23 @@
 **/test/logs
 *.pyc
+*.so
 default_categories.json
 

src/sas/models/resolution.py

@@ -11 +11 @@
 MINIMUM_RESOLUTION = 1e-8
 
-class Resolution1D(object):
+class Resolution(object):
     """
     Abstract base class defining a 1D resolution function.
@@ -32 +32 @@
 
 
-class Perfect1D(Resolution1D):
+class Perfect1D(Resolution):
     """
     Resolution function to use when there is no actual resolution smearing
@@ -45 +45 @@
 
 
-class Pinhole1D(Resolution1D):
+class Pinhole1D(Resolution):
     r"""
     Pinhole aperture with q-dependent gaussian resolution.
@@ -77 +77 @@
 
 
-class Slit1D(Resolution1D):
+class Slit1D(Resolution):
     """
     Slit aperture with a complicated resolution function.
@@ -92 +92 @@
     The *weight_matrix* is computed by :func:`slit1d_resolution`
     """
-    def __init__(self, q, width, height, q_calc=None):
-        # TODO: maybe issue warnings rather than raising errors
-        if not np.isscalar(width):
-            if np.any(np.diff(width) > 0.0):
-                raise ValueError("Slit resolution requires fixed width slits")
-            width = width[0]
-        if not np.isscalar(height):
-            if np.any(np.diff(height) > 0.0):
-                raise ValueError("Slit resolution requires fixed height slits")
-            height = height[0]
-
+    def __init__(self, q, width, height=0., q_calc=None):
         # Remember what width/height was used even though we won't need them
         # after the weight matrix is constructed
         self.width, self.height = width, height
+
+        # Allow independent resolution on each point even though it is not
+        # needed in practice.
+        if np.isscalar(width):
+            width = np.ones(len(q))*width
+        else:
+            width = np.asarray(width)
+        if np.isscalar(height):
+            height = np.ones(len(q))*height
+        else:
+            height = np.asarray(height)
 
         self.q = q.flatten()
@@ -169 +170 @@
             \int_{-\Delta q_v}^{\Delta q_v} I(u) du
     """
-    if width == 0.0 and height == 0.0:
-        #print "condition zero"
-        return 1
-
+    #np.set_printoptions(precision=6, linewidth=10000)
+
+    # The current algorithm is a midpoint rectangle rule.
     q_edges = bin_edges(q_calc) # Note: requires q > 0
     q_edges[q_edges<0.0] = 0.0 # clip edges below zero
-
-    #np.set_printoptions(linewidth=10000)
-    if width <= 100.0 * height or height == 0:
-        # The current algorithm is a midpoint rectangle rule.  In the test case,
-        # neither trapezoid nor Simpson's rule improved the accuracy.
-        #print "condition h", q_edges.shape, q.shape, q_calc.shape
-        weights = np.zeros((len(q), len(q_calc)), 'd')
-        for i, qi in enumerate(q):
-            weights[i, :] = np.diff(q_to_u(q_edges, qi))
-        weights /= width
-        weights = weights.T
-    else:
-        #print "condition w"
-        # Make q_calc into a row vector, and q into a column vector
-        q, q_calc = q[None,:], q_calc[:,None]
-        in_x = (q_calc >= q-width) * (q_calc <= q+width)
-        weights = np.diff(q_edges)[:,None] * in_x
-
-    return weights
+    weights = np.zeros((len(q), len(q_calc)), 'd')
+
+    for i, (qi, w, h) in enumerate(zip(q, width, height)):
+        if w == 0. and h == 0.:
+            # Perfect resolution, so return the theory value directly.
+            # Note: assumes that q is a subset of q_calc.  If qi need not be
+            # in q_calc, then we can do a weighted interpolation by looking
+            # up qi in q_calc, then weighting the result by the relative
+            # distance to the neighbouring points.
+            weights[i, :] = (q_calc == qi)
+        elif h == 0 or w > 100.0 * h:
+            # Convert bin edges from q to u
+            u_limit = np.sqrt(qi**2 + w**2)
+            u_edges = q_edges**2 - qi**2
+            u_edges[q_edges < qi] = 0.
+            u_edges[q_edges > u_limit] = u_limit**2 - qi**2
+            weights[i, :] = np.diff(np.sqrt(u_edges))/w
+            #print "i, qi",i,qi,qi+width
+            #print q_calc
+            #print weights[i,:]
+        elif w == 0:
+            in_x = (q_calc >= qi-h) * (q_calc <= qi+h)
+            weights[i,:] = in_x * np.diff(q_edges) / (2*h)
+
+    return weights.T
 
 
@@ -212 +218 @@
     function.
     """
-    height # keep lint happy
     q_min, q_max = np.min(q), np.max(np.sqrt(q**2 + width**2))
     return geometric_extrapolation(q, q_min, q_max)
@@ -233 +238 @@
         ])
     return edges
-
-
-def q_to_u(q, q0):
-    """
-    Convert *q* values to *u* values for the integral computed at *q0*.
-    """
-    # array([value])**2 - value**2 is not always zero
-    qpsq = q**2 - q0**2
-    qpsq[qpsq<0] = 0
-    return sqrt(qpsq)
 
 
@@ -297 +292 @@
     *points_per_decade* sets the ratio between consecutive steps such
     that there will be $n$ points used for every factor of 10 increase
-    in $q$.
+    in *q*.
 
     If *points_per_decade* is not given, it will be estimated as follows.
@@ -316 +311 @@
     Substituting:
 
-    .. math::
-
         n_\text{extend} = (n-1) (\log q_\text{max} - \log q_n)
-            / (\log q_n - \log q_1)
+            / (\log q_n - log q_1)
     """
     q = np.sort(q)
@@ -359 +352 @@
 
 
-def romberg_slit_1d(q, delta_qv, form, pars):
+def romberg_slit_1d(q, width, height, form, pars):
     """
     Romberg integration for slit resolution.
@@ -366 +359 @@
     that make it slow to evaluate but give it good accuracy.
     """
-    from scipy.integrate import romberg
+    from scipy.integrate import romberg, dblquad
 
     if any(k not in form.info['defaults'] for k in pars.keys()):
@@ -374 +367 @@
                          (", ".join(sorted(extra)), ", ".join(sorted(keys))))
 
-    _fn = lambda u, q0: eval_form(sqrt(q0**2 + u**2), form, pars)
-    r = [romberg(_fn, 0, delta_qv[0], args=(qi,),
-                 divmax=100, vec_func=True, tol=0, rtol=1e-8)
-         for qi in q]
+    if np.isscalar(width):
+        width = [width]*len(q)
+    if np.isscalar(height):
+        height = [height]*len(q)
+    _int_w = lambda w, qi: eval_form(sqrt(qi**2 + w**2), form, pars)
+    _int_h = lambda h, qi: eval_form(abs(qi+h), form, pars)
+    _int_wh = lambda w, h, qi: eval_form(sqrt((qi+h)**2 + w**2), form, pars)
+    result = np.empty(len(q))
+    for i, (qi, w, h) in enumerate(zip(q, width, height)):
+        if h == 0.:
+            r = romberg(_int_w, 0, w, args=(qi,),
+                        divmax=100, vec_func=True, tol=0, rtol=1e-8)
+            result[i] = r/w
+        elif w == 0.:
+            r = romberg(_int_h, -h, h, args=(qi,),
+                        divmax=100, vec_func=True, tol=0, rtol=1e-8)
+            result[i] = r/(2*h)
+        else:
+            r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: width,
+                             args=(qi,))
+            result[i] = r/(w*2*h)
+
     # r should be [float, ...], but it is [array([float]), array([float]),...]
-    return np.asarray(r).flatten()/delta_qv[0]
+    return result
 
 
@@ -520 +531 @@
 
     def compare(self, q, output, answer, tolerance):
-        err = (output - answer)/answer
-        idx = abs(err) >= tolerance
-        problem = zip(q[idx], output[idx], answer[idx], err[idx])
-        print "\n".join(str(v) for v in problem)
+        #err = (output - answer)/answer
+        #idx = abs(err) >= tolerance
+        #problem = zip(q[idx], output[idx], answer[idx], err[idx])
+        #print "\n".join(str(v) for v in problem)
         np.testing.assert_allclose(output, answer, rtol=tolerance)
 
@@ -609 +620 @@
         data = np.loadtxt(data_string.split('\n')).T
         q, delta_qv, _, answer = data
-        answer = romberg_slit_1d(q, delta_qv, self.model, pars)
+        answer = romberg_slit_1d(q, delta_qv, 0., self.model, pars)
         q_calc = slit_extend_q(interpolate(q, 2*np.pi/radius/20),
-                               delta_qv[0], delta_qv[0])
+                               delta_qv[0], 0.)
         resolution = Slit1D(q, width=delta_qv, height=0, q_calc=q_calc)
         output = self.Iq_sphere(pars, resolution)
@@ -629 +640 @@
         form = load_model('ellipsoid', dtype='double')
         q = np.logspace(log10(4e-5),log10(2.5e-2), 68)
-        delta_qv = [0.117]
-        resolution = Slit1D(q, width=delta_qv, height=0)
-        answer = romberg_slit_1d(q, delta_qv, form, pars)
+        width, height = 0.117, 0.
+        resolution = Slit1D(q, width=width, height=height)
+        answer = romberg_slit_1d(q, width, height, form, pars)
         output = resolution.apply(eval_form(resolution.q_calc, form, pars))
         # TODO: 10% is too much error; use better algorithm
+        #print np.max(abs(answer-output)/answer)
         self.compare(q, output, answer, 0.1)
 
@@ -860 +872 @@
 
 def _eval_demo_1d(resolution, title):
+    import sys
     from sasmodels import core
-    from sasmodels.models import cylinder
-    ## or alternatively:
-    # cylinder = core.load_model_definition('cylinder')
-    model = core.load_model(cylinder)
+    name = sys.argv[1] if len(sys.argv) > 1 else 'cylinder'
+
+    if name == 'cylinder':
+        pars = {'length':210, 'radius':500}
+    elif name == 'teubner_strey':
+        pars = {'a2':0.003, 'c1':-1e4, 'c2':1e10, 'background':0.312643}
+    elif name == 'sphere' or name == 'spherepy':
+        pars = TEST_PARS_SLIT_SPHERE
+    elif name == 'ellipsoid':
+        pars = {
+            'scale':0.05,
+            'rpolar':500, 'requatorial':15000,
+            'sld':6, 'solvent_sld': 1,
+            }
+    else:
+        pars = {}
+    defn = core.load_model_definition(name)
+    model = core.load_model(defn)
 
     kernel = core.make_kernel(model, [resolution.q_calc])
-    theory = core.call_kernel(kernel, {'length':210, 'radius':500})
+    theory = core.call_kernel(kernel, pars)
     Iq = resolution.apply(theory)
+
+    if isinstance(resolution, Slit1D):
+        width, height = resolution.width, resolution.height
+        Iq_romb = romberg_slit_1d(resolution.q, width, height, model, pars)
+    else:
+        dq = resolution.q_width
+        Iq_romb = romberg_pinhole_1d(resolution.q, dq, model, pars)
 
     import matplotlib.pyplot as plt
     plt.loglog(resolution.q_calc, theory, label='unsmeared')
     plt.loglog(resolution.q, Iq, label='smeared', hold=True)
+    plt.loglog(resolution.q, Iq_romb, label='romberg smeared', hold=True)
     plt.legend()
     plt.title(title)
@@ -879 +914 @@
 
 def demo_pinhole_1d():
-    q = np.logspace(-3, -1, 400)
+    q = np.logspace(-4, np.log10(0.2), 400)
     q_width = 0.1*q
     resolution = Pinhole1D(q, q_width)
@@ -885 +920 @@
 
 def demo_slit_1d():
-    q = np.logspace(-3, -1, 400)
-    qx_width = 0.005
-    qy_width = 0.0
-    resolution = Slit1D(q, qx_width, qy_width)
-    _eval_demo_1d(resolution, title="0.005 Qx Slit Resolution")
+    q = np.logspace(-4, np.log10(0.2), 400)
+    w = h = 0.
+    #w = 0.005
+    w = 0.0277790
+    #h = 0.00277790
+    resolution = Slit1D(q, w, h)
+    _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w,h))
 
 def demo():
@@ -895 +932 @@
     plt.subplot(121)
     demo_pinhole_1d()
+    #plt.yscale('linear')
     plt.subplot(122)
     demo_slit_1d()
+    #plt.yscale('linear')
     plt.show()
 

test/sasfit/test/utest_fit_line.py

@@ -11 +11 @@
 from sas.models.LineModel import LineModel
 from sas.models.Constant import Constant
+
+from bumps import fitters
+try:
+    from bumps.options import FIT_CONFIG
+    def set_fitter(alg, opts):
+        FIT_CONFIG.selected_id = alg
+        FIT_CONFIG.values[alg].update(opts, monitors=[])
+except:
+    # CRUFT: Bumps changed its handling of fit options around 0.7.5.6
+    def set_fitter(alg, opts):
+        #print "fitting",alg,opts
+        #print "options",fitters.FIT_OPTIONS[alg].__dict__
+        fitters.FIT_DEFAULT = alg
+        fitters.FIT_OPTIONS[alg].options.update(opts, monitors=[])
+
 
 class testFitModule(unittest.TestCase):
@@ -67 +82 @@
 
     def fit_bumps(self, alg, **opts):
-        #Importing the Fit module
-        from bumps import fitters
-        fitters.FIT_DEFAULT = alg
-        fitters.FIT_OPTIONS[alg].options.update(opts)
-        fitters.FIT_OPTIONS[alg].options.update(monitors=[])
-        #print "fitting",alg,opts
-        #kprint "options",fitters.FIT_OPTIONS[alg].__dict__
+        set_fitter(alg, opts)
         self.fit_single(isdream=(alg=='dream'))
 

test/sasfit/test/utest_fit_smeared.py

@@ -13 +13 @@
 from sas.models.CylinderModel import CylinderModel
 from sas.models.SphereModel import SphereModel
+
+
+from bumps import fitters
+try:
+    from bumps.options import FIT_CONFIG
+    def set_fitter(alg, opts):
+        FIT_CONFIG.selected_id= alg
+        FIT_CONFIG.values[alg].update(opts, monitors=[])
+except:
+    # CRUFT: Bumps changed its handling of fit options around 0.7.5.6
+    def set_fitter(alg, opts):
+        #print "fitting",alg,opts
+        #print "options",fitters.FIT_OPTIONS[alg].__dict__
+        fitters.FIT_DEFAULT = alg
+        fitters.FIT_OPTIONS[alg].options.update(opts, monitors=[])
+
 
 class testFitModule(unittest.TestCase):
@@ -57 +73 @@
             Cylinder fit with dispersion
         """
-        alg = 'lm'
-        from bumps import fitters
-        fitters.FIT_DEFAULT = alg
-        #fitters.FIT_OPTIONS[alg].options.update(opts)
-        fitters.FIT_OPTIONS[alg].options.update(monitors=[])
+        set_fitter('lm', {})
         self._dispersion(fitter = Fit())