Index: compare.sh
===================================================================
--- compare.sh	(revision 5753e4ebc8ad092cc172a8a013acd92a6a19e86a)
+++ compare.sh	(revision ec7e360eb9b8ec7bb8f5006ed71720a9fa7d593b)
@@ -12,3 +12,3 @@
 #SAS_OPENMP=1; export SAS_OPENMP
 
-python -m sasmodels.compare $*
+${PYTHON:-python} -m sasmodels.compare $*
Index: sasmodels/bumps_model.py
===================================================================
--- sasmodels/bumps_model.py	(revision 9404dd3f599fcbe72790aae8a6d176e64f292ea6)
+++ sasmodels/bumps_model.py	(revision ec7e360eb9b8ec7bb8f5006ed71720a9fa7d593b)
@@ -12,11 +12,8 @@
 """
 
-import datetime
 import warnings
 
 import numpy as np
 
-from . import sesans
-from . import weights
 from .data import plot_theory
 from .direct_model import DataMixin
@@ -31,37 +28,48 @@
     return experiment
 
+def create_parameters(model_info, **kwargs):
+    # lazy import; this allows the doc builder and nosetests to run even
+    # when bumps is not on the path.
+    from bumps.names import Parameter
+
+    partype = model_info['partype']
+
+    pars = {}
+    for p in model_info['parameters']:
+        name, default, limits = p[0], p[2], p[3]
+        value = kwargs.pop(name, default)
+        pars[name] = Parameter.default(value, name=name, limits=limits)
+    for name in partype['pd-2d']:
+        for xpart, xdefault, xlimits in [
+            ('_pd', 0., limits),
+            ('_pd_n', 35., (0, 1000)),
+            ('_pd_nsigma', 3., (0, 10)),
+        ]:
+            xname = name + xpart
+            xvalue = kwargs.pop(xname, xdefault)
+            pars[xname] = Parameter.default(xvalue, name=xname, limits=xlimits)
+
+    pd_types = {}
+    for name in partype['pd-2d']:
+        xname = name + '_pd_type'
+        xvalue = kwargs.pop(xname, 'gaussian')
+        pd_types[xname] = xvalue
+
+    if kwargs:  # args not corresponding to parameters
+        raise TypeError("unexpected parameters: %s"
+                        % (", ".join(sorted(kwargs.keys()))))
+
+    return pars, pd_types
 
 class Model(object):
-    def __init__(self, model, **kw):
-        # lazy import; this allows the doc builder and nosetests to run even
-        # when bumps is not on the path.
-        from bumps.names import Parameter
-
+    def __init__(self, model, **kwargs):
         self._sasmodel = model
-        partype = model.info['partype']
-
-        pars = []
-        for p in model.info['parameters']:
-            name, default, limits = p[0], p[2], p[3]
-            value = kw.pop(name, default)
-            setattr(self, name, Parameter.default(value, name=name, limits=limits))
-            pars.append(name)
-        for name in partype['pd-2d']:
-            for xpart, xdefault, xlimits in [
-                ('_pd', 0, limits),
-                ('_pd_n', 35, (0, 1000)),
-                ('_pd_nsigma', 3, (0, 10)),
-                ('_pd_type', 'gaussian', None),
-                ]:
-                xname = name + xpart
-                xvalue = kw.pop(xname, xdefault)
-                if xlimits is not None:
-                    xvalue = Parameter.default(xvalue, name=xname, limits=xlimits)
-                    pars.append(xname)
-                setattr(self, xname, xvalue)
-        self._parameter_names = pars
-        if kw:
-            raise TypeError("unexpected parameters: %s"
-                            % (", ".join(sorted(kw.keys()))))
+        pars, pd_types = create_parameters(model.info, **kwargs)
+        for k,v in pars.items():
+            setattr(self, k, v)
+        for k,v in pd_types.items():
+            setattr(self, k, v)
+        self._parameter_names = list(pars.keys())
+        self._pd_type_names = list(pd_types.keys())
 
     def parameters(self):
@@ -71,4 +79,8 @@
         return dict((k, getattr(self, k)) for k in self._parameter_names)
 
+    def state(self):
+        pars = dict((k, getattr(self, k).value) for k in self._parameter_names)
+        pars.update((k, getattr(self, k)) for k in self._pd_type_names)
+        return pars
 
 class Experiment(DataMixin):
@@ -113,25 +125,6 @@
     def theory(self):
         if 'theory' not in self._cache:
-            pars = dict((k, v.value) for k,v in self.model.parameters().items())
+            pars = self.model.state()
             self._cache['theory'] = self._calc_theory(pars, cutoff=self.cutoff)
-            """
-            if self._fn is None:
-                q_input = self.model.kernel.make_input(self._kernel_inputs)
-                self._fn = self.model.kernel(q_input)
-
-            fixed_pars = [getattr(self.model, p).value for p in self._fn.fixed_pars]
-            pd_pars = [self._get_weights(p) for p in self._fn.pd_pars]
-            #print(fixed_pars,pd_pars)
-            Iq_calc = self._fn(fixed_pars, pd_pars, self.cutoff)
-            #self._theory[:] = self._fn.eval(pars, pd_pars)
-            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)
-                self._cache['theory'] = result
-            else:
-                Iq = self.resolution.apply(Iq_calc)
-                self._cache['theory'] = Iq
-            """
         return self._cache['theory']
 
@@ -162,19 +155,4 @@
         pass
 
-    def remove_get_weights(self, name):
-        """
-        Get parameter dispersion weights
-        """
-        info = self.model.kernel.info
-        relative = name in info['partype']['pd-rel']
-        limits = info['limits'][name]
-        disperser, value, npts, width, nsigma = [
-            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, relative)
-        return value, weight / np.sum(weight)
-
     def __getstate__(self):
         # Can't pickle gpu functions, so instead make them lazy
Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision f4f3919e224d120c63a5e16bba1709c37088b522)
+++ sasmodels/compare.py	(revision ec7e360eb9b8ec7bb8f5006ed71720a9fa7d593b)
@@ -66,11 +66,86 @@
 
 
-def sasview_model(model_definition, **pars):
-    """
-    Load a sasview model given the model name.
-    """
+def parameter_range(p, v):
+    """
+    Choose a parameter range based on parameter name and initial value.
+    """
+    if p.endswith('_pd_n'):
+        return [0, 100]
+    elif p.endswith('_pd_nsigma'):
+        return [0, 5]
+    elif p.endswith('_pd_type'):
+        return v
+    elif any(s in p for s in ('theta','phi','psi')):
+        # orientation in [-180,180], orientation pd in [0,45]
+        if p.endswith('_pd'):
+            return [0,45]
+        else:
+            return [-180, 180]
+    elif 'sld' in p:
+        return [-0.5, 10]
+    elif p.endswith('_pd'):
+        return [0, 1]
+    elif p in ['background', 'scale']:
+        return [0, 1e3]
+    else:
+        return [0, (2*v if v>0 else 1)]
+
+def _randomize_one(p, v):
+    """
+    Randomizing parameter.
+    """
+    if any(p.endswith(s) for s in ('_pd_n','_pd_nsigma','_pd_type')):
+        return v
+    else:
+        return np.random.uniform(*parameter_range(p, v))
+
+def randomize_pars(pars, seed=None):
+    np.random.seed(seed)
+    # Note: the sort guarantees order `of calls to random number generator
+    pars = dict((p,_randomize_one(p,v))
+                for p,v in sorted(pars.items()))
+    return pars
+
+def constrain_pars(model_definition, pars):
+    """
+    Restrict parameters to valid values.
+    """
+    name = model_definition.name
+    if name == 'capped_cylinder' and pars['cap_radius'] < pars['radius']:
+        pars['radius'],pars['cap_radius'] = pars['cap_radius'],pars['radius']
+    if name == 'barbell' and pars['bell_radius'] < pars['radius']:
+        pars['radius'],pars['bell_radius'] = pars['bell_radius'],pars['radius']
+
+    # Limit guinier to an Rg such that Iq > 1e-30 (single precision cutoff)
+    if name == 'guinier':
+        #q_max = 0.2  # mid q maximum
+        q_max = 1.0  # high q maximum
+        rg_max = np.sqrt(90*np.log(10) + 3*np.log(pars['scale']))/q_max
+        pars['rg'] = min(pars['rg'],rg_max)
+
+def parlist(pars):
+    return "\n".join("%s: %s"%(p,v) for p,v in sorted(pars.items()))
+
+def suppress_pd(pars):
+    """
+    Suppress theta_pd for now until the normalization is resolved.
+
+    May also suppress complete polydispersity of the model to test
+    models more quickly.
+    """
+    pars = pars.copy()
+    for p in pars:
+        if p.endswith("_pd"): pars[p] = 0
+    return pars
+
+def eval_sasview(model_definition, data):
+    # importing sas here so that the error message will be that sas failed to
+    # import rather than the more obscure smear_selection not imported error
+    import sas
+    from sas.models.qsmearing import smear_selection
+
     # convert model parameters from sasmodel form to sasview form
     #print("old",sorted(pars.items()))
-    modelname, pars = revert_model(model_definition, pars)
+    modelname, pars = revert_model(model_definition, {})
     #print("new",sorted(pars.items()))
     sas = __import__('sas.models.'+modelname)
@@ -79,240 +154,124 @@
         raise ValueError("could not find model %r in sas.models"%modelname)
     model = ModelClass()
-
-    for k,v in pars.items():
-        parts = k.split('.')  # polydispersity components
-        if len(parts) == 2:
-            model.dispersion[parts[0]][parts[1]] = v
+    smearer = smear_selection(data, model=model)
+
+    if hasattr(data, 'qx_data'):
+        q = np.sqrt(data.qx_data**2 + data.qy_data**2)
+        index = ((~data.mask) & (~np.isnan(data.data))
+                 & (q >= data.qmin) & (q <= data.qmax))
+        if smearer is not None:
+            smearer.model = model  # because smear_selection has a bug
+            smearer.accuracy = data.accuracy
+            smearer.set_index(index)
+            theory = lambda: smearer.get_value()
         else:
-            model.setParam(k, v)
-    return model
-
-def randomize(p, v):
-    """
-    Randomizing parameter.
-
-    Guess the parameter type from name.
-    """
-    if any(p.endswith(s) for s in ('_pd_n','_pd_nsigma','_pd_type')):
-        return v
-    elif any(s in p for s in ('theta','phi','psi')):
-        # orientation in [-180,180], orientation pd in [0,45]
-        if p.endswith('_pd'):
-            return 45*np.random.rand()
-        else:
-            return 360*np.random.rand() - 180
-    elif 'sld' in p:
-        # sld in in [-0.5,10]
-        return 10.5*np.random.rand() - 0.5
-    elif p.endswith('_pd'):
-        # length pd in [0,1]
-        return np.random.rand()
-    else:
-        # values from 0 to 2*x for all other parameters
-        return 2*np.random.rand()*(v if v != 0 else 1)
-
-def randomize_model(pars, seed=None):
-    if seed is None:
-        seed = np.random.randint(1e9)
-    np.random.seed(seed)
-    # Note: the sort guarantees order of calls to random number generator
-    pars = dict((p,randomize(p,v)) for p,v in sorted(pars.items()))
-
-    return pars, seed
-
-def constrain_pars(model_definition, pars):
-    """
-    Restrict parameters to valid values.
-    """
-    name = model_definition.name
-    if name == 'capped_cylinder' and pars['cap_radius'] < pars['radius']:
-        pars['radius'],pars['cap_radius'] = pars['cap_radius'],pars['radius']
-    if name == 'barbell' and pars['bell_radius'] < pars['radius']:
-        pars['radius'],pars['bell_radius'] = pars['bell_radius'],pars['radius']
-
-    # Limit guinier to an Rg such that Iq > 1e-30 (single precision cutoff)
-    if name == 'guinier':
-        #q_max = 0.2  # mid q maximum
-        q_max = 1.0  # high q maximum
-        rg_max = np.sqrt(90*np.log(10) + 3*np.log(pars['scale']))/q_max
-        pars['rg'] = min(pars['rg'],rg_max)
-
-def parlist(pars):
-    return "\n".join("%s: %s"%(p,v) for p,v in sorted(pars.items()))
-
-def suppress_pd(pars):
-    """
-    Suppress theta_pd for now until the normalization is resolved.
-
-    May also suppress complete polydispersity of the model to test
-    models more quickly.
-    """
-    pars = pars.copy()
-    for p in pars:
-        if p.endswith("_pd"): pars[p] = 0
-    return pars
-
-def eval_sasview(model_definition, pars, data, Nevals=1):
-    # importing sas here so that the error message will be that sas failed to
-    # import rather than the more obscure smear_selection not imported error
-    import sas
-    from sas.models.qsmearing import smear_selection
-    model = sasview_model(model_definition, **pars)
-    smearer = smear_selection(data, model=model)
-    value = None  # silence the linter
-    toc = tic()
-    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)
-            index = ((~data.mask) & (~np.isnan(data.data))
-                     & (q >= data.qmin) & (q <= data.qmax))
-            if smearer is not None:
-                smearer.model = model  # because smear_selection has a bug
-                smearer.accuracy = data.accuracy
-                smearer.set_index(index)
-                value = smearer.get_value()
+            theory = lambda: model.evalDistribution([data.qx_data[index], data.qy_data[index]])
+    elif smearer is not None:
+        theory = lambda: smearer(model.evalDistribution(data.x))
+    else:
+        theory = lambda: model.evalDistribution(data.x)
+
+    def calculator(**pars):
+        # paying for parameter conversion each time to keep life simple, if not fast
+        _, pars = revert_model(model_definition, pars)
+        for k,v in pars.items():
+            parts = k.split('.')  # polydispersity components
+            if len(parts) == 2:
+                model.dispersion[parts[0]][parts[1]] = v
             else:
-                value = model.evalDistribution([data.qx_data[index], data.qy_data[index]])
-        else:
-            value = model.evalDistribution(data.x)
-            if smearer is not None:
-                value = smearer(value)
-    average_time = toc()*1000./Nevals
-    return value, average_time
-
-def eval_opencl(model_definition, pars, data, dtype='single', Nevals=1,
-                cutoff=0., fast=False):
+                model.setParam(k, v)
+        return theory()
+
+    calculator.engine = "sasview"
+    return calculator
+
+DTYPE_MAP = {
+    'half': '16',
+    'fast': 'fast',
+    'single': '32',
+    'double': '64',
+    'quad': '128',
+    'f16': '16',
+    'f32': '32',
+    'f64': '64',
+    'longdouble': '128',
+}
+def eval_opencl(model_definition, data, dtype='single', cutoff=0.):
     try:
-        model = core.load_model(model_definition, dtype=dtype,
-                                platform="ocl", fast=fast)
+        model = core.load_model(model_definition, dtype=dtype, platform="ocl")
     except Exception as exc:
         print(exc)
         print("... trying again with single precision")
-        model = core.load_model(model_definition, dtype='single',
-                                platform="ocl", fast=fast)
+        dtype = 'single'
+        model = core.load_model(model_definition, dtype=dtype, platform="ocl")
     calculator = DirectModel(data, model, cutoff=cutoff)
-    # Run the calculator once before starting the timing loop
+    calculator.engine = "OCL%s"%DTYPE_MAP[dtype]
+    return calculator
+
+def eval_ctypes(model_definition, data, dtype='double', cutoff=0.):
+    if dtype=='quad':
+        dtype = 'longdouble'
+    model = core.load_model(model_definition, dtype=dtype, platform="dll")
+    calculator = DirectModel(data, model, cutoff=cutoff)
+    calculator.engine = "OMP%s"%DTYPE_MAP[dtype]
+    return calculator
+
+def time_calculation(calculator, pars, Nevals=1):
+    # initialize the code so time is more accurate
     value = calculator(**suppress_pd(pars))
-    # Now run the timing loop
     toc = tic()
-    for _ in range(max(Nevals, 1)):  # force at least one eval
+    for _ in range(max(Nevals, 1)):  # make sure there is at least one eval
         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")
-    calculator = DirectModel(data, model, cutoff=cutoff)
-    # Run the calculator once before starting the timing loop
-    value = calculator(**suppress_pd(pars))
-    # Now run the timing loop
-    toc = tic()
-    for _ in range(max(Nevals, 1)):  # force at least one eval
-        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:
-        data = empty_data2D(np.linspace(-qmax, qmax, Nq), resolution=resolution)
-        data.accuracy = accuracy
+def make_data(opts):
+    qmax, nq, res = opts['qmax'], opts['nq'], opts['res']
+    if opts['is2d']:
+        data = empty_data2D(np.linspace(-qmax, qmax, nq), resolution=res)
+        data.accuracy = opts['accuracy']
         set_beam_stop(data, 0.004)
         index = ~data.mask
     else:
-        if view == 'log':
+        if opts['view'] == 'log':
             qmax = math.log10(qmax)
-            q = np.logspace(qmax-3, qmax, Nq)
+            q = np.logspace(qmax-3, qmax, nq)
         else:
-            q = np.linspace(0.001*qmax, qmax, Nq)
-        data = empty_data1D(q, resolution=resolution)
+            q = np.linspace(0.001*qmax, qmax, nq)
+        data = empty_data1D(q, resolution=res)
         index = slice(None, None)
     return data, index
 
-def compare(name, pars, Ncomp, Nbase, opts, set_pars):
-    model_definition = core.load_model_definition(name)
-
-    view = ('linear' if '-linear' in opts
-            else 'log' if '-log' in opts
-            else 'q4' if '-q4' in opts
-            else 'log')
-
-    opt_values = dict(split
-                      for s in opts for split in ((s.split('='),))
-                      if len(split) == 2)
-    # Sort out data
-    qmax = (10.0 if '-exq' in opts
-            else 1.0 if '-highq' in opts
-            else 0.2 if '-midq' in opts
-            else 0.05)
-    Nq = int(opt_values.get('-Nq', '128'))
-    res = float(opt_values.get('-res', '0'))
-    accuracy = opt_values.get('-accuracy', 'Low')
-    is2D = "-2d" in opts
-    data, index = make_data(qmax, is2D, Nq, res, accuracy, view=view)
-
-
-    # modelling accuracy is determined by dtype and cutoff
-    dtype = ('longdouble' if '-quad' in opts
-             else 'double' if '-double' in opts
-             else 'half' if '-half' in opts
-             else 'single')
-    cutoff = float(opt_values.get('-cutoff','1e-5'))
-    fast = "-fast" in opts and dtype is 'single'
-
-    # randomize parameters
-    #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(pars, seed=seed)
-        print("Randomize using -random=%i"%seed)
-    pars.update(set_pars)  # set value after random to control value
-    constrain_pars(model_definition, pars)
-    constrain_new_to_old(model_definition, pars)
-
-    # parameter selection
-    if '-mono' in opts:
-        pars = suppress_pd(pars)
-    if '-pars' in opts:
-        print("pars "+str(parlist(pars)))
+def make_engine(model_definition, data, dtype, cutoff):
+    if dtype == 'sasview':
+        return eval_sasview(model_definition, data)
+    elif dtype.endswith('!'):
+        return eval_ctypes(model_definition, data, dtype=dtype[:-1],
+                           cutoff=cutoff)
+    else:
+        return eval_opencl(model_definition, data, dtype=dtype,
+                           cutoff=cutoff)
+
+def compare(opts):
+    Nbase, Ncomp = opts['N1'], opts['N2']
+    pars = opts['pars']
+    data = opts['data']
 
     # Base calculation
-    if 0:
-        from sasmodels.models import sphere as target
-        base_name = target.name
-        base, base_time = eval_ctypes(target, pars, data,
-                dtype='longdouble', cutoff=0., Nevals=Ncomp)
-    elif Nbase > 0 and "-ctypes" in opts and "-sasview" in opts:
+    if Nbase > 0:
+        base = opts['engines'][0]
         try:
-            base, base_time = eval_sasview(model_definition, pars, data, Ncomp)
-            base_name = "sasview"
-            #print("base/sasview", (base-pars['background'])/(comp-pars['background']))
-            print("sasview t=%.1f ms, intensity=%.0f"%(base_time, sum(base)))
-            #print("sasview",comp)
+            base_value, base_time = time_calculation(base, pars, Nbase)
+            print("%s t=%.1f ms, intensity=%.0f"%(base.engine, base_time, sum(base_value)))
         except ImportError:
             traceback.print_exc()
             Nbase = 0
-    elif Nbase > 0:
-        base, base_time = eval_opencl(model_definition, pars, data,
-                dtype=dtype, cutoff=cutoff, Nevals=Nbase, fast=fast)
-        base_name = "ocl"
-        print("opencl t=%.1f ms, intensity=%.0f"%(base_time, sum(base)))
-        #print("base " + base)
-        #print(max(base), min(base))
 
     # Comparison calculation
-    if Ncomp > 0 and "-ctypes" in opts:
-        comp, comp_time = eval_ctypes(model_definition, pars, data,
-                dtype=dtype, cutoff=cutoff, Nevals=Ncomp)
-        comp_name = "ctypes"
-        print("ctypes t=%.1f ms, intensity=%.0f"%(comp_time, sum(comp)))
-    elif Ncomp > 0:
+    if Ncomp > 0:
+        comp = opts['engines'][1]
         try:
-            comp, comp_time = eval_sasview(model_definition, pars, data, Ncomp)
-            comp_name = "sasview"
-            #print("base/sasview", (base-pars['background'])/(comp-pars['background']))
-            print("sasview t=%.1f ms, intensity=%.0f"%(comp_time, sum(comp)))
-            #print("sasview",comp)
+            comp_value, comp_time = time_calculation(comp, pars, Ncomp)
+            print("%s t=%.1f ms, intensity=%.0f"%(comp.engine, comp_time, sum(comp_value)))
         except ImportError:
             traceback.print_exc()
@@ -322,25 +281,24 @@
     if Nbase > 0 and Ncomp > 0:
         #print("speedup %.2g"%(comp_time/base_time))
-        #print("max |base/comp|", max(abs(base/comp)), "%.15g"%max(abs(base)), "%.15g"%max(abs(comp)))
-        #comp *= max(base/comp)
-        resid = (base - comp)
-        relerr = resid/comp
-        #bad = (relerr>1e-4)
-        #print(relerr[bad],comp[bad],base[bad],data.qx_data[bad],data.qy_data[bad])
-        _print_stats("|%s-%s|"%(base_name,comp_name)+(" "*(3+len(comp_name))), resid)
-        _print_stats("|(%s-%s)/%s|"%(base_name,comp_name,comp_name), relerr)
+        #print("max |base/comp|", max(abs(base_value/comp_value)), "%.15g"%max(abs(base_value)), "%.15g"%max(abs(comp_value)))
+        #comp *= max(base_value/comp_value)
+        resid = (base_value - comp_value)
+        relerr = resid/comp_value
+        _print_stats("|%s - %s|"%(base.engine,comp.engine)+(" "*(3+len(comp.engine))), resid)
+        _print_stats("|(%s - %s) / %s|"%(base.engine,comp.engine,comp.engine), relerr)
 
     # Plot if requested
-    if '-noplot' in opts: return
+    if not opts['plot'] and not opts['explore']: return
+    view = opts['view']
     import matplotlib.pyplot as plt
+    if Nbase > 0:
+        if Ncomp > 0: plt.subplot(131)
+        plot_theory(data, base_value, view=view, plot_data=False)
+        plt.title("%s t=%.1f ms"%(base.engine, base_time))
+        #cbar_title = "log I"
     if Ncomp > 0:
-        if Nbase > 0: plt.subplot(131)
-        plot_theory(data, comp, view=view, plot_data=False)
-        plt.title("%s t=%.1f ms"%(comp_name,comp_time))
-        #cbar_title = "log I"
-    if Nbase > 0:
-        if Ncomp > 0: plt.subplot(132)
-        plot_theory(data, base, view=view, plot_data=False)
-        plt.title("%s t=%.1f ms"%(base_name,base_time))
+        if Nbase > 0: plt.subplot(132)
+        plot_theory(data, comp_value, view=view, plot_data=False)
+        plt.title("%s t=%.1f ms"%(comp.engine,comp_time))
         #cbar_title = "log I"
     if Ncomp > 0 and Nbase > 0:
@@ -363,9 +321,10 @@
         v[v==0] = 0.5*np.min(np.abs(v[v!=0]))
         plt.hist(np.log10(np.abs(v)), normed=1, bins=50);
-        plt.xlabel('log10(err), err = | F(q) single - F(q) double| / | F(q) double |');
+        plt.xlabel('log10(err), err = |(%s - %s) / %s|'%(base.engine, comp.engine, comp.engine));
         plt.ylabel('P(err)')
-        plt.title('Comparison of single and double precision models for %s'%name)
-
-    plt.show()
+        plt.title('Distribution of relative error between calculation engines')
+
+    if not opts['explore']:
+        plt.show()
 
 def _print_stats(label, err):
@@ -387,33 +346,39 @@
 #
 USAGE="""
-usage: compare.py model [Nopencl] [Nsasview] [options...] [key=val]
+usage: compare.py model N1 N2 [options...] [key=val]
 
 Compare the speed and value for a model between the SasView original and the
-OpenCL rewrite.
+sasmodels rewrite.
 
 model is the name of the model to compare (see below).
-Nopencl is the number of times to run the OpenCL model (default=5)
-Nsasview is the number of times to run the Sasview model (default=1)
+N1 is the number of times to run sasmodels (default=1).
+N2 is the number times to run sasview (default=1).
 
 Options (* for default):
 
     -plot*/-noplot plots or suppress the plot of the model
-    -half/-single*/-double/-quad/-fast sets the calculation precision
     -lowq*/-midq/-highq/-exq use q values up to 0.05, 0.2, 1.0, 10.0
-    -Nq=128 sets the number of Q points in the data set
+    -nq=128 sets the number of Q points in the data set
     -1d*/-2d computes 1d or 2d data
     -preset*/-random[=seed] preset or random parameters
     -mono/-poly* force monodisperse/polydisperse
-    -ctypes/-sasview* selects gpu:cpu, gpu:sasview, or sasview:cpu if both
     -cutoff=1e-5* cutoff value for including a point in polydispersity
     -pars/-nopars* prints the parameter set or not
     -abs/-rel* plot relative or absolute error
-    -linear/-log/-q4 intensity scaling
+    -linear/-log*/-q4 intensity scaling
     -hist/-nohist* plot histogram of relative error
     -res=0 sets the resolution width dQ/Q if calculating with resolution
     -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh
-
-Key=value pairs allow you to set specific values to any of the model
-parameters.
+    -edit starts the parameter explorer
+
+Any two calculation engines can be selected for comparison:
+
+    -single/-double/-half/-fast sets an OpenCL calculation engine
+    -single!/-double!/-quad! sets an OpenMP calculation engine
+    -sasview sets the sasview calculation engine
+
+The default is -single -sasview
+
+Key=value pairs allow you to set specific values for the model parameters.
 
 Available models:
@@ -423,18 +388,19 @@
 NAME_OPTIONS = set([
     'plot', 'noplot',
-    'half', 'single', 'double', 'quad', 'fast',
+    'half', 'fast', 'single', 'double',
+    'single!', 'double!', 'quad!', 'sasview',
     'lowq', 'midq', 'highq', 'exq',
     '2d', '1d',
     'preset', 'random',
     'poly', 'mono',
-    'sasview', 'ctypes',
     'nopars', 'pars',
     'rel', 'abs',
     'linear', 'log', 'q4',
     'hist', 'nohist',
+    'edit',
     ])
 VALUE_OPTIONS = [
     # Note: random is both a name option and a value option
-    'cutoff', 'random', 'Nq', 'res', 'accuracy',
+    'cutoff', 'random', 'nq', 'res', 'accuracy',
     ]
 
@@ -453,11 +419,22 @@
 def get_demo_pars(model_definition):
     info = generate.make_info(model_definition)
+    # Get the default values for the parameters
     pars = dict((p[0],p[2]) for p in info['parameters'])
+
+    # Fill in default values for the polydispersity parameters
+    for p in info['parameters']:
+        if p[4] in ('volume', 'orientation'):
+            pars[p[0]+'_pd'] = 0.0
+            pars[p[0]+'_pd_n'] = 0
+            pars[p[0]+'_pd_nsigma'] = 3.0
+            pars[p[0]+'_pd_type'] = "gaussian"
+
+    # Plug in values given in demo
     pars.update(info['demo'])
     return pars
 
-def main():
-    opts = [arg for arg in sys.argv[1:] if arg.startswith('-')]
-    popts = [arg for arg in sys.argv[1:] if not arg.startswith('-') and '=' in arg]
+def parse_opts():
+    flags = [arg for arg in sys.argv[1:] if arg.startswith('-')]
+    values = [arg for arg in sys.argv[1:] if not arg.startswith('-') and '=' in arg]
     args = [arg for arg in sys.argv[1:] if not arg.startswith('-') and '=' not in arg]
     models = "\n    ".join("%-15s"%v for v in MODELS)
@@ -472,5 +449,5 @@
         print("expected parameters: model Nopencl Nsasview")
 
-    invalid = [o[1:] for o in opts
+    invalid = [o[1:] for o in flags
                if o[1:] not in NAME_OPTIONS
                   and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)]
@@ -479,31 +456,202 @@
         sys.exit(1)
 
+
+    # Interpret the flags
+    opts = {
+        'plot'      : True,
+        'view'      : 'log',
+        'is2d'      : False,
+        'qmax'      : 0.05,
+        'nq'        : 128,
+        'res'       : 0.0,
+        'accuracy'  : 'Low',
+        'cutoff'    : 1e-5,
+        'seed'      : -1,  # default to preset
+        'mono'      : False,
+        'show_pars' : False,
+        'show_hist' : False,
+        'rel_err'   : True,
+        'explore'   : False,
+    }
+    engines = []
+    for arg in flags:
+        if arg == '-noplot':    opts['plot'] = False
+        elif arg == '-plot':    opts['plot'] = True
+        elif arg == '-linear':  opts['view'] = 'linear'
+        elif arg == '-log':     opts['view'] = 'log'
+        elif arg == '-q4':      opts['view'] = 'q4'
+        elif arg == '-1d':      opts['is2d'] = False
+        elif arg == '-2d':      opts['is2d'] = True
+        elif arg == '-exq':     opts['qmax'] = 10.0
+        elif arg == '-highq':   opts['qmax'] = 1.0
+        elif arg == '-midq':    opts['qmax'] = 0.2
+        elif arg == '-loq':     opts['qmax'] = 0.05
+        elif arg.startswith('-nq='):       opts['nq'] = int(arg[4:])
+        elif arg.startswith('-res='):      opts['res'] = float(arg[5:])
+        elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:]
+        elif arg.startswith('-cutoff='):   opts['cutoff'] = float(arg[8:])
+        elif arg.startswith('-random='):   opts['seed'] = int(arg[8:])
+        elif arg == '-random':  opts['seed'] = np.random.randint(1e6)
+        elif arg == '-preset':  opts['seed'] = -1
+        elif arg == '-mono':    opts['mono'] = True
+        elif arg == '-poly':    opts['mono'] = False
+        elif arg == '-pars':    opts['show_pars'] = True
+        elif arg == '-nopars':  opts['show_pars'] = False
+        elif arg == '-hist':    opts['show_hist'] = True
+        elif arg == '-nohist':  opts['show_hist'] = False
+        elif arg == '-rel':     opts['rel_err'] = True
+        elif arg == '-abs':     opts['rel_err'] = False
+        elif arg == '-half':    engines.append(arg[1:])
+        elif arg == '-fast':    engines.append(arg[1:])
+        elif arg == '-single':  engines.append(arg[1:])
+        elif arg == '-double':  engines.append(arg[1:])
+        elif arg == '-single!': engines.append(arg[1:])
+        elif arg == '-double!': engines.append(arg[1:])
+        elif arg == '-quad!':   engines.append(arg[1:])
+        elif arg == '-sasview': engines.append(arg[1:])
+        elif arg == '-edit':    opts['explore'] = True
+
+    if len(engines) == 0:
+        engines.extend(['single','sasview'])
+    elif len(engines) == 1:
+        if engines[0][0] != 'sasview':
+            engines.append('sasview')
+        else:
+            engines.append('single')
+    elif len(engines) > 2:
+        del engines[2:]
+
+    name = args[0]
+    model_definition = core.load_model_definition(name)
+
+    N1 = int(args[1]) if len(args) > 1 else 1
+    N2 = int(args[2]) if len(args) > 2 else 1
+
     # Get demo parameters from model definition, or use default parameters
     # if model does not define demo parameters
-    name = args[0]
-    model_definition = core.load_model_definition(name)
     pars = get_demo_pars(model_definition)
 
-    Ncomp = int(args[1]) if len(args) > 1 else 5
-    Nbase = int(args[2]) if len(args) > 2 else 1
-
-    # Fill in default polydispersity parameters
-    pds = set(p.split('_pd')[0] for p in pars if p.endswith('_pd'))
-    for p in pds:
-        if p+"_pd_nsigma" not in pars: pars[p+"_pd_nsigma"] = 3
-        if p+"_pd_type" not in pars: pars[p+"_pd_type"] = "gaussian"
-
     # Fill in parameters given on the command line
-    set_pars = {}
-    for arg in popts:
+    presets = {}
+    for arg in values:
         k,v = arg.split('=',1)
         if k not in pars:
-            # extract base name without distribution
+            # extract base name without polydispersity info
             s = set(p.split('_pd')[0] for p in pars)
             print("%r invalid; parameters are: %s"%(k,", ".join(sorted(s))))
             sys.exit(1)
-        set_pars[k] = float(v) if not v.endswith('type') else v
-
-    compare(name, pars, Ncomp, Nbase, opts, set_pars)
+        presets[k] = float(v) if not k.endswith('type') else v
+
+    # randomize parameters
+    #pars.update(set_pars)  # set value before random to control range
+    if opts['seed'] > -1:
+        pars = randomize_pars(pars, seed=opts['seed'])
+        print("Randomize using -random=%i"%opts['seed'])
+    pars.update(presets)  # set value after random to control value
+    constrain_pars(model_definition, pars)
+    constrain_new_to_old(model_definition, pars)
+    if opts['mono']:
+        pars = suppress_pd(pars)
+    if opts['show_pars']:
+        print("pars " + str(parlist(pars)))
+
+    # Create the computational engines
+    data, _index = make_data(opts)
+    if N1:
+        base = make_engine(model_definition, data, engines[0], opts['cutoff'])
+    else:
+        base = None
+    if N2:
+        comp = make_engine(model_definition, data, engines[1], opts['cutoff'])
+    else:
+        comp = None
+
+    # Remember it all
+    opts.update({
+        'name'      : name,
+        'def'       : model_definition,
+        'N1'        : N1,
+        'N2'        : N2,
+        'presets'   : presets,
+        'pars'      : pars,
+        'data'      : data,
+        'engines'   : [base, comp],
+    })
+
+    return opts
+
+def main():
+    opts = parse_opts()
+    if opts['explore']:
+        explore(opts)
+    else:
+        compare(opts)
+
+def explore(opts):
+    import wx
+    from bumps.names import FitProblem
+    from bumps.gui.app_frame import AppFrame
+
+    problem = FitProblem(Explore(opts))
+    isMac = "cocoa" in wx.version()
+    app = wx.App()
+    frame = AppFrame(parent=None, title="explore")
+    if not isMac: frame.Show()
+    frame.panel.set_model(model=problem)
+    frame.panel.Layout()
+    frame.panel.aui.Split(0, wx.TOP)
+    if isMac: frame.Show()
+    app.MainLoop()
+
+class Explore(object):
+    """
+    Return a bumps wrapper for a SAS model comparison.
+    """
+    def __init__(self, opts):
+        from bumps.cli import config_matplotlib
+        import bumps_model
+        config_matplotlib()
+        self.opts = opts
+        info = generate.make_info(opts['def'])
+        pars, pd_types = bumps_model.create_parameters(info, **opts['pars'])
+        if not opts['is2d']:
+            active = [base + ext
+                      for base in info['partype']['pd-1d']
+                      for ext in ['','_pd','_pd_n','_pd_nsigma']]
+            active.extend(info['partype']['fixed-1d'])
+            for k in active:
+                v = pars[k]
+                v.range(*parameter_range(k, v.value))
+        else:
+            for k, v in self.pars.items():
+                v.range(*parameter_range(k, v.value))
+
+        self.pars = pars
+        self.pd_types = pd_types
+
+    def numpoints(self):
+        """
+        Return the number of points
+        """
+        return len(self.pars) + 1  # so dof is 1
+
+    def parameters(self):
+        """
+        Return a dictionary of parameters
+        """
+        return self.pars
+
+    def nllf(self):
+        return 0.  # No nllf
+
+    def plot(self, view='log'):
+        """
+        Plot the data and residuals.
+        """
+        pars = dict((k, v.value) for k,v in self.pars.items())
+        pars.update(self.pd_types)
+        self.opts['pars'] = pars
+        compare(self.opts)
+
 
 if __name__ == "__main__":
Index: sasmodels/compare_many.py
===================================================================
--- sasmodels/compare_many.py	(revision f4f3919e224d120c63a5e16bba1709c37088b522)
+++ sasmodels/compare_many.py	(revision ec7e360eb9b8ec7bb8f5006ed71720a9fa7d593b)
@@ -7,7 +7,8 @@
 from . import core
 from .kernelcl import environment
-from .compare import (MODELS, randomize_model, suppress_pd, eval_sasview,
+from .compare import (MODELS, randomize_pars, suppress_pd, eval_sasview,
                       eval_opencl, eval_ctypes, make_data, get_demo_pars,
-                      columnize, constrain_pars, constrain_new_to_old)
+                      columnize, constrain_pars, constrain_new_to_old,
+                      make_engine)
 
 def calc_stats(target, value, index):
@@ -34,6 +35,16 @@
     print(','.join('"%s"'%c for c in columns))
 
+PRECISION = {
+    'fast': 1e-3,
+    'half': 1e-3,
+    'single': 5e-5,
+    'double': 5e-14,
+    'single!': 5e-5,
+    'double!': 5e-14,
+    'quad!': 5e-18,
+    'sasview': 5e-14,
+}
 def compare_instance(name, data, index, N=1, mono=True, cutoff=1e-5,
-                     precision='double'):
+                     base='sasview', comp='double'):
     model_definition = core.load_model_definition(name)
     pars = get_demo_pars(model_definition)
@@ -47,7 +58,7 @@
     # to allow them to update values in the current scope since nonlocal
     # declarations are not available in python 2.7.
-    def try_model(fn, *args, **kw):
+    def try_model(fn, pars):
         try:
-            result, _ = fn(model_definition, pars_i, data, *args, **kw)
+            result = fn(**pars)
         except KeyboardInterrupt:
             raise
@@ -60,10 +71,16 @@
                 result = np.NaN*data.x
         return result
-    def check_model(label, target, value, acceptable):
-        stats = calc_stats(target, value, index)
-        columns.extend(stats)
-        labels.append('GPU single')
+    def check_model(pars):
+        base_value = try_model(calc_base, pars)
+        comp_value = try_model(calc_comp, pars)
+        stats = calc_stats(base_value, comp_value, index)
         max_diff[0] = max(max_diff[0], stats[0])
-        good[0] = good[0] and (stats[0] < acceptable)
+        good[0] = good[0] and (stats[0] < expected)
+        return list(stats)
+
+
+    calc_base = make_engine(model_definition, data, base, cutoff)
+    calc_comp = make_engine(model_definition, data, comp, cutoff)
+    expected = max(PRECISION[base], PRECISION[comp])
 
     num_good = 0
@@ -72,5 +89,6 @@
     for k in range(N):
         print("%s %d"%(name, k))
-        pars_i, seed = randomize_model(pars)
+        seed = np.random.randint(1e6)
+        pars_i = randomize_pars(pars, seed)
         constrain_pars(model_definition, pars_i)
         constrain_new_to_old(model_definition, pars_i)
@@ -79,30 +97,8 @@
 
         good = [True]
-        labels = []
-        columns = []
-        target = try_model(eval_sasview)
-        #target = try_model(eval_ctypes, dtype='double', cutoff=0.)
-        #target = try_model(eval_ctypes, dtype='longdouble', cutoff=0.)
-        if precision == 'single':
-            value = try_model(eval_opencl, dtype='single', cutoff=cutoff)
-            check_model('GPU single', target, value, 5e-5)
-            single_value = value  # remember for single/double comparison
-        elif precision == 'double':
-            if environment().has_type('double'):
-                label = 'GPU double'
-                value = try_model(eval_opencl, dtype='double', cutoff=cutoff)
-            else:
-                label = 'CPU double'
-                value = try_model(eval_ctypes, dtype='double', cutoff=cutoff)
-            check_model(label, target, value, 5e-14)
-            double_value = value  # remember for single/double comparison
-        elif precision == 'quad':
-            value = try_model(eval_opencl, dtype='longdouble', cutoff=cutoff)
-            check_model('CPU quad', target, value, 5e-14)
-        if 0:
-            check_model('single/double', double_value, single_value, 5e-5)
-
+        columns = check_model(pars_i)
         columns += [v for _,v in sorted(pars_i.items())]
         if first:
+            labels = [" vs. ".join((calc_base.engine, calc_comp.engine))]
             print_column_headers(pars_i, labels)
             first = False
@@ -110,5 +106,5 @@
             num_good += 1
         else:
-            print(("%d,"%seed)+','.join("%g"%v for v in columns))
+            print(("%d,"%seed)+','.join("%s"%v for v in columns))
     print('"good","%d/%d","max diff",%g'%(num_good, N, max_diff[0]))
 
@@ -144,5 +140,6 @@
 is set in compare.py defaults for each model.
 
-PRECISION is the floating point precision to use for comparisons.
+PRECISION is the floating point precision to use for comparisons.  If two
+precisions are given, then compare one to the other, ignoring sasview.
 
 Available models:
@@ -151,5 +148,5 @@
 
 def main():
-    if len(sys.argv) != 6:
+    if len(sys.argv) not in (6,7):
         print_help()
         sys.exit(1)
@@ -167,5 +164,6 @@
         mono = sys.argv[4] == 'mono'
         cutoff = float(sys.argv[4]) if not mono else 0
-        precision = sys.argv[5]
+        base = sys.argv[5]
+        comp = sys.argv[6] if len(sys.argv) > 6 else "sasview"
     except:
         traceback.print_exc()
@@ -173,9 +171,10 @@
         sys.exit(1)
 
-    data, index = make_data(qmax=1.0, is2D=is2D, Nq=Nq)
+    data, index = make_data({'qmax':1.0, 'is2d':is2D, 'nq':Nq, 'res':0.,
+                              'accuracy': 'Low', 'view':'log'})
     model_list = [model] if model != "all" else MODELS
     for model in model_list:
         compare_instance(model, data, index, N=count, mono=mono,
-                         cutoff=cutoff, precision=precision)
+                         cutoff=cutoff, base=base, comp=comp)
 
 if __name__ == "__main__":
