Changeset 786117e in sasmodels

Timestamp:

Nov 12, 2015 11:43:33 AM (8 years ago)

Author:

krzywon

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

77ad412

Parents:

eff3d66 (diff), dbde9f8 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'master' of ​https://github.com/SasView/sasmodels.git

Files:

• compare.py (modified) (6 diffs)
• doc/genapi.py (modified) (1 diff)
• sasmodels/exception.py (modified) (3 diffs)
• sasmodels/kernelcl.py (modified) (1 diff)
• sasmodels/kerneldll.py (modified) (1 diff)
• sasmodels/model_test.py (modified) (4 diffs)
• sasmodels/models/gaussian_peak.py (modified) (3 diffs)
• sasmodels/models/spherepy.py (modified) (2 diffs)
• sasmodels/models/teubner_strey.py (added)
• sasmodels/resolution.py (modified) (15 diffs)
• sasmodels/models/guinier.py (added)
• sasmodels/models/img/image061.jpg (added)
• sasmodels/models/img/image062.jpg (added)
• sasmodels/models/img/image074.jpg (added)
• sasmodels/models/img/image179.jpg (added)
• sasmodels/models/img/pearl_fig.jpg (added)
• sasmodels/models/img/pearl_plot.jpg (added)

compare.py

@@ -100 +100 @@
     model = sasview_model(name, **pars)
     smearer = smear_selection(data, model=model)
+    value = None  # silence the linter
     toc = tic()
-    for _ in range(Nevals):
+    for _ in range(max(Nevals, 1)):  # make sure there is at least one eval
         if hasattr(data, 'qx_data'):
             q = np.sqrt(data.qx_data**2 + data.qy_data**2)
@@ -120 +121 @@
     return value, average_time
 
-def eval_opencl(model_definition, pars, data, dtype='single', Nevals=1, cutoff=0):
+def eval_opencl(model_definition, pars, data, dtype='single', Nevals=1, cutoff=0.):
     try:
         model = core.load_model(model_definition, dtype=dtype, platform="ocl")
@@ -128 +129 @@
         model = core.load_model(model_definition, dtype='single', platform="ocl")
     problem = Experiment(data, Model(model, **pars), cutoff=cutoff)
+    value = None  # silence the linter
     toc = tic()
-    for _ in range(Nevals):
+    for _ in range(max(Nevals, 1)):  # force at least one eval
         #pars['scale'] = np.random.rand()
         problem.update()
@@ -136 +138 @@
     return value, average_time
 
-def eval_ctypes(model_definition, pars, data, dtype='double', Nevals=1, cutoff=0):
+def eval_ctypes(model_definition, pars, data, dtype='double', Nevals=1, cutoff=0.):
     model = core.load_model(model_definition, dtype=dtype, platform="dll")
     problem = Experiment(data, Model(model, **pars), cutoff=cutoff)
+    value = None  # silence the linter
     toc = tic()
-    for _ in range(Nevals):
+    for _ in range(max(Nevals, 1)):  # force at least one eval
         problem.update()
         value = problem.theory()
@@ -224 +227 @@
         #bad = (relerr>1e-4)
         #print relerr[bad],cpu[bad],ocl[bad],data.qx_data[bad],data.qy_data[bad]
-        def stats(label,err):
-            sorted_err = np.sort(abs(err))
-            p50 = int((len(err)-1)*0.50)
-            p98 = int((len(err)-1)*0.98)
-            data = [
-                "max:%.3e"%sorted_err[-1],
-                "median:%.3e"%sorted_err[p50],
-                "98%%:%.3e"%sorted_err[p98],
-                "rms:%.3e"%np.sqrt(np.mean(err**2)),
-                "zero-offset:%+.3e"%np.mean(err),
-                ]
-            print label,"  ".join(data)
-        stats("|ocl-%s|"%comp+(" "*(3+len(comp))), resid)
-        stats("|(ocl-%s)/%s|"%(comp,comp), relerr)
+        _print_stats("|ocl-%s|"%comp+(" "*(3+len(comp))), resid)
+        _print_stats("|(ocl-%s)/%s|"%(comp,comp), relerr)
 
     # Plot if requested
@@ -277 +268 @@
     plt.show()
 
+def _print_stats(label, err):
+    sorted_err = np.sort(abs(err))
+    p50 = int((len(err)-1)*0.50)
+    p98 = int((len(err)-1)*0.98)
+    data = [
+        "max:%.3e"%sorted_err[-1],
+        "median:%.3e"%sorted_err[p50],
+        "98%%:%.3e"%sorted_err[p98],
+        "rms:%.3e"%np.sqrt(np.mean(err**2)),
+        "zero-offset:%+.3e"%np.mean(err),
+        ]
+    print label,"  ".join(data)
+
+
+
 # ===========================================================================
 #

doc/genapi.py

@@ -57 +57 @@
 modules=[
     #('__init__', 'Top level namespace'),
+    #('alignment', 'GPU data alignment [unused]'),
     ('bumps_model', 'Bumps interface'),
+    ('convert', 'Sasview to sasmodel converter'),
+    ('core', 'Model access'),
+    ('direct_model', 'Simple interface'),
+    ('exception', 'Annotate exceptions'),
+    ('generate', 'Model parser'),
+    ('kernelcl', 'OpenCL model evaluator'),
+    ('kerneldll', 'Ctypes model evaluator'),
+    ('kernelpy', 'Python model evaluator'),
+    ('model_test', 'Unit test support'),
+    ('resolution', '1-D resolution functions'),
+    ('resolution2d', '2-D resolution functions'),
     ('sasview_model', 'Sasview interface'),
+    ('sesans', 'SESANS model evaluator'),
     ('weights', 'Distribution functions'),
-    ('gen', 'Model generator'),
-    ('convert', 'Model parameter converter'),
-    ('gpu', 'OpenCL model loader'),
-    ('dll', 'Ctypes model loader'),
     #('transition', 'Model stepper for automatic model selection'),
 ]

sasmodels/exception.py

@@ -2 +2 @@
 Utility to add annotations to python exceptions.
 """
+
+# Platform cruft: WindowsError is only defined on Windows.
+try:
+    WindowsError
+except NameError:
+    class WindowsError(Exception):
+        pass
 
 def annotate_exception(exc, msg):
@@ -19 +26 @@
         KeyError: "hello while accessing 'D'"
     """
+    # Can't extend WindowsError exceptions; instead raise a new exception.
+    # TODO: try to incorporate current stack trace in the raised exception
+    if isinstance(exc, WindowsError):
+        raise OSError(str(exc)+msg)
+
     args = exc.args
     if not args:
@@ -28 +40 @@
         except:
             exc.args = (" ".join((str(exc),msg)),)
-

sasmodels/kernelcl.py

@@ -39 +39 @@
 
 from . import generate
-from .kernelpy import PyModel
 
 F64_DEFS = """\

sasmodels/kerneldll.py

@@ -148 +148 @@
         try:
             self.dll = ct.CDLL(self.dllpath)
-        except WindowsError, exc:
-            # TODO: update annotate_exception so it can handle WindowsError
-            raise RuntimeError(str(exc)+" while loading "+self.dllpath)
         except Exception, exc:
             annotate_exception(exc, "while loading "+self.dllpath)

sasmodels/model_test.py

@@ -68 +68 @@
         model_definition = load_model_definition(model_name)
 
-        smoke_tests = [
-            [{},0.1,None],
-            [{},(0.1,0.1),None],
-            [{},'ER',None],
-            [{},'VR',None],
-            ]
-        tests = smoke_tests + getattr(model_definition, 'tests', [])
-
-        if tests: # in case there are no smoke tests...
-            #print '------'
-            #print 'found tests in', model_name
-            #print '------'
-
-            # if ispy then use the dll loader to call pykernel
-            # don't try to call cl kernel since it will not be
-            # available in some environmentes.
-            ispy = callable(getattr(model_definition,'Iq', None))
-
-            # test using opencl if desired
-            if not ispy and ('opencl' in loaders and HAVE_OPENCL):
+        #print '------'
+        #print 'found tests in', model_name
+        #print '------'
+
+        # if ispy then use the dll loader to call pykernel
+        # don't try to call cl kernel since it will not be
+        # available in some environmentes.
+        is_py = callable(getattr(model_definition,'Iq', None))
+
+        if is_py:  # kernel implemented in python
+            test_name = "Model: %s, Kernel: python"%model_name
+            test_method_name = "test_%s_python" % model_name
+            test = ModelTestCase(test_name, model_definition,
+                                 test_method_name,
+                                 platform="dll",  # so that
+                                 dtype="double")
+            suite.addTest(test)
+        else:   # kernel implemented in C
+            # test using opencl if desired and available
+            if 'opencl' in loaders and HAVE_OPENCL:
                 test_name = "Model: %s, Kernel: OpenCL"%model_name
-                test_method = "test_%s_opencl" % model_name
+                test_method_name = "test_%s_opencl" % model_name
                 test = ModelTestCase(test_name, model_definition,
-                                     tests, test_method,
+                                     test_method_name,
                                      platform="ocl", dtype='single')
                 #print "defining", test_name
@@ -97 +97 @@
 
             # test using dll if desired
-            if ispy or 'dll' in loaders:
+            if 'dll' in loaders:
                 test_name = "Model: %s, Kernel: dll"%model_name
-                test_method = "test_%s_dll" % model_name
+                test_method_name = "test_%s_dll" % model_name
                 test = ModelTestCase(test_name, model_definition,
-                                     tests, test_method,
-                                     platform="dll", dtype="double")
+                                     test_method_name,
+                                     platform="dll",
+                                     dtype="double")
                 suite.addTest(test)
 
@@ -110 +111 @@
 def _hide_model_case_from_nosetests():
     class ModelTestCase(unittest.TestCase):
-        def __init__(self, test_name, definition, tests, test_method,
+        def __init__(self, test_name, definition, test_method_name,
                      platform, dtype):
             self.test_name = test_name
             self.definition = definition
-            self.tests = tests
             self.platform = platform
             self.dtype = dtype
 
-            setattr(self, test_method, self._runTest)
-            unittest.TestCase.__init__(self, test_method)
+            setattr(self, test_method_name, self._runTest)
+            unittest.TestCase.__init__(self, test_method_name)
 
         def _runTest(self):
+            smoke_tests = [
+                [{},0.1,None],
+                [{},(0.1,0.1),None],
+                [{},'ER',None],
+                [{},'VR',None],
+                ]
+
+            tests = getattr(self.definition, 'tests', [])
             try:
                 model = load_model(self.definition, dtype=self.dtype,
                                    platform=self.platform)
-                for test in self.tests:
+                for test in smoke_tests + tests:
                     self._run_one_test(model, test)
+
+                if not tests and self.platform == "dll":
+                    ## Uncomment the following to make forgetting the test
+                    ## values an error.  Only do so for the "dll" tests
+                    ## to reduce noise from both opencl and dll, and because
+                    ## python kernels us
+                    #raise Exception("No test cases provided")
+                    pass
 
             except Exception,exc:
@@ -199 +215 @@
         loaders = ['opencl', 'dll']
     if not models:
-        print >>sys.stderr, "usage: python -m sasmodels.model_test [opencl|dll|opencl_and_dll] model1 model2 ..."
-        print >>sys.stderr, "if model1 is 'all', then all except the remaining models will be tested"
+        print >>sys.stderr, """\
+usage:
+  python -m sasmodels.model_test [opencl|dll|opencl_and_dll] model1 model2 ...
+
+If model1 is 'all', then all except the remaining models will be tested.
+If no compute target is specified, then models will be tested with both opencl
+and dll; the compute target is ignored for pure python models."""
+
         return 1
 

sasmodels/models/gaussian_peak.py

@@ -1 +1 @@
 r"""
+
+Definition
+----------
+
 This model describes a Gaussian shaped peak on a flat background
 
-.. image:: img/image198.PNG
+.. math::
 
-with the peak having height of *I0* centered at *q0* and having a standard deviation of *B*. The FWHM (full-width
-half-maximum) is 2.354 B.
+    I(q) = (\text{scale}) \exp\left[ -\tfrac12 (q-q_0)^2 / \sigma^2 \right]
+        + \text{background}
 
-For 2D data: The 2D scattering intensity is calculated in the same way as 1D, where the *q* vector is defined as
+with the peak having height of *scale* centered at $q_0$ and having a standard
+deviation of $\sigma$. The FWHM (full-width half-maximum) is $2.354 \sigma$.
 
-.. image:: img/image040.gif
+For 2D data, scattering intensity is calculated in the same way as 1D,
+where the $q$ vector is defined as
 
-.. image:: img/image199.jpg
+.. math::
 
-*Figure. 1D plot using the default values (w/500 data points).*
+    q = \sqrt{q_x^2 + q_y^2}
 
-REFERENCE
+
+.. image:: img/gaussian_peak_1d.jpg
+
+    1D plot using the default values (w/500 data points).
+
+Reference
 ---------
 
@@ -27 +38 @@
 description = """
     Model describes a Gaussian shaped peak including a flat background
-    Provide F(q) = scale*exp( -1/2 *[(q-q0)/B]^2 )+ background
+    Provide F(q) = scale*exp( -1/2 *[(q-q0)/sigma]^2 )+ background
 """
 category = "shape-independent"
@@ -33 +44 @@
 #             ["name", "units", default, [lower, upper], "type","description"],
 parameters = [["q0", "1/Ang", 0.05, [-inf, inf], "", "Peak position"],
-              ["sigma", "1/Ang", 0.005, [-inf, inf], "",
+              ["sigma", "1/Ang", 0.005, [0, inf], "",
                "Peak width (standard deviation)"],
              ]

sasmodels/models/spherepy.py

@@ -91 +91 @@
     qr = q * radius
     sn, cn = sin(qr), cos(qr)
-    # FOR VECTORIZED VERSION, UNCOMMENT THE NEXT TWO LINES
+    ## The natural expression for the bessel function is the following:
+    ##     bes = 3 * (sn-qr*cn)/qr**3 if qr>0 else 1
+    ## however, to support vector q values we need to handle the conditional
+    ## as a vector, which we do by first evaluating the full expression
+    ## everywhere, then fixing it up where it is broken.   We should probably
+    ## set numpy to ignore the 0/0 error before we do though...
     bes = 3 * (sn - qr * cn) / qr ** 3 # may be 0/0 but we fix that next line
     bes[qr == 0] = 1
-    # FOR NON VECTORIZED VERSION, UNCOMMENT THE NEXT LINE
-    #bes = 3 * (sn-qr*cn)/qr**3 if qr>0 else 1
     fq = bes * (sld - solvent_sld) * form_volume(radius)
     return 1.0e-4 * fq ** 2
-# FOR VECTORIZED VERSION, UNCOMMENT THE NEXT LINE
-Iq.vectorized = True
+Iq.vectorized = True  # Iq accepts an array of Q values
 
 def Iqxy(qx, qy, sld, solvent_sld, radius):
     return Iq(sqrt(qx ** 2 + qy ** 2), sld, solvent_sld, radius)
-Iqxy.vectorized = True
+Iqxy.vectorized = True  # Iqxy accepts arrays of Qx, Qy values
 
 def sesans(z, sld, solvent_sld, radius):
@@ -119 +121 @@
     g[low] = sqrt(1 - dlow2 / 4.) * (1 + dlow2 / 8.) + dlow2 / 2.*(1 - dlow2 / 16.) * log(dlow / (2. + sqrt(4. - dlow2)))
     return g
-sesans.vectorized = True
+sesans.vectorized = True  # sesans accepts and array of z values
 
 def ER(radius):

sasmodels/resolution.py

@@ -74 +74 @@
 
     def apply(self, theory):
-        print "q calc", self.q_calc
-        print "weights", self.weight_matrix.shape
         return apply_resolution_matrix(self.weight_matrix, theory)
 
@@ -94 +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()
@@ -171 +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
 
 
@@ -214 +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)
@@ -235 +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)
 
 
@@ -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()