-                      r346bc88
+                      rdbde9f8
     def apply(self, theory):
-        print "q calc", self.q_calc
-        print "weights", self.weight_matrix.shape
         return apply_resolution_matrix(self.weight_matrix, theory)
 …
     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()
 …
             \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
 …
     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)
 …
         ])
     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)
 …
 def romberg_slit_1d(q, delta_qv, form, pars):
+def romberg_slit_1d(q, width, height, form, pars):
     """
     Romberg integration for slit resolution.
 …
     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()):
 …
                          (", ".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
 …
     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)
 …
         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)
 …
         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)
 …
 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)
 …
 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)
 …
 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():
 …
     plt.subplot(121)
     demo_pinhole_1d()
+    #plt.yscale('linear')
     plt.subplot(122)
     demo_slit_1d()
+    #plt.yscale('linear')
     plt.show()

Note: See TracChangeset for help on using the changeset viewer.

SasView

Changeset dbde9f8 in sasmodels

Legend:

sasmodels/resolution.py

Download in other formats: