Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision 630156b6317631be1fced3257e8862a36cf47eaf)
+++ sasmodels/compare.py	(revision 0bdddc2288f8fc5bc5e416af84e456e811c3b733)
@@ -264,5 +264,5 @@
 
 
-def _randomize_one(model_info, p, v):
+def _randomize_one(model_info, name, value):
     # type: (ModelInfo, str, float) -> float
     # type: (ModelInfo, str, str) -> str
@@ -270,22 +270,50 @@
     Randomize a single parameter.
     """
-    if any(p.endswith(s) for s in ('_pd', '_pd_n', '_pd_nsigma', '_pd_type')):
-        return v
-
-    # Find the parameter definition
-    for par in model_info.parameters.call_parameters:
-        if par.name == p:
-            break
-    else:
-        raise ValueError("unknown parameter %r"%p)
-
+    if name.endswith('_pd'):
+        par = model_info.parameters[name[:-3]]
+        if par.type == 'orientation':
+            # Let oriention variation peak around 13 degrees; 95% < 42 degrees
+            return 180*np.random.beta(2.5, 20)
+        else:
+            # Let polydispersity peak around 15%; 95% < 0.4; max=100%
+            return np.random.beta(1.5, 7)
+
+    if name.endswith('_pd_n'):
+        # let pd be selected globally rather than per parameter
+        return 0
+
+    if name.endswith('_pd_type'):
+        # Don't mess with distribution type for now
+        return 'gaussian'
+
+    if name.endswith('_pd_nsigma'):
+        # type-dependent value; for gaussian use 3.
+        return 3.
+
+    if name == 'background':
+        return np.random.uniform(0, 1)
+
+    if name == 'scale':
+        return 10**np.random.uniform(-5,0)
+
+    par = model_info.parameters[name]
     if len(par.limits) > 2:  # choice list
         return np.random.randint(len(par.limits))
 
-    limits = par.limits
-    if not np.isfinite(limits).all():
-        low, high = parameter_range(p, v)
-        limits = (max(limits[0], low), min(limits[1], high))
-
+    if np.isfinite(par.limits).all():
+        return np.random.uniform(*par.limits)
+
+    if par.type == 'sld':
+        # Range of neutron SLDs
+        return np.random.uniform(-0.5, 12)
+
+    if par.type == 'volume':
+        if ('length' in par.name or
+                'radius' in par.name or
+                'thick' in par.name):
+            return 10**np.random.uniform(2,4)
+
+    low, high = parameter_range(par.name, value)
+    limits = (max(par.limits[0], low), min(par.limits[1], high))
     return np.random.uniform(*limits)
 
@@ -301,8 +329,10 @@
     :func:`constrain_pars` needs to be called afterward..
     """
-    with push_seed(seed):
-        # Note: the sort guarantees order `of calls to random number generator
-        random_pars = dict((p, _randomize_one(model_info, p, v))
-                           for p, v in sorted(pars.items()))
+    # Note: the sort guarantees order of calls to random number generator
+    random_pars = dict((p, _randomize_one(model_info, p, v))
+                       for p, v in sorted(pars.items()))
+    if model_info.random is not None:
+        random_pars.update(model_info.random())
+
     return random_pars
 
@@ -662,7 +692,13 @@
     parameters.
     """
-    result = run_models(opts, verbose=True)
-    if opts['plot']:  # Note: never called from explore
-        plot_models(opts, result, limits=limits)
+    limits = np.Inf, -np.Inf
+    for k in range(opts['sets']):
+        opts['pars'] = parse_pars(opts)
+        result = run_models(opts, verbose=True)
+        if opts['plot']:
+            limits = plot_models(opts, result, limits=limits, setnum=k)
+    if opts['plot']:
+        import matplotlib.pyplot as plt
+        plt.show()
 
 def run_models(opts, verbose=False):
@@ -743,5 +779,5 @@
 
 
-def plot_models(opts, result, limits=None):
+def plot_models(opts, result, limits=(np.Inf, -np.Inf), setnum=0):
     # type: (Dict[str, Any], Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float]
     base_value, comp_value= result['base_value'], result['comp_value']
@@ -758,13 +794,12 @@
     view = opts['view']
     import matplotlib.pyplot as plt
-    if limits is None and not use_data:
-        vmin, vmax = np.Inf, -np.Inf
-        if have_base:
-            vmin = min(vmin, base_value.min())
-            vmax = max(vmax, base_value.max())
-        if have_comp:
-            vmin = min(vmin, comp_value.min())
-            vmax = max(vmax, comp_value.max())
-        limits = vmin, vmax
+    vmin, vmax = limits
+    if have_base:
+        vmin = min(vmin, base_value.min())
+        vmax = max(vmax, base_value.max())
+    if have_comp:
+        vmin = min(vmin, comp_value.min())
+        vmax = max(vmax, comp_value.max())
+    limits = vmin, vmax
 
     if have_base:
@@ -813,10 +848,5 @@
         plt.title('Distribution of relative error between calculation engines')
 
-    if not opts['explore']:
-        plt.show()
-
     return limits
-
-
 
 
@@ -841,5 +871,5 @@
 VALUE_OPTIONS = [
     # Note: random is both a name option and a value option
-    'cutoff', 'random', 'nq', 'res', 'accuracy', 'title', 'data',
+    'cutoff', 'random', 'nq', 'res', 'accuracy', 'title', 'data', 'sets'
     ]
 
@@ -914,5 +944,5 @@
               if not arg.startswith('-') and '=' in arg]
     positional_args = [arg for arg in argv
-            if not arg.startswith('-') and '=' not in arg]
+                       if not arg.startswith('-') and '=' not in arg]
     models = "\n    ".join("%-15s"%v for v in MODELS)
     if len(positional_args) == 0:
@@ -959,4 +989,5 @@
         'title'     : None,
         'datafile'  : None,
+        'sets'      : 1,
     }
     engines = []
@@ -976,4 +1007,5 @@
         elif arg.startswith('-nq='):       opts['nq'] = int(arg[4:])
         elif arg.startswith('-res='):      opts['res'] = float(arg[5:])
+        elif arg.startswith('-sets='):     opts['sets'] = int(arg[6:])
         elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:]
         elif arg.startswith('-cutoff='):   opts['cutoff'] = float(arg[8:])
@@ -1008,4 +1040,8 @@
     # pylint: enable=bad-whitespace
 
+    # Force random if more than one set
+    if opts['sets'] > 1 and opts['seed'] < 0:
+        opts['seed'] = np.random.randint(1000000)
+
     if MODEL_SPLIT in name:
         name, name2 = name.split(MODEL_SPLIT, 2)
@@ -1020,75 +1056,6 @@
         return None
 
-    # Get demo parameters from model definition, or use default parameters
-    # if model does not define demo parameters
-    pars = get_pars(model_info, opts['use_demo'])
-    pars2 = get_pars(model_info2, opts['use_demo'])
-    pars2.update((k, v) for k, v in pars.items() if k in pars2)
-    # randomize parameters
-    #pars.update(set_pars)  # set value before random to control range
-    if opts['seed'] > -1:
-        pars = randomize_pars(model_info, pars, seed=opts['seed'])
-        if model_info != model_info2:
-            pars2 = randomize_pars(model_info2, pars2, seed=opts['seed'])
-            # Share values for parameters with the same name
-            for k, v in pars.items():
-                if k in pars2:
-                    pars2[k] = v
-        else:
-            pars2 = pars.copy()
-        constrain_pars(model_info, pars)
-        constrain_pars(model_info2, pars2)
-        print("Randomize using -random=%i"%opts['seed'])
-    if opts['mono']:
-        pars = suppress_pd(pars)
-        pars2 = suppress_pd(pars2)
-    if not opts['magnetic']:
-        pars = suppress_magnetism(pars)
-        pars2 = suppress_magnetism(pars2)
-
-    # Fill in parameters given on the command line
-    presets = {}
-    presets2 = {}
-    for arg in values:
-        k, v = arg.split('=', 1)
-        if k not in pars and k not in pars2:
-            # 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))))
-            return None
-        v1, v2 = v.split(MODEL_SPLIT, 2) if MODEL_SPLIT in v else (v,v)
-        if v1 and k in pars:
-            presets[k] = float(v1) if isnumber(v1) else v1
-        if v2 and k in pars2:
-            presets2[k] = float(v2) if isnumber(v2) else v2
-
-    # If pd given on the command line, default pd_n to 35
-    for k, v in list(presets.items()):
-        if k.endswith('_pd'):
-            presets.setdefault(k+'_n', 35.)
-    for k, v in list(presets2.items()):
-        if k.endswith('_pd'):
-            presets2.setdefault(k+'_n', 35.)
-
-    # Evaluate preset parameter expressions
-    context = MATH.copy()
-    context['np'] = np
-    context.update(pars)
-    context.update((k,v) for k,v in presets.items() if isinstance(v, float))
-    for k, v in presets.items():
-        if not isinstance(v, float) and not k.endswith('_type'):
-            presets[k] = eval(v, context)
-    context.update(presets)
-    context.update((k,v) for k,v in presets2.items() if isinstance(v, float))
-    for k, v in presets2.items():
-        if not isinstance(v, float) and not k.endswith('_type'):
-            presets2[k] = eval(v, context)
-
-    # update parameters with presets
-    pars.update(presets)  # set value after random to control value
-    pars2.update(presets2)  # set value after random to control value
-    #import pprint; pprint.pprint(model_info)
-
-    same_model = name == name2 and pars == pars
+    # TODO: check if presets are different when deciding if models are same
+    same_model = name == name2
     if len(engines) == 0:
         if same_model:
@@ -1106,19 +1073,4 @@
         del engines[2:]
 
-    use_sasview = any(engine == 'sasview' and count > 0
-                      for engine, count in zip(engines, [n1, n2]))
-    if use_sasview:
-        constrain_new_to_old(model_info, pars)
-        constrain_new_to_old(model_info2, pars2)
-
-    if opts['show_pars']:
-        if not same_model:
-            print("==== %s ====="%model_info.name)
-            print(str(parlist(model_info, pars, opts['is2d'])))
-            print("==== %s ====="%model_info2.name)
-            print(str(parlist(model_info2, pars2, opts['is2d'])))
-        else:
-            print(str(parlist(model_info, pars, opts['is2d'])))
-
     # Create the computational engines
     if opts['datafile'] is not None:
@@ -1142,11 +1094,93 @@
         'def'       : [model_info, model_info2],
         'count'     : [n1, n2],
-        'presets'   : [presets, presets2],
-        'pars'      : [pars, pars2],
         'engines'   : [base, comp],
+        'values'    : values,
     })
     # pylint: enable=bad-whitespace
 
     return opts
+
+def parse_pars(opts):
+    model_info, model_info2 = opts['def']
+
+    # Get demo parameters from model definition, or use default parameters
+    # if model does not define demo parameters
+    pars = get_pars(model_info, opts['use_demo'])
+    pars2 = get_pars(model_info2, opts['use_demo'])
+    pars2.update((k, v) for k, v in pars.items() if k in pars2)
+    # randomize parameters
+    #pars.update(set_pars)  # set value before random to control range
+    if opts['seed'] > -1:
+        pars = randomize_pars(model_info, pars)
+        if model_info != model_info2:
+            pars2 = randomize_pars(model_info2, pars2)
+            # Share values for parameters with the same name
+            for k, v in pars.items():
+                if k in pars2:
+                    pars2[k] = v
+        else:
+            pars2 = pars.copy()
+        constrain_pars(model_info, pars)
+        constrain_pars(model_info2, pars2)
+    if opts['mono']:
+        pars = suppress_pd(pars)
+        pars2 = suppress_pd(pars2)
+    if not opts['magnetic']:
+        pars = suppress_magnetism(pars)
+        pars2 = suppress_magnetism(pars2)
+
+    # Fill in parameters given on the command line
+    presets = {}
+    presets2 = {}
+    for arg in opts['values']:
+        k, v = arg.split('=', 1)
+        if k not in pars and k not in pars2:
+            # 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))))
+            return None
+        v1, v2 = v.split(MODEL_SPLIT, 2) if MODEL_SPLIT in v else (v,v)
+        if v1 and k in pars:
+            presets[k] = float(v1) if isnumber(v1) else v1
+        if v2 and k in pars2:
+            presets2[k] = float(v2) if isnumber(v2) else v2
+
+    # If pd given on the command line, default pd_n to 35
+    for k, v in list(presets.items()):
+        if k.endswith('_pd'):
+            presets.setdefault(k+'_n', 35.)
+    for k, v in list(presets2.items()):
+        if k.endswith('_pd'):
+            presets2.setdefault(k+'_n', 35.)
+
+    # Evaluate preset parameter expressions
+    context = MATH.copy()
+    context['np'] = np
+    context.update(pars)
+    context.update((k, v) for k, v in presets.items() if isinstance(v, float))
+    for k, v in presets.items():
+        if not isinstance(v, float) and not k.endswith('_type'):
+            presets[k] = eval(v, context)
+    context.update(presets)
+    context.update((k, v) for k, v in presets2.items() if isinstance(v, float))
+    for k, v in presets2.items():
+        if not isinstance(v, float) and not k.endswith('_type'):
+            presets2[k] = eval(v, context)
+
+    # update parameters with presets
+    pars.update(presets)  # set value after random to control value
+    pars2.update(presets2)  # set value after random to control value
+    #import pprint; pprint.pprint(model_info)
+
+    if opts['show_pars']:
+        if model_info.name != model_info2.name or pars != pars2:
+            print("==== %s ====="%model_info.name)
+            print(str(parlist(model_info, pars, opts['is2d'])))
+            print("==== %s ====="%model_info2.name)
+            print(str(parlist(model_info2, pars2, opts['is2d'])))
+        else:
+            print(str(parlist(model_info, pars, opts['is2d'])))
+
+    return pars, pars2
 
 def show_docs(opts):
@@ -1180,6 +1214,7 @@
     model = Explore(opts)
     problem = FitProblem(model)
-    frame = AppFrame(parent=None, title="explore", size=(1000,700))
-    if not is_mac: frame.Show()
+    frame = AppFrame(parent=None, title="explore", size=(1000, 700))
+    if not is_mac:
+        frame.Show()
     frame.panel.set_model(model=problem)
     frame.panel.Layout()
@@ -1191,5 +1226,6 @@
     if is_mac: frame.Show()
     # If running withing an app, start the main loop
-    if app: app.MainLoop()
+    if app:
+        app.MainLoop()
 
 class Explore(object):
@@ -1206,4 +1242,5 @@
         config_matplotlib()
         self.opts = opts
+        opts['pars'] = list(opts['pars'])
         p1, p2 = opts['pars']
         m1, m2 = opts['def']
@@ -1226,5 +1263,5 @@
         self.starting_values = dict((k, v.value) for k, v in pars.items())
         self.pd_types = pd_types
-        self.limits = None
+        self.limits = np.Inf, -np.Inf
 
     def revert_values(self):
@@ -1282,8 +1319,12 @@
     """
     opts = parse_opts(argv)
+    if opts['seed'] > -1:
+        print("Randomize using -random=%i"%opts['seed'])
+        np.random.seed(opts['seed'])
     if opts is not None:
         if opts['html']:
             show_docs(opts)
         elif opts['explore']:
+            opts['pars'] = parse_pars(opts)
             explore(opts)
         else:
Index: sasmodels/modelinfo.py
===================================================================
--- sasmodels/modelinfo.py	(revision 724257c92ca5f608f1dc18211ee76a0f7d2fe6f1)
+++ sasmodels/modelinfo.py	(revision 0bdddc2288f8fc5bc5e416af84e456e811c3b733)
@@ -467,4 +467,13 @@
                          if p.polydisperse and p.type not in ('orientation', 'magnetic'))
         self.pd_2d = set(p.name for p in self.call_parameters if p.polydisperse)
+
+    def __getitem__(self, key):
+        # Find the parameter definition
+        for par in self.call_parameters:
+            if par.name == key:
+                break
+        else:
+            raise KeyError("unknown parameter %r"%key)
+        return par
 
     def _set_vector_lengths(self):
@@ -754,4 +763,5 @@
     info.opencl = getattr(kernel_module, 'opencl', not callable(info.Iq))
     info.single = getattr(kernel_module, 'single', not callable(info.Iq))
+    info.random = getattr(kernel_module, 'random', None)
 
     # multiplicity info
Index: sasmodels/models/adsorbed_layer.py
===================================================================
--- sasmodels/models/adsorbed_layer.py	(revision b0c42719ac9096c35137bf644db59e1c5915b06b)
+++ sasmodels/models/adsorbed_layer.py	(revision 0bdddc2288f8fc5bc5e416af84e456e811c3b733)
@@ -94,4 +94,17 @@
 Iq.vectorized = True  # Iq accepts an array of q values
 
+def random():
+    # only care about the value of second_moment:
+    #    curve = scale * e**(-second_moment^2 q^2)/q^2
+    #    scale = 6 pi/100 (contrast/density*absorbed_amount)^2 * Vf/radius
+    # the remaining parameters can be randomly generated from zero to
+    # twice the default value.
+    import numpy as np
+    pars = dict(
+        scale=1,
+        second_moment=10**np.random.uniform(1, 3),
+    )
+    return pars
+
 # unit test values taken from SasView 3.1.2
 tests = [
Index: sasmodels/models/bcc_paracrystal.py
===================================================================
--- sasmodels/models/bcc_paracrystal.py	(revision 69e1afc57270f67b2a0d7157a6d263fc2db98cc0)
+++ sasmodels/models/bcc_paracrystal.py	(revision 0bdddc2288f8fc5bc5e416af84e456e811c3b733)
@@ -81,5 +81,5 @@
 .. figure:: img/parallelepiped_angle_definition.png
 
-    Orientation of the crystal with respect to the scattering plane, when 
+    Orientation of the crystal with respect to the scattering plane, when
     $\theta = \phi = 0$ the $c$ axis is along the beam direction (the $z$ axis).
 
@@ -131,4 +131,24 @@
 source = ["lib/sas_3j1x_x.c", "lib/gauss150.c", "lib/sphere_form.c", "bcc_paracrystal.c"]
 
+def random():
+    import numpy as np
+    # Define lattice spacing as a multiple of the particle radius
+    # using the formulat a = 4 r/sqrt(3).  Systems which are ordered
+    # are probably mostly filled, so use a distribution which goes from
+    # zero to one, but leaving 90% of them within 80% of the
+    # maximum bcc packing.  Lattice distortion values are empirically
+    # useful between 0.01 and 0.7.  Use an exponential distribution
+    # in this range 'cuz its easy.
+    dnn_fraction = np.random.beta(a=10, b=1)
+    pars = dict(
+        #sld=1, sld_solvent=0, scale=1, background=1e-32,
+        radius=10**np.random.uniform(1.3, 4),
+        d_factor=10**np.random.uniform(-2, -0.7),  # sigma_d in 0.01-0.7
+    )
+    pars['dnn'] = pars['radius']*4/np.sqrt(3)/dnn_fraction
+    #pars['scale'] = 1/(0.68*dnn_fraction**3)  # bcc packing fraction is 0.68
+    pars['scale'] = 1
+    return pars
+
 # parameters for demo
 demo = dict(
Index: sasmodels/models/broad_peak.py
===================================================================
--- sasmodels/models/broad_peak.py	(revision 43fe34b36fb4d68e5cf894e939f88cc312c730be)
+++ sasmodels/models/broad_peak.py	(revision 0bdddc2288f8fc5bc5e416af84e456e811c3b733)
@@ -94,4 +94,20 @@
 Iq.vectorized = True  # Iq accepts an array of q values
 
+def random():
+    import numpy as np
+    pars = dict(
+        scale=1,
+        porod_scale=10**np.random.uniform(-8, -5),
+        porod_exp=np.random.uniform(1, 6),
+        lorentz_scale=10**np.random.uniform(0.3, 6),
+        lorentz_length=10**np.random.uniform(0, 2),
+        peak_pos=10**np.random.uniform(-3, -1),
+        lorentz_exp=np.random.uniform(1, 4),
+    )
+    pars['lorentz_length'] /= pars['peak_pos']
+    pars['lorentz_scale'] *= pars['porod_scale'] / pars['peak_pos']**pars['porod_exp']
+    #pars['porod_scale'] = 0.
+    return pars
+
 demo = dict(scale=1, background=0,
             porod_scale=1.0e-05, porod_exp=3,
