Index: compare.py
===================================================================
--- compare.py	(revision 29fc2a3903d566dd506b21caa40f5ff7ef0e9d17)
+++ compare.py	(revision 7cf2cfd8695ecfd2684525ab15c14bc00f23616a)
@@ -6,4 +6,5 @@
 from os.path import basename, dirname, join as joinpath
 import glob
+import datetime
 
 import numpy as np
@@ -13,7 +14,8 @@
 
 
-from sasmodels.bumps_model import Model, Experiment, plot_theory, tic
 from sasmodels import core
 from sasmodels import kerneldll
+from sasmodels.data import plot_theory, empty_data1D, empty_data2D
+from sasmodels.direct_model import DirectModel
 from sasmodels.convert import revert_model
 kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True
@@ -22,4 +24,42 @@
 MODELS = [basename(f)[:-3]
           for f in sorted(glob.glob(joinpath(ROOT,"sasmodels","models","[a-zA-Z]*.py")))]
+
+# CRUFT python 2.6
+if not hasattr(datetime.timedelta, 'total_seconds'):
+    def delay(dt):
+        """Return number date-time delta as number seconds"""
+        return dt.days * 86400 + dt.seconds + 1e-6 * dt.microseconds
+else:
+    def delay(dt):
+        """Return number date-time delta as number seconds"""
+        return dt.total_seconds()
+
+
+def tic():
+    """
+    Timer function.
+
+    Use "toc=tic()" to start the clock and "toc()" to measure
+    a time interval.
+    """
+    then = datetime.datetime.now()
+    return lambda: delay(datetime.datetime.now() - then)
+
+
+def set_beam_stop(data, radius, outer=None):
+    """
+    Add a beam stop of the given *radius*.  If *outer*, make an annulus.
+
+    Note: this function does not use the sasview package
+    """
+    if hasattr(data, 'qx_data'):
+        q = np.sqrt(data.qx_data**2 + data.qy_data**2)
+        data.mask = (q < radius)
+        if outer is not None:
+            data.mask |= (q >= outer)
+    else:
+        data.mask = (data.x < radius)
+        if outer is not None:
+            data.mask |= (data.x >= outer)
 
 
@@ -99,7 +139,7 @@
         if p.endswith("_pd"): pars[p] = 0
 
-def eval_sasview(name, pars, data, Nevals=1):
+def eval_sasview(model_definition, pars, data, Nevals=1):
     from sas.models.qsmearing import smear_selection
-    model = sasview_model(name, **pars)
+    model = sasview_model(model_definition, **pars)
     smearer = smear_selection(data, model=model)
     value = None  # silence the linter
@@ -131,28 +171,26 @@
         print "... trying again with single precision"
         model = core.load_model(model_definition, dtype='single', platform="ocl")
-    problem = Experiment(data, Model(model, **pars), cutoff=cutoff)
+    calculator = DirectModel(data, model, cutoff=cutoff)
     value = None  # silence the linter
     toc = tic()
     for _ in range(max(Nevals, 1)):  # force at least one eval
-        #pars['scale'] = np.random.rand()
-        problem.update()
-        value = problem.theory()
+        value = calculator(**pars)
     average_time = toc()*1000./Nevals
     return value, average_time
 
+
 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)
+    calculator = DirectModel(data, model, cutoff=cutoff)
     value = None  # silence the linter
     toc = tic()
     for _ in range(max(Nevals, 1)):  # force at least one eval
-        problem.update()
-        value = problem.theory()
+        value = calculator(**pars)
     average_time = toc()*1000./Nevals
     return value, average_time
 
+
 def make_data(qmax, is2D, Nq=128, resolution=0.0, accuracy='Low', view='log'):
     if is2D:
-        from sasmodels.bumps_model import empty_data2D, set_beam_stop
         data = empty_data2D(np.linspace(-qmax, qmax, Nq), resolution=resolution)
         data.accuracy = accuracy
@@ -160,5 +198,4 @@
         index = ~data.mask
     else:
-        from sasmodels.bumps_model import empty_data1D
         if view == 'log':
             qmax = math.log10(qmax)
@@ -190,9 +227,10 @@
 
     # randomize parameters
-    pars.update(set_pars)
+    #pars.update(set_pars)  # set value before random to control range
     if '-random' in opts or '-random' in opt_values:
         seed = int(opt_values['-random']) if '-random' in opt_values else None
         pars, seed = randomize_model(name, pars, seed=seed)
         print "Randomize using -random=%i"%seed
+    pars.update(set_pars)  # set value after random to control value
 
     # parameter selection
@@ -217,7 +255,11 @@
         print "ctypes t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu))
     elif Ncpu > 0:
-        cpu, cpu_time = eval_sasview(model_definition, pars, data, Ncpu)
-        comp = "sasview"
-        print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu))
+        try:
+            cpu, cpu_time = eval_sasview(model_definition, pars, data, Ncpu)
+            comp = "sasview"
+            #print "ocl/sasview", (ocl-pars['background'])/(cpu-pars['background'])
+            print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu))
+        except ImportError:
+            Ncpu = 0
 
     # Compare, but only if computing both forms
@@ -238,12 +280,12 @@
     if Ncpu > 0:
         if Nocl > 0: plt.subplot(131)
-        plot_theory(data, cpu, view=view)
+        plot_theory(data, cpu, view=view, plot_data=False)
         plt.title("%s t=%.1f ms"%(comp,cpu_time))
-        cbar_title = "log I"
+        #cbar_title = "log I"
     if Nocl > 0:
         if Ncpu > 0: plt.subplot(132)
-        plot_theory(data, ocl, view=view)
+        plot_theory(data, ocl, view=view, plot_data=False)
         plt.title("opencl t=%.1f ms"%ocl_time)
-        cbar_title = "log I"
+        #cbar_title = "log I"
     if Ncpu > 0 and Nocl > 0:
         plt.subplot(133)
@@ -253,10 +295,10 @@
             err,errstr,errview = abs(relerr), "rel err", "log"
         #err,errstr = ocl/cpu,"ratio"
-        plot_theory(data, err, view=errview)
+        plot_theory(data, None, resid=err, view=errview, plot_data=False)
         plt.title("max %s = %.3g"%(errstr, max(abs(err))))
-        cbar_title = errstr if errview=="linear" else "log "+errstr
-    if is2D:
-        h = plt.colorbar()
-        h.ax.set_title(cbar_title)
+        #cbar_title = errstr if errview=="linear" else "log "+errstr
+    #if is2D:
+    #    h = plt.colorbar()
+    #    h.ax.set_title(cbar_title)
 
     if Ncpu > 0 and Nocl > 0 and '-hist' in opts:
@@ -320,7 +362,6 @@
 
 Available models:
-
-    %s
 """
+
 
 NAME_OPTIONS = set([
@@ -342,4 +383,16 @@
     ]
 
+def columnize(L, indent="", width=79):
+    column_width = max(len(w) for w in L) + 1
+    num_columns = (width - len(indent)) // column_width
+    num_rows = len(L) // num_columns
+    L = L + [""] * (num_rows*num_columns - len(L))
+    columns = [L[k*num_rows:(k+1)*num_rows] for k in range(num_columns)]
+    lines = [" ".join("%-*s"%(column_width, entry) for entry in row)
+             for row in zip(*columns)]
+    output = indent + ("\n"+indent).join(lines)
+    return output
+
+
 def get_demo_pars(name):
     import sasmodels.models
@@ -355,5 +408,6 @@
     models = "\n    ".join("%-15s"%v for v in MODELS)
     if len(args) == 0:
-        print(USAGE%models)
+        print(USAGE)
+        print(columnize(MODELS, indent="  "))
         sys.exit(1)
     if args[0] not in MODELS:
Index: mpare_all.sh
===================================================================
--- compare_all.sh	(revision 5134b2c3f75004ca3a032a5e4ca3c5dcc3fe9156)
+++ 	(revision )
@@ -1,15 +1,0 @@
-#!/bin/sh
-
-SASVIEW=$(ls -d ../sasview/build/lib.*)
-PYTHONPATH=../bumps:../periodictable:$SASVIEW
-export PYTHONPATH
-
-echo PYTHONPATH=$PYTHONPATH
-set -x
-
-for dim in 1d100 2d32; do
-    ./compare_many.py all 100 $dim 0 > ${dim}_poly_0.csv
-    ./compare_many.py all 100 $dim 1e-5 > ${dim}_poly_1e-5.csv
-    ./compare_many.py all 100 $dim 1e-3 > ${dim}_poly_1e-3.csv
-    ./compare_many.py all 10000 $dim mono > ${dim}_mono.csv
-done
Index: compare_many.py
===================================================================
--- compare_many.py	(revision ab559434e60c5d32fa4b8815251b1af5f697ed22)
+++ compare_many.py	(revision 7cf2cfd8695ecfd2684525ab15c14bc00f23616a)
@@ -1,11 +1,13 @@
 #!/usr/bin/env python
-
 import sys
+import traceback
 
 import numpy as np
 
+from sasmodels import core
 from sasmodels.kernelcl import environment
 from compare import (MODELS, randomize_model, suppress_pd, eval_sasview,
-                     eval_opencl, eval_ctypes, make_data, get_demo_pars)
+                     eval_opencl, eval_ctypes, make_data, get_demo_pars,
+                     columnize)
 
 def get_stats(target, value, index):
@@ -13,5 +15,5 @@
     relerr = resid/target[index]
     srel = np.argsort(relerr)
-    p90 = int(len(relerr)*0.90)
+    #p90 = int(len(relerr)*0.90)
     p95 = int(len(relerr)*0.95)
     maxrel = np.max(relerr)
@@ -27,4 +29,5 @@
         groups.append(p)
         groups.extend(['']*(len(stats)-1))
+    groups.append("Parameters")
     columns = ['Seed'] + stats*len(parts) +  list(sorted(pars.keys()))
     print(','.join('"%s"'%c for c in groups))
@@ -32,8 +35,24 @@
 
 def compare_instance(name, data, index, N=1, mono=True, cutoff=1e-5):
+    model_definition = core.load_model_definition(name)
     pars = get_demo_pars(name)
     header = '\n"Model","%s","Count","%d"'%(name, N)
     if not mono: header += ',"Cutoff",%g'%(cutoff,)
     print(header)
+
+    # Stuff the failure flag into a mutable object so we can update it from
+    # within the nested function.  Note that the nested function uses "pars"
+    # which is dynamically scoped, not lexically scoped in this context.  That
+    # is, pars is replaced each time in the loop, so don't assume that it is
+    # the default values defined above.
+    def trymodel(fn, *args, **kw):
+        try:
+            result, _ = fn(model_definition, pars, data, *args, **kw)
+        except:
+            result = np.NaN
+            traceback.print_exc()
+        return result
+
+    num_good = 0
     first = True
     for _ in range(N):
@@ -41,46 +60,66 @@
         if mono: suppress_pd(pars)
 
-        target, _ = eval_sasview(name, pars, data)
+        # Force parameter constraints on a per-model basis.
+        if name in ('teubner_strey','broad_peak'):
+            pars['scale'] = 1.0
+        #if name == 'parallelepiped':
+        #    pars['a_side'],pars['b_side'],pars['c_side'] = \
+        #        sorted([pars['a_side'],pars['b_side'],pars['c_side']])
 
-        env = environment()
-        gpu_single_value,_ = eval_opencl(name, pars, data, dtype='single', cutoff=cutoff)
-        gpu_single = get_stats(target, gpu_single_value, index)
-        if env.has_double:
-            gpu_double_value,_ = eval_opencl(name, pars, data, dtype='double', cutoff=cutoff)
-            gpu_double = get_stats(target, gpu_double_value, index)
+
+        good = True
+        labels = []
+        columns = []
+        if 1:
+            sasview_value = trymodel(eval_sasview)
+        if 0:
+            gpu_single_value = trymodel(eval_opencl, dtype='single', cutoff=cutoff)
+            stats = get_stats(sasview_value, gpu_single_value, index)
+            columns.extend(stats)
+            labels.append('GPU single')
+            good = good and (stats[0] < 1e-14)
+        if 0 and environment().has_double:
+            gpu_double_value = trymodel(eval_opencl, dtype='double', cutoff=cutoff)
+            stats = get_stats(sasview_value, gpu_double_value, index)
+            columns.extend(stats)
+            labels.append('GPU double')
+            good = good and (stats[0] < 1e-14)
+        if 1:
+            cpu_double_value = trymodel(eval_ctypes, dtype='double', cutoff=cutoff)
+            stats = get_stats(sasview_value, cpu_double_value, index)
+            columns.extend(stats)
+            labels.append('CPU double')
+            good = good and (stats[0] < 1e-14)
+        if 0:
+            stats = get_stats(cpu_double_value, gpu_single_value, index)
+            columns.extend(stats)
+            labels.append('single/double')
+            good = good and (stats[0] < 1e-14)
+
+        columns += [v for _,v in sorted(pars.items())]
+        if first:
+            print_column_headers(pars, labels)
+            first = False
+        if good:
+            num_good += 1
         else:
-            gpu_double = [0]*len(gpu_single)
-        cpu_double_value,_ =  eval_ctypes(name, pars, data, dtype='double', cutoff=cutoff)
-        cpu_double = get_stats(target, cpu_double_value, index)
-        single_double = get_stats(cpu_double_value, gpu_single_value, index)
+            print(("%d,"%seed)+','.join("%g"%v for v in columns))
+    print '"%d/%d good"'%(num_good, N)
 
-        values = (list(gpu_single) + list(gpu_double) + list(cpu_double)
-                  + list(single_double) + [v for _,v in sorted(pars.items())])
-        if gpu_single[0] > 5e-5:
-            if first:
-                print_column_headers(pars,'GPU single|GPU double|CPU double|single/double'.split('|'))
-                first = False
-            print(("%d,"%seed)+','.join("%g"%v for v in values))
 
-def main():
-    try:
-        model = sys.argv[1]
-        assert (model in MODELS) or (model == "all")
-        count = int(sys.argv[2])
-        is2D = sys.argv[3].startswith('2d')
-        assert sys.argv[3][1] == 'd'
-        Nq = int(sys.argv[3][2:])
-        mono = sys.argv[4] == 'mono'
-        cutoff = float(sys.argv[4]) if not mono else 0
-    except:
-        import traceback; traceback.print_exc()
-        models = "\n    ".join("%-7s: %s"%(k,v.__name__.replace('_',' '))
-                               for k,v in sorted(MODELS.items()))
-        print("""\
-usage: compare_many.py MODEL COUNT (1dNQ|2dNQ) (CUTOFF|mono)
+def print_usage():
+    print "usage: compare_many.py MODEL COUNT (1dNQ|2dNQ) (CUTOFF|mono)"
 
-MODEL is the model name of the model, which is one of:
-    %s
-or "all" for all the models in alphabetical order.
+
+def print_models():
+    print(columnize(MODELS, indent="  "))
+
+
+def print_help():
+    print_usage()
+    print("""\
+
+MODEL is the model name of the model or "all" for all the models
+in alphabetical order.
 
 COUNT is the number of randomly generated parameter sets to try. A value
@@ -98,5 +137,28 @@
 choice of polydisperse parameters, and the number of points in the distribution
 is set in compare.py defaults for each model.
-"""%(models,))
+
+Available models:
+""")
+    print_models()
+
+def main():
+    if len(sys.argv) == 1:
+        print_help()
+        sys.exit(1)
+
+    model = sys.argv[1]
+    if not (model in MODELS) and (model != "all"):
+        print 'Bad model %s.  Use "all" or one of:'
+        print_models()
+        sys.exit(1)
+    try:
+        count = int(sys.argv[2])
+        is2D = sys.argv[3].startswith('2d')
+        assert sys.argv[3][1] == 'd'
+        Nq = int(sys.argv[3][2:])
+        mono = sys.argv[4] == 'mono'
+        cutoff = float(sys.argv[4]) if not mono else 0
+    except:
+        print_usage()
         sys.exit(1)
 
Index: compare_many.sh
===================================================================
--- compare_many.sh	(revision 7cf2cfd8695ecfd2684525ab15c14bc00f23616a)
+++ compare_many.sh	(revision 7cf2cfd8695ecfd2684525ab15c14bc00f23616a)
@@ -0,0 +1,17 @@
+#!/bin/sh
+
+SASVIEW=$(ls -d ../sasview/build/lib.*)
+PYTHONPATH=../bumps:../periodictable:$SASVIEW
+export PYTHONPATH
+
+#echo PYTHONPATH=$PYTHONPATH
+
+./compare_many.py $*
+
+#for dim in 1d100 2d32; do
+#    ./compare_many.py all 100 $dim 0 > ${dim}_poly_0.csv
+#    ./compare_many.py all 100 $dim 1e-5 > ${dim}_poly_1e-5.csv
+#    ./compare_many.py all 100 $dim 1e-3 > ${dim}_poly_1e-3.csv
+#    ./compare_many.py all 100 $dim mono > ${dim}_mono.csv
+#done
+
Index: example/fit.py
===================================================================
--- example/fit.py	(revision 346bc886df88ad4f5c0e3441573497ef436dbf18)
+++ example/fit.py	(revision 7cf2cfd8695ecfd2684525ab15c14bc00f23616a)
@@ -5,14 +5,15 @@
 from bumps.names import *
 from sasmodels.core import load_model
-from sasmodels import bumps_model as sas
+from sasmodels.bumps_model import Model, Experiment
+from sasmodels.data import load_data, set_beam_stop, set_top
 
 """ IMPORT THE DATA USED """
-radial_data = sas.load_data('DEC07267.DAT')
-sas.set_beam_stop(radial_data, 0.00669, outer=0.025)
-sas.set_top(radial_data, -.0185)
+radial_data = load_data('DEC07267.DAT')
+set_beam_stop(radial_data, 0.00669, outer=0.025)
+set_top(radial_data, -.0185)
 
-tan_data = sas.load_data('DEC07266.DAT')
-sas.set_beam_stop(tan_data, 0.00669, outer=0.025)
-sas.set_top(tan_data, -.0185)
+tan_data = load_data('DEC07266.DAT')
+set_beam_stop(tan_data, 0.00669, outer=0.025)
+set_top(tan_data, -.0185)
 #sas.set_half(tan_data, 'right')
 
@@ -28,5 +29,5 @@
 
 if name == "ellipsoid":
-    model = sas.Model(kernel,
+    model = Model(kernel,
         scale=0.08,
         rpolar=15, requatorial=800,
@@ -51,5 +52,5 @@
 
 elif name == "lamellar":
-    model = sas.Model(kernel,
+    model = Model(kernel,
         scale=0.08,
         thickness=19.2946,
@@ -87,5 +88,5 @@
         theta_pd=10, theta_pd_n=50, theta_pd_nsigma=3,
         phi_pd=0, phi_pd_n=10, phi_pd_nsigma=3)
-    model = sas.Model(kernel, **pars)
+    model = Model(kernel, **pars)
 
     # SET THE FITTING PARAMETERS
@@ -102,5 +103,5 @@
 
 elif name == "core_shell_cylinder":
-    model = sas.Model(kernel,
+    model = Model(kernel,
         scale= .031, radius=19.5, thickness=30, length=22,
         core_sld=7.105, shell_sld=.291, solvent_sld=7.105,
@@ -129,5 +130,5 @@
 
 elif name == "capped_cylinder":
-    model = sas.Model(kernel,
+    model = Model(kernel,
         scale=.08, radius=20, cap_radius=40, length=400,
         sld_capcyl=1, sld_solv=6.3,
@@ -144,5 +145,5 @@
 
 elif name == "triaxial_ellipsoid":
-    model = sas.Model(kernel,
+    model = Model(kernel,
         scale=0.08, req_minor=15, req_major=20, rpolar=500,
         sldEll=7.105, solvent_sld=.291,
@@ -170,10 +171,10 @@
 
 model.cutoff = cutoff
-M = sas.Experiment(data=data, model=model)
+M = Experiment(data=data, model=model)
 if section == "both":
-   tan_model = sas.Model(model.kernel, **model.parameters())
+   tan_model = Model(model.kernel, **model.parameters())
    tan_model.phi = model.phi - 90
    tan_model.cutoff = cutoff
-   tan_M = sas.Experiment(data=tan_data, model=tan_model)
+   tan_M = Experiment(data=tan_data, model=tan_model)
    problem = FitProblem([M, tan_M])
 else:
Index: sasmodels/bumps_model.py
===================================================================
--- sasmodels/bumps_model.py	(revision 5d80bbff8d089e442a2552ccfa8ff6bbd9dd4702)
+++ sasmodels/bumps_model.py	(revision 7cf2cfd8695ecfd2684525ab15c14bc00f23616a)
@@ -10,47 +10,17 @@
 how far the polydispersity integral extends.
 
-A variety of helper functions are provided:
-
-    :func:`load_data` loads a sasview data file.
-
-    :func:`empty_data1D` creates an empty dataset, which is useful for plotting
-    a theory function before the data is measured.
-
-    :func:`empty_data2D` creates an empty 2D dataset.
-
-    :func:`set_beam_stop` masks the beam stop from the data.
-
-    :func:`set_half` selects the right or left half of the data, which can
-    be useful for shear measurements which have not been properly corrected
-    for path length and reflections.
-
-    :func:`set_top` cuts the top part off the data.
-
-    :func:`plot_data` plots the data file.
-
-    :func:`plot_theory` plots a calculated result from the model.
-
 """
 
 import datetime
 import warnings
-import traceback
 
 import numpy as np
 
+from bumps.names import Parameter
+
 from . import sesans
-from .resolution import Perfect1D, Pinhole1D, Slit1D
-from .resolution2d import Pinhole2D
-
-# CRUFT python 2.6
-if not hasattr(datetime.timedelta, 'total_seconds'):
-    def delay(dt):
-        """Return number date-time delta as number seconds"""
-        return dt.days * 86400 + dt.seconds + 1e-6 * dt.microseconds
-else:
-    def delay(dt):
-        """Return number date-time delta as number seconds"""
-        return dt.total_seconds()
-
+from . import weights
+from .data import plot_theory
+from .direct_model import DataMixin
 
 # CRUFT: old style bumps wrapper which doesn't separate data and model
@@ -64,290 +34,11 @@
 
 
-
-def tic():
-    """
-    Timer function.
-
-    Use "toc=tic()" to start the clock and "toc()" to measure
-    a time interval.
-    """
-    then = datetime.datetime.now()
-    return lambda: delay(datetime.datetime.now() - then)
-
-
-def load_data(filename):
-    """
-    Load data using a sasview loader.
-    """
-    from sas.dataloader.loader import Loader
-    loader = Loader()
-    data = loader.load(filename)
-    if data is None:
-        raise IOError("Data %r could not be loaded" % filename)
-    return data
-
-def plot_data(data, view='log'):
-    """
-    Plot data loaded by the sasview loader.
-    """
-    if hasattr(data, 'qx_data'):
-        _plot_2d_signal(data, data.data, view=view)
-    else:
-        # Note: kind of weird using the _plot_result1D to plot just the
-        # data, but it handles the masking and graph markup already, so
-        # do not repeat.
-        _plot_result1D(data, None, None, view)
-
-def plot_theory(data, theory, view='log'):
-    if hasattr(data, 'qx_data'):
-        _plot_2d_signal(data, theory, view=view)
-    else:
-        _plot_result1D(data, theory, None, view, include_data=False)
-
-
-def empty_data1D(q, resolution=0.05):
-    """
-    Create empty 1D data using the given *q* as the x value.
-
-    *resolution* dq/q defaults to 5%.
-    """
-
-    from sas.dataloader.data_info import Data1D
-
-    Iq = 100 * np.ones_like(q)
-    dIq = np.sqrt(Iq)
-    data = Data1D(q, Iq, dx=resolution * q, dy=dIq)
-    data.filename = "fake data"
-    data.qmin, data.qmax = q.min(), q.max()
-    data.mask = np.zeros(len(Iq), dtype='bool')
-    return data
-
-
-def empty_data2D(qx, qy=None, resolution=0.05):
-    """
-    Create empty 2D data using the given mesh.
-
-    If *qy* is missing, create a square mesh with *qy=qx*.
-
-    *resolution* dq/q defaults to 5%.
-    """
-    from sas.dataloader.data_info import Data2D, Detector
-
-    if qy is None:
-        qy = qx
-    Qx, Qy = np.meshgrid(qx, qy)
-    Qx, Qy = Qx.flatten(), Qy.flatten()
-    Iq = 100 * np.ones_like(Qx)
-    dIq = np.sqrt(Iq)
-    mask = np.ones(len(Iq), dtype='bool')
-
-    data = Data2D()
-    data.filename = "fake data"
-    data.qx_data = Qx
-    data.qy_data = Qy
-    data.data = Iq
-    data.err_data = dIq
-    data.mask = mask
-    data.qmin = 1e-16
-    data.qmax = np.inf
-
-    # 5% dQ/Q resolution
-    if resolution != 0:
-        # https://www.ncnr.nist.gov/staff/hammouda/distance_learning/chapter_15.pdf
-        # Should have an additional constant which depends on distances and
-        # radii of the aperture, pixel dimensions and wavelength spread
-        # Instead, assume radial dQ/Q is constant, and perpendicular matches
-        # radial (which instead it should be inverse).
-        Q = np.sqrt(Qx**2 + Qy**2)
-        data.dqx_data = resolution * Q
-        data.dqy_data = resolution * Q
-
-    detector = Detector()
-    detector.pixel_size.x = 5 # mm
-    detector.pixel_size.y = 5 # mm
-    detector.distance = 4 # m
-    data.detector.append(detector)
-    data.xbins = qx
-    data.ybins = qy
-    data.source.wavelength = 5 # angstroms
-    data.source.wavelength_unit = "A"
-    data.Q_unit = "1/A"
-    data.I_unit = "1/cm"
-    data.q_data = np.sqrt(Qx ** 2 + Qy ** 2)
-    data.xaxis("Q_x", "A^{-1}")
-    data.yaxis("Q_y", "A^{-1}")
-    data.zaxis("Intensity", r"\text{cm}^{-1}")
-    return data
-
-
-def set_beam_stop(data, radius, outer=None):
-    """
-    Add a beam stop of the given *radius*.  If *outer*, make an annulus.
-    """
-    from sas.dataloader.manipulations import Ringcut
-    if hasattr(data, 'qx_data'):
-        data.mask = Ringcut(0, radius)(data)
-        if outer is not None:
-            data.mask += Ringcut(outer, np.inf)(data)
-    else:
-        data.mask = (data.x >= radius)
-        if outer is not None:
-            data.mask &= (data.x < outer)
-
-
-def set_half(data, half):
-    """
-    Select half of the data, either "right" or "left".
-    """
-    from sas.dataloader.manipulations import Boxcut
-    if half == 'right':
-        data.mask += \
-            Boxcut(x_min=-np.inf, x_max=0.0, y_min=-np.inf, y_max=np.inf)(data)
-    if half == 'left':
-        data.mask += \
-            Boxcut(x_min=0.0, x_max=np.inf, y_min=-np.inf, y_max=np.inf)(data)
-
-
-def set_top(data, cutoff):
-    """
-    Chop the top off the data, above *cutoff*.
-    """
-    from sas.dataloader.manipulations import Boxcut
-    data.mask += \
-        Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=cutoff)(data)
-
-def protect(fn):
-    def wrapper(*args, **kw):
-        try: 
-            return fn(*args, **kw)
-        except:
-            traceback.print_exc()
-    return wrapper
-
-@protect
-def _plot_result1D(data, theory, resid, view, include_data=True):
-    """
-    Plot the data and residuals for 1D data.
-    """
-    import matplotlib.pyplot as plt
-    from numpy.ma import masked_array, masked
-    #print "not a number",sum(np.isnan(data.y))
-    #data.y[data.y<0.05] = 0.5
-    mdata = masked_array(data.y, data.mask)
-    mdata[~np.isfinite(mdata)] = masked
-    if view is 'log':
-        mdata[mdata <= 0] = masked
-
-    scale = data.x**4 if view == 'q4' else 1.0
-    if resid is not None:
-        plt.subplot(121)
-
-    positive = False
-    if include_data:
-        plt.errorbar(data.x, scale*mdata, yerr=data.dy, fmt='.')
-        positive = positive or (mdata>0).any()
-    if theory is not None:
-        mtheory = masked_array(theory, mdata.mask)
-        plt.plot(data.x, scale*mtheory, '-', hold=True)
-        positive = positive or (mtheory>0).any()
-    plt.xscale(view)
-    plt.yscale('linear' if view == 'q4' or not positive else view)
-    plt.xlabel('Q')
-    plt.ylabel('I(Q)')
-    if resid is not None:
-        mresid = masked_array(resid, mdata.mask)
-        plt.subplot(122)
-        plt.plot(data.x, mresid, 'x')
-        plt.ylabel('residuals')
-        plt.xlabel('Q')
-        plt.xscale(view)
-
-# pylint: disable=unused-argument
-@protect
-def _plot_sesans(data, theory, resid, view):
-    import matplotlib.pyplot as plt
-    plt.subplot(121)
-    plt.errorbar(data.x, data.y, yerr=data.dy)
-    plt.plot(data.x, theory, '-', hold=True)
-    plt.xlabel('spin echo length (nm)')
-    plt.ylabel('polarization (P/P0)')
-    plt.subplot(122)
-    plt.plot(data.x, resid, 'x')
-    plt.xlabel('spin echo length (nm)')
-    plt.ylabel('residuals (P/P0)')
-
-@protect
-def _plot_result2D(data, theory, resid, view):
-    """
-    Plot the data and residuals for 2D data.
-    """
-    import matplotlib.pyplot as plt
-    target = data.data[~data.mask]
-    if view == 'log':
-        vmin = min(target[target>0].min(), theory[theory>0].min())
-        vmax = max(target.max(), theory.max())
-    else:
-        vmin = min(target.min(), theory.min())
-        vmax = max(target.max(), theory.max())
-    #print vmin, vmax
-    plt.subplot(131)
-    _plot_2d_signal(data, target, view=view, vmin=vmin, vmax=vmax)
-    plt.title('data')
-    plt.colorbar()
-    plt.subplot(132)
-    _plot_2d_signal(data, theory, view=view, vmin=vmin, vmax=vmax)
-    plt.title('theory')
-    plt.colorbar()
-    plt.subplot(133)
-    _plot_2d_signal(data, resid, view='linear')
-    plt.title('residuals')
-    plt.colorbar()
-
-@protect
-def _plot_2d_signal(data, signal, vmin=None, vmax=None, view='log'):
-    """
-    Plot the target value for the data.  This could be the data itself,
-    the theory calculation, or the residuals.
-
-    *scale* can be 'log' for log scale data, or 'linear'.
-    """
-    import matplotlib.pyplot as plt
-    from numpy.ma import masked_array
-
-    image = np.zeros_like(data.qx_data)
-    image[~data.mask] = signal
-    valid = np.isfinite(image)
-    if view == 'log':
-        valid[valid] = (image[valid] > 0)
-        image[valid] = np.log10(image[valid])
-    elif view == 'q4':
-        image[valid] *= (data.qx_data[valid]**2+data.qy_data[valid]**2)**2
-    image[~valid | data.mask] = 0
-    #plottable = Iq
-    plottable = masked_array(image, ~valid | data.mask)
-    xmin, xmax = min(data.qx_data), max(data.qx_data)
-    ymin, ymax = min(data.qy_data), max(data.qy_data)
-    # TODO: fix vmin, vmax so it is shared for theory/resid
-    vmin = vmax = None
-    try:
-        if vmin is None: vmin = image[valid & ~data.mask].min()
-        if vmax is None: vmax = image[valid & ~data.mask].max()
-    except:
-        vmin, vmax = 0, 1
-    #print vmin,vmax
-    plt.imshow(plottable.reshape(128, 128),
-               interpolation='nearest', aspect=1, origin='upper',
-               extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax)
-
-
 class Model(object):
-    def __init__(self, kernel, **kw):
-        from bumps.names import Parameter
-
-        self.kernel = kernel
-        partype = kernel.info['partype']
+    def __init__(self, model, **kw):
+        self._sasmodel = model
+        partype = model.info['partype']
 
         pars = []
-        for p in kernel.info['parameters']:
+        for p in model.info['parameters']:
             name, default, limits = p[0], p[2], p[3]
             value = kw.pop(name, default)
@@ -378,5 +69,6 @@
         return dict((k, getattr(self, k)) for k in self._parameter_names)
 
-class Experiment(object):
+
+class Experiment(DataMixin):
     """
     Return a bumps wrapper for a SAS model.
@@ -397,71 +89,7 @@
 
         # remember inputs so we can inspect from outside
-        self.data = data
         self.model = model
         self.cutoff = cutoff
-        if hasattr(data, 'lam'):
-            self.data_type = 'sesans'
-        elif hasattr(data, 'qx_data'):
-            self.data_type = 'Iqxy'
-        else:
-            self.data_type = 'Iq'
-
-        # interpret data
-        partype = model.kernel.info['partype']
-        if self.data_type == 'sesans':
-            q = sesans.make_q(data.sample.zacceptance, data.Rmax)
-            self.index = slice(None, None)
-            self.Iq = data.y
-            self.dIq = data.dy
-            #self._theory = np.zeros_like(q)
-            q_vectors = [q]
-        elif self.data_type == 'Iqxy':
-            q = np.sqrt(data.qx_data**2 + data.qy_data**2)
-            qmin = getattr(data, 'qmin', 1e-16)
-            qmax = getattr(data, 'qmax', np.inf)
-            accuracy = getattr(data, 'accuracy', 'Low')
-            self.index = (~data.mask) & (~np.isnan(data.data)) \
-                         & (q >= qmin) & (q <= qmax)
-            self.Iq = data.data[self.index]
-            self.dIq = data.err_data[self.index]
-            self.resolution = Pinhole2D(data=data, index=self.index,
-                                        nsigma=3.0, accuracy=accuracy)
-            #self._theory = np.zeros_like(self.Iq)
-            if not partype['orientation'] and not partype['magnetic']:
-                raise ValueError("not 2D without orientation or magnetic parameters")
-                #qx,qy = self.resolution.q_calc
-                #q_vectors = [np.sqrt(qx**2 + qy**2)]
-            else:
-                q_vectors = self.resolution.q_calc
-        elif self.data_type == 'Iq':
-            self.index = (data.x >= data.qmin) & (data.x <= data.qmax) & ~np.isnan(data.y)
-            self.Iq = data.y[self.index]
-            self.dIq = data.dy[self.index]
-            if getattr(data, 'dx', None) is not None:
-                q, dq = data.x[self.index], data.dx[self.index]
-                if (dq>0).any():
-                    self.resolution = Pinhole1D(q, dq)
-                else:
-                    self.resolution = Perfect1D(q)
-            elif (getattr(data, 'dxl', None) is not None and
-                  getattr(data, 'dxw', None) is not None):
-                q = data.x[self.index]
-                width = data.dxh[self.index]  # Note: dx
-                self.resolution = Slit1D(data.x[self.index],
-                                         width=data.dxh[self.index],
-                                         height=data.dxw[self.index])
-            else:
-                self.resolution = Perfect1D(data.x[self.index])
-
-            #self._theory = np.zeros_like(self.Iq)
-            q_vectors = [self.resolution.q_calc]
-        else:
-            raise ValueError("Unknown data type") # never gets here
-
-        # Remember function inputs so we can delay loading the function and
-        # so we can save/restore state
-        self._fn_inputs = [v for v in q_vectors]
-        self._fn = None
-
+        self._interpret_data(data, model._sasmodel)
         self.update()
 
@@ -483,6 +111,9 @@
     def theory(self):
         if 'theory' not in self._cache:
+            pars = dict((k, v.value) for k,v in self.model.parameters().items())
+            self._cache['theory'] = self._calc_theory(pars, cutoff=self.cutoff)
+            """
             if self._fn is None:
-                q_input = self.model.kernel.make_input(self._fn_inputs)
+                q_input = self.model.kernel.make_input(self._kernel_inputs)
                 self._fn = self.model.kernel(q_input)
 
@@ -495,9 +126,10 @@
                 result = sesans.hankel(self.data.x, self.data.lam * 1e-9,
                                        self.data.sample.thickness / 10,
-                                       self._fn_inputs[0], Iq_calc)
+                                       self._kernel_inputs[0], Iq_calc)
                 self._cache['theory'] = result
             else:
                 Iq = self.resolution.apply(Iq_calc)
                 self._cache['theory'] = Iq
+            """
         return self._cache['theory']
 
@@ -518,48 +150,27 @@
         Plot the data and residuals.
         """
-        data, theory, resid = self.data, self.theory(), self.residuals()
-        if self.data_type == 'Iq':
-            _plot_result1D(data, theory, resid, view)
-        elif self.data_type == 'Iqxy':
-            _plot_result2D(data, theory, resid, view)
-        elif self.data_type == 'sesans':
-            _plot_sesans(data, theory, resid, view)
-        else:
-            raise ValueError("Unknown data type")
+        data, theory, resid = self._data, self.theory(), self.residuals()
+        plot_theory(data, theory, resid, view)
 
     def simulate_data(self, noise=None):
-        theory = self.theory()
-        if noise is not None:
-            self.dIq = theory*noise*0.01
-        dy = self.dIq
-        y = theory + np.random.randn(*dy.shape) * dy
-        self.Iq = y
-        if self.data_type == 'Iq':
-            self.data.dy[self.index] = dy
-            self.data.y[self.index] = y
-        elif self.data_type == 'Iqxy':
-            self.data.data[self.index] = y
-        elif self.data_type == 'sesans':
-            self.data.y[self.index] = y
-        else:
-            raise ValueError("Unknown model")
+        Iq = self.theory()
+        self._set_data(Iq, noise)
 
     def save(self, basename):
         pass
 
-    def _get_weights(self, par):
+    def remove_get_weights(self, name):
         """
         Get parameter dispersion weights
         """
-        from . import weights
-
-        relative = self.model.kernel.info['partype']['pd-rel']
-        limits = self.model.kernel.info['limits']
+        info = self.model.kernel.info
+        relative = name in info['partype']['pd-rel']
+        limits = info['limits'][name]
         disperser, value, npts, width, nsigma = [
-            getattr(self.model, par + ext)
+            getattr(self.model, name + ext)
             for ext in ('_pd_type', '', '_pd_n', '_pd', '_pd_nsigma')]
         value, weight = weights.get_weights(
             disperser, int(npts.value), width.value, nsigma.value,
-            value.value, limits[par], par in relative)
+            value.value, limits, relative)
         return value, weight / np.sum(weight)
 
@@ -567,5 +178,5 @@
         # Can't pickle gpu functions, so instead make them lazy
         state = self.__dict__.copy()
-        state['_fn'] = None
+        state['_kernel'] = None
         return state
 
@@ -573,13 +184,2 @@
         # pylint: disable=attribute-defined-outside-init
         self.__dict__ = state
-
-
-def demo():
-    data = load_data('DEC07086.DAT')
-    set_beam_stop(data, 0.004)
-    plot_data(data)
-    import matplotlib.pyplot as plt; plt.show()
-
-
-if __name__ == "__main__":
-    demo()
Index: sasmodels/core.py
===================================================================
--- sasmodels/core.py	(revision bcd3aa36447d1e388f9ed30473e1874c3ae584d6)
+++ sasmodels/core.py	(revision 7cf2cfd8695ecfd2684525ab15c14bc00f23616a)
@@ -123,5 +123,5 @@
     """
     relative = name in info['partype']['pd-rel']
-    limits = info['limits']
+    limits = info['limits'][name]
     disperser = pars.get(name+'_pd_type', 'gaussian')
     value = pars.get(name, info['defaults'][name])
@@ -130,7 +130,6 @@
     nsigma = pars.get(name+'_pd_nsigma', 3.0)
     value,weight = weights.get_weights(
-        disperser, npts, width, nsigma,
-        value, limits[name], relative)
-    return value,weight/np.sum(weight)
+        disperser, npts, width, nsigma, value, limits, relative)
+    return value, weight / np.sum(weight)
 
 def dispersion_mesh(pars):
Index: sasmodels/direct_model.py
===================================================================
--- sasmodels/direct_model.py	(revision aa4946b68cdcd8686f2c8386cafe8a9e446f0ed2)
+++ sasmodels/direct_model.py	(revision 7cf2cfd8695ecfd2684525ab15c14bc00f23616a)
@@ -5,21 +5,138 @@
 from .core import load_model_definition, load_model, make_kernel
 from .core import call_kernel, call_ER, call_VR
+from . import sesans
+from . import resolution
+from . import resolution2d
 
-class DirectModel:
-    def __init__(self, name, q_vectors=None, dtype='single'):
-        self.model_definition = load_model_definition(name)
-        self.model = load_model(self.model_definition, dtype=dtype)
-        if q_vectors is not None:
-            q_vectors = [np.ascontiguousarray(q,dtype=dtype) for q in q_vectors]
-            self.kernel = make_kernel(self.model, q_vectors)
+class DataMixin(object):
+    """
+    DataMixin captures the common aspects of evaluating a SAS model for a
+    particular data set, including calculating Iq and evaluating the
+    resolution function.  It is used in particular by :class:`DirectModel`,
+    which evaluates a SAS model parameters as key word arguments to the
+    calculator method, and by :class:`bumps_model.Experiment`, which wraps the
+    model and data for use with the Bumps fitting engine.  It is not
+    currently used by :class:`sasview_model.SasviewModel` since this will
+    require a number of changes to SasView before we can do it.
+    """
+    def _interpret_data(self, data, model):
+        self._data = data
+        self._model = model
+
+        # interpret data
+        if hasattr(data, 'lam'):
+            self.data_type = 'sesans'
+        elif hasattr(data, 'qx_data'):
+            self.data_type = 'Iqxy'
+        else:
+            self.data_type = 'Iq'
+
+        partype = model.info['partype']
+
+        if self.data_type == 'sesans':
+            q = sesans.make_q(data.sample.zacceptance, data.Rmax)
+            self.index = slice(None, None)
+            if data.y is not None:
+                self.Iq = data.y
+                self.dIq = data.dy
+            #self._theory = np.zeros_like(q)
+            q_vectors = [q]
+        elif self.data_type == 'Iqxy':
+            if not partype['orientation'] and not partype['magnetic']:
+                raise ValueError("not 2D without orientation or magnetic parameters")
+            q = np.sqrt(data.qx_data**2 + data.qy_data**2)
+            qmin = getattr(data, 'qmin', 1e-16)
+            qmax = getattr(data, 'qmax', np.inf)
+            accuracy = getattr(data, 'accuracy', 'Low')
+            self.index = ~data.mask & (q >= qmin) & (q <= qmax)
+            if data.data is not None:
+                self.index &= ~np.isnan(data.data)
+                self.Iq = data.data[self.index]
+                self.dIq = data.err_data[self.index]
+            self.resolution = resolution2d.Pinhole2D(data=data, index=self.index,
+                                                     nsigma=3.0, accuracy=accuracy)
+            #self._theory = np.zeros_like(self.Iq)
+            q_vectors = self.resolution.q_calc
+        elif self.data_type == 'Iq':
+            self.index = (data.x >= data.qmin) & (data.x <= data.qmax)
+            if data.y is not None:
+                self.index &= ~np.isnan(data.y)
+                self.Iq = data.y[self.index]
+                self.dIq = data.dy[self.index]
+            if getattr(data, 'dx', None) is not None:
+                q, dq = data.x[self.index], data.dx[self.index]
+                if (dq>0).any():
+                    self.resolution = resolution.Pinhole1D(q, dq)
+                else:
+                    self.resolution = resolution.Perfect1D(q)
+            elif (getattr(data, 'dxl', None) is not None and
+                          getattr(data, 'dxw', None) is not None):
+                self.resolution = resolution.Slit1D(data.x[self.index],
+                                                    width=data.dxh[self.index],
+                                                    height=data.dxw[self.index])
+            else:
+                self.resolution = resolution.Perfect1D(data.x[self.index])
+
+            #self._theory = np.zeros_like(self.Iq)
+            q_vectors = [self.resolution.q_calc]
+        else:
+            raise ValueError("Unknown data type") # never gets here
+
+        # Remember function inputs so we can delay loading the function and
+        # so we can save/restore state
+        self._kernel_inputs = [v for v in q_vectors]
+        self._kernel = None
+
+    def _set_data(self, Iq, noise=None):
+        if noise is not None:
+            self.dIq = Iq*noise*0.01
+        dy = self.dIq
+        y = Iq + np.random.randn(*dy.shape) * dy
+        self.Iq = y
+        if self.data_type == 'Iq':
+            self._data.dy[self.index] = dy
+            self._data.y[self.index] = y
+        elif self.data_type == 'Iqxy':
+            self._data.data[self.index] = y
+        elif self.data_type == 'sesans':
+            self._data.y[self.index] = y
+        else:
+            raise ValueError("Unknown model")
+
+    def _calc_theory(self, pars, cutoff=0.0):
+        if self._kernel is None:
+            q_input = self._model.make_input(self._kernel_inputs)
+            self._kernel = self._model(q_input)
+
+        Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff)
+        if self.data_type == 'sesans':
+            result = sesans.hankel(self._data.x, self._data.lam * 1e-9,
+                                   self._data.sample.thickness / 10,
+                                   self._kernel_inputs[0], Iq_calc)
+        else:
+            result = self.resolution.apply(Iq_calc)
+        return result
+
+
+class DirectModel(DataMixin):
+    def __init__(self, data, model, cutoff=1e-5):
+        self.model = model
+        self.cutoff = cutoff
+        self._interpret_data(data, model)
+        self.kernel = make_kernel(self.model, self._kernel_inputs)
     def __call__(self, **pars):
-        return call_kernel(self.kernel, pars)
+        return self._calc_theory(pars, cutoff=self.cutoff)
     def ER(self, **pars):
         return call_ER(self.model.info, pars)
     def VR(self, **pars):
         return call_VR(self.model.info, pars)
+    def simulate_data(self, noise=None, **pars):
+        Iq = self.__call__(**pars)
+        self._set_data(Iq, noise=noise)
 
 def demo():
     import sys
+    from .data import empty_data1D, empty_data2D
+
     if len(sys.argv) < 3:
         print "usage: python -m sasmodels.direct_model modelname (q|qx,qy) par=val ..."
@@ -27,27 +144,33 @@
     model_name = sys.argv[1]
     call = sys.argv[2].upper()
-    if call in ("ER","VR"):
-        q_vectors = None
-    else:
-        values = [float(v) for v in sys.argv[2].split(',')]
+    if call not in ("ER","VR"):
+        try:
+            values = [float(v) for v in call.split(',')]
+        except:
+            values = []
         if len(values) == 1:
-            q = values[0]
-            q_vectors = [[q]]
+            q, = values
+            data = empty_data1D([q])
         elif len(values) == 2:
             qx,qy = values
-            q_vectors = [[qx],[qy]]
+            data = empty_data2D([qx],[qy])
         else:
             print "use q or qx,qy or ER or VR"
             sys.exit(1)
-    model = DirectModel(model_name, q_vectors)
+    else:
+        data = empty_data1D([0.001])  # Data not used in ER/VR
+
+    model_definition = load_model_definition(model_name)
+    model = load_model(model_definition, dtype='single')
+    calculator = DirectModel(data, model)
     pars = dict((k,float(v))
                 for pair in sys.argv[3:]
                 for k,v in [pair.split('=')])
     if call == "ER":
-        print model.ER(**pars)
+        print calculator.ER(**pars)
     elif call == "VR":
-        print model.VR(**pars)
+        print calculator.VR(**pars)
     else:
-        Iq = model(**pars)
+        Iq = calculator(**pars)
         print Iq[0]
 
